Rossby modes in slowly rotating stars: depth dependence in distorted polytropes with uniform rotation

Large-scale Rossby waves have recently been discovered from measurements of horizontal surface and near-surface solar flows (L\"optien at al. 2018). We are interested in understanding why only the sectoral modes are seen in the observations and also in modelling the radial structure of the observed modes. To do so, we characterise here the radial eigenfunctions of r modes for slowly-rotating polytropes in uniform rotation. We find that for free-surface boundary conditions on a spheroid of non-vanishing surface density, r modes can only exist for $\ell=m$ spherical harmonics in the inviscid case, and we compute their depth dependence and frequencies to leading order. For quasi-adiabatic stratification the sectoral modes with no radial nodes are the only modes which are almost toroidal and the depth dependence of the corresponding horizontal motion scales as $r^m$. For all r modes except the zero radial order sectoral ones, non-adiabatic stratification plays a crucial role in the radial force balance. The lack of quasi-toroidal solutions when stratification is close to neutral, except for the sectoral modes without nodes in radius, follows from the statement that the system needs to be in both horizontal and radial force balance. In the absence of super- or subadiabatic stratification and viscosity, both the horizontal and radial force balances independently determine the pressure perturbation. The only quasi-toroidal cases in which the two determinations of the pressure perturbation are consistent are the special cases where $\ell=m$, and the horizontal displacement scales with $r^m$.


Introduction
Rossby waves, large-scale waves of radial vorticity with retrograde phase speed, have recently been discovered from measurements of horizontal surface and near-surface solar flows (Löptien et al. 2018, and confirmed by Liang et al. 2019). The clearly observed waves have frequencies near those of sectoral traditional Rossby waves in a uniformly rotating fluid system (e. g. Longuet-Higgins 1964), corresponding to sectoral spherical harmonics of azimuthal order 3 ≤ m ≤ 15. Löptien et al. found that the amplitudes of these Rossby waves do not depend strongly on depth down to 21 Mm below the photosphere, but could not further characterise the radial dependence of the eigenfunctions. Assuming the motion is incompressible, they argued that viscous damping is the reason why they observe only sectoral Rossby modes in the Sun. Here we are interested in understanding, for a more realistic stellar stratification, why only the sectoral modes are seen in the observations and also, in modelling the radial structure of the observed modes.
We restrict our attention to the Rossby waves discussed in Löptien et al. which have a dispersion relation close to that of traditional (non-magnetic) Rossby waves. We do not discuss magnetic Rossby waves (e.g. Zaqarashvili et al. 2010;McIntosh et al. 2017;Dikpati et al. 2018).
The restoring force for traditional Rossby waves is the Coriolis force. Rossby waves have been studied extensively in the geophysical context (see, for example, the textbook by Vallis 2006), with special interest in their horizontal motion and with applications to the Earth's atmosphere and oceans (Rossby 1939;Dickinson 1978), but also the atmospheres of Jupiter and Venus (e.g. Allison 1990;Covey & Schubert 1982;Nara et al. 2019).
In the stellar context, waves analogous to planetary Rossby waves are known as r modes (or quasi-toroidal modes, see e.g. Papaloizou & Pringle 1978;Unno et al. 1989). They have been considered in the photosphere of the Sun (starting with the speculative work of Plaskett 1966), as well as near the base of the convection zone (see the series of papers starting with Gilman 1969) where the stratification is assumed to be subadiabatic. Wolff & Blizard (1986) studied the properties of the r modes in the convective zone of the Sun, but did not predict any restriction on the existence of non-sectoral modes for uniform rotation.
In general, it can be shown that for low-frequency nonradial oscillations of a rotating star, the spheroidal components associated to the spherical harmonics of degree couple with the toroidal components of adjacent degrees ±1. Furthermore, these toroidal components ± 1 themselves couple with the spheroidal components associated with and ± 2. Thus, without significant simplifications, a nonradial oscillation mode in a rotating star is given by an infinite sum of terms proportional to spherical harmonics with different degrees for a given azimuthal index m (Zahn 1966;Berthomieu et al. 1978). In numerical analysis, a truncation of the series is inevitable. Previous works have opted Article number, page 1 of 10 arXiv:2003.05276v1 [astro-ph.SR] 11 Mar 2020 A&A proofs: manuscript no. _arXiv for drastic truncation, retaining only the first two terms (see e.g. Lee & Saio 1987), which may affect the results significantly.
Two approaches have been considered in order to remove this difficulty. One of these is the so-called traditional approximation, in which the horizontal component of the rotation vector is neglected (Lee & Saio 1989, 1997. Then the Coriolis force associated with radial motion and the radial component of the Coriolis force associated with horizontal motion are both neglected. Alternatively, the solution can be sought using asymptotic expansion relative to a small parameter proportional to the rotation frequency (Provost et al. 1981;Smeyers et al. 1981). The former approximation is valid locally in regions of the star where both the rotation frequency and the pulsation angular frequency in the corotating frame are significantly smaller than the Brunt-Väisälä frequency (Lee & Saio 1989;Townsend 2003), whereas the latter only requires slow stellar rotation to be valid.
In this paper, we are interested in the low-degree modes observed in a slowly rotating star, so we will consider the framework of Provost et al. (1981, hereafter P81), concentrating on the cases where the stratification is close to, but not exactly, adiabatic. Our results are not merely an addition to those of P81, but they provide better insight into the nature of r modes, and arguably amend those of P81. Assuming that the stellar interior is inviscid and the motions are adiabatic, we show that the = m mode with no radial nodes is the only almost toroidal Rossby mode which can be present, for uniform rotation. The corresponding eigenfunctions scale as r m .

Digest of Provost et al. analysis
This papers follows the formalism developed in P81, with minor changes in the notation. In this section, we summarise the main points of their method that we think are important for understanding our results. P81 considers a rotating star, in a corotating reference frame of basis vectors (ê r ,ê θ ,ê φ ), with the origin at the star's centre of mass and spherical coordinates (r, ϑ, ϕ), where r is the radial distance to the origin, ϑ is the polar angle, ϕ is the azimuthal angle. Rotation is assumed to be uniform, with angular frequency Ω = Ω ê z parallel tô e z = cos ϑê r −sin ϑê θ . The densities are in units of ρ = M /R 3 , the pressures are in units of p = GM 2 /R 4 , the times are in units of τ = (GM /R 3 ) −1/2 and the lengths are in units of R , where R and M are the radius and mass of the star and G the universal gravitational constant. The dimensionless angular frequency of the star is then denoted as ε = Ω τ . (1) For the Sun 1 , ε = Ω τ = 4.5×10 −3 . In a uniformly rotating star, the isobaric and isopsynic surfaces coincide with the level surfaces of constant total potential (gravitational and centrifugal). This is know as the Poincaré-Wavre theorem and holds whatever the equation of state of the gas. For slow rotators (ε 1), those surfaces can be expressed through a distortion term of the order of ε 2 and a function α determined by the internal structure. These surface levels can be used to implicitly define a set of curvilinear coordinates (x, θ, φ), 1 The subscript is replaced by the symbol to denote solar quantities throughout the paper.
where the new coordinate x is constant on surfaces of constant density and pressure. The surface of the star is an isobaric surface x = constant and we choose the normalisation so that x = 1 at the surface. With respect to the new variable x, the equilibrium pressure p(x) and density ρ(x) are then independent of θ and φ. The oscillations are treated as a small perturbation around this static equilibrium state. To study linear modes of oscillations, the temporal and longitudinal structure of all perturbed quantities and the displacement ξ are assumed to be proportional to 2 e i(mφ−σt) .
The equations governing the small amplitude, periodic, adiabatic oscillations of a uniformly rotating star are obtained by writing the linearised equations for the conservation of angular momentum, mass, and energy: ρ + ∇ · (ρξ) = 0, where p and ρ are the Eulerian perturbations of pressure and density, and γ = ∂ ln p ∂ ln ρ ad is the first adiabatic exponent. Here we have neglected the perturbation of the gravitational potential (Cowling's approximation) and have assumed the flows are adiabatic and the viscosity is negligible. Cowling's approximation was shown to be justified in most cases by P81. Saio (1982) argues that the adiabatic assumption is justified in stellar radiative zones except near the boundaries. The assumption that the flow is inviscid is also probably justifiable in radiative zones, however in stellar convection zones nonadiabatic mixing by turbulent convective motions and turbulent viscosity are likely to be relevant on timescales shorter than or comparable to the rotation period of the star. Notwithstanding, we chose to follow P81 in considering the inviscid and adiabatic case because this reveals interesting basic physics. In addition, we expect that modes which rely for their existence on viscosity or non-adiabatic processes will decay faster than modes which exist in the adiabatic, inviscid case. Equations (6)-(8) form a system of partial differential equations that requires appropriate boundary conditions to constitute a well posed boundary value problem. We will consider the boundary conditions at the centre and the surface of the rotating star. At the centre, the displacement must remain finite. The conservation of momentum across the non-spherical surface of the star requires that the Lagrangian pressure perturbation δp = 0 at the surface.
2.1. Series expansion in terms of ε P81 considers solutions for ξ = (ξ r , ξ θ , ξ φ ), p , ρ , σ, that are solutions of system Eqs. (6) -(8) for cases where the rotation rate is small. This motivates an expansion of the form Our methodology is identical to P81 but our notation is slightly different. P81 use the indices j = 0, 1, ... only for nonnull terms in the final expansion: for example, in their notation σ = Ωσ Provost 0 +Ω 3 (σ 0 σ 1 ) Provost (their equation 2). In our notation we index all terms. We thus write σ = σ 0 +εσ 1 +ε 2 σ 2 +ε 3 σ 3 +.... Our notation is closer to that of Smeyers et al. (1981). As is implicit in the notation of P81, symmetry arguments lead to the conclusion that the coefficients of the even powers of ε in the expansion of σ must be 0, so this difference with respect to P81 is purely one of notation. Similarly an inspection of Eqs. (6) -(8) reveals that the even and odd coefficients in the expansion of ξ θ , ξ φ , ξ r , p and ρ decouple. The solutions for the odd terms in the expansion of these quantities are trivial if we know the solution for the even powers (they correspond to the solutions for the expansion keeping only the even terms where all quantities except σ are multiplied by ε).
P81 is concerned with the generalisation of modes which are purely toroidal for non-rotating stars. The displacement vector of these toroidal modes satisfies ∇ · ξ = 0 and ξ r = 0. They also have σ 2 = 0. In the case where rotation is present, these toroidal modes become non-trivial and develop characteristics similar to Rossby waves in the Earth's atmosphere and oceans. They are often referred to as r modes, after the seminal work of Papaloizou & Pringle (1978), or quasi-toroidal mode according to the nomenclature of P81. They are called quasi-toroidal because, to zeroth-order, they have the same properties as the toroidal modes (∇·ξ (0) = 0, ξ (0) r = 0 and σ 0 = 0). To summarise, the quasi-toroidal modes are non-radial modes of low-frequency, whose radial displacement is small compared to their horizontal motion.
P81 performs the equivalent of a fourth-order expansion in terms of the small parameter ε in Eqs. (6) -(8). At the zeroth order, the system (6) -(8) only retains zeroth-order quantities and reduces to: which are the non-rotating toroidal modes. At higher-order, we have to consider the value of the Ledoux discriminant which is a measure of convective instability as will be discussed further in Sec. 3. The case A = 0 will be discussed later. If A 0, the second-order approximation of the conservation of momentum (ê r ,ê θ , andê φ components) yields where is the unperturbed gravity. The radial part of the conservation of momentum Eq. (12) involves the component of the displacement that is normal to isopotential surfaces, ξ (2) x = ξ (2) r − 2α cos θ sin θξ (0) θ . Equation (12) is P81's (7a). Combining our Eqs. (13) and (14) gives their (6a), which only involves zerothorder quantities.
Equations (13), (14) and (10) express the conservation of total vertical angular momentum and of mass, to zero order. By elimination of p (2) , they form a classical Legendre equation. This allows us to find exact solutions for σ 1 and the angular dependence of ξ (0) θ and ξ (0) φ in terms of , m, and the associated Legendre polynomials. The radial dependence of the eigenfunctions is not determined at this order. The solutions are Rossby waves on surfaces of constant x, as in Eq. (9a) of P81: Here we are interested in retrograde Rossby waves so we are only going to consider m > 0. At the fourth-order level of the approximation, the θ and φ components of the conservation of momentum can be combined to eliminate the pressure perturbation terms p (4) and p (2) , yielding Eq. (7c) from P81. Finally, the system of equations is closed by taking the second-order expansion of the continuity equation, yielding Provost et al.'s Eq. (7d). After some manipulation, the aforementioned closed system of equations can be reduced to a single differential equation for the amplitude C (0) ,m (x) of the horizontal displacement, Eq. (11) of P81, where λ 1 , λ 2 , λ 3 and λ 4 are functions of m and and are defined in P81. When A 0, for x ∈ (0, 1) and with appropriate boundary conditions, Eq. (19) is a Sturm-Liouville eigenvalue problem.

Boundary condition at the centre
Near the centre x ∼ 0, the quantities that appear in Eq. (19) are in order of magnitudes as follows where the subscript c denotes values at the centre. Therefore, by asymptotic analysis of Eq. (19) in the vicinity of x = 0, we have Using the Frobenius method (Bender & Orszag 1999), we look for a solution of the form C (0) ,m (x) = ∞ j=0 a j x j+β in the vicinity of x = 0. According to (25), the coefficient of the lowest power of x must satisfy which has solutions β = and β = − − 1. The only solution that ensures the regularity of and

Boundary condition at the surface
At the free surface the Lagrangian pressure perturbation must vanish. This is a very stringent condition for the existence of non-sectoral Rossby waves in the inviscid case, and its implications depend on whether ρ/A vanishes or not at the surface.
2.3.1. Boundary conditions for a star where the surface has ρ/A = 0 As a model of a star P81 considers a complete polytrope 3 characterised by a polytropic index n. The stellar surface for such models is defined by ρ = p = 0. In general, for a complete polytrope or otherwise, when ρ/A vanishes at the surface, δp = 0 is met at the surface as long as C (0) ,m and its derivative are regular. This is the boundary condition that was used in P81.
For the particular case of a polytropic stellar model, the behaviour of the quantities that appear in (19) can be approximated by Taylor expansion ofθ, the solution to the Lane-Emden equation, near x = 1 as follows A ∼ − n − n + 1 γ 3 Polytropes will be described in more details in Sec. 3.
By asymptotic analysis of Eq. (19) in the vicinity of x = 1, we have where We remind the reader that λ 1 , λ 2 and λ 3 are simple, but tediously long, rational functions of and m. They are given explicitly in P81.
Thus the condition of regularity for Here and from now on, the notation in bracket with subscript S means that we take the value at the surface.

Boundary conditions when ρ/A 0 at the surface
Truncated polytropes have been used as models of stars which include an atmosphere (e.g. Hendry 1993; Bogdan & Cally 1995). In these models, the polytrope is truncated at some location, with non-zero pressure and density, that represents the stellar surface. In the limit that the density of the atmosphere is small compared to the surface density, the boundary condition corresponds to a free surface (with vanishing Lagrangian pressure perturbation at the surface). Such truncated polytropes with a free surface have previously been used in the study of helioseismic acoustic waves (e.g. Bogdan & Cally 1995). The analysis for Rossby waves proceeds differently according to whether = m or not.
Cases where = m. When = m, the free-surface boundary yields one condition In the cases where {ρ/A} S 0, the boundary condition reduces to Cases where m. When m, the free-surface boundary condition involve two associated Legendre polynomials P m −1 and P m +1 , whose coefficients must both vanish and this thus yields In general, ρ/A does not vanish at the surface of stars, Eqs. (40) and (41)  ,m (x) = 0. Consequently, for non-vanishing density at the surface, there is no non-trivial solution when m. The sectoral mode is the only quasi-toroidal mode that can satisfy the free-surface boundary condition in the inviscid case. Let us note that the non-sectoral modes found by Wolff & Blizard (1986), using P81's derivations applied to the Sun, are obtained by imposing dC (0) ,m dx = 0 at x = 0.999, which is inconsistent with a free-surface boundary condition. We remark as well that the no-penetration boundary condition is ξ (2) so that Wolff & Blizard (1986) appears to be also inconsistent with the no-penetration boundary condition.

Comparison of boundary conditions of polytrope and truncated polytrope.
It is interesting to compare the boundary conditions for the complete and truncated polytropes. For = m these are Eqs. 37 and 39 respectively. The quantity Q(1) is given by Eq. 36 evaluated at the surface and depends on m. For the special cases where = m, γ = 5/3 and n = 3/2, we find that Q(1)/(n + 1) = m. The values of Q/(n + 1) evaluated at the surface of the full polytrope (x = 1) are given in Table 1 for some of the other cases studied in this paper with = m = 3. In Table 1 we see that Q/(n + 1) is approximately equal to m near n = 1.5 (it is exactly m at n = 1.5). A consequence of this is that in the neighbourhood of n = 1.5, the boundary condition for the truncated polytrope is equivalent to that for the complete polytrope. This indicates that the results have some robustness to the details of the model. For m the situation is more complicated, and no solutions exist for the truncated polytrope. This is also the case for the complete polytrope when n = 1.5, (see Appendix I in P81). For other values of n, the solution for the complete polytrope of course exists, and is as described in P81.

Results
To solve the boundary-value problem for (19), we use an implementation of a fourth-order collocation algorithm based on control of residuals, provided by the function solve_bvp from the scipy.integrate python module. We follow P81 in considering polytropes, characterised by a polytropic index n. In the stellar context, a polytrope is a gas spheroid in gravitational equilibrium, where the pressure is related to the density by the relation where K is a constant of proportionality. By definition of the free surface, a complete polytrope has p = 0 and ρ = 0 at the surface. For all polytropes, the Ledoux discriminant A is a monotonic function of x and it is everywhere positive when n < 1.5 (superadiabatic stratification). The polytrope with index n = 1.5 has the interesting property that then A = 0 everywhere in the star (adiabatic stratification). Finally, for n > 1.5, A is everywhere negative (sub-adiabatic stratification. The Ledoux discriminant is the argument of the criterion for convection which develops or not whether A > 0 or A < 0 (Ledoux & Walraven 1958). So we shall also refer to the convectively unstable case when n < 1.5, neutrally stratified case when n = 1.5, and radiative case when n > 1.5. The shape function α(x) in Eq.
(2) can be derived for the distorted polytrope (Chandrasekhar 1933), and is obtained by solving the Lane-Emden equation modified for rotation. We obtain the full polytrope by setting the surface x = 1 at the first zero of the Lane-Emden function, and the truncated polytrope by setting the surface at x = 0.999.
To test our solver, we computed the eigenvalues for n = 1 and n = 3 given in Tables 1 and 2 of P81 to the same decimal place accuracy, for a complete polytrope using the boundary conditions (27), (28) and (37). We also reproduced the radial dependence of C (0) ,m (x) for the same ( , m, k) as in their figures 1, 2 and 3.
We also computed the eigen-solutions for the truncated polytropes when = m, using the boundary conditions (27), (28) and (39). As is to be expected from Table 1, for sectoral modes, the eigen-solutions for the full polytrope and the same polytrope truncated very close to the surface are nearly identical close to n = 1.5, and the differences even for n = 1 and n = 3 are small and not distinguishable in the figures that follow. Hence we only show solutions for a polytrope truncated at x = 0.999.
For cases when m, we refer the reader to P81 for the solutions in a complete polytrope and again note that there are no such solutions satisfying the free-surface boundary conditions for truncated polytropes.

Non-adiabatic stratification is essential to the normal force balance for all sectoral r-modes except those with no radial nodes
To get a better understanding of the effect of stratification on the r modes, we also solve the problem for 1 ≤ n ≤ 3. As an example, Figure 1 shows the eigenfrequencies found for = 3, and m = 3 as a function of the polytropic index n for several values of the radial order k. We see that for all modes, with the notable exception of the sectoral mode with k = 0, the eigenfrequencies become increasingly large as n gets close to 1.5 (where A = 0). Let us stress here that the derivation of Eq. (19) is obtained by a singular perturbation method. Thus the asymptotic expansion (9) does not necessarily converge, and the solution is valid when |σ 3 | σ 1 /ε 2 , which means here ε 2 |σ 3 |/(2πτ ) 10 5 pHz for = m = 3 (see Fig. 1).
Solving Eq. (19) allows the derivation of all the perturbed quantities of our problem. Figure 2 shows the radial force bal-A&A proofs: manuscript no. _arXiv  Fig. 1. Third-order term in the frequency expansion σ 3 ε 3 as a function of the polytropic index n, for the sectoral r modes = m = 3 and several radial orders k. The corresponding frequency term in the incompressible case, which has only = m, k = 0 solutions, is given by a blue dotted line (see text). The frequency is given in pHz, using the relevant solar quantities, and is displayed in symmetric log scaling, with linear scaling between ±6 pHz. The quantity τ is our unit of time and corresponds to the Sun's dynamical time scale.
ance associated to r modes given by Eq. (12) for the polytrope n = 1.49 and for = m = 3. In the momentum equation, the Coriolis term is balanced by a non-trivial combination For the modes with k 0 (Fig. 2, left) the non-adiabatic stratification plays an essential role in the radial force balance. This turns out to be true for all the modes (i.e. also those with m) except the = m, k = 0 modes (this will be discussed in Section 4). For the case with = m and k = 0, the radial force balance is essentially between the Coriolis term and p 1/γ ∂ ∂x p (2) p 1/γ (Fig. 2, right) -the term involving the non-adiabatic stratification plays essentially no role. This is a special property of the = m, k = 0, r modes.
3.2. The depth dependence of the sectoral modes of zero radial order is x m for quasi-adiabatic stratification Figure 3 shows the radial structure C (0) ,m (x) for ( , m, k) = (3, 3, 0) for different values of the polytropic index n. As shown by P81, when A = 0, the only non-trivial modes are the sectoral modes with = m. The solution must then have the form C (0) m,m (x) = x m , which has no radial nodes. However, in this case there is no finite value of σ 3 that can give a finite radial displacement at the surface of the complete polytrope, since in this case the solution is divergent at the surface (P81). This problem does not exist for truncated polytropes.
In the limit of n → 1.5, as we can see in Fig. 3, we find that the depth dependence of the sectoral mode with k = 0 is very close to x m . Since this mode is the only one allowed to exist in the case A = 0, we find that there is no discontinuity of solutions for this mode, and the solutions slowly depart from an x m dependence as the stratification departs from neutral.
The solution near n = 1.5 (where A = 0) has a radial dependence proportional to x m . This is also the form of the incompressible Rossby wave (Bryan 1889;Provost et al. 1981), which has no dependence on the stratification, as will be discussed in Section 4.

Symmetries of the eigenfunctions about n = 1.5
Figs. 4 and 5 shows eigenfunctions for the (3, 3, 1) and (3, 3, 0) modes, both for n = 1.49 and n = 1.51. These values of n were chosen because the background is convectively stable in one case and unstable in the other. Here, instead of the radial component of the displacement, we plot the more relevant component of the displacement ξ (2) x that is normal to isopotential surfaces. The solution for the case with k 0 (Fig. 4) shows that most of the physical quantities vary only slightly across n = 1.5 except for ξ x , which flips sign. The reason for this flip in sign can be inferred from Eq. (12), and will be discussed in Section 4. Conversely, in the case = m and k = 0 (Fig. 5) all the quantities vary only slightly across n = 1.5 and in particular the sign of ξ x does not change.
This difference results from the fact that the radial force balance in the cases with m or k 0 the super-or subadiabaticity plays an essential role and ξ x changes sign with A. By way of contrast, the sectoral mode of zero radial order does not depend essentially on the adiabaticity and ξ x varies smoothly with n.

Discussion
The solutions with = m and k = 0 are qualitatively different from all other solutions. We showed the mode with ( , m, k) = (3, 3, 0) as an example. In the neighbourhood of n = 1.5 (corresponding to A ∼ 0), all values of , m, and k admit solutions in a complete distorted polytrope, as discussed above. Only the sectoral modes ( = m) would be admissible if the density at the surface did not vanish. For the cases with = m and k = 0, the term Aξ (2) r remains small and the solution is not sensitive to A; the solution does not depend on non-adiabatic stratification. All other solutions are baroclinic modes and the term Aξ (2) r is approximately constant in the neighbourhood of n = 1.5, with |ξ r | being arbitrarily large in the limit n → 1.5 (A → 0) and flipping sign at n = 1.5. Similar reasoning explains the behaviour shown in Fig. 1 for the eigenfrequencies.
The radial displacements associated with solutions where the non-adiabatic stratification is essential (i.e., where the contribution to the radial force balance from Agξ r is substantial), go to infinity as n approaches 1.5 (where A → 0). These solutions cease to be quasi-toroidal. Hence the only solutions which are valid 5 near n = 1.5 are those with = m and k = 0.
The lack of quasi-toroidal solutions as n approaches 1.5, except when = m and k = 0, is not a consequence of the chosen expansion. It follows from the statement that the system needs to be in both horizontal and radial force balance. In the absence of super-or subadiabatic stratification and viscosity, both the horizontal and radial force balances independently determine the pressure perturbation. The only case in which the two determinations of the pressure perturbation are consistent and quasitoroidal are the special cases where = m, k = 0 and the horizontal displacement scales with x m . It is here that we make contact with the example Rossby waves in an incompressible, unstratified spherical shell discussed by Löptien et al. (2018).
The existence of non-trivial solutions only in the = m case with an r m radial dependence is also that found by P81 in the incompressible case with arbitrary stratification. This is understood by considering that the horizontal force balance sets the horizontal structure of the pressure perturbation independently of A or incompressibilitity. In the incompressible case (or in the case with A = 0) and for quasi-toroidal motion, the pressure perturbation must alone balance the radial component of the Coriolis force, and this is only possible in the case = m, and results in an r m dependence, and has no dependence on stratification.
The expression of the third-order term σ 3 of the frequency expansion for the = m mode in the incompressible case for non-spherical shapes have been obtained by Bryan (1889) and P81, Fig. 1 shows these values of σ 3 for m = 3 as a function of n as a blue dotted line 6 ; they are the same order of magnitude as the eigenfrequency associated to the sectoral mode of zero radial order. This is to be expected since here, we are seeking solutions that are quasi-toroidal, which means that the flow is divergence free to zeroth-order. All these solutions are approximately the same, although the problem is formulated differently in each case (compressible vs incompressible, polytrope vs neutrally stable stratification). In the incompressible case or the neutrally stratified case, the radial force balance is independent of the normal displacement. In this case, the problem admits only solutions that have = m and k = 0 and an x m dependence. For the quasi-toroidal modes, the solutions are again approximately the same because in this case also the toroidal components of the motion are the same (spherical harmonics characterised by , m). This in turns determines the pressure perturbation, which then must balance the radial component of the Coriolis force if the motions are to remain quasi-toroidal. This only happens for = m and k = 0. For m or k 0, the modes will develop substantial radial velocity whenever the stratification is close to, but not strictly, neutral.
4.1. Non-spherical geometry is important away from n = 1.5 We also found the eigenvalues and eigenfunctions for the case where we artificially set α = 0, which then treats the star as a perfect sphere. Consistent with P81, we found the changes to σ 3 were substantial relative to σ 3 (i.e. can have the wrong sign), but since |σ 3 |ε 3 is small, this amounts to a change of the order of 10 pHz for the solar rotation rate. More important are the changes in the radial structure of the eigenfunctions, which can be seen by comparing Fig. 6 with the equivalent in Fig. 3. This would suggest that the distortion of the background star due to rotation should be included when determining the radial structure of the eigenfunctions, especially for non-adiabatic stratifications. Fig. 4. Meridional cuts (x,θ) showing the leading order terms, for solar rotation, of the latitudinal σ 1 ξ θ and azimuthal σ 1 ξ φ flow velocity, the relative normal displacement ξ x /H p , the relative pressure p /p and density ρ /ρ perturbations (from left to right respectively, as labelled in the colour bar) for the ( , m, k) = (3, 3, 1) mode and for n = 1.49 (top) and n = 1.51 (bottom). H p is the pressure height scale. The radial dependence of all the quantities has been normalised consistently so that at phase 0, a dimensionless velocity of 1 on this scale corresponds to a velocity of 1 m/s on the Sun.

Conclusions
In this paper, we used polytropes to understand the fundamental properties of quasi-toroidal modes for slowly and uniformly rotating stars.
The sectoral r modes of zero radial order are qualitatively different from the other r modes in that they do not rely on nonadiabatic stratification to balance the radial component of the Coriolis force. This is critical when the stratification is close to neutral (for polytropes, in the neighbourhood of n = 1.5). In this neighbourhood, the modes with m or k 0 all involve large radial displacements and are no longer nearly toroidal. The = m, k = 0 modes retain a small radial displacement through the star (except possibly at the surface if the surface has p = ρ = T = 0 as pointed out by P81). The depth dependence of the horizontal displacement is close to x m , as it is the only so-lution allowing both horizontal and vertical force balance, in the absence of viscosity and the lack of a buoyant contribution from non-adiabatic stratification.
Consequently, in the case of the Sun, we speculate that only the = m, k = 0 quasi-toroidal modes can exist in the convection zone, which is very close to adiabatically stratified (corresponding approximately to the polytropic index n = 1.49). They are presumably the modes that are observed at the solar surface. Figure 7 shows the kinetic energy density associated with those modes on a meridional cut. In Fig. 8, we plot the corresponding radial dependence at the equator. These figures suggest that solar Rossby waves have diagnostic potential because different modes have different radial and latitudinal distribution of the kinetic energy density. It can be noticed that the kinetic energy density of the modes peaks in the interior (not at the surface as p modes would), and that the m = 3 mode has a kinetic energy density that peaks near x = 0.75. However the results of this study cannot be directly applied to the solar case, since the Sun is not a uniformly rotating inviscid polytrope. In particular, latitudinal differential rotation will modify the latitudinal eigenfunctions (Gizon et al. 2019, in prep.). Also, the inclusion of radial differential rotation has not been considered here. Furthermore, the Sun has a convectively stable radiative interior, beneath its convective envelope that requires a careful consideration of the matching conditions.