The r-modes of slowly rotating, stratified neutron stars

The only r-modes that exist in a globally barotropic, rotating, Newtonian star are the fundamental $l = |m|$ solutions, where $l$ and $m$ are the indices of the spherical harmonic $Y_l^m$ that describe the mode's angular dependence. This is in stark contrast to a stellar model that is non-barotropic throughout its interior, which hosts all the $l \geq |m|$ perturbations including radial overtones. In reality, neutron stars are stratified with locally barotropic regions. Therefore, we explore how stratification alters a star's ability to support r-modes. We consider the globally stratified case and examine the behaviour of the modes as the star gets close to barotropicity. In this limit, we find that all but the fundamental $l = |m|$ perturbations change character and become generic inertial modes. Restricting the analysis to $l = |m|$ perturbations, we develop the r-mode equations in order to consider stellar models that exhibit local barotropicity. Our results for such models show that the r-mode overtones diverge and join the inertial modes. In order to see which r-modes persist and retain their character in realistic neutron stars, these calculations will need to be brought into full general relativity.


INTRODUCTION
Stellar interiors host a rich spectrum of oscillation modes (Cox 1980;Unno et al. 1989).These perturbations are sensitive to the characteristics of the star.Indeed, each ingredient of physics -density, stratification, rotation, magnetic field etc. -corresponds (more or less) directly to a unique set of oscillation modes.It would seem, however, that not all oscillation modes are created equal.An interesting family of modes that arise due to stellar rotation are the r-modes.These Coriolis-dominated perturbations possess the remarkable property of being generically unstable in a perfect-fluid star due to the emission of gravitational radiation (Andersson 1998;Friedman & Morsink 1998).This result inspired a wide body of literature (see reviews Andersson & Kokkotas 2001;Andersson 2003), including suggestions that the r-mode instability may limit the rotation rates of newly born pulsars (Lindblom et al. 1998;Andersson et al. 1999a) and more mature accreting systems (Bildsten 1998;Andersson et al. 1999b).The r-modes are also candidates for gravitational-wave observations, with a recent search focused on the glitching neutron star PSR J0537-6910 (Abbott et al. 2021).However, so far no gravitational-wave signatures consistent with an r-mode instability have been seen.
It was Cowling (1941), in his work on non-rotating polytropes, who provided the first classification of modes according to the physics dominating their behaviour.The simplest stellar model is that of a spherical, incompressible fluid.Such a star will only have one family of oscillation modes, the fundamental f -modes.The f -modes are distinguished by having no nodes in their radial eigenfunctions and inducing large density perturbations in the star.Should one allow for compressibility of the stellar fluid, by including an equation of state, then the p-modes arise.These are high-frequency acoustic waves restored by the pressure of the fluid, also associated with large perturbations in the density.Suppose the star becomes stratified such that the matter is no longer barotropic. 1 Then g-modes will appear, restored by gravity gradients.The g-modes have low frequencies and small density perturbations.
These three families of fluid oscillations, the f -, p-and g-modes, all belong to the class of polar modes.A perturbation is polar (spheroidal or even parity) if it varies like a spherical harmonic    under a parity transformation.In a spherically symmetric, fluid star, this is the only class of modes that exists.Furthermore, the perturbations of a spherical star are especially simple as each mode is associated with a single    .The class of axial (toroidal or odd parity) modes, which transform opposite under parity to polar modes, require some level of anisotropy in order to exist.
When the star begins to rotate, a new family of modes appears.These are known as the inertial modes and they are restored by the Coriolis force (Bryan 1889;Greenspan 1964;Lindblom & Ipser 1999;Rieutord 2001).Even in the non-rotating limit, the inertial modes generically involve a mix of polar and axial perturbations with definite parity and couple multiple spherical harmonics (Lee et al. 1992;Lockitch & Friedman 1999;Yoshida & Lee 2000b).Formally, these modes inhabit a zero-frequency subspace on the spherical star, being stationary convective fluid currents, and attain non-zero oscillation frequencies at first order in the star's angular frequency.Among these modes exists a special subclass, the purely axial r-modes.
The r-modes were first studied in the context of astrophysics by Papaloizou & Pringle (1978), who named them for their similarity to the Rossby waves of terrestrial meteorology (see Zaqarashvili et al. 2021 for a recent review).Contrary to a typical inertial mode, an r-mode is an axial perturbation associated with a single    on the spherical star.Moreover in Newtonian gravity, their leading-order frequency in a slow-rotation expansion can be determined analytically to be in the rotating frame of the star with angular frequency Ω.Since the frequency of a mode as measured by an inertial observer outside the star is simply  0 − Ω, the r-modes are retrograde in the frame of the star, but prograde in the inertial frame at all rates of rotation.Therefore, they satisfy the well-known Chandrasekhar-Friedman-Schutz instability criterion and are generically unstable to perturbations driven by gravitational radiation (Chandrasekhar 1970;Friedman & Schutz 1978;Andersson 1998;Friedman & Morsink 1998).
To date, most of the r-mode calculations have involved Newtonian stellar models that are globally barotropic (Lindblom et al. 1998;Lockitch & Friedman 1999;Yoshida & Lee 2000b) or globally non-barotropic (Provost et al. 1981;Smeyers et al. 1981;Saio 1982;Andersson et al. 1999a;Yoshida & Lee 2000a).The barotropic star is especially simple: there exist only the fundamental  = || r-modes in the star and the leading-order axial eigenfunctions can be obtained analytically, independent of the equation of state.However, the non-barotropic case is more complicated.One must work to beyond leading order in rotation, where the star departs from spherical symmetry (Papaloizou & Pringle 1978;Provost et al. 1981;Saio 1982).The equations then admit a standard eigenvalue problem for the  || r-mode solutions.The reason for the concentration on Newtonian stars has been in part due to the challenges in calculating the r-modes in general relativity.
Our particular focus is on neutron stars, which are highly relativistic bodies.In order to involve a realistic description for the nuclear matter, we need to formulate the r-mode problem beyond Newtonian gravity.In relativity, there are no longer any purely axial inertial modes in barotropes (except for stationary dipole perturbations; Lockitch 1999).In this direction, there have been a number of calculations of the relativistic inertial modes (Lockitch et al. 2000(Lockitch et al. , 2003;;Ruoff et al. 2003), including physically motivated equations of state (Idrisy et al. 2015).The relativistic inertial modes may or may not be a reasonable approximation of the problem.In reality, neutron stars are stratified due to varying chemical composition (Reisenegger & Goldreich 1992;Andersson & Pnigouras 2019).The r-modes exist in non-barotropic, relativistic stars.However, the relativistic perturbation equations imply a continuous spectrum (Kojima 1998;Beyer & Kokkotas 1999).This is surprising, since the r-modes have well-defined frequencies in Newtonian gravity.It is unclear whether the continuous spectrum is physical or an artefact of some simplifying assumptions.Adopting the latter view, there have been efforts to regularise the problem (Lockitch et al. 2004;Pons et al. 2005), as well as studies using the Cowling approximation (Kojima & Hosonuma 1999;Kraav et al. 2022a,b).
In this paper, our goal is to study the role of stratification for the r-modes of Newtonian stars.In particular, we want to move beyond global assumptions about the matter and consider the more realistic case where the fluid may be locally barotropic.This is expected to be the case for neutron stars; their high-density cores will likely have composition gradients, whereas the outer layers will be barotropic (since matter at low densities is composed of single nuclei).Our hope is that, in understanding the Newtonian problem, we may make progress towards calculating the relativistic r-modes.
This paper is organised as follows.We begin in Section 2 with a brief discussion on constructing slowly rotating stellar models, which will form the background for the r-mode oscillations.We move on to Section 3 to describe the perturbation formalism for rotating stars and present the r-mode equations.We develop these equations into an eigenvalue problem in Section 4, assuming that the star is globally non-barotropic.These expressions are used to calculate the r-modes as the star gets close to the barotropic limit.In Section 5, seeking to consider more realistic stellar models, we derive the  = || system of equations, where the matter may be locally barotropic.We implement some neutron-star equations of state for the perturbations with a polytropic background and compute the oscillations.Finally, we summarise and suggest future directions in Section 6.

SLOWLY ROTATING BACKGROUND
We will examine the r-modes of a slowly rotating, Newtonian star.In principle, this limits the rates of rotation our analysis is accurate to.However, this will enable us to explore the character of the modes and, in practice, many stars fall comfortably within the slow-rotation regime.
The structure of a uniformly rotating star is a solution to the following system of equations: with an equation of state where is the effective potential -the sum of the gravitational Φ and centrifugal potentials - is the mass density and  is the pressure.Here, we use spherical polar coordinates (, , ), where the star rotates with angular velocity Ω  about the  = 0 axis,   is a Legendre polynomial of degree  and ∇  is the covariant derivative.We assume the star to be rotating slowly such that where  is the total mass of the star and  is the radius of its corresponding spherical configuration. 2 In this context, a perturbative approach to solving equations ( 2) is permitted, where the departure from sphericity is small. 3n a spherical star, the surfaces of constant density, pressure and gravitational potential coincide and depend solely on .This is not so when the star begins to rotate, as the centrifugal force spoils the symmetry with respect to the polar angle .Therefore, in calculating rotating configurations, it is convenient to adjust the coordinate system in the following way.We use the Clairaut-Legendre expansion (see, e.g., Section 5.2 of Tassoul 1978) and introduce a new coordinate  that corresponds to the isobaric surfaces such that (, ) = () (and thus also corresponds to the isopycnic surfaces).By the Euler equation (2a), these surfaces will also coincide with the surfaces of constant Ψ.Since the star rotates slowly, we may define where  =  ( 2 ) characterises the deviation from spherical symmetry.(That the expansion only involves even powers of  is a consequence of the rotational symmetry.)This coordinate change is associated with the metric tensor    defined by the line element We note that the metric differs from the usual spherical polar coordinates at  ( 2 ) and the coordinate basis is no longer orthogonal.We can see from (3) that the centrifugal potential only has contributions from the  = 0, 2 Legendre polynomials.Hence, the problem decouples into these two sectors and we must have The  = 2 sector characterises the shape of the star and reduces to Clairaut's equation where  0 =  (1) denotes the mass distribution of the non-rotating configuration.We also note  0 =  (1),  0 =  (1) and Φ 0 =  (1) as the mass density, pressure and gravitational potential of the spherical star, respectively.An examination of the behaviour at the centre shows that  2 approaches a constant.The second boundary condition comes from matching the interior solution for Φ to the exterior solution that decays as 1/ 3 .Therefore, we find Equipped with the non-spherical shape of the slowly rotating background to second order in rotation, we go on to formulate the perturbation problem.

THE PERTURBATIONS
Oscillation modes are harmonic solutions to the following linearised equations of motion (expressed in the rotating frame): with an equation of state for the perturbations Δ where  and Δ denote the Eulerian and Lagrangian variations of a quantity, respectively,   is the Lagrangian displacement vector of the fluid elements and Γ 1 = Γ 1 () is the adiabatic index of the perturbations.The stratification of the fluid enters the description through the linearised equation of state (10d).As a starting assumption, it is common to choose Γ 1 = const.This leads to qualitative insight, but one has to be careful in drawing quantitative conclusions.In general, when Γ 1 differs from the adiabatic index of the background the perturbations obey a different equation of state and thus the star is non-barotropic.
Since the equilibrium configuration is axisymmetric with respect to the azimuthal coordinate , each mode will have a definite order  and we may assume the following form for the Lagrangian displacement in our (, , ) coordinate basis: where  is the angular frequency of the mode and    (, ) is a spherical harmonic.In the Ω → 0 limit, the displacement vector (12) will tend towards the familiar vector decomposition in (, , ) coordinates (see, e.g., Lockitch & Friedman 1999), where   and   are the polar functions and   is the axial function.Indeed, since a mode's parity does not change as it varies continuously along a sequence of equilibrium configurations, beginning with a spherical star and with increasing rotation, it is appropriate to identify (  ,   ) and   with the polar and axial perturbations, respectively, as they correspond to these classes on a spherical star. 4The scalar perturbations are consequently decomposed as

Stratification on the spherical star
Although we want to formulate the r-mode problem on a rotating star, stratification already plays an important role on the non-rotating star.It is instructive to consider this case briefly.
A barotropic (Γ 1 = Γ), spherically symmetric star only admits f -and p-modes.It can be shown that there are time-independent solutions to the Ω = 0 perturbation equations (10) with vanishing perturbed mass densities that are purely polar or axial in nature (Lockitch & Friedman 1999) (see also Lockitch et al. 2000 for the corresponding result in relativity).These stationary currents are associated with the polar g-modes and the axial r-modes and they reside in the zero-frequency subspace on such a star.If the star spins up, these stationary currents will become oscillatory, being restored by the Coriolis force.In general, these polar and axial perturbations will mix, forming the inertial modes.However, there will persist perturbations that are purely axial at zeroth order in rotation, associated with a definite (, ).These solutions are the fundamental  = || r-modes.
Suppose we now consider a stratified (Γ 1 ≠ Γ), spherical star.Alongside the f -and p-modes, it will host g-modes with non-zero frequencies, supported by the buoyancy.The only trivial solutions that exist in this case are axial.Therefore, when the star rotates, the only inertial modes that appear are the r-modes.Because the axial perturbations do not have stationary polar currents to mix with, a non-barotropic, rotating star has the complete set of  || r-modes including radial overtones.
In between these two extremal cases, there is a third regime where the stratification and rotation are of the same order of magnitude.This is a form of weak (but non-zero) stratification relative to the rotation.We discuss this scenario in more detail in Andersson & Gittins (2022).
Neutron stars are, in general, non-barotropic since the chemical composition changes throughout the interior.This is illustrated in Fig. 1 for the two realistic equations of state, BSk19 and BSk21 (Fantina et al. 2013;Potekhin et al. 2013).The exact chemical composition of neutron-star interiors is at present unknown and is related to the nuclear reactions going on under the surface (Reisenegger & Goldreich 1992;Andersson & Pnigouras 2019).However, we understand low-density nuclear matter quite well and expect the outer layers to be barotropic.(This is shown in Fig. 1 for both equations of state where Γ 1 = Γ at low densities.)

Slow rotation
Moving beyond the spherical star, in a slow-rotation expansion, an r-mode has the ordering5 To uniquely identify the r-modes from the perturbation equations ( 10), it is sufficient to look for solutions of the form ( 14).By rotational symmetry, we can assume the series expansion for the mode frequency (Papaloizou & Pringle 1978;Smeyers et al. 1981) where  0 =  () and  2 =  ( 3 ).Clearly, the validity of this expansion relies on | 2 / 0 | =  ( 2 ) 1. (This is a feature we will pay close attention to later.)The remaining terms we need only calculate to leading order.
In addition to the fact that the r-modes are purely axial, they are special among the inertial modes in that their frequency at leading order is analytic and independent of the equation of state.This can be seen from the angular components of the linearised Euler equation (10b).At  ( 2 ), we find the simple expression There are three solutions to equation ( 16): (i) the equation permits a zero-frequency solution with  0 = 0; (ii) the axial functions are   = 0 for all  ; or (iii) a single   survives with frequency (1).The physically interesting case that corresponds to an r-mode is solution (iii).We note that there are no axisymmetric ( = 0) r-modes and all the perturbations travel retrograde to the rotation of the star.However, at this order, the axial function   is undetermined.To calculate these functions, we need to develop the perturbation equations ( 10) into an eigenvalue problem for the frequency correction  2 .
The equations we need were first derived by Saio (1982).We will focus on the simple adiabatic case appropriate for (cold) neutron stars.In this context, the oscillations will be normal modes with manifestly real frequencies, since there are no dissipative effects.As is typically useful in the numerical computation of oscillation modes, we define the following dimensionless variables (Unno et al. 1989): where  = Φ 0 / is the gravitational acceleration.Rotation couples spherical harmonics with different values of  .In order to see this, one may use the standard recurrence relations where It is useful to note the orthogonality of the spherical harmonics where Ω is the element of solid angle.Using these identities, the linearised equations ( 10) with the ordering ( 14) provide the following system of differential equations: with the algebraic relations where is the Schwarzschild discriminant (useful for characterising the stratification of the fluid) and we have defined the dimensionless terms Equations ( 21) describe an r-mode and there are a few remarks worth making: (i) the leading-order axial perturbation, associated with the spherical harmonic    , sources  ± 1 polar perturbations at  ( 2 ) (as well as  ± 2 axial perturbations at this order, as we show in Appendix A); (ii) since   = 0 for  = ||, the  − 1 polar terms must vanish and thus the lowest degree mode that exists is  = || [this has been assumed in the decomposition of equations ( 12) and ( 13)]; and (iii) the only rotational quantity from the background that comes into the mode equations is the non-spherical shape  2 .At this point, no assumptions have been made about the state of the matter.
In the globally barotropic case Γ 1 = Γ, where  = 0, the linearised equations ( 21) admit a simple solution.We can combine (21c) with ( 21g) and ( 21d) with (21h) to obtain two equations that an r-mode must satisfy, Equations ( 24) only provide a consistent, non-trivial solution when  = || and  6, || ∝  ||−1 .This is the fundamental r-mode that we alluded to previously.Thus, there are no overtones when the star is barotropic and there are no  > || solutions.The higher order features, including the other eigenfunctions and eigenvalue ω2 , are determined from the remaining  + 1 equations ( 21).

GLOBALLY NON-BAROTROPIC STARS
We want to express the r-mode equations ( 21) in a form that is suitable for integration.As it stands, we are unable to determine  5,±1 , since we only have one algebraic relation (21i) for these two functions that appear in equations ( 21a) and (21b).To begin with, we will assume the star to be stably stratified throughout its interior such that Γ 1 ≠ Γ. 6 As we have noted above, this assumption is inappropriate for realistic stars with some degree of local barotropicity, but we can learn something about how the r-mode solutions behave as one gets close to the barotropic limit.
Since we assume that  ≠ 0, this allows us to eliminate  1,−1 from the system.Here, we note that ( 27) implies that  2,+1 = 0 in barotropic regions, for  > ||.Next, we can differentiate ( 27) and combine the result with the differential equations ( 21a) and (21b) to obtain Along with equation (21i), we are in a position to decouple  5,±1 .
We now have enough information to remove the variables  1,−1 ,  2,−1 ,  5,±1 and  6, from equations ( 21).We end up with where In order to solve the eigenvalue problem (29), we must provide boundary conditions that constrain the solutions.At the stellar centre, the functions must be well behaved and regular.As  → 0, the regular solutions are given by (see Section 18.1 of Unno et al. 1989) Inserting (31) into the perturbation equations ( 29), we find the following boundary conditions at the centre of the star: At the surface, the Lagrangian variation of the pressure must vanish Δ/ 0 = 0 and the perturbed gravitational potential must match smoothly to the exterior solution, which decays as Φ  ∝ 1/  +1 for a given multipole  .Thus, for  = , we must have Equations ( 29) constitute six first-order ordinary differential equations, supplemented by three boundary conditions at the centre (32) and three at the surface (33).The boundary conditions will only be satisfied for an eigenvalue ω2 .
Table 1.The r-mode eigenfrequencies ω2 / 3 of the  = 1 polytrope with Γ 1 = 5/3.We solve this eigenvalue problem using the following numerical approach (our technique is similar to that used by Lindblom & Detweiler 1983).Since the system of equations ( 29) is linear (by construction, as we are using first-order perturbation theory), we may express it as with a matrix Q and abstract vector field Y = ( 1,+1 ,  2,+1 ,  3,+1 ,  3,−1 ,  4,+1 ,  4,−1 ).This is a sixth-order system of linear equations, so there will exist six linearly independent solutions for a given (, , ω2 ).However, only for specific values of ω2 will the linearly independent solutions combine to satisfy all the boundary conditions; these are the eigenfrequencies.At the centre of the star, there are three linearly independent vectors Y( → 0) that satisfy boundary conditions (32).We select three such initial vectors and integrate them using equation ( 34) to a point 0 <  0 <  in the star.This generates three solutions Y 1 , Y 2 , Y 3 defined on the domain 0 <   0 , each of which satisfy the central boundary conditions (32).In a similar fashion, we also produce three linearly independent solutions Y 4 , Y 5 , Y 6 for the region  0   out of the surface boundary conditions (33).Therefore, we obtain the general solution with constants  1 , . . .,  6 .Hence, ω2 is an eigenvalue if and only if This matching can be written as the matrix equation A • x = 0, where A depends non-linearly on ω2 .For non-trivial eigenvalues, the matrix A is singular.Thus, we look for values of ω2 , for a given (, ), such that det A = 0.
Once the eigenfrequency ω2 is determined, the eigenfunctions may be calculated.Since we are calculating normal modes, there is a free amplitude in the functions.We choose to normalise the modes at the surface by  6, () = 1. (37) We implement this normalisation by replacing a row in A in favour of this condition, Ã • x = b, where Ã becomes a non-singular matrix and b is a known column vector.This concludes the discussion of our numerical method.
For our results, we assume a polytropic equation of state for the equilibrium star where  and  are the polytropic constant and index, respectively.We consider  = 1 to approximate a neutron star and obtain the shape  2 from the results of Chandrasekhar & Lebovitz (1962).
As a consistency check of our computational technique, we first consider Γ 1 = 5/3 (also employed by Provost et al. 1981;Saio 1982). 7The results are summarised by Table 1, where the overtones, with nodes in their displacement eigenfunctions, are denoted by .We can compare the  = 3 eigenfrequencies with Provost et al. (1981), who calculated within the Cowling approximation.Our results are compatible within the estimated errors of the approximation.Additionally, we considered the non-analytic  = 3 polytrope, numerically solving (8) for the shape, to compare with Saio (1982), finding excellent agreement.
The system of equations ( 29) is constrained to stellar models that have Γ 1 ≠ Γ. (Indeed, Saio 1982 was aware of this limitation and so only considered  = || r-modes for more realistic stellar models in his calculation.)Although we are unable to consider barotropic stars, we are in a position to explore what happens to the r-mode solutions as the stellar model tends towards barotropicity.Motivated by Fig. 1, we will consider the range 2 < Γ 1 2.25 to approximate  = 1 neutron stars.In Figs.2-4, we show how the eigenfrequencies vary in the Γ 1 → Γ limit.All but the fundamental ( = 0)  = || eigenfrequencies diverge (and the divergences are worse for the higher overtones).That is, the frequency correction ω2 grows exponentially as the star becomes more barotropic, showing that care must be taken with the frequency expansion (15) in this limit.If ω2 becomes comparable in magnitude to ω0 , the frequency will no longer satisfy the Euler equation ( 16) at leading order, which in turn will result in a breakdown in the assumed ordering (14).
This dependence on  is important to note since, from the outset, we merely assumed slow rotation according to (4).Now, it is evident that, as we approach barotropicity, the non-fundamental  = || r-modes only have support at even slower rotation rates.This is an additional constraint on the solutions.In general, if | ω2 /( 2 ω0 )| 1, then there is an additional dependence on the rotation  in addition to slow rotation.Here, we see the competition between the rotation  of the star and its stratification, parametrised by the adiabatic indices Γ and Γ 1 , in supporting r-mode oscillations.For rotations of  ∼ 0.1 (appropriate for rapidly rotating neutron stars), the frequency corrections of the non-fundamental  = || solutions become ω2 ∼  () in the Γ 1 → Γ limit.Hence, the leading-order frequencies are no longer simply given by equation (1).
In Figs. 5 and 6, we present the eigenfunctions for the  = 0 and  = 1 (, ) = (2, 2) r-modes, respectively.As we witnessed for the eigenfrequencies, the eigenfunctions of the  = 0 r-mode in Fig. 5 retain their character and are well behaved in the barotropic limit.However, the  = 1 overtone has markedly different behaviour, shown in Fig. 6.As Γ 1 → Γ, the polar displacement eigenfunctions  1,+1 and  5,+1 The eigenfrequencies ω2 of the  = 2 r-modes of the  = 1 polytrope with varying Γ 1 .Note that only the (, ) = (2, 0) solution remains well behaved in the barotropic limit.The other solutions diverge, implying that ω2 is promoted to the same order as the leading-order frequency ω0 and the frequency expansion is spoiled.
The eigenfrequencies ω2 of the (, ) = (1, 1) r-modes of the  = 1 polytrope with varying Γ 1 .Note that only the  = 0 solution remains well behaved in the barotropic limit.The other solutions diverge, implying that ω2 is promoted to the same order as the leading-order frequency ω0 and the frequency expansion is spoiled.diverge, again presenting an issue with the assumed ordering ( 14) for moderate rotation rates .We note that we saw the same qualitative features in all the r-modes we calculated, although the divergences had varying levels of severity (as was the case for the eigenfrequencies).
Indeed, for  ∼ 0.1,  1,+1 and  5,+1 get promoted to  (1) perturbations like  6, such that the displacement at leading order becomes a mix of polar and axial functions.We show this feature analytically in Andersson & Gittins (2022).The scalar perturbations remain at  ( 2 ).This, along with the divergence of the eigenfrequency ω2 ∼  (), meets the definition of a generic inertial mode.Therefore, the  1,  = || perturbations and the  > || perturbations become similar in character to inertial modes as the matter becomes barotropic.

ALLOWING FOR LOCAL BAROTROPICITY
As we have discussed, equations (29) assume the star to be globally non-barotropic such that Γ 1 ≠ Γ.In reality, neutron stars will have regions where the matter is barotropic (see Fig. 1).We want to explore the extent to which we can drop this assumption and move towards more realistic matter models.We begin with equations ( 21), where no assumption about the matter has been made.
One problem immediately rears its ugly head.Without using the algebraic relation ( 27), which has divergent behaviour as Γ 1 → Γ, one is  unable to decouple the functions  5,±1 that appear in equations ( 21a), ( 21b) and (21i).But, given what we found in Section 4, this may not be too surprising.As a star goes from being globally non-barotropic to globally barotropic, all but the fundamental  = || r-mode solutions diverge in such a way that they change character and become like generic inertial modes.Should this result hold if the star becomes locally barotropic (even, say, at a point), then it stands to reason that the perturbation equations constructed from the assumed ordering ( 14) will not admit any solutions, since no such modes would exist.This seems to be the situation we find ourselves in and we will provide further evidence for this in this section.
As one might expect, the problem simplifies in the  = || case.We know that a barotropic star only has one solution for a given  = || [the fundamental r-mode; see equations ( 24)].Therefore, we can examine whether a star that is locally barotropic supports  = || r-modes with overtones  1.
These modes have no  − 1 couplings.Equations ( 21g) and (21i) enable us to remove  5, ||+1 and  6, || from the system of equations.Hence, we obtain the following system for the  = || r-modes: The behaviour of the eigenfunctions at the centre is given by (31) and the  + 1 surface conditions from (33) remain the same.
For our calculation, we use the ratio Γ 1 /Γ from the BSk19 and BSk21 nuclear-matter models shown in Fig. 1, where Γ is obtained from the stellar background.These equations of state depend on the number density of baryons  b in the neutron star and therefore introduce dimensionality to the problem.We assume our  = 1 stellar model to have  = 1.4 M and  = 10 km.The eigenfrequencies we calculate are listed in Table 2.We first note that the eigenfrequencies with  1 are larger in magnitude for the BSk19 model than BSk21.This is related to the fact that BSk19 is more weakly stratified and is consistent with our results above with Γ 1 = const.Some of the eigenfrequencies of the BSk19 and BSk21 equations of state do not exhibit particularly strong divergences.As we expect, the  = 0 r-modes have reasonable values of ω2 .However, we find that all the modes with  1 that we considered have issues with their eigenfunctions.As representative examples, we show the eigenfunctions of the  = 0 and  = 1 (, ) = (||, 3) r-modes in Figs.7 and 8, respectively.These two modes do not have strong divergences in ω2 (see Table 2).(Although, the  = 1 solution will begin to breakdown at spins of  2 ∼ 0.1.)The eigenfunctions of the  = 0 r-mode in Fig. 7 are perfectly well behaved and seem to be relatively insensitive to the linearised equation of state.This is in contrast to the  = 1 solution, shown in Fig. 8, which has divergent behaviour in the polar displacement functions  1, ||+1 and  5, ||+1 .Thus, violating the assumed ordering (14) at reasonable rates of rotation.To summarise, we do not find any well-behaved solutions with  1; all these solutions have divergences in the  1, ||+1 and  5, ||+1 eigenfunctions and many also have divergences in their eigenfrequency ω2 .
Complementary to what we saw for the Γ 1 ≠ Γ calculation in Section 4, our results for the BSk19 and BSk21 models show that most of the r-mode solutions diverge when the star hosts barotropic regions.The only solution that retains its character is the fundamental  = || r-mode.

CONCLUSIONS
We have considered the role stratification plays in supporting r-mode oscillations on slowly rotating, Newtonian stars.We focused on stellar models approximating neutron stars.However, our qualitative results should be the same for other stars that are locally barotropic.
In using the linearised equations derived by Saio (1982), which are only valid for stellar models that are globally non-barotropic Γ 1 ≠ Γ, we found that the majority of solutions exhibited divergences in the barotropic limit Γ 1 → Γ.These divergences occur for the  1,  = || r-modes and the  > || r-modes and manifest themselves in two ways.(1) The frequency correction ω2 becomes comparable in magnitude to the linear term ω0 such that ω2 / ω0 ∼  (1) for modest rates of rotation .This spoils the usual identification of the r-mode frequency at leading order.(2) The polar terms in the displacement vector increase to  (1); the same order as the axial term.Thus, the displacement becomes a mixture of polar and axial functions at zeroth order in rotation.These two divergences lead to the solutions becoming generic inertial modes.However, as one may expect, none of these features are seen in the fundamental ( = 0)  = || r-mode solutions, which remain well behaved in the non-stratified limit.
Moving beyond the globally non-barotropic approximation, we considered stars that have varying stratification with barotropic regions.The perturbation equations cannot be expressed as a standard eigenvalue problem for general  || r-modes if the star is barotropic at some point.However, the situation becomes tractable for  = || modes.We calculated the  = || r-modes of an  = 1 polytrope with the perturbations described by the BSk19 and BSk21 nuclear-matter equations of state.Numerically, we obtained solutions in addition to the fundamental r-mode with  1.Although these numerical solutions exist, they also present divergent behaviour in ω2 and the polar displacement functions.Formally, these results confirm the expectation that for stratified stars there will be a critical rotation rate above which the fluid only supports the fundamental  = || r-modes.In addition, we have shown that the same is true for stars that are barotropic in a local region.The remaining Coriolis-driven perturbations must join the general inertial-mode family.Furthermore, depending on how stratified the equation of state is, the rotation at which the solutions change character can be very modest indeed.This implies that rapidly rotating neutron stars, which are of interest for gravitational-wave observations, may only host the fundamental  = || r-modes.However, in order to verify which modes persist at fast rotations, we must address the issues with the relativistic problem.Only then can we construct realistic neutron-star models using nuclear-matter equations of state.
This work has some natural extensions.Our analysis was limited to examine perturbations that have the r-mode ordering in a slow-rotation expansion.Clearly, the divergences we have witnessed spoil this ordering and the equations we use are not strictly valid when they arise.It would therefore be interesting to examine the behaviour of the divergent r-modes using machinery appropriate for the inertial modes.It is notable that the vast majority of the work on inertial modes has focused on barotropic stars.Moreover, efforts in the non-barotropic context tend to assume Γ 1 = const (Lee & Saio 1987;Yoshida & Lee 2000a), with Villain et al. (2005) and Kraav et al. (2022b) as exceptions.Further work should consider whether realistic stratification has important consequences for the r-modes of relativistic stars.

Figure 1 .
Figure1.The relative difference between the adiabatic indices of the equilibrium Γ and the perturbations Γ 1 against the baryon-number density  b for two nuclear-matter equations of state, BSk19 and BSk21.These models show that, although neutron stars are predominantly non-barotropic, they are expected to have barotropic layers.In particular, both equations of state are barotropic towards the surface and BSk19 is also close to barotropic near  b = 0.2 fm −3 , 1 fm −3 .SeeAndersson & Gittins (2022) for a description of how Γ 1 is defined.