Formulating the r-mode Problem for Slowly Rotating Neutron Stars

We revisit the problem of inertial r-modes in stratified stars, drawing on a more precise description of the composition stratification in a mature neutron star. The results highlight issues with the traditional approach to the problem, leading us to rethink the computational strategy for the r-modes of nonbarotropic neutron stars. We outline two strategies for dealing with the problem. For moderate to slowly rotating neutron stars the only viable alternative may be to approach the problem numerically from the outset, while a meaningful slow-rotation calculation can be carried out for the fastest known spinning stars (which may be close to being driven unstable by the emission of gravitational waves). We demonstrate that the latter approach leads to a problem close, but not identical, to that for barotropic inertial modes. We also suggest that these reformulations of the problem likely resolve the long-standing problem of singular behavior associated with a corotation point in rotating relativistic neutron stars. This issue needs to be resolved in order to guide future gravitational-wave searches.

The inertial r-modes of a spinning neutron star have attracted a fair amount of attention.Theorists have explored the precise nature of the r-modes and how they depend on the complex physics of a neutron star interior, while observers have tried to establish the presence of the predicted r-mode signature in observational data.Much of this interest stems from the discovery (now a quarter of a century ago!) that the r-modes may be driven unstable by the emission of gravitational waves (Andersson 1998;Friedman & Morsink 1998;Lindblom et al. 1998;Andersson et al. 1999a).The r-mode instability may limit the spin-up of accreting neutron stars in low-mass x-ray binaries (Andersson et al. 1999b), providing a natural explanation for the apparent absence of sub-millisecond radio pulsars.The mechanism may also lead to the emission of detectable gravitational waves from newly born neutron stars (Owen et al. 1998), mainly through the current multipoles associated with the induced fluid motion.This has motivated a sequence of observational gravitational-wave papers (Abadie et al. 2010;Aasi et al. 2015;Fesik & Papa 2020;Abbott et al. 2021aAbbott et al. ,b, 2022a,b;,b;Covas et al. 2022), so far mainly setting upper limits on the attainable r-mode amplitude.There have also been tantalizing hints of r-mode oscillations in the the x-ray emission from two fast spinning, accreting neutron stars (Strohmayer & Mahmoodifar 2014a,b), see also (Lee 2014;Andersson et al. 2014), but these results are far from conclusive.
The need to understand the impact of different aspects of neutron star physics-ranging from the main dissipation channels (shear and bulk viscosity), the state of matter (superfluid mutual friction and the interfaces with the elastic crust) and the role of the star's magnetic field-led to a flurry of activity following the original instability discovery.This work is summarised in early review articles, see Andersson & Kokkotas (2001) and Andersson (2003), with more recent additions like Lasky (2015) and Haskell & Schwenzer (2021) providing a mature perspective.While many aspects of the problem are fairly well understood some vexing issues remain.Arguably, the most important of these issues relates to the r-modes in relativity.
In order to make robust predictions for (say) the r-mode frequency in a neutron star we need to involve a realistic matter equation of state.This, in turn, requires a general relativistic mode calculation.This problem has not-in our view-yet been solved in a satisfactory fashion.Let us explain.The most important contributions to the discussion, from the initial relativistic inertial-mode calculations by Lockitch et al. (2000Lockitch et al. ( , 2003) ) and Ruoff et al. (2003) through to the more recent work for real equations of state by Idrisy et al. (2015), assumes that the matter is barotropic.However, this is not expected to be a realistic assumption.Instead, as established by Reisenegger & Goldreich (1992), the stratification associated with internal composition gradients makes the problem non-barotropic.This seems to complicate the r-mode calculation.In fact, the problem appears to become singular leading to a continuous spectrum (Kojima 1998;Kojima & Hosonuma 1999;Beyer & Kokkotas 1999).The implications of this are not well understood, but it seems reasonable to argue that the continuous spectrum arises because of simplifying assumptions introduced in the analysis.Adopting this view, the question becomes how we can regularise the problem.While different strategies have been proposed, see Lockitch et al. (2004); Yoshida & Lee (2002); Pons et al. (2005) and the recent effort from Kraav et al. (2021Kraav et al. ( , 2022)), it is probably fair to suggest that the issue has not yet been resolved-at least not completely.This motivates us to return to the problem.
Our aim is to formulate the r-mode problem for neutron stars stratified by composition gradients from first principles.This forces us to consider how the non-barotropic aspects arise and how this affects the fluid perturbation equations.
The key point is that the matter composition may be considered frozen provided the relevant nuclear reactions are slow compared to the dynamics.The argument was already outlined by Andersson & Pnigouras (2019), although in that case the main focus was on the composition g-modes.Here we take one step further by framing the discussion in the context of a realistic matter model.This leads to an important-in hindsight probably obvious-insight.While realistic neutron star models are likely to have non-barotropic high-density cores they will always have barotropic-basically because the matter is composed of single nuclei at lower densities-outer layers.This has two immediate repercussions.First, the standard assumption of a constant adiabatic index (Γ 1 ) for the perturbations is never appropriate.Second, none of the existing r-mode calculations actually solve the problem we should be considering.We have have to rethink how we approach the problem.
Ultimately, the implications of our new perspective on the problem will depend on the extent to which numerical mode results differ from existing ones.This problem will not be solved here.In this first paper we are mainly interested in the qualitative aspects and so we introduce a number of approximations that sacrifice accuracy for clarity.We want to make the key points as transparent as possible in order to motivate renewed effort in several directions.First and foremost, we need relativistic mode calculations for true equations of state in order to inform future gravitationalwave searches.It may well be that the current models are "good enough" but-as we are arguing for changes to the formulation of the problem-this is by no means guaranteed.Second, we need to revisit the (technically challenging) problem of nonlinear mode-coupling and saturation of the r-mode instability (Schenk et al. 2001;Arras et al. 2003;Brink et al. 2004aBrink et al. ,b, 2005;;Bondarescu et al. 2007;Bondarescu & Wasserman 2013;Bondarescu et al. 2009).Available mode-coupling results relate to barotropic models and hence may change when we add realism to the discussion.Third, the precise nature of the r-modes impact on the gravito-magnetic interaction and the dynamical tides in spinning neutron star binaries (Flanagan & Racine 2007).Again, this is a problem that has so far only been explored for barotropic stellar models (Xu & Lai 2017;Poisson 2020;Poisson & Buisson 2020;Ma et al. 2021).
It appears that there is quite a lot of work to get on with, so let us get started.

LOW-FREQUENCY OSCILLATIONS OF A ROTATING STAR
In general, the oscillation modes of a rotating star belongs to one of two categories: modes that exist (have a finite eigenfrequency) already in a non-rotating star and modes that owe their existence to the rotation (the Coriolis force).Our main interest here is in the latter class, collectively referred to as inertial modes, but the story will not be complete unless we also touch upon the former.The argument at the centre of our discussion also impacts on the gravity g-modes, which are present if the star is stably stratified, either in terms of an entropy gradient or a varying composition.This is a non-trivial issue.Work on main-sequence stars demonstrates that the impact of stratification may vary, with global oscillation modes having different character in different parts of a star (Lee & Saio 1987).An illustrative example concerns slowly pulsating B stars (Lee 2006).
Work on the oscillations of rotating stars is complicated by the fact that the Coriolis force breaks the spherical symmetry of the non-rotating case, leading to a coupling of the angular harmonics traditionally used to represent the modes.This coupling is particularly significant for modes with frequency comparable to or smaller than the star's rotation frequency.Higher frequency modes may be approached "perturbatively", adding rotational corrections to the modes of a non-rotating star, but the low-frequency problem is intricate.This issue is not at all new.It was recognised already in the seminal work on rotating ellipoids by Bryan (1889) (for a modern version of the calculation, see Lindblom & Ipser (1999)).It is also well known from work on waves in shallow ocean/atmospheres, mainly focused on weather and climate studies.The r-modes-the main focus of our attention-are in fact analogous to the Rossby waves from the shallow water problem (see, for example, Longuet-Higgins (1968); Zaqarashvili et al. (2021)).
For rotating stars, the r-modes were first introduced by Papaloizou & Pringle (1978).Their work was followed by a number of detailed studies in the early 1980s (Provost et al. 1981;Saio 1982Saio , 1981;;Smeyers & Martens 1983).The nature of the problem was laid out in detail, relaxing the Cowling approximation (i.e.allowing for perturbations of the gravitational potential), by Smeyers et al. (1981).This body of work establishes the main features of the r-modes.They are represented by retrograde waves in the co-rotating frame of the star, but appear prograde in the inertial frame.It is precisely this character that makes them unstable to the emission of gravitational waves (Andersson 1998;Friedman & Morsink 1998).
The problem is complicated by the fact that the modes are degenerate in barotropic stars, belonging to the broader class of inertial modes (with frequency proportional to the star's angular frequency) (Lockitch & Friedman 1999).This changes when the stellar fluid is stably stratified.Stratification couples the radial layers in the star, breaking the degeneracy of the modes and allowing for a richer spectrum of r-modes (including radial overtones) (Provost et al. 1981;Saio 1982).In the context of the r-mode instability, the work by Yoshida & Lee (2000a) is particularly notable although they consider either models with stable stratification throughout the interior or models that are fully convective.However-and this is important-the arguments they put forward for the stratification relate to entropy gradients, which are only expected to be relevant for newly born neutron stars.
We will set up the problem in such a way that our discussion connects with the inertial-mode analysis of Lockitch & Friedman (1999).This makes sense because one of the questions we want to address involves how the non-barotropic r-modes morph into barotropic inertial modes as the stratification weakens.
Starting from the velocity perturbations (δv i ), we consider the perturbed Euler equation in the rotating frame of the star.We then have (in a coordinate basis, making due distinction between co-and contra-variant components, and using δ to indicate Eulerian variations) where p is the pressure, ρ is the mass density and Φ the gravitaitonal potential-along with the continuity equation and the Poisson equation for the perturbed gravitational potential In the static limit (where time variations vanish and the rotation frequency Ω i → 0) the equations decouple into two sets.First we have along with (3) and, secondly, ∇ i (ρδv i ) = 0 . (5) At this point, Lockitch & Friedman (1999) note that the first two equations represent perturbations that take us to a neighbouring equilibrium star.The argument for this is straightforward-a static perturbation of the equation for hydrostatic equilibrium takes us to a new configuration with pressure p = p + δp.In a slowly rotating star, this solution would pick up rotational corrections at order Ω, which suggests a solution such that [δρ, δp, δΦ] = O(1) and δv i = O(Ω) . ( This problem is equivalent to considering the dynamics of the original configuration (albeit for a slight different central density).The dynamical aspects of the problem are contained in the second set of perturbations, represented by equation ( 2) for which we would have This is the assumption that leads to the inertial modes.Alternatively, as we can always multiply the linearised equations by a constant, we may normalise the Lagrangian displacement associated with the perturbation, in the rotation frame simply given by such that This is the convention we assume in the following.It is important to appreciate that, regardless of the chosen normalisation, we cannot (completely) determine the density perturbations etc without accounting for the change in shape due to the centrifugal force (which also enters at order Ω 2 ).This obviously complicates the analysis and it would be natural to turn to numerics.However, in order to understand the nature of the problem we wish to proceed analytically (as far as we can).We consider the numerical problem in a companion paper (Gittins & Andersson 2022).
In order to make progress analytically, we will study the perturbations to second order in Ω/Ω 0 , where However, in the interests of clarity, we will keep the background spherical (Lee & Saio 1987).This is a useful simplification as it removes some of the rotational multipole couplings from the problem.It also makes "sense" since the centrifugal force does not introduce additional oscillation modes.In parts of the discussion we will also, again for clarity, make use of the Cowling approximation (neglect the perturbed gravitational potential, δΦ).Both assumptions are relaxed in the companion numerical work.

The frozen composition argument
It is easy to see how the assumption of non-barotropic perturbations upsets the inertial-mode logic.The usual argument introduces the adiabatic index Γ 1 in such a way that (with ∆ representing Lagrangian variations) This immediately leads to a conflict with the assumed ordering for the rotating-star perturbations (Lockitch & Friedman 1999).Since the background is spherical we must either have (introducing both the adiabatic index Γ and the speed of sound c 2 s for the background configuration) which would represent a barotropic model, or the ordering of the solution must change in such a way that ξ r = O(Ω 2 ) to balance the density and pressure perturbations.This then leads to the different starting assumption for non-barotropic models suggested by Lockitch & Friedman (1999) (and eventually brings us to the vexing issue for relativistic r-modes, see (Lockitch et al. 2004)).However, for realistic neutron star physics, the argument turns out to be a bit more subtle.
In order to explain the issue we need to explore the physics that give rise to the stratification in a mature neutron star (internal composition gradients) in the first place, leading to Γ 1 = Γ.The argument draws heavily on the discussion of g-modes and reactions from Andersson & Pnigouras (2019).Essentially, once we account for out-of equilibrium nuclear reactions1 , the perturbed proton fraction x p evolves according to (for small deviations from equilibrium-i.e. in the so-called sub-thermal limit (Alford & Harris 2018)) Here ∆β represents the deviation from beta-equilibrium and γ encodes the (dominant) reaction rate.Considering β a function of ρ and the proton fraction x p , and assuming that the star is non-rotating (so that v i = 0, which is also true in the rotating frame), we have which, once we consider (13), becomes That is, we have with The coefficients A and B are time independent, as they are evaluated for the equilibrium background, so if we work in the frequency domain (as we typically do when we consider stellar oscillations) then we have a harmonic time dependence e iωt and it follows that This relation is commonly taken as the starting point for discussions of bulk viscosity, see Schmitt & Shternin (2018) for a recent review of this issue.Here we want to make a slightly different emphasis.
Let us consider the timescales involved.Noting that A needs to be negative in order for the system to relax towards equilibrium, we introduce a characteristic reaction time as Then we see that, if the reactions are fast compared to the dynamics (associated with a timescale ∼ 1/ω) then |t R ω| 1 and we have ∆β ≈ 0. In effect, the fluid elements reach equilibrium before executing an oscillation.The fluid remains in beta-equilibrium and hence the perturbations are (effectively) barotropic.
As a ballpark estimate of the relevant timescale, we draw on the discussion by Haensel et al. (2002) and assume for the modified Urca reactions (ignoring density dependence as we only need a rough idea here).This estimate suggests that-for all modes/rotating rates we may conceivably be interested in-mature neutron star matter will not be in the fast-reaction regime.We have to consider the slow-reaction (stratified) problem.
In the limit of slow reactions we have |t R ω| 1 and we can Taylor expand (18) to get we have Comparing to (11), we have an expression for Γ 1 in terms of thermodynamical derivatives: where the last equality-demonstrated by Andersson & Pnigouras (2019)-holds as this limit represents frozen composition (∆x p = 0).In this case we need to consider the impact of composition stratification on the fluid dynamics.Moving on, introducing the gravitational acceleration (for a spherical star) we see that (11) leads to It is important to note that the composition of matter impacts on both terms on the right-hand side of this relation.
In the following, when we consider the impact of composition stratification on the oscillations of a slowly rotating star, it is convenient to consider the Brunt-Väisälä frequency, given by This relation illustrates some of the subtleties we need to consider.For example, if Γ and Γ 1 are taken to be constant (as is commonly assumed) then N 2 must diverge as we approach the surface of the star, where p/ρ → 0. The only way to avoid this problem is if the fluid becomes barotropic in the low-density limit, as we then have Γ 1 → Γ as ρ → 0. Also, it is important to keep in mind that the assumption that Γ 1 is constant is different from holding N 2 fixed.However, as we will soon see, neither assumption is realistic.
The equation of state relation allows us to remove the density perturbation from the discussion.We then have Evidently, and quite intuitively, stratification does not affect incompressible flows, for which ∆ρ = 0 =⇒ ∆p = 0.For later convenience, if we introduce the density scale height we may remove one of the three background quantities from the discussion.

The impact of stratification
As already stated, we will work with the displacement rather than the velocity perturbations in the following.We may then write the perturbed continuity equation as (ignoring the rotational deformation of the background star, as advertised, and making use of the equation of state relation (28)) Clearly, it makes sense to consider the dimensionless quantity An example of this quantity, for the BSk19 and BSk21 equations of state (Fantina et al. 2013;Potekhin et al. 2013) (two models with very different proton fraction profiles), is provided in Figure 1.The results show that ˆ 2 typically varies by about an order of magnitude throughout the star's core, reaches a peak in the low-density region and then drops sharply towards the surface.This behaviour will guide the discussion in the following.
Next consider the Euler equations (in a rotating frame), which take the form  2013)).The stellar models are obtained by integrating the relativistic equations of stellar structure with appropriate generalisations for the Brunt-Väisälä frequency N and gravitational acceleration g.The horizontal lines indicate the value of2 = Ω 2 /Ω 2 0 for the two spin frequencies 100 Hz and 716 Hz, the latter representing the fastest known spinning neutron star (Hessels et al. 2006).
For an axisymmetric rotating star, the ϕ-component becomes (noting that all perturbations behave as e imϕ ) In essence, given that the inertial modes we are interested in have frequency ω ∼ Ω, we must have δp ∼ O(Ω 2 ) in order to have ξ i ∼ O(1).The radial component of the Euler equation then tells us that we must also have δρ ∼ O(Ω 2 ) and the equation of state relation (28) then leads to For a given stratification, expressed in terms of N 2 , this constrains the slow-rotation ordering of the radial displacement.This accords with the point we made earlier (Lockitch & Friedman 1999).It is, however, useful to make this argument more precise.For inertial modes of uniformly rotating stars it is natural to use the slow-rotation expansion with = Ω/Ω 0 as the small parameter, noting that the Kepler break-up frequency corresponds to roughly Given the previous discussion we know that we should have δp = 2 δ p using the tilde to indicate an O(1) quantity (notation we will adopt in the following), and δρ = 2 δ ρ, as well.If, in addition, we focus on inertial modes we know that ω ∼ Ω and it is evident from (33) that we do not have to consider terms of order ; we can take 2 as the small rotation parameter 2 .Hence, we expand the mode frequency as The equation of state relation (28) then becomes Let us see what we learn from this.Noting that the Kepler frequency corresponds to 1 while the results from Figure 1 show that ˆ 2 < 1 as well, we may consider a double expansion in and ˆ .Of course, since ˆ varies throughout the star this would have to be a local argument.The first few terms of such an expansion would be and it is easy to see that we need to consider several different cases.First, suppose ˆ as would be the case if we combine the results in Figure 1 with a fairly slowly spinning star.Effectively, this situation corresponds to taking ˆ ∼ O( 0 ) and-formally letting ξ r 0 + ˆ 2 ξr 2 =⇒ ξ r 0 -it follows that ξ r 0 = 0 , at order 0 , (40) This is the usual stratified r-mode ordering (see, for example, Provost et al. (1981)), but it is not clear that this is the case we should consider.The argument would only apply to slow and perhaps moderately fast spinning neutron stars (say, the 100 Hz case illustrated in Figure 1), but not the fastest observed systems.Instead, the results in Figure 1 suggest that, for stars spinning with ∼ 0.3 (about half the Kepler rate), roughly corresponding to the fastest known pulsars with spin frequency close to 700 Hz, we should instead consider ˆ ∼ O( ), in which case we may ignore terms of order 2 ˆ 2 in the expansion.This then leads to This is different from what we usually assume, yet seems a case we ought to consider.Perhaps significantly, the relation suggests that we may have ξ r 0 = 0, which changes the mode structure.The model is also interesting as it limits to the barotropic case in regions where ˆ 2 → 0.
Finally, from Figure 1 it is clear that there will always be always a low-density region where ˆ O( 2 ).If the ˆ 2 are smaller than (say) 4 then we may ignore the stratification and consider the barotropic result If these assumptions hold throughout the star, then we must end up with the inertial modes from Lockitch & Friedman (1999).Globally, this is unlikely to be the relevant case but Figure 1 suggests that we always have to to consider the region close to the surface as barotropic.This impacts on the surface boundary condition and may, in turn, also affect the modes.
Clearly, the profile for ˆ is fixed for any given stellar model, while can be varied (up to the Kepler limit, which means that we should have 0.6).From Figure 1 it is easy to see that all the suggested orderings may apply locally in a neutron star core, making the formulation of consistent model tricky.The main conclusion is that we have to consider all results obtained with the standard "constant-Γ 1 " prescription as unrealistic.In fact, if we take the results in Figure 1 at face value-and there is no reason why we should not-then we have to reconsider our strategy for the fastest spinning stratified neutron stars.For the core of these stars, relation (42) should apply, in which case the perturbation problem is closer to-but not exactly the same as-that for inertial modes (Lockitch & Friedman 1999).To what extent this affects the mode frequencies remains to be established.

FORMULATING THE MODE PROBLEM
A mode solution to the perturbation problem must satisfy the perturbation equations (obviously) and relevant boundary conditions (typically, regularity at the centre of the star and the vanishing of the Lagrangian variation of the pressure at the star's surface).It is well established that the equations allow for oscillation modes of different character.Moreover, the problem gets richer as more detailed neutron star physics is considered.Somewhat simplistically, each aspect of physics added to the model-matter composition, rotation, superfluidity, electromagnetism, elasticity...brings new families of modes into play.This makes the general problem complex.
In order to build useful intuition, it is natural to focus on particular aspects.If we are mainly interested in the oscillations of a rotating star, the natural starting point would be to work out how the rotation impacts on modes that exist already in a non-rotating star.For a mode with frequency ω n in the non-rotating star, the strategy isat least in principle-fairly straightforward, although the mode calculation may get quite involved as the rotation couples different multipole contributions.Still, for slow to moderate rotation rates, this problem can be dealt with perturbatively as long as Ω ω n .This should be the case for the fundamental mode of the star, which has frequency of order the Kepler breakup frequency, and the (even higher frequency) pressure p-modes.We will not consider those problems here.We will also not consider the rotational corrections to the gravity g-modes, a slightly more subtle issue given that the high-overtone g-modes are expected to have very low frequencies in a non-rotating star.Hence, for these modes the Coriolis force may dominate over the buoyancy already at fairly low rotation rates and they would then become part of the problem we are considering.An example of this behaviour can be found in the study by Yoshida & Lee (2000b), where it is shown how inertial modes are strongly modified when the buoyancy force becomes comparable to, or stronger than the Coriolis force.
In this exploratory analysis our main focus is on the qualitative nature of the low-frequency modes of a rotating star.In this spirit, we rely on simpliying assumptions.In particular, following Lee & Saio (1997), we ignore the change in shape of the background star associated with the centrifugal force.This assumption is not expected to affect the qualitative nature of the problem.
As we have already seen, the discussion necessarily gets somewhat involved and some of the issues we need to consider are subtle.In general, we need to consider both the impact of rotation on oscillation modes that exist already in a non-rotating star and modes that are brought into existence when we consider the impact of the Coriolis force.Given this, it makes sense to work out the general perturbation equations to first slow-rotation order.This involves making choices already at the outset.
There are three common strategies for investigating the oscillations of rotating stars (Unno et al. 1989).The first, and formally most elegant, approach expresses the rotational corrections to a given mode as a sum over all the modes of the corresponding non-rotating star (which form a suitable complete basis as long as we ignore dissipation).The second option builds on an explicit expansion in angular harmonics while the third involves time evolving the perturbation equations.The last strategy has the advantage that one can readily deal with fast spinning stars, for which the algebra of the other approaches becomes daunting, but it also has the drawback that one loses track of the fine details of the problem (you get what you get from the simulation, depending on the chosen initial data).Examples of work in this direction can be found in Jones et al. (2002); Passamonti et al. (2009) for Newtonian models and Gaertig & Kokkotas (2009); Gaertig et al. (2011); Gaertig & Kokkotas (2011); Krüger et al. (2021) for efforts in relativity.In the following we will carry out an expansion in harmonics.This approach has the advantage that it highlights the nature of the fluid motion.
Assuming that the oscillation modes-with label n and frequency ω n -are associated with a displacement vector (where the hat indicates a quantity that is independent of t and we adopt the Lockitch & Friedman (1999) sign convention) we have (in a coordinate basis with e i r = ∂r/∂r etc) where we refer to W l and V l as polar perturbations, while U l is axial (and noting that the m-multipoles decouple for an axisymmetric system, like a rotating star).For a given multipole l, these perturbations have different parity.This follows since the equilibrium state of a rotating star is invariant under the parity transformation (defined by a reflection through the origin, θ → π − θ and ϕ → ϕ + π), the linear perturbations have definite parity for this transformation.Alternatively, the different modes are sometimes described as even and odd, see for example, Lee & Saio (1987).
Along with the decomposition of the displacement, all scalar perturbations are expanded in spherical harmonics.That is, we have with (dropping the hats on the individual l-multipole components to keep the equations that follow as tidy as possible) and similar for all other scalar quantities.We also know that the rotating equilibrium remains spherical to linear order in Ω so all background quantities depend only on r as long as we consider the first order slow-rotation corrections.
Working to this order of approximation, let us summarise the equations we need.

The perturbation equations
First, it follows from (45) that This result is important because it shows that-up to order Ω-only the polar contributions to a given mode [W l , V l ] contribute to the density perturbation δρ l .This is brought out by the continuity equation, which leads to (changing l → j for consistency with the recurrence relations to be derived in the following) Notably-as long as we ignore the rotational deformation-this equation does not involve the coupling of different multipoles.
Turning to the perturbed Euler equations, in the rotating frame we have leading to the radial component the θ component and the ϕ component One possible strategy would be to work with these equations as they are, deal with the fact that different multipoles couple head on and solve the problem numerically (see, for example, Lee (2006)).However, this may not be the most "transparent" option as it obscures the nature of the different mode solutions.As we will see later, it follows from the angular components ( 52) and ( 53) that we have to consider coupling between the l-components and the ones for l ± 2. This leads to a larger set of equations to solve so it makes to ask if this coupling can be avoided.It turns out that it cannot, but we can find a somewhat more "intuitive" set of equations to solve.
From the Euler equations it is easy to see that it makes sense to introduce (not to be confused with the axial amplitude U l ).This variable is used both in the classic (dimensionless) formulation from Unno et al. (1989) and the two-potential formalism used by, for example, Lindblom & Ipser (1999).
Making use of the standard recurrence relation where the radial Euler equations leads to Multiplying this by Y m j , integrating over the angles and executing the sum over l we have the recurrence relation Next we consider the combination (the radial component of the vorticity equation) where we have made use of Legendre's equation to simplify the result.
Using the previous recurrence relation ( 55), along with we arrive at the recurrence relation Keeping Legendre's equation in mind, it may also be useful to consider the combination (Lee & Strohmayer 1996;Glampedakis & Andersson 2006) which leads to For inertial modes, it is notable that the radial vorticity equation ( 62), links only variables that have a leading order contribution.Another such relation follows from the θ component of the vorticity equation This equation is identical to equation ( 39) from Lockitch and Friedman, apart from the right-hand side which only vanishes for barotropic stars.In general, the equation is a bit messy as it couples radial derivatives of all displacement components and it also involves the l ± 2 multipoles.However, we will find the equation useful for non-barotropic stars satisfying the traditional r-mode slow-rotation ordering as it simplifies considerably in that case.Finally, we also need the perturbed Poisson equation for the gravitational potential.However, we will make the Cowling approximation in our explicit examples (set δΦ j = 0) so will not give the equation here.

THE "TRADITIONAL" R-MODES
As already advertised, we will focus on low-frequency modes such that Ω ω 0 .This includes high-order gravity g-modes and inertial modes.We have already seen that the perturbed Euler equations then imply that we must have [δρ l , δp l , δΦ l ] ∼ O(Ω 2 ).Moreover, the continuity equation requires the polar components W l and V l to be of the same order, and we already know that if we consider a strongly stratified star (with ˆ ) then we must have [W l , V l ] ∼ O(Ω 2 ).We are then left to consider if it is possible to combine these assumptions with U l ∼ O( 1).The answer to this question is affirmative, but it follows from (62) that we must then have That is, for given values of l = l (say) and m, we may have U l = 0 as long as the leading order mode frequency is given by These are the r-modes (Papaloizou & Pringle 1978;Provost et al. 1981;Saio 1982).In addition to the leading order displacement they will have polar components as well as other axial multipoles U l =l , but these enter at O(Ω 2 ).We will deliberate on these contributions in the following.With our conventions, the pattern speed of a mode is −ω/m.Therefore, all r-modes travel in the same direction across the star (retrograde with respect spin, in the rotating frame).We also note that there are no axisymmetric r-modes; we must have m = 0.
Focusing on the r-mode problem, we consider the perturbations for mode solutions such that ω n ∼ O(Ω) and As before, we use tildes to identify terms that enter at order Ω 2 , i.e.
with δ pl ∼ O(1) by definition.In addition, we need to keep track of the rotational correction to the frequency, so recall (37) from which it is worth noting that assumed slow-rotation ordering is only valid as long as With these assumptions, the continuity equation ( 49) relates order Ω 2 quantities, and we have So far, the different relations only involve single-multipole polar components.This changes when we turn to the perturbed Euler equations.
Let us first consider the radial component of the vorticity equation ( 62).We know already that, at leading order we may have a single axial contribution U l = 0 as long as the frequency is given by (67).However, at this point, we cannot determine the axial eigenfunction.Essentially, the leading order r-mode solutions are degenerate.To break this degeneracy, we need to go to higher orders, keeping in mind that the polar contributions enter at order Ω 2 .Equation ( 62) provides a recurrence relation involving these multipole contributions.At order Ω 3 we have (with This is the only relation we get that involves the leading order eigenfunction and the frequency correction ω2 .However, for j = l + 2 we get from ( 62) From these relations, the pattern is clear.At order Ω 2 , the r-mode solution involves a number of multipoles.The question then becomes, does this sequence truncate?To answer this question, first note that Q m = 0, which helps establish the lowest order term in the series.We have three options.First, we may have l = m in (72).In this case, the (leading order) U l =m term corresponds to the lowest order multipole in the solution.In the language of Lockitch & Friedman (1999), the mode is axial-led.This case corresponds to the traditional l = m r-mode (Papaloizou & Pringle 1978).Another option would be to have l = m + 2 in (74).This would also lead to an axial-led mode, but now the lowest multipole is given by Ũl −2 ∼ O(Ω 2 ).A third option follows by setting l = m − 1 in (74) which then decouples and from ( 72) we arrive at a polar-led mode with the lowest multipole contributions given by [ Wl −1 , Ṽl −1 ] (again at order Ω 2 ).The main lesson here is that the nature of the r-modes is quite similar to that of the general inertial modes discussed by Lockitch & Friedman (1999).Each mode has several multipole contributions, but the U l term is elevated above the other contributions in the slow-rotation expansion.The close relation between the two problems may not been very clearly explained in the existing literature.It is, however, important for what follows.
In order to complete the formulation of the problem, we will use the other vorticity equation ( 65).With the ordering we have, at order Ω 2 , this reduces to From this relation, we infer two relations involving the leading order U l term.First, with j = l + 1 we have Second, with j = l − 1 we get Finally, we have the divergence equation ( 64), which leads to For j = l + 1 we have In essence, if we want to determine the leading order eigenfunction and the frequency correction ω2 , we need to solve a coupled system for U l and Wl ±1 .The other contributions to the mode-solution (like Ũl ±2 ) can be calculated as a second step.
Finally, the mode-solution must satisfy the condition that the Lagrangian perturbation in the pressure vanishes at the surface.That is, we require At this point, we may return to the question of whether the multipole sum truncates for the r-modes.First, the relation ( 75) also tells us, for j = l + 3 and j = l − 3, respectively, that we must have (as long as N 2 = 0) Second, in ( 79) and ( 80), we can use j = l ± 2 to show that we must have Third, with the assumed r-mode ordering, the radial component of the Euler equations (58) leads to so, for j = l + 1 and making use of (79) we have It also follows that Similarly, with j = l − 1 we get and we also have Combining the results, we see that we must have Wl ±3 = δ pl ±3 = δ Φl ±3 = 0. Finally, the continuity equation leads to Ṽl ±3 = 0, while δ ρl ±3 = 0 follows from the equation of state relation.Is essence, all polar l ± 3 multipole contributions will vanish.In turn, this means that a general r-mode must truncate with the Ũl ±2 terms obtained from ( 73) and ( 74).This accords with the discussion in Smeyers et al. (1981) and Smeyers & Martens (1983).

The l = m modes
Having written down the equations we need to solve to determine the frequency correction ω2 and the multipole structure of an r-mode to order Ω 2 -notably without any simplifying assumptions other than neglecting the rotational change in shape of the star-we have a decision to make.Do we want to consider a model that is as "realistic" as possible-which will require a numerical solution-or are we more focused on the formal structure of the mode solution?The initial answer is quite simple.As we are not including the rotational shape change it is natural to focus on the qualitative nature of the solution.This leads us to the question of which further simplifying assumptions we may consider.
As already advertised, we are now going to make the Cowling approximation.That is, we assume that δ Φl = 0.For the problem at hand, this is pragmatic (as we are focusing on qualitative aspects) and reasonable (as we do not have to solve the Poisson equation for the perturbed gravitational potential).We want to keep the problem tractable enough that we may proceed to solve it by analytic means.We also know from available numerical results that the r-modes are determined with reasonable precision within this approximation (at least in the context of Newtonian gravity).In the Cowling approximation, we have In the l = m case, we have Q l =m = 0 which means that we only need to consider the coupling between the leading order U l and the polar l + 1 contributions.(The axial second order contribution Ũl +2 can be calculated at a second stage.)The set of equations to consider now are: (i) the continuity equation ( 71) (ii) the differential equation ( 76) along with (iii) the algebraic relation (72 (iv) the relation ( 79) and the surface boundary condition (81) (obviously).It is easy to see that we end up with two coupled first-order equations for U l and Wl +1 .
It is instructive to introduce U l = r l +1 Ūl and rewrite (91) as It follows immediately that, in the barotropic limit (when N 2 → 0 for a fixed Ω), we must have The only alternative would be for to remain finite in the barotropic limit.However, this would violate the assumed slow-rotation ordering for the r-mode solution.That this happens should not be a surprise given the general discussion in Section 2. The result is simply an illustration of the fact that we need to make different assumption in barotropic regions of the star.
In general, we need to solve (94) along with the continuity equation ( 90), which becomes The two equations (plus the boundary conditions) constitute a Sturm-Liouville problem (Provost et al. 1981) so we expect to have an infinite set of eigenvalues (see Saio (1982) and Gittins & Andersson (2022) for indicative results for the r-mode overtones).However, as we have seen, the problem changes in barotropic limit.The overtones disappear and we are left with a single r-mode, represented by (95).
As a simple example of the single remaining r-mode, we may consider an incompressible barotropic star, for which ρ = constant so c 2 s → ∞ and we are left with with Ūl = A = constant.This leads to and we also have Finally, the surface boundary condition becomes so we arrive at ω2 = − 8m (l + 1) 4 , when l = m.We briefly compare this result to available results from the literature in Appendix A. In summary, our arguments clearly illustrate the known fact that the nature of the l = m r-mode problem changes as stratification weakens.The evidence is clear.We have to approach the N 2 → 0 limit with care.In fact, for neutron stars the problem is particularly intricate.Taking the results in Figure 1 at face value, we always have to assume the region close to the surface of the star to be barotropic, while the high-density region may not be.This, in turn, means that the assumed slow-rotation ordering associated with the non-barotropic r-mode must break (as ˆ close to the star's surface).In effect, the formulation of the problem-as we have presented it-is not consistent.This presents a technical challenge as the solution needs to smoothly join the stratified region where the stratified assumptions hold with a barotropic region where the equation becomes those associated with the general inertial modes.As far as we are aware, this problem has not been considered (at least not for neutron stars), although the required strategy-basically abandoning the slow-rotation ordering for the perturbation-has been developed and employed to good effect for main sequence stars (Lee 2006).

The l = m modes
Let us now turn to the l = m r-modes.In general, we then need to consider the continuity equations for the polar l ± 1 contributions.In particular, we need to solve the differential equations ( 76) and (77): and It is easy to see that we run into trouble in the barotropic limit.Combining ( 104) and ( 105) we have an algebraic relation: or, alternatively, The first relation shows that, if Wl ±1 remain finite then we must have U l → 0 in the barotropic limit.This would be incompatible with (104).We get a hint of the resolution to the problem from the alternative version, which suggests that if we insist that U l remains O(1) when N 2 → 0 then Wl −1 must diverge.
The unavoidable conclusion is that the l = m r-modes cannot exist in the barotropic limit-as expected from the arguments by Lockitch & Friedman (1999).If N 2 → 0 for a fixed rotation rate Ω, then the assumed r-mode ordering must break.In fact, if N 2 = 0 at some point in the star we have a problem.Given the available equations there does not seem to be a way to avoid dividing by N 2 so the problem will be singular.
In summary, while the stratified problem can be solved for l = m r-modes (Saio 1982), the solution does not apply for realistic neutron star models, see Gittins & Andersson (2022) for related numerical results.If we want to consider the actual problem, then we have to rethink our strategy.This again suggests that we may need to abandon the slow-rotation ordering a tackle the general problem numerically (as in the body of work by Lee and collaborators (Lee & Saio 1987;Lee & Baraffe 1995;Lee & Saio 1997;Yoshida & Lee 2000a,b;Lee 2006)).

A PHYSICALLY MOTIVATED ALTERNATIVE
Based on the stratification results from Figure 1, it makes sense-for the fastest spinning stars-to consider the stratification to be second order in the slow rotation expansion.We then have (42) which leads to (at order 2 ) where W j ∼ O(1).This suggests that we change the assumed ordering in such a way that W j → W j + 2 Wj and similar for V j .The axial displacement remains as before.With the polar displacement components present already at leading order, the problem is close to that for a general inertial mode.With these assumptions, the leading order continuity equation requires 1 r 2 ∂ r (rρW j ) − j(j + 1) The radial vorticity equation leads to, at order : It is also convenient to use the algebraic relation from the divergence equation, which at order 2 provides : ω 0 [j(j + 1)ω 0 − 2m] V j − 2mω 0 W j − j(j + 1)δ Ũj − 2ω 0 [(j − 1)(j + 1)Q j U j−1 + j(j + 2)Q j+1 U j+1 ] = 0 .
Finally, in this case it seems natural (given that the horizontal vorticity equation involves derivatives of all three displacement components to use the radial Euler equation, which leads to, at order 2 : − ω 2 0 W j + 2mω 0 V j + r∂ r δ Ũj + 2ω 0 [(j − 1)Q j U j−1 − (j + 2)Q j+1 Finally, in the Cowling approximation we have (89) and we also need ( 65), which at order 2 leads to ω 0 2 1 − Q 2 j − Q 2 j+1 r∂ r W j + mω 0 W j − ω 0 mω 0 + 2 (j + 1)Q 2 j − jQ 2 j+1 r∂ r V j + 2m 2 V j − ω 0 [(ω 0 (j − 1) − 2m) r∂ r U j−1 + 2m(j − 1)U j−1 ] Q j + ω 0 {[ω 0 (j + 2) + 2m] r∂ r U j+1 + 2m(j + 2)U j+1 } Q j+1 Let us focus on the problem we would have to solve in order to identify a solution "close to" the traditional r-mode.That is, we are looking for modes such that U l ∼ O(1) with l = m and with a frequency given by ( 37) and ( 67).With the usual ordering for stratified stars, this would include the r-mode overtones.The only difference here is that we are no longer (necessarily) assuming that the polar displacement contributions enter at higher slow-rotation order.The equations that involve U l are then, first of all, the leading order relation which, for a mode with the usual leading order r-mode frequency (67), reduces to This can be combined with the leading order continuity equation to give ∂ r (ρW l +1 ) + (l + 2)ρW l +1 = 0 , which leads to, with A constant, r l +2 ρW l +1 = A =⇒ W l +1 = A ρr l +2 .
(117) This is clearly problematic as the solution diverges at the centre of the star (and at the surface as well, if ρ → 0 as r → R).The only way to avoid trouble is to have the trivial solution, A = 0, and move on to the equations for the higher order terms.If we do this then we immediately see that the problem is identical to the one we already solved for barotropic stars.There will be a single r-mode for each l = m.This tells us that the frequencies of the r-mode overtones can no longer be given by (37).As expected-and in accordance with the results from Figure 4 of Yoshida & Lee (2000a)-they have to change character.Next, for l = m it is easy to see that the equations we are now considering still lead to a singular problem unless the polar components W j and V j are O(1).This is as expected: We need to consider solutions close to the barotropic inertial modes.
Finally, for the modes of the fastest spinning neutron stars we see that the problem (to leading order) is very close to the inertial-mode problem as formulated by Lockitch & Friedman (1999).The only difference is the right-hand side of (113).In essence, the problem we need to consider if we want to establish the astrophysical role of the gravitationalwave driven r-mode instability is close, but not identical, to the inertial-mode problem.As far as we are aware, this problem has not been stated despite the numerous discussions of the r-mode instability in the literature.This problem clearly needs further attention and our intention is to approach it numerically (also accounting for the rotational shape corrections, following the strategy outlined in Gittins & Andersson (2022)) in the near future.

CONCLUSIONS AND OUTLOOK
We have revisited the problem of inertial r-modes in stratified neutron stars.Our motivation for this was twofold.First, we wanted to add realism to the discussion by introducing a more precise description of the composition stratification in a mature neutron star.Our analysis of the problem highlights issues with the traditional approach to the problem.In order to account for the expected variation of the internal composition stratification with density, we need to rethink the computational strategy for determining the r-modes.There appears to be two strategies for

Figure 1 .
Figure1.The dimensionless quantity ˆ 2 inside two M = 1.4 M neutron stars described by the BSk19 (left panel) and BSk21 (right panel) equations of state fromFantina et al. (2013).The two models are suitably indicative because the dependence of the proton fraction with density is very different in the two cases (see Figure4fromPotekhin et al. (2013)).The stellar models are obtained by integrating the relativistic equations of stellar structure with appropriate generalisations for the Brunt-Väisälä frequency N and gravitational acceleration g.The horizontal lines indicate the value of 2 = Ω 2 /Ω 2 0 for the two spin frequencies 100 Hz and 716 Hz, the latter representing the fastest known spinning neutron star(Hessels et al. 2006).