© The Authors 2026

1. Introduction

Kelvin waves are a class of waves that contain a component of surface gravity waves. However, their definition has varied over time. The first reference to Kelvin’s work on waves is likely found in Cowling (1941), which is dedicated to the non-radial oscillations of a polytropic star. In Kelvin’s original work (Thomson 1863), waves are defined as oscillations of a self-gravitating sphere made of an incompressible fluid (see Lamb 1932, §262), but Cowling did not connect his set of waves to those of Kelvin. In fact, Chandrasekhar & Lebovitz (1963) appear to have been the first to explicitly introduce Kelvin waves into the astrophysical literature. Later, Hurley et al. (1966) performed a full numerical analysis of non-radial oscillations of polytropes and provided the first non-perturbative approach to Kelvin-mode frequencies that fully included compressibility, thus moving a step beyond the works of Chandrasekhar & Lebovitz (1963) and Chandrasekhar (1964). The nomenclature in the astrophysical literature has changed, and these Kelvin modes are now referred to as f-modes (fundamental gravity modes; Cowling (1941)). Robe & Brandt (1966) established the link between Kelvin modes and f-modes. Since then, Cowling’s terminology has prevailed.

Also in 1966, a new Kelvin wave appeared in the literature with the work of Matsuno (1966), who studied the oscillations of the Earth’s atmosphere and oceans. These Kelvin waves refer to two families: coastal Kelvin waves and equatorial Kelvin waves. They are also related to Kelvin’s work (Thomson 1880), which is dedicated to the study of horizontal oscillations of a thin fluid layer in a rotating frame. Their study was a follow-up of Laplace’s seminal work on tides (Laplace 1799), which introduced the shallow water approximation (see Sect. 2). Matsuno (1966) revisited the system originally set out by Laplace and introduced a less restrictive assumption than that used by Kelvin, namely the β-plane approximation, which includes a linear variation of the projected rotation vector as a function of latitude. This approximation remains a standard framework in Earth sciences. This new class of Kelvin waves exists only when a background rotation is present. It can be viewed as a special class of gravito-inertial waves¹.

Two years after Matsuno’s study, Wallace & Kousky (1968) reported the first observation of equatorial Kelvin waves in the Earth’s stratosphere. Kelvin waves play a major role in terrestrial fluid dynamics, notably through their impact on phenomena such as El Niño (Cushman-Roisin 1994). However, while Kelvin waves have been extensively studied in oceanography and atmospheric sciences, their role in astrophysical contexts remains relatively unexplored and is usually treated as a secondary question. A first in-depth discussion was provided by Townsend (2003), which relied on the traditional approximation of rotation (TAR). Later, while analysing the light curves of α Ophiuchi (Rasalhague) obtained with the Microvariability and Oscillations of STars (MOST) satellite (Walker et al. 2003), Monnier et al. (2010) suggested the possible detection of equatorial Kelvin waves in that star. More recently, Takata et al. (2020) introduced the ‘internal equatorial Kelvin waves’, referring to the role of buoyancy in their dynamics. However, the specificity of equatorial Kelvin waves was revealed by the seminal work of Delplace et al. (2017) who underline their topological nature. These properties and their consequences in stellar physics are studied in very recent works such as Perez et al. (2022, 2025) and Leclerc et al. (2022, 2024). For completeness, we should mention that Kelvin waves, in their diverse definitions, have also been studied for their potential efficiency in the emission of gravitational radiation by neutron stars (Andersson 2021).

The stability of equatorial Kelvin waves is of great interest given their known properties and their robustness arising from their topological nature, particularly when rotation and gravity combine. Indeed, their prograde nature and equatorial confinement make them interesting candidates for launching matter into orbit around rapidly rotating stars, provided that they can be destabilised by differential rotation. Hence, they may be important contributors to the Be phenomenon, which is now associated with near-critically rotating B stars (Porter & Rivinius 2003). The purpose of this paper is to explore this possibility using the simplified setup of an incompressible fluid in a rotating spherical shell.

The paper is organised as follows. To start from a well-posed problem, we first investigate how equatorial Kelvin waves transform when moving from the shallow water approximation (Sect. 2) to a thick fluid spherical shell mimicking a stellar envelope (Sect. 3). We then impose a shellular differential rotation on this fluid layer, including viscosity, and investigate the stability of the waves (Sect. 4). We conclude the paper with the conclusions and outlook.

2. The shallow water system

To investigate Kelvin modes in a fluid inside a spherical shell rotating at constant angular velocity Ω = Ωe_z, we begin with the shallow water system, where these modes have been studied.

The shallow water system considers a thin layer of fluid of thickness H₀, which is very small compared to the wavelength of the modes. In this case, the equations can be simplified by neglecting the radial component of the velocity. The radial momentum equation ensures predominant hydrostatic balance in the vertical direction.

These two simplifications imply the removal of the tangential component of the rotation vector, which means that the shallow water system falls within the framework of the traditional approximation of rotation². Furthermore, the pressure variable can be expressed as a function of a surface elevation variable ζ_*. The linearised shallow water equations thus read as (but see also Cushman-Roisin 1994)

$$\begin{array}{c}\hfill \{\begin{array}{cc}& {\partial }_t\mathit{v}+2\mathrm{\Omega }{\mathit{e}}_{\mathit{z}}\times \mathit{v}=-g\mathbf{\nabla }{\zeta }_{\ast }\hfill \\ \hfill & {\partial }_t{\zeta }_{\ast }+H_0\mathbf{\nabla }·\mathit{v}=0\hfill \end{array},\end{array}$$(1)

for momentum and mass conservation, respectively. The Coriolis angular frequency is given by 2Ω, and g is the effective gravity, defined as the sum of the gravitational and centrifugal accelerations. The velocity field is v = v_θe_θ + v_ϕe_ϕ.

Matsuno (1966) found solutions of the shallow water system (1) using the β-plane approximation. He identified four types of gravito-inertial modes: Poincaré modes, Rossby modes, Kelvin modes, and Yanai modes, which we now briefly discuss in the context of spherical geometry.

Since perturbations develop over an axisymmetric background, we can write their components as

$$\begin{array}{c}\hfill (v_{\theta },v_{\varphi },{\zeta }_{\ast })=(v_{\theta }(\theta ),v_{\varphi }(\theta ),{\zeta }_{\ast }(\theta ))e^{i{\omega }_{\ast }t+im\varphi },\end{array}$$(2)

where θ is the co-latitude and ϕ is the longitude. Next, we introduce the variable μ = cos θ and define the following non-dimensional quantities:

$$\begin{array}{c}\hfill w_{\theta }=\frac{v_{\theta }sin\theta }{\mathrm{\Omega }{H}_0},w_{\varphi }=\frac{v_{\varphi }sin\theta }{\mathrm{\Omega }{H}_0},\zeta =\frac{{\zeta }_{\ast }R}{H_0^2}.\end{array}$$(3)

Introducing $\omega =\frac{{\omega }_{\ast }}{\mathrm{\Omega }}$, we rewrite (1) as

$$\begin{array}{c}\hfill \{\begin{array}{cc}& i\omega w_{\theta }=2\mu w_{\varphi }+\mathrm{\Gamma }(1-{\mu }^2){\partial }_{\mu }\zeta \hfill \\ \hfill & i\omega w_{\varphi }=-2\mu w_{\theta }-im\mathrm{\Gamma }\zeta \hfill \\ \hfill & i\omega \zeta ={\partial }_{\mu }{w}_{\theta }-\frac{\mathit{im}}{1-{\mu }^2}{w}_{\varphi }\hfill \end{array},\end{array}$$(4)

with

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{g{H}_0}{{\left(\mathrm{\Omega }R\right)}^2}.\end{array}$$(5)

System (4) thus depends on a single non-dimensional parameter, Γ, which can also be expressed as Γ = (c/c_eq)², where $c=\sqrt{g{H}_0}$ is the phase and group velocities of surface gravity waves in a non-rotating plane layer, and c_eq is the background equatorial velocity.

2.1. Rossby modes

Rossby modes, also called planetary waves, are a subclass of inertial modes whose restoring force is the Coriolis force. They are characterised by their low frequency (ω ≪ 1) and their 2D nature, which arises naturally in the shallow water approximation.

Unlike gravity modes, Rossby waves do not require a free surface to exist. They are also supported in a fluid shell with a rigid, non-deformable upper boundary (Longuet-Higgins 1964). Thus, we assume ζ → 0 but keep Γζ finite to preserve the pressure perturbation associated with the velocity field. Hence, the third equation of (4) reduces to

$$\begin{array}{c}\hfill {\partial }_{\mu }{w}_{\theta }-\frac{\mathit{im}}{1-{\mu }^2}{w}_{\varphi }=0,\end{array}$$

which implies a divergence-free velocity field. This motivated the introduction of a stream function χ, through which we express the velocity field as v = ∇ × (χe_r). Thus,

$$\begin{array}{c}\hfill w_{\theta }=im\chi ,w_{\varphi }=(1-{\mu }^2){\partial }_{\mu }\chi .\end{array}$$

We can then rewrite system (4) in terms of the stream function as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_h\chi =-\frac{2m}{\omega }\chi ,{\mathrm{\Delta }}_h=(1-{\mu }^2){\partial }_{\mu }^2-2\mu {\partial }_{\mu }-\frac{m^2}{1-{\mu }^2},\end{array}$$

where Δ_h is the horizontal Laplacian, whose eigenfunctions are the spherical harmonics. Setting χ = χ_m^ℓY_ℓ^m, we deduce the dispersion relation of Rossby modes, namely

$$\begin{array}{c}\hfill \omega =\frac{2m}{\ell (\ell +1)}\mathrm{with}\ell \ge |m|.\end{array}$$(6)

This relation is well-known (Longuet-Higgins 1964; Rieutord 2015). These waves are retrograde and their frequency is bounded |ω|≤1. We also note that these waves are independent of Γ.

2.2. Poincaré modes

In the absence of rotation, an incompressible fluid with a free surface supports surface gravity waves. These waves occupy the high-frequency band, and their dispersion relation depends on the fluid depth. However, when rotation is present, the Coriolis force modifies their dynamics. In the shallow water approximation combined with the so-called f-plane approximation, their dispersion relation is ω² = c²k² + f², where k is the wave number and f = 2Ω sin θ is the projected Coriolis frequency at latitude θ. These waves are commonly known as Poincaré waves (Cushman-Roisin 1994; Delplace et al. 2017).

When the domain is the whole surface of the sphere, Poincaré modes are rotation-modified surface gravity modes, where both gravity and the Coriolis force act as restoring forces. Unlike planetary waves, their frequency is unbounded.

In this high-frequency regime, the dominant solutions correspond to pure surface gravity modes with a zeroth-order frequency given by

$$\begin{array}{c}\hfill {\omega }_0=\pm \sqrt{\mathrm{\Gamma }\ell (\ell +1)},\end{array}$$

which is associated with the spherical harmonics Y_ℓ^m as the eigenfunction. Introducing a first-order correction due to rotation, the frequency becomes

$$\begin{array}{c}\hfill \omega =\pm \sqrt{\mathrm{\Gamma }\ell (\ell +1)}+\frac{m}{\ell (\ell +1)},\end{array}$$(7)

or, with dimensional quantities

$$\begin{array}{c}\hfill {\omega }_{\ast }=\pm \frac{\sqrt{g{H}_0}}{R}\sqrt{\ell (\ell +1)}+\frac{m\mathrm{\Omega }}{\ell (\ell +1)}.\end{array}$$(8)

The corresponding eigenfunctions are Hough functions (Wang et al. 2016), which reduce to spherical harmonics in the high frequency limit. We give the derivation of (7) in Appendix A.

2.3. Equatorial modes

The spectral space associated with equations (4) has topological properties that are related to the breaking of the time-reversal symmetry. This leads to two families of equatorially trapped prograde modes, namely the Kelvin and Yanai modes³ (Delplace et al. 2017). Kelvin modes are equatorially symmetric, whereas Yanai modes are anti-symmetric. Because of their topological origin, these two series of modes are robust to variations in the background physics. Hence, we expect their presence if the background is no longer as simple as in the shallow water model. In the following, we extend the shallow water model to a differentially rotating thick layer. However, as a preliminary step, we first derived the dispersion relation of these modes. To this end, since these modes are equatorially trapped, we assume μ² ≪ 1, so that (4) can be simplified as follows:

$$\begin{array}{c}\hfill \{\begin{array}{cc}& i\omega w_{\theta }=2\mu w_{\varphi }+\mathrm{\Gamma }{\partial }_{\mu }\zeta \hfill \\ \hfill & i\omega w_{\varphi }=-2\mu w_{\theta }-im\mathrm{\Gamma }\zeta \hfill \\ \hfill & i\omega \zeta ={\partial }_{\mu }{w}_{\theta }-im{w}_{\varphi }\hfill \end{array},\end{array}$$(9)

which are the β-plane equations (e.g. Matsuno 1966; Cushman-Roisin 1994).

2.3.1. Kelvin modes

Among equatorial solutions, Matsuno (1966) identified the Kelvin mode as a particular solution of the shallow water system in which the latitudinal velocity vanishes. Thus, setting w_θ = 0, we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}& 0=2\mu w_{\varphi }+\mathrm{\Gamma }{\partial }_{\mu }\zeta \hfill \\ \hfill & i\omega w_{\varphi }=-im\mathrm{\Gamma }\zeta \hfill \\ \hfill & i\omega \zeta =-im{w}_{\varphi }\hfill \end{array},\end{array}$$(10)

from which we derive the dispersion relation

$$\begin{array}{c}\hfill \omega =-m\sqrt{\mathrm{\Gamma }}.\end{array}$$(11)

We chose the sign to ensure the non-divergence of the corresponding eigenfunction,

$$\begin{array}{c}\hfill \zeta ={\zeta }_0{e}^{-\frac{{\mu }^2}{\sqrt{\mathrm{\Gamma }}}}.\end{array}$$(12)

The eigenfunction becomes equatorially trapped as Γ decreases, i.e. as the rotation rate increases. Equation (11) shows that Kelvin waves are prograde and regularly spaced in frequency. Equation (12) shows that the latitudinal width of these equatorial waves scales as Γ^1/4.

2.3.2. Yanai modes

The second equatorially trapped wave in the shallow water system is the Yanai mode. It is anti-symmetric with respect to the equator for w_ϕ and ζ, and has a non-zero meridional velocity at the equator (Zeitlin 2018). We retrieve it by imposing

$$\begin{array}{c}\hfill \zeta =\frac{2i\mu }{\omega \sqrt{\mathrm{\Gamma }}+m\mathrm{\Gamma }}{w}_{\theta }.\end{array}$$(13)

This leads to a Gaussian solution for w_θ, namely

$$\begin{array}{c}\hfill w_{\theta }=w_0{e}^{-\frac{{\mu }^2}{\sqrt{\mathrm{\Gamma }}}},\end{array}$$(14)

which is symmetric with respect to the equator, as expected. The dispersion relation of these modes reads

$$\begin{array}{c}\hfill \omega =-\frac{m\sqrt{\mathrm{\Gamma }}}{2}\pm \frac{1}{2}\sqrt{m^2\mathrm{\Gamma }+8\sqrt{\mathrm{\Gamma }}}.\end{array}$$(15)

2.4. Numerical solutions

To further illustrate and visualise the eigenspectrum of the shallow water system, we next solved (4) numerically. To this end, we discretised (4) on the Gauss-Lobatto collocation grid associated with Chebyshev polynomials (Fornberg 1998). This spectral discretization ensures rapid convergence of the numerical solutions. We then solved the resulting eigenvalue problem with the classical QZ algorithm (e.g. Valdettaro et al. 2007; Chatelin 2012).

The value of Γ plays a crucial role in the shape of gravito-inertial modes, especially for equatorial modes. On the one hand, Γ determines their equatorial confinement, as shown in (12) and (14). The smaller Γ, the more confined the Kelvin and Yanai modes are at the equator. On the other hand, equatorial modes are characterised by frequencies that asymptotically approach those of other gravito-inertial modes. Specifically, Kelvin modes connect Poincaré modes at high wave numbers, while Yanai modes transit between Rossby modes and Poincaré modes. Figure 1 illustrates this. If Γ is too large, the dispersion relation of gravito-inertial modes shows no distinction for equatorial modes, revealing only Poincaré and Rossby modes⁴. For equatorial modes to stand out among other modes, Γ must be sufficiently small, typically less than 10⁻¹.

Fig. 1. Dispersion relation of gravito-inertial modes in the shallow water system for Γ = 0.01. Left: Eigenfrequency of modes symmetric with respect to the equator as a function of the azimuthal wave number m. Right : Same as the left panel but for antisymmetric modes. Red dots show Poincaré modes, blue dots show Rossby modes, and green dots show Kelvin (left) and Yanai (right) modes.

2.5. Planetary and stellar context

The shallow water model demonstrates the existence of Kelvin and Yanai modes when the parameter Γ is small compared to unity. On Earth, the kilometric thickness of the atmosphere or the oceans makes Γ on the order of 10⁻¹, a value at which equatorial modes are expected to be prominent. This has been confirmed through years of oceanic observations. For a rotating main-sequence star, no such thin layers exist and surface waves can spread in depth. However, because our study considers Be stars, whose rotation is near critical, the effective equatorial surface gravity is weak. We therefore expect small values of Γ even in a deep layer. Assuming H₀ = R/2, we obtain

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{g_{\mathrm{eq}}^{\mathrm{eff}}}{2{\mathrm{\Omega }}^{2}R}=\frac{1}{2}(\frac{{\mathrm{\Omega }}_k^2}{{\mathrm{\Omega }}^2}-1),\end{array}$$(16)

where $g_{\mathrm{eq}}^{\mathrm{eff}}$ is the effective gravity at the equator and ${\mathrm{\Omega }}_k=\sqrt{GM/R^3}$ is the Keplerian angular velocity at the same location. For a star rotating at 90% of the critical angular velocity, Γ ∼ 0.1, which is similar to the Earth value. Rapidly rotating intermediate-mass stars evolve at nearly constant angular momentum (Gagnier et al. 2019) and slowly enough to remain in a quasi-steady state. As shown by Mombarg et al. (2024), the ratio Ω²/Ω_k² increases to unity during the main sequence. Hence, equatorially confined Kelvin or Yanai modes should be expected when the star reaches near-critical rotation.

The foregoing considerations prompted us to examine the role of shell thickness on the properties of Kelvin waves. We wanted to determine whether the low value of the Γ parameter was still a sufficient condition for the existence of equatorially confined Kelvin modes.

3. Extension to the thick shell

To appreciate the effects of thickening the fluid layer, we considered a spherical shell of fluid extending from r = η to r = 1, where the outer surface of the shell is free to support gravity waves. Simultaneously, we introduced viscosity to numerically regularise the problem and avoid singularities of inertial modes in the spherical shell. Indeed, as shown in Rieutord & Valdettaro (1997) and Rieutord et al. (2001), the Poincaré equation, which governs the inertial modes of an inviscid rotating fluid, is spatially hyperbolic and has mainly singular solutions, wrapped around attractors of characteristics generated by boundary conditions. Viscosity smooths these singularities into oscillating shear layers (Rieutord et al. 2002).

Moreover, we ignored self-gravity, hence working under Cowling’s approximation. Surface gravity waves combined with self-gravity are unstable at sufficiently high rotation rates and allow the bifurcation of the axisymmetric MacLaurin spheroid to triaxial spheroids when rotation is high enough (Chandrasekhar 1969). We neglected these complex matters because our incompressible fluid layer serves as a simple model for describing Kelvin waves in stellar envelopes, which exhibit density variations, entropy stratification, and at least differential rotation. As shown in Appendix B, our uniformly rotating fluid shell is always stable. Hence, we considered the ratio

$$\begin{array}{c}\hfill \gamma =g/{\mathrm{\Omega }}^{2}R,\end{array}$$(17)

namely the surface gravity divided by the centrifugal acceleration, and allowed it to be greater or less than unity to explore the physical or mathematical properties of the global modes. As previously observed, γ < 1 is not physically absurd if g refers to the effective gravity in the equatorial region of a rapidly rotating star.

3.1. Mathematical formulation

Consider a fluid particle in a rapidly rotating fluid. Its velocity can be expressed as Ω × r + v, with v small compared to Ω × r. The particle dynamics is characterised by a small Rossby number and is governed by a linear system. Using Ω⁻¹ as the timescale and the outer radius of the shell R as the length scale, the linearised momentum and mass conservation equations in an inertial frame read as

$$\begin{array}{c}\hfill \begin{array}{cc}& (\lambda +im)\mathit{u}+2{\mathit{e}}_{\mathit{z}}\times \mathit{u}=-\mathbf{\nabla }p+E\mathbf{\Delta }u,\hfill \\ \hfill & \mathbf{\nabla }·\mathit{u}=0.\hfill \\ \hfill \end{array}\end{array}$$(18)

We assumed velocity and pressure perturbations of the form

$$\begin{array}{c}\hfill \mathit{u}=\mathit{u}(r,\theta )e^{\lambda t+im\varphi },p=p(r,\theta )e^{\lambda t+im\varphi },\end{array}$$(19)

where λ = τ + iω is the complex eigenvalue. We also introduced the Ekman number,

$$\begin{array}{c}\hfill E=\frac{\nu }{\mathrm{\Omega }{R}^2},\end{array}$$(20)

as a measure of the kinematic viscosity ν.

At the inner boundary r = η, stress-free conditions are imposed, namely,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill u_r(\eta )={r\frac{\partial }{\partial r}\frac{u_{\theta }}{r}|}_{r=\eta }={r\frac{\partial }{\partial r}\frac{u_{\varphi }}{r}|}_{r=\eta }=0.\end{array}\end{array}$$(21)

On the outer free surface, the kinematic boundary condition implies

$$\begin{array}{c}\hfill u_r(1)=(\lambda +im)\zeta ,\end{array}$$(22)

where ζ is the dimensionless elevation of the surface. Since the surface experiences no stress, the dynamical boundary condition imposes

$$\begin{array}{cc}\hfill p(1)-2E{\left(\frac{\partial u_r}{\partial r}\right)}_{r=1}& =\gamma \zeta ,\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill c_{r\theta }{|}_{r=1}& =0,\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill c_{r\varphi }{|}_{r=1}& =0.\hfill \end{array}$$(25)

Here, [c] is the shear tensor with c_ij = ∂_iu_j + ∂_ju_i and we define γ by (17). We can eliminate the surface elevation, ζ, between the kinematic condition (22) and the pressure condition (23). In doing so, we note that as γ → ∞, the combined boundary condition becomes equivalent to u_r = 0, i.e., a rigid boundary, as expected. We also note that the parameters Γ and γ are related by

$$\begin{array}{c}\hfill \mathrm{\Gamma }=(1-\eta )\gamma ,\end{array}$$(26)

since H₀ = (1 − η)R.

3.2. Numerical method

We solved system (18), together with boundary conditions (21-25) numerically using spectral methods. We expanded the pressure field in the usual spherical harmonics as

$$\begin{array}{c}\hfill p(r,\theta ,\phi )={\displaystyle \sum _{\ell }}{p}_m^{\ell }(r)Y_{\ell }^m,\end{array}$$(27)

while we expanded the velocity field on vectorial spherical harmonics as

$$\begin{array}{c}\hfill \mathit{u}(r,\theta ,\phi )={\displaystyle \sum _{\ell =0}^{L_{max}}}{\displaystyle \sum _{m=-l}^{+l}}{u}_m^{\ell }(r){\mathit{R}}_{\ell }^m+v_m^{\ell }(r){\mathit{S}}_{\ell }^m+w_m^{\ell }(r){\mathit{T}}_{\ell }^m,\end{array}$$

with

$$\begin{array}{c}\hfill {\mathit{R}}_{\ell }^m=Y_{\ell }^m(\theta ,\phi ){\mathit{e}}_r,{\mathit{S}}_{\ell }^m=\mathbf{\nabla }{Y}_{\ell }^m,{\mathit{T}}_{\ell }^m=\mathbf{\nabla }\times {\mathit{R}}_{\ell }^m,\end{array}$$

where gradients are taken on the unit sphere. Here, L_max denotes the truncation order of the spherical harmonics expansion.

We then sampled the radial functions p_m^ℓ, u_m^ℓ, v_m^ℓ, and w_m^ℓ on N_r grid points of the Gauss-Lobatto collocation grid associated with Chebyshev polynomials. We determined eigenvalues λ either with the QZ-method for a global calculation or with the incomplete Arnoldi-Chebyshev method for computing specific eigenvalues (e.g. Valdettaro et al. 2007).

3.3. Properties of Kelvin waves in a thick layer

Kelvin waves do not exist without rotation, but they share several properties with surface gravity waves, which they actually are. They are non-dispersive waves in the shallow water limit (see Eq. (11)) and propagate at velocity $c=\sqrt{g{H}_0}$, just like surface gravity waves in this limit.

To identify Kelvin waves in the thick spherical shell, we tracked their eigenvalues in the complex plane as the core radius of the shell, η, decreases. Starting at η ≲ 1, we easily identify the equatorial Kelvin wave because of its shallow water frequency. Figure 2 shows that the Kelvin wave can be continuously traced as the core size decreases until the core disappears. In this limit, it joins the spectrum of eigenmodes of the full sphere, which we can derive analytically (e.g. Appendix C and Bryan 1889). In this case, the dispersion relation of the Kelvin waves reads $\omega =1-\sqrt{1+m\mathrm{\Gamma }}$ (note that in the full sphere case Γ = γ).

Fig. 2. Evolution of the eigenfrequency of the m = 1 Kelvin mode with the size of the core at γ = 1. In the shallow water and full sphere limits, the value given by their corresponding analytic expression, (11) and (C.10) are recovered, respectively.

In a thick layer, Kelvin waves also behave similarly to surface gravity waves since, for instance, their amplitude decreases faster with depth the shorter their wavelength, as illustrated in Fig. 3. The m = 10 Kelvin mode (Fig. 3 top) shows a smooth, nodeless structure, as in normal surface gravity waves. However, this Kelvin mode lies outside the inertial frequency band in the corotating frame (i.e. |ω|> 2). This is not the case for the m = 3 Kelvin mode (Fig. 3 bottom), which shows features of inertial modes, such as the emission of a shear layer by the critical latitude singularity on the inner boundary (see for instance He et al. 2022). The excitation of this shear layer increases the damping rate of the mode. In the present example, we note that the damping rate of the m = 3 Kelvin mode is larger than that of the m = 10 equivalent (see caption of Fig. 3).

Fig. 3. Top: Distribution of kinetic energy on a log10 scale in a meridional section of the spherical shell for the m = 10 Kelvin mode at λ = −1.31 × 10−5 − 12.31i. Bottom: Same as the top panel but for the m = 3 Kelvin mode at λ = −5.13 × 10−4  −  3.99i. In both cases η = 0.5, γ = 1, E = 10−7, Nr = 200, and Lmax = 200.

As with surface gravity waves, Kelvin waves become dispersive as the shell thickens. Figure 4 illustrates this, where we clearly see the increase of phase velocity with increasing thickness of the layer. When the core radius vanishes, that is, when the sphere is full, an analytic expression of the Kelvin modes frequency can be derived (see Appendix C). Figure 4 also shows the continuous relation between the shallow water (η ≃ 1) frequency and the full sphere (η = 0) frequency. We note that the m = 10 Kelvin mode reaches its full-sphere frequency even with a core radius η = 0.7. This is explained by the fact that for such a high m, the eigenfunction has its amplitude mainly close to the upper boundary and almost does not ‘feel’ the presence of the core.

Fig. 4. Evolution of the phase velocity of Kelvin modes with various m in the co-rotating frame as a function of core size η for γ = 1, E = 10−3, and Nr = Lmax = 20. The dashed red line shows the phase velocity of the m = 10 Kelvin mode in the full-sphere case.

Another view of the dispersive effect of the finite thickness of the fluid layer is the frequency difference between consecutive Kelvin modes. In the shallow water case this is a constant quantity that can be used to identify a series of Kelvin modes in photometric data (e.g. Monnier et al. 2010). In Fig. 5 we show how this frequency difference evolves with the size of the core and that it is not the same for low azimuthal wave numbers (i.e. when m ≲ 10).

Fig. 5. Frequency difference Δω = ωm + 1 − ωm between two consecutive Kelvin modes as a function of core radius. Red dots show the values for the full sphere case as derived from (C.10), while the dashed black line shows the asymptotic shallow-water case. Parameters are γ = 1, E = 10−3 with Nr = Lmax = 20.

As expressed in (12), within the limits of the shallow water approximation, the eigenfunction of Kelvin modes is a Gaussian centred on the equator. When the layer thickens, the surface shape of Kelvin modes gradually transits towards the sin^mθ-shape, which characterises the surface eigenfunction of the Kelvin mode in the full sphere. We illustrate this in Fig. 6. In the shallow water regime, the meridian profile of the eigenfunction does not depend on the azimuthal wave number m and is only a function of the parameter Γ. When the spherical shell thickens, the latitudinal extension of the Kelvin wave depends on m and somehow weakens its dependence on Γ. When the core disappears, the latitude dependence is purely that of the associate Legendre polynomial P_m^m(cos θ), that is, in sin^mθ.

Fig. 6. Radial velocity profile at the spherical-shell surface shown as a function of colatitude θ for the m = 1-Kelvin wave at various core sizes. The parameters are γ = 1, E = 10−3, Nr = 30, and Lmax = 60.

4. Influence of differential rotation

In the previous section we showed that Kelvin waves remain identifiable in a thick layer with a latitudinal spread that increases with the layer thickness similarly to the dispersion in phase velocity. We next investigated whether these surface waves could be destabilised if the background flow differs from a solid-body rotation. We aimed to determine whether the ubiquitous differential rotation of stellar envelopes could destabilise Kelvin waves. We neglected Yanai waves because we expected them to be more stable than Kelvin waves. Indeed, sampling their eigenvalue shows that the damping rate of Yanai modes is generally larger than that of Kelvin modes. Yanai modes are anti-symmetric with respect to the equator; for similar m, their global wave number is higher than that of Kelvin modes, making them presumably more damped. This is consistent with other examples in fluid mechanics where anti-symmetric modes are more stable than symmetric ones (e.g. in thermal convection, Zhang & Busse 1987; Chandrasekhar 1961).

Regarding differential rotation in early-type stars, Espinosa Lara & Rieutord (2013) showed that it results from radiative envelope baroclinicity. The resulting shear is latitudinal and radial, but the latter is usually stronger. For simplicity, we restricted our modelling of differential rotation to a shellular profile, thus depending only on the radial coordinate. We next investigated the stability conditions of Kelvin modes over this type of flows.

4.1. Equations of motion

Perturbations evolving over a shellular differential rotation satisfy the following momentum and mass conservation non-dimensional equations, written in an inertial frame:

$$\begin{array}{c}\hfill \begin{array}{cc}& (\lambda +im\mathrm{\Omega }(r))\mathit{u}+2\mathrm{\Omega }(r){\mathit{e}}_z\times \mathit{u}\hfill \\ \hfill & +r{u}_r{\partial }_r\mathrm{\Omega }sin\theta {\mathit{e}}_{\varphi }=-\mathbf{\nabla }p+E\mathbf{\Delta }\mathit{u},\hfill \\ \hfill & \mathbf{\nabla }·\mathit{u}=0,\hfill \\ \hfill \end{array}\end{array}$$(28)

where we scaled the angular velocity by its surface value. We used the outer radius of the shell as the length scale. We assumed velocity and pressure perturbations to be proportional to exp(λt + imφ).

We chose the background shellular rotation to be

$$\begin{array}{c}\hfill \mathrm{\Omega }(r)=1+({\mathrm{\Omega }}_{\eta }-1){\left(\frac{1-r}{1-\eta }\right)}^2,\end{array}$$(29)

which satisfies Ω(1) = 1, and we introduce the new parameter Ω_η = Ω(η), namely the rotation rate at the core boundary. This parameter controls the strength of the differential rotation.

Profile (29) also verifies ∂_rΩ(r = 1) = 0, so that it exerts no viscous stress at the surface, as required. This profile thus combines the simplicity of a polynomial radial dependence and the right behaviour near the surface. To mimic envelope differential rotation in more realistic models (e.g. Espinosa Lara & Rieutord 2013), we considered Ω_η > 1, hence decreasing rotation rates with radius.

The strength of the differential rotation Ω_η should not be too high, as we did not want our set-up to be unstable with respect to centrifugal (Taylor-Couette) instability (e.g. Drazin & Reid 1981). This instability involves axisymmetric perturbations and redistributes the fluid’s angular momentum. Since it develops on a dynamical time scale, it is unlikely to be active in a star. In fact, we assume that profile (29) represents a large-scale differential rotation as driven by baroclinicity and possibly including small-scale turbulence generated by shear instabilities. In Appendix D we show that if

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\eta }\le 1+8{(1-\eta )}^2,\end{array}$$(30)

the differential rotation is stable with respect to the centrifugal instability. We therefore limited the strength of the differential rotation with (30).

Concerning the boundary conditions for the perturbations, we kept the same ones as in the uniform rotation case, but modified condition (25) to account for the variations of Ω, such that

$$\begin{array}{c}\hfill c_{r\varphi }(1)+sin\theta \mathrm{\Omega }″(1)\zeta =0,\end{array}$$(31)

where ζ is the radial elevation of the surface as introduced in (3).

4.2. Destabilisation of Kelvin modes

As demonstrated in Appendix B, the gravito-inertial modes of our system are stable in the case of solid-body rotation. As illustrated in Fig. 7, this is no longer the case when a shellular differential rotation is present. We plot the growth rate τ = ℜe(λ) of three Kelvin modes as a function of the Ekman number E (the non-dimensional kinematic viscosity) or the strength of the differential rotation Ω_η. Figure 7 (top) clearly shows that for a given differential rotation, there is a critical Ekman number below which a Kelvin mode becomes unstable, and the higher the azimuthal wave number m, the lower the critical Ekman number (as expected). Conversely, for a given Ekman number, there is a critical Ω_η beyond which a Kelvin mode becomes unstable.

Fig. 7. Top: Growth rate τ = ℜe(λ) of selected Kelvin modes as a function of the Ekman number E, for Ωη = 2, Γ = 0.01, and η = 0.18. Bottom: Same as the top panel but as a function of differential rotation parametrized by Ωη, with E = 10−5. The numerical resolution is Nr = Lmax = 100.

4.3. Origin of the instability

The growth rate of the modes of our system can be expressed using the momentum equation in (28). To do so, we took the dot product of this equation with the complex conjugate of the velocity field, u^*, and integrated the equation over the fluid volume. After some rearrangement of the terms, we obtain:

$$\begin{array}{cc}\hfill A_0\mathfrak{R}e(\lambda )& =\stackrel{\mathrm{I}}{\overbrace{-E{\displaystyle {\int }_{(V)}}{c}_{\mathit{ij}}{c}_{\mathit{ij}}^{\ast }dV}}+\stackrel{\mathrm{II}}{\overbrace{E{\displaystyle {\int }_{(S)}}\mathfrak{R}e(u_{\phi }^{\ast }{c}_{r\phi })dS}}\hfill \\ \hfill & \stackrel{\mathrm{III}}{\overbrace{-{\displaystyle {\int }_{(V)}}rsin\theta {\mathrm{\Omega }}^{\prime }\mathfrak{R}e(u_{\phi }^{\ast }{u}_r)dV}},\hfill \end{array}$$(32)

where

$$\begin{array}{c}\hfill A_0={\displaystyle ({\int }_V{|u|}^{2}dV+\frac{1}{\gamma }{\int }_{{\partial V|}_{r=1}}{|p-2E{u}_r^{\prime }|}^{2}dS)}>0\end{array}$$(33)

is a positive-definite term. The first integral is the bulk viscous dissipation, which is always negative, as expected (we used Einstein implicit summation on repeated indices). The second term is a surface integral with no obvious sign. Using the surface boundary condition (31), we note that it relates to the differential rotation via Ω″(1) and to the phase difference between u_r and u_φ through the kinematic boundary condition (22). We numerically investigated its role in the instability and found that it is always negative, thus playing a stabilising role, although we could not prove that this is always the case. The third term is responsible for the instability of Kelvin modes. Because Ω′(r) < 0 over the volume, modes can become unstable when ℜe(u_ru_φ^*) is positive or, equivalently, when the phase difference between u_r and u_φ is less than π/2 over most of the volume.

To better characterise this instability, we investigated the dependence of τ = ℜe(λ) as a function of the differential rotation parameter Ω_η for various Ekman numbers. We recall that as Ω_η increases, the core rotation rate increases as the differential rotation. Figure 8 shows the variation of τ with Ω_η for three Ekman numbers. The shape of the curve remains similar for the Ekman numbers considered. We clearly see that for sufficiently low E, τ is positive when Ω_η is large enough; however, remarkably, τ becomes negative again when Ω_η is too large. Since τ is the sum of three integrals, its behaviour can be explained by the dependence of these integrals on Ω_η. Figure 9 shows the variations of the bulk and surface viscous integrals with Ω_η. The viscous contribution varies monotonically with Ω_η, the damping effect being more important when E and Ω_η increase. Figure 10 shows that the variations of τ come from the coupling integral, which first increases with Ω_η and then decreases if Ω_η ≳ 5. Figures 11 and 12 summarise this behaviour of the growth and damping rate. There, a Kelvin wave becomes unstable below a critical Ekman number, but this is conditioned by the strength of the differential rotation, which must satisfy Ω_m < Ω_η < Ω_M, where the bounds Ω_m and Ω_M depend on the Ekman number.

Fig. 8. Evolution of τ, the real part of the eigenvalue, for the m = 5 Kelvin mode as a function of Ωη, for three Ekman numbers. The parameters are γ = 2, Nr = 100, and Lmax = 60.

Fig. 9. Evolution of the sum of the viscous integrals (I) and (II) in (32) as a function of Ωη for the m = 5 Kelvin mode for three Ekman numbers. The parameters are γ = 2, Nr = 100, and Lmax = 60.

Fig. 10. Evolution of the coupling integral (III) in equation (32) as a function of Ωη for the m = 5 Kelvin mode for three Ekman numbers. The parameters are γ = 2, Nr = 100, and Lmax = 60.

Fig. 11. Evolution of the eigenvalue of the m = 5 Kelvin mode in the complex plane. The parameters are γ = 2, η = 0.18, E = 1 × 10−4, Nr = 100, and Lmax = 100.

Fig. 12. Stability diagram of the m = 5 Kelvin mode at η = 0.18 and γ = 2. The purple region denotes the unstable domain in the explored parameter space and forms a single connected region. The instability does not exist beyond a critical Ekman number E ≈ 5 × 10−4.

The foregoing situation resembles shear flow instabilities where a critical layer plays a crucial role. The critical layer is the location where the phase speed of the wave matches that of the background flow. In our case, the critical layer lies on a sphere of radius r_c such that

$$\begin{array}{c}\hfill \mathrm{\Omega }(r_c)=-\frac{\omega }{m}.\end{array}$$(34)

Figure 13 illustrates the wave and critical layer interactions. We first note that as differential rotation strengthens, the critical layer moves towards the surface. Indeed, as shown in Fig. 11, the wave frequency does not vary much when Ω_η increases by a factor of ∼4, hence condition (34) is satisfied with an approximately constant Ω(r_c), requiring a larger r_c as Ω_η increases (since Ω is a decreasing function of r).

Fig. 13. Radial profiles of the coupling term in the equatorial plane θ = π/2, at longitude ϕ = 0, for an m = 5-Kelvin mode and for three differential rotations. The parameters are γ = 2, E = 1 × 10−5, η = 0.18, Nr = 100, and Lmax = 100. Inertial frame eigenvalues are λ = −2.968 × 10−3 − 7.372i for Ωη = 2 implying rcri = 0.435; λ = 1.231 × 10−2 − 7.591i for Ωη = 5 implying and rcri = 0.704. λ = −4.083 × 10−2 − 7.929i, for Ωη = 8 implying rcri = 0.762.

Figure 13 also shows that at high differential rotation (here at Ω_η = 8), the internal part of the Kelvin wave is trapped around the critical layer. This part of the perturbation displays a shear layer, whose width scales as E^0.4 in this example. This shear layer is reminiscent of what occurs in plane-parallel shear flows: in an inviscid plane-parallel shear flow, the critical layer is the location of a singularity in the perturbations, which is regularised by viscosity into a thin shear layer (Drazin & Reid 1981; Charru 2011). Numerical exploration of Kelvin waves shows that, at least for some parameters, their instability can also disappear when the Ekman number is sufficiently low, as illustrated in Fig. 14. This property can be understood qualitatively from the scaling of the shear-layer width, which makes the coupling integral (III in Eq. (32)) decay more rapidly with E than the viscous integrals. Here we encounter the double role of a critical layer, which can be either destabilising or stabilising, as shown in other contexts (e.g. Riedinger & Gilbert 2014, in ocean dynamics). We did not analyse this instability as it is beyond the scope of this paper, and this constant-density model remains far from the physical conditions met in stars.

Fig. 14. Evolution of the growth rate τ of the m = 1 Kelvin mode as a function of Ekman number E. The frequency is ω ≃ −1.030303. The parameters are η = 0.35, γ−1 = 16.25, Ωη = 1.065 Nr = 200, and Lmax = 200.

5. Conclusions

We explored the stability of equatorial Kelvin waves propagating over an incompressible, differentially rotating fluid layer contained in a spherical shell. These waves have presumably been detected in the fast rotating star Rasalhague (Monnier et al. 2010) and could also represent a mechanism by which Be stars eject matter in their equatorial plane. Since equatorial Kelvin waves are also special waves of the stellar wave zoo, owing to their topological origin (Delplace et al. 2017), their linear stability is of interest. We did not consider Yanai waves, the anti-symmetric counterpart of Kevin waves, since they are expected to be more stable, as suggested by the cases we computed. Our simplified model shows that even if Kelvin waves have been studied in the shallow water framework (for terrestrial applications), they still exist when the fluid layer is no longer a thin spherical shell. The key difference is that their equatorial confinement is no longer well pronounced unless the azimuthal wave number is large. We also note that at low azimuthal wave number, their frequency lies in the inertial frequency band. In this case, these waves exhibit features of inertial waves, such as shear layers emitted by the critical-latitude singularity at the inner core boundary. Such shear layers provide an additional source of viscous dissipation.

With this model, we also show that Kelvin waves persist in the presence of a shellular (radial) differential rotation. This differential rotation mimics the shear imposed by the baroclinicity in a stellar radiative envelope. Moreover, we show that Kelvin waves become unstable if the viscosity is not too high and differential rotation is strong enough. However, this instability disappears when either the differential rotation is too strong or the Ekman number is very low. We find that the critical layer associated with each Kelvin wave presumably plays a destabilising or stabilising role. As in plane-parallel shear flows, viscosity is a key parameter. The shear layers that appear when Kelvin waves lie in the inertial frequency band are not sufficiently dissipative to prevent the rise of the instability.

The foregoing results suggest that the instability of Kelvin waves will persist when we relax the simplification of a constant density fluid. In this case, surface gravity waves are replaced by the so-called f-modes, which are also topological waves at low wave numbers (Le Saux et al. 2025). The investigation of these waves with a compressible fluid, which more realistically represents a stellar envelope, is a natural follow-up of this work.

Acknowledgments

We would like to thank the referee very much for constructive remarks, which helped us improve the original manuscript. We are also very grateful to Lorenzo Valdettaro for his help in an early phase of the project and to Armand Leclerc for enlightening discussions on topological waves. The research leading to these results has received funding from the European Research Council (ERC) under the Horizon Europe programme (Synergy Grant agreement N°101071505: 4D-STAR). While partially funded by the European Union, views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. Computations have been possible thanks to HPC resources from CALMIP supercomputing centre (Grant 2025-P0107).

References

Appendix A: Tidal Laplace Equation and Poincaré Modes

As in the classical β-plane approximation, it is possible to combine the full set of shallow water equations (4) into a single equation, known as the Tidal Laplace Equation (e.g. Longuet-Higgins 1968):

$$\begin{array}{c}\hfill \begin{array}{cc}& \frac{\partial }{\partial \mu }\left(\frac{1-{\mu }^2}{1-4q^2{\mu }^2}\frac{\partial \zeta }{\partial \mu }\right)-\frac{2mq(1+q^2{\mu }^2)}{{(1-4q^2{\mu }^2)}^2}\zeta \hfill \\ \hfill & -\frac{m^2}{(1-{\mu }^2)(1-4q^2{\mu }^2)}\zeta =-\frac{1}{q^2\mathrm{\Gamma }}\zeta \hfill \end{array}\end{array}$$(A.1)

where we have introduced the spin parameter q = 1/ω. In the limit of large frequency (q → 0), this equation simplifies and transforms to the classical spherical Laplacian eigenvalue problem:

$$\begin{array}{c}\hfill \frac{\partial }{\partial \mu }\left((1-{\mu }^2)\frac{\partial \zeta }{\partial \mu }\right)-\frac{m^2}{1-{\mu }^2}\zeta =-(\frac{1}{q^2\mathrm{\Gamma }}-2mq)\zeta \end{array}$$(A.2)

which leads to eigenvalues solutions of:

$$\begin{array}{c}\hfill \frac{{\omega }^2}{\mathrm{\Gamma }}-\frac{2m}{\omega }=\ell (\ell +1)\end{array}$$(A.3)

This cubic equation in ω can be solved perturbatively in the limit of large frequencies by introducing a small correction to the zeroth-order solution, which reads

$$\begin{array}{c}\hfill {\omega }_0=\pm \sqrt{\mathrm{\Gamma }\ell (\ell +1)}.\end{array}$$(A.4)

Introducing the perturbation δω ≪ ω₀ and expanding each term to first order:

$$\begin{array}{cc}\hfill {\omega }^2& ={\omega }_0^2+2{\omega }_0\delta \omega +\mathcal{O}(\delta {\omega }^2)\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill \frac{1}{\omega }& =\frac{1}{{\omega }_0}-\frac{\delta \omega }{{\omega }_0^2}+\mathcal{O}(\delta {\omega }^2)\hfill \end{array}$$(A.6)

we finally get:

$$\begin{array}{c}\hfill \omega =\pm \sqrt{\mathrm{\Gamma }\ell (\ell +1)}+\frac{m}{\ell (\ell +1)}+\mathcal{O}({\ell }^{-4})\end{array}$$(A.7)

This expression is the first-order-corrected Poincaré mode frequency in the limit of large ℓ.

Appendix B: The stability of the uniformly rotating viscous spherical shell

Ignoring self-gravity, perturbations of a uniformly rotating fluid in a spherical shell verify (18) together with boundary conditions (22) and (23). Taking the dot product of the momentum equation with u^* (the complex conjugate of u), integrating over the volume of the spherical shell and taking the real part of the equation gives:

$$\begin{array}{cc}& \mathfrak{R}e(\lambda ){\displaystyle {\int }_{(V)}}{|\mathit{u}|}^{2}dV\hfill \\ \hfill & =-\mathfrak{R}e{\displaystyle {\int }_{(S)}}(p-2E{u}_r^{\prime })u_r^{\ast }dS-\frac{E}{2}{\displaystyle {\int }_{(V)}}{|c_{\mathit{ij}}|}^{2}dV\hfill \end{array}$$(B.1)

where c_ij are the components of the shear tensor. Combining boundary conditions (22) and (23) gives:

$$\begin{array}{c}\hfill \gamma u_r=(\lambda +im)(p-2E{u}_r^{\prime })\mathrm{at}r=1\end{array}$$(B.2)

where ′ indicates the radial derivative. This allows us to rewrite (B.1) as

$$\begin{array}{cc}\hfill \mathfrak{R}e(\lambda )[{\displaystyle {\int }_{(V)}}{|\mathit{u}|}^{2}dV+\frac{1}{\gamma }{\displaystyle {\int }_{(S)}}{|p-2E{u}_r^{\prime }|}^{2}dS]& \hfill \\ \hfill =-\frac{E}{2}{\displaystyle {\int }_{(V)}}{|c_{\mathit{ij}}|}^{2}dV\end{array}$$(B.3)

which shows that ℜe(λ) < 0 when the Ekman number E is non-zero, hence when the fluid is viscous.

Appendix C: The full sphere case

Starting from boundary conditions (22) and (23) in the inviscid limit and in the corotating frame, one obtains

$$\begin{array}{c}\hfill \mathit{v}·\mathit{n}=i\omega \frac{p}{\gamma }\end{array}$$(C.1)

in cylindrical coordinates (s, z, φ) on the unit sphere (r = 1),

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathit{n}=s{\mathit{e}}_{\mathit{s}}+z{\mathit{e}}_{\mathit{z}}\hfill \\ \hfill & \text{with,}\hfill \\ \hfill & v_s=-\frac{1}{4-{\omega }^2}(i\omega \frac{\partial P}{\partial s}+\frac{2}{s}\frac{\partial P}{\partial \phi })\hfill & \hfill v_z=-\frac{1}{i\omega }\frac{\partial P}{\partial z}\end{array}\end{array}$$(C.2)

Equation (C.1) can thus be rewritten as

$$\begin{array}{c}\hfill s\frac{\partial P}{\partial s}+(\frac{2m}{\omega }+\frac{4-{\omega }^2}{\gamma })P-\frac{4-{\omega }^2}{{\omega }^2}z\frac{\partial P}{\partial z}=0\end{array}$$(C.3)

By introducing the same change of variables as in Rieutord (2015), one finally obtains the following dispersion relation:

$$\begin{array}{c}\hfill [2m+\frac{(4-{\omega }^2)\omega }{\gamma }]P_l^m\left(\frac{\omega }{2}\right)=\frac{4-{\omega }^2}{2}{P_l^m}^{\prime }\left(\frac{\omega }{2}\right).\end{array}$$(C.4)

C.1. Kelvin mode

We consider the sectoral mode m = ℓ > 0. From the definition of the associated Legendre polynomials,

$$\begin{array}{c}\hfill P_{\ell }^m={(-1)}^m{(1-{\omega }^2)}^{m/2}\frac{d^m{P}_{\ell }(\omega )}{d{\omega }^m},\end{array}$$(C.5)

we deduce,

$$\begin{array}{c}\hfill \frac{d{P}_m^m(\omega )}{d\omega }={(-1)}^m(-\omega m){(1-{\omega }^2)}^{m/2-1}\frac{d^m{P}_m(\omega )}{d{\omega }^m}\end{array}$$(C.6)

So, (C.4) yields for m = ℓ,

$$\begin{array}{c}\hfill (2m+\omega \frac{4-{\omega }^2}{\gamma })=-\omega m.\end{array}$$(C.7)

Noting that ω = −2 is a solution, we rewrite (C.7) as

$$\begin{array}{c}\hfill (\omega +2)({\omega }^2-2\omega -m\gamma )=0.\end{array}$$(C.8)

The two other roots are

$$\begin{array}{c}\hfill {\omega }_{1,2}=1\pm \sqrt{1+m\gamma },\end{array}$$(C.9)

and we identify the Kelvin mode as the prograde mode, ω < 0, namely

$$\begin{array}{c}\hfill \omega =1-\sqrt{1+m\gamma }\end{array}$$(C.10)

The other root is that of the first retrograde Poincaré mode.

Appendix D: Centrifugal stability

The flow is stable with respect to centrifugal instability if and only if the specific angular momentum L = s²Ω increases with the distance s to the rotation axis, namely if

$$\begin{array}{c}\hfill {\partial }_{s}L>0,\end{array}$$(D.1)

with s = r sin θ. Since

$$\begin{array}{c}\hfill \mathrm{\Omega }(r)=1+K{(1-r)}^2\mathrm{with}K=\frac{{\mathrm{\Omega }}_{\eta }-1}{{(1-\eta )}^2},\end{array}$$(D.2)

Condition (D.1) implies

$$\begin{array}{c}\hfill 1+K(1-r)(1-r(1+{sin}^2\theta ))>0\forall r,\theta \end{array}$$(D.3)

which is true if this second order equation has no root for r or that

$$\begin{array}{c}\hfill 4+\frac{4}{K}>\frac{{(2+{sin}^2\theta )}^2}{1+{sin}^2\theta }\end{array}$$(D.4)

for all θ. Since the RHS is a monotonic increasing function of sin²θ, it reaches a maximum at θ = π/2. Hence, inequality (D.1) is always satisfied if $4+\frac{4}{K}>9/2$, which implies (30).

All Figures

	Fig. 1. Dispersion relation of gravito-inertial modes in the shallow water system for Γ = 0.01. Left: Eigenfrequency of modes symmetric with respect to the equator as a function of the azimuthal wave number m. Right : Same as the left panel but for antisymmetric modes. Red dots show Poincaré modes, blue dots show Rossby modes, and green dots show Kelvin (left) and Yanai (right) modes.
In the text

	Fig. 2. Evolution of the eigenfrequency of the m = 1 Kelvin mode with the size of the core at γ = 1. In the shallow water and full sphere limits, the value given by their corresponding analytic expression, (11) and (C.10) are recovered, respectively.
In the text

	Fig. 3. Top: Distribution of kinetic energy on a log10 scale in a meridional section of the spherical shell for the m = 10 Kelvin mode at λ = −1.31 × 10−5 − 12.31i. Bottom: Same as the top panel but for the m = 3 Kelvin mode at λ = −5.13 × 10−4  −  3.99i. In both cases η = 0.5, γ = 1, E = 10−7, Nr = 200, and Lmax = 200.
In the text

	Fig. 4. Evolution of the phase velocity of Kelvin modes with various m in the co-rotating frame as a function of core size η for γ = 1, E = 10−3, and Nr = Lmax = 20. The dashed red line shows the phase velocity of the m = 10 Kelvin mode in the full-sphere case.
In the text

	Fig. 5. Frequency difference Δω = ωm + 1 − ωm between two consecutive Kelvin modes as a function of core radius. Red dots show the values for the full sphere case as derived from (C.10), while the dashed black line shows the asymptotic shallow-water case. Parameters are γ = 1, E = 10−3 with Nr = Lmax = 20.
In the text

	Fig. 6. Radial velocity profile at the spherical-shell surface shown as a function of colatitude θ for the m = 1-Kelvin wave at various core sizes. The parameters are γ = 1, E = 10−3, Nr = 30, and Lmax = 60.
In the text

	Fig. 7. Top: Growth rate τ = ℜe(λ) of selected Kelvin modes as a function of the Ekman number E, for Ωη = 2, Γ = 0.01, and η = 0.18. Bottom: Same as the top panel but as a function of differential rotation parametrized by Ωη, with E = 10−5. The numerical resolution is Nr = Lmax = 100.
In the text

	Fig. 8. Evolution of τ, the real part of the eigenvalue, for the m = 5 Kelvin mode as a function of Ωη, for three Ekman numbers. The parameters are γ = 2, Nr = 100, and Lmax = 60.
In the text

	Fig. 9. Evolution of the sum of the viscous integrals (I) and (II) in (32) as a function of Ωη for the m = 5 Kelvin mode for three Ekman numbers. The parameters are γ = 2, Nr = 100, and Lmax = 60.
In the text

	Fig. 10. Evolution of the coupling integral (III) in equation (32) as a function of Ωη for the m = 5 Kelvin mode for three Ekman numbers. The parameters are γ = 2, Nr = 100, and Lmax = 60.
In the text

	Fig. 11. Evolution of the eigenvalue of the m = 5 Kelvin mode in the complex plane. The parameters are γ = 2, η = 0.18, E = 1 × 10−4, Nr = 100, and Lmax = 100.
In the text

	Fig. 12. Stability diagram of the m = 5 Kelvin mode at η = 0.18 and γ = 2. The purple region denotes the unstable domain in the explored parameter space and forms a single connected region. The instability does not exist beyond a critical Ekman number E ≈ 5 × 10−4.
In the text

Fig. 13. Radial profiles of the coupling term in the equatorial plane θ = π/2, at longitude ϕ = 0, for an m = 5-Kelvin mode and for three differential rotations. The parameters are γ = 2, E = 1 × 10−5, η = 0.18, Nr = 100, and Lmax = 100. Inertial frame eigenvalues are λ = −2.968 × 10−3 − 7.372i for Ωη = 2 implying rcri = 0.435; λ = 1.231 × 10−2 − 7.591i for Ωη = 5 implying and rcri = 0.704. λ = −4.083 × 10−2 − 7.929i, for Ωη = 8 implying rcri = 0.762.

In the text

	Fig. 14. Evolution of the growth rate τ of the m = 1 Kelvin mode as a function of Ekman number E. The frequency is ω ≃ −1.030303. The parameters are η = 0.35, γ−1 = 16.25, Ωη = 1.065 Nr = 200, and Lmax = 200.
In the text