# On the stability of viscous Riemann ellipsoids

Joris Labarbe

► To cite this version:

Joris Labarbe. On the stability of viscous Riemann ellipsoids. 2026. hal-05495396v2

HAL Id: hal-05495396

https://hal.science/hal-05495396v2

Preprint submitted on 5 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
ON THE STABILITY OF VISCOUS RIEMANN ELLIPSOIDS

JORIS LABARBE

Université Côte d'Azur, CNRS UMR 7351, Laboratoire J. A. Dieudonné, Parc Valrose,
06108 Nice cedex 2, France

ABSTRACT. The present study investigates the linear stability of Riemann ellipsoids in both the inviscid limit and in the presence of weak viscosity. In the inviscid regime, we derive a generalised Poincaré equation governing small fluid oscillations and construct a family of polynomial solutions that extends the classical results of Cartan to flows with a uniform strain field. This formulation provides an analytic dispersion relation for three-dimensional ellipsoidal disturbances and remains computationally efficient at arbitrary harmonic degree, in contrast to the virial tensor method or to short-wavelength (WKB) approximations. The viscous effects are incorporated through a boundary-layer analysis based on Prandtl's theory, leading to first-order viscous corrections to the inviscid spectrum and allowing a systematic investigation of viscosity-driven instabilities. Stability diagrams are presented over the space of admissible Riemann ellipsoids, illustrating the roles of rotation, internal strain, and diffusion, with implications for rotating shear flows in geophysical and astrophysical contexts.

# 1. INTRODUCTION

Celestial mechanics and the laws of planetary motion rank among the oldest and most
prominent areas of research in classical physics. Historically, the formulation of Newtonian
mechanics together with the universal law of gravitation provided a unified and quantitative
framework within which such problems could be addressed rigorously. Since then, sustained
interest has developed in the study of the equilibrium and stability of self-gravitating fluids,
motivated both by fundamental theoretical considerations and by applications in astrophysics
and geophysics. In his treatise of fluxions, Maclaurin (1742) identified the only stable ax-
isymmetric family of equilibrium figures corresponding to rigidly rotating, homogeneous fluid
bodies in a gravitational field. These configurations, now known as Maclaurin spheroids,
were subsequently investigated by numerous authors as simplified models of isolated rotating
stars (or planets) subject to the Coriolis effect and centrifugal acceleration. Over a century
later, Meyer (1842) and Liouville (1851), employing distinct analytical approaches, demon-
strated that the Maclaurin sequence undergoes a bifurcation at a finite eccentricity towards
a non-axisymmetric equilibrium, namely the Jacobi ellipsoid (Jacobi, 1839). These results
were placed on a rigorous footing by Poincaré (1885), who interpreted the departure from the
initial state within the modern framework of dynamical systems theory and bifurcation anal-
ysis, thereby identifying one of the earliest paradigmatic example of spontaneous symmetry
breaking in continuum mechanics. In the same monograph, Poincaré derived a differential
equation governing the small oscillations of a rotating, incompressible, self-gravitating, invis-
cid, and strainless mass of fluid. This formulation stimulated the work of Bryan (1889), who

<footnote>E-mail address: joris.LABARBE@univ-cotedazur.fr.</footnote>
analysed the normal modes of Maclaurin spheroids through solutions of what is now referred
to as the Poincaré equation. Several decades later, Cartan (1922) revisited the problem and
rigorously identified the unique family of polynomial solutions satisfying this equation. In the
absence of self-gravitational effects, Greenspan (1965) showed that the associated Poincaré
operator is self-adjoint and possesses a bounded real spectrum. He further raised the question
of the completeness of the corresponding inertial modes, a fundamental issue for the modal
decomposition of rotating flows, which was resolved affirmatively only much later for arbitrary
bounded ellipsoidal domains (Backus & Rieutord, 2017).

As shown by Dedekind (1861) in a seminal contribution, to every equilibrium figure whose
velocity field is generated by a linear mapping of the position vector there exists an adjoint
configuration whose velocity field is determined by the transpose of this linear map. Re-
markably, the original and adjoint configurations possess the same ellipsoidal shape and both
satisfy the conditions of gravitational equilibrium. Owing to the linear dependence of the ve-
locity field on the spatial coordinates, these flows are characterised, in the rotating frame, by
a uniform rate of strain and a uniform vorticity. Building on Dedekind’s theorem, Riemann
(1860) provided a complete classification of all admissible equilibrium configurations with re-
spect to Dedekind’s assumptions. In particular, Riemann identified three distinct families of
non-axisymmetric ellipsoids obtained by superposing solid-body rotation with internal shear
motions of uniform vorticity. Among these, the so-called S-type ellipsoids — whose vorticity
vector is aligned with the rotation axis and therefore with the Coriolis force — are of par-
ticular relevance for physical applications. Despite their potential importance for modelling
internal differential motions in planets and stars, and their role as the simplest exact solu-
tions incorporating both rotation and strain, these equilibria were largely overlooked for many
decades. A systematic linear stability analysis of Riemann ellipsoids was eventually under-
taken by Chandrasekhar (1965, 1966), who employed the virial tensor method developed and
applied extensively in related contexts (Chandrasekhar, 1969). Nevertheless, a comprehensive
linear stability theory for Riemann ellipsoids, valid for harmonic perturbations of arbitrary
order and comparable to that established for Maclaurin spheroids (Comins, 1979), has yet to
be developed.

In the specific context of rigidly rotating spheroids, Riemann (1860) identified the onset of
a dynamical instability at a finite eccentricity, corresponding to disk-like configurations. The
associated critical threshold describes a Hamilton–Hopf bifurcation, which can be interpreted
within Krein theory as the collision of two eigenvalues with opposite Krein signatures, or
energy signs (Labarbe, 2021). Such instabilities are nowadays recognised as generic features of
non-dissipative systems with specific symmetries, and similar phenomena arise well beyond the
setting of Maclaurin spheroids, appearing in a broad class of out-of-equilibrium Hamiltonian
systems (Kirillov, 2021). On the basis of Dirichlet’s principle, Riemann (1860) further argued
that his non-homogeneous ellipsoids are stable for almost all admissible shapes and rotation
rates. This conclusion was subsequently challenged by Chandrasekhar (1965, 1966), who
demonstrated that a wide range of self-gravitating configurations are unstable with respect
to second- and third-order harmonic perturbations, thereby contradicting Riemann’s original
claim. These instability results were later corroborated by Lebovitz and Lifschitz (1996) using
a short-wavelength, geometrical-optics (WKB) approximation of the linearised equations of
motion. It is now understood that the discrepancy between the conclusions of Riemann and
those of Chandrasekhar arises from an incorrect application of Dirichlet’s principle, since the
governing equations of motion were not casted into an admissible Lagrangian form (Morrison,
Lebovitz, & Biello, 2009).
Instability is an inherent feature of inviscid fluid systems with elliptic streamlines. From a physical standpoint, such instability arises from resonant interactions between the normal modes of oscillation and an imposed planar strain field (Kerswell, 2002). Riemann ellipsoids may therefore be regarded as among the earliest examples of this phenomenon, although their connection to elliptic instability remained naturally unrecognised for many decades. Subsequent studies demonstrated that a uniform strain field can excite Kelvin modes in the form of helical waves, in axisymmetric vortices (Moore & Saffman, 1975) as well as in unbounded vortex flows (Bayly, 1986; Pierrehumbert, 1986). These works established that elliptic instability generates local enstrophy along the vortex axis through vortex stretching, thereby highlighting its fundamental three-dimensional character. Given the central role of enstrophy production in turbulence, this mechanism has also been extensively investigated in the context of shear-flow transition (Bayly, Orszag, & Herbert, 1988). When viscosity is incorporated into the stability analysis, it introduces a damping contribution that weakens growth rates without suppressing the underlying instability mechanism (Landman & Saffman, 1987). However, because dissipation alters the global energy balance, it can lead to significant modifications of the nonlinear saturation and long-term dynamics relative to the conservative case.

Energy dissipation, although often regarded as a stabilising phenomenon, can in certain circumstances act as a destabilising mechanism. This counterintuitive behaviour was first identified in the classical study of Thompson and Tait (1867), in the context of Maclaurin spheroids endowed with arbitrarily small viscosity. Nearly a century later, Roberts and Stewarton (1963) rigorously established the existence of this secular instability using asymptotic methods and boundary-layer theory. They actually showed that sectoral modes of rigidly rotating stars exhibit exponential growth once the eccentricity exceeds a critical value, which coincides with the bifurcation point from the Maclaurin to the Jacobi sequence identified by Meyer (1842) and Liouville (1851). This secular instability is now understood as a prototypical example of dissipation-induced instability and is commonly attributed to the destruction of gyroscopic stabilisation by diffusion (Bloch, Krishnaprasad, Marsden, & Ratiu, 1994; Kirillov, 2021; Labarbe, 2021). With regard to this phenomenon, and as already noted by Lebovitz and Lifschitz (1996), a systematic treatment of viscosity-induced instabilities in the context of Riemann ellipsoids remains conspicuously absent from the literature.

In a further seminal contribution, Chandrasekhar (1970) demonstrated the existence of dissipation-induced instabilities in Maclaurin spheroids arising from diffusive effects associated with gravitational radiation reaction, within a post-Newtonian framework. One of the principal implications of this result is that sufficiently massive, isolated rotating stars may emit gravitational waves as a consequence of the Chandrasekhar–Friedman–Schutz (CFS) instability Chandrasekhar (1970); Schutz and Friedman (1978a, 1978b). This mechanism has since become a promising objective in gravitational-wave astrophysics, particularly following the direct detection of gravitational waves from compact binary mergers by ground-based interferometers (Abott et al., 2016). Subsequently, Comins (1979) generalised these results to harmonic perturbations of arbitrary order and demonstrated that the secular modes destabilised by viscosity are distinct from those driven unstable by gravitational radiation reaction. To date, no comparable comprehensive analysis has been carried out for Riemann ellipsoids, beyond the inviscid investigations of Chandrasekhar (1965, 1966) and Lebovitz and Lifschitz (1996).
The present work aims to fill the aforementioned gaps by developing a unified and tractable linear stability theory for Riemann ellipsoids, valid both in the inviscid limit and in the presence of weak viscosity. In the inviscid setting, we introduce a generalised Poincaré equation governing small-amplitude oscillations about the equilibrium state and construct an associated family of polynomial solutions that encompasses, as a special case, the classical solutions obtained by Cartan. This formulation yields an explicit analytic dispersion relation for inviscid Riemann ellipsoids and allows stability properties to be computed efficiently for ellipsoidal harmonics of arbitrary order. In contrast to the virial tensor method introduced by Chandrasekhar (1969), the present approach remains computationally efficient at high harmonic degree and provides direct spectral access to the full mode structure. Moreover, unlike the short-wavelength (WKB) approximation employed by Lebovitz and Lifschitz (1996), the analysis applies uniformly across all wavelengths. In the presence of viscosity, we derive a corresponding analytic dispersion relation by means of a boundary-layer analysis based on Prandtl's theory (Landau & Lifshitz, 1987), yielding first-order viscous corrections to the inviscid spectrum and enabling a systematic investigation of viscosity-driven effects.

The structure of the article is as follows. In Section 2, we introduce the governing equations of motion together with the relevant boundary conditions, and describe the inviscid equilibrium state as well as the linearised equations about this steady solution. Section 3 presents the derivation of the equation governing the normal modes of oscillations and its general polynomial solution. In Section 4, we develop the small-viscosity approximation and obtain the first-order viscous corrections to the dispersion relation. Selected stability results across the relevant parameter space are then presented. Finally, Section 5 discusses the physical interpretation of these findings in the broader contexts of geophysics and astrophysics, highlighting their implications for rotating fluid bodies with internal motion.

## 2. FORMULATION OF THE PROBLEM

### 2.1. Equations of motion and equilibrium state.
Riemann ellipsoids (of S-type family, as exclusively restricted in this work) are equilibrium configurations of self-gravitating, homogeneous and incompressible masses of fluid. They are characterised by a combination of solid-body rotation and uniform internal vorticity, both aligned with a fixed vertical axis. The rigid rotation of the ellipsoid is described by the angular velocity $\Omega$, gravitational effects are modelled through the Newtonian potential $\Phi$, and the constant vorticity in the co-rotating frame is denoted by $\zeta$. The assumption of uniform vorticity implies that the associated velocity field depends linearly on the spatial coordinates, according to Dedekind's theorem for adjoint ellipsoidal configurations (Dedekind, 1861). A schematic representation of an arbitrary Riemann ellipsoid is depicted in Figure 1.

The fluid motion is described either in Cartesian coordinates $(x, y, z)$, with respect to an orthonormal basis $(e_x, e_y, e_z)$, or in ellipsoidal coordinates, better suited to the geometry of the equilibrium figure. The latter coordinate system, together with a selection of geometrical relations, is presented in Appendix A. Without loss of generality, we assume the semi-axes of the ellipsoid to respect $0 < a_3 \le a_2 \le a_1 < \infty$. We denote the characteristic time frequency by $(\pi\rho G)^{1/2}$ and the unit length by $a_1$, the largest semi-axis of the ellipsoid. Thus, the zero level set of the function

$$ (1) \qquad F(x, y, z) = x^2 + \frac{y^2}{\Gamma^2} + \frac{z^2}{\Xi^2} - 1, $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Sketch of a rigidly rotating and self-gravitating Riemann ellipsoid with uniform strain at its interior. In blue are depicted the elliptic streamlines.</figcaption></figure>

denotes the surface of the Riemann ellipsoid. In the latter, $\Gamma := a_2/a_1$ and $\Xi := a_3/a_1$ are the semi-axes aspect ratios.

The governing equations for the velocity and pressure fields are given by the incompressible Navier-Stokes equations. When written in dimensionless form and in the non-inertial frame of reference, they read

$$ (2a) \qquad \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} + 2\mathrm{Ro}(\mathbf{e}_z \times \mathbf{u}) + \mathrm{Ro}^2 \mathbf{e}_z \times (\mathbf{e}_z \times \mathbf{r}) = \nabla (\Phi - p) + \mathrm{RoEk} \nabla^2 \mathbf{u}, $$

$$ (2b) \qquad \nabla \cdot \mathbf{u} = 0, $$

where $\mathbf{u}(\mathbf{r}, t)$ denotes the velocity field, $p(\mathbf{r}, t)$ the pressure and $\mathbf{r} = x\mathbf{e}_x + y\mathbf{e}_y + z\mathbf{e}_z$ the Cartesian position vector. The Newtonian gravitational potential $\Phi$ is given by

$$ (3) \qquad \Phi(\mathbf{r}, t) = \frac{1}{\pi} \int_{\mathcal{D}_t} \frac{d\mathbf{r}'}{|\mathbf{r} - \mathbf{r}'|}, $$

with the integration performed over the time-dependent fluid domain $\mathcal{D}_t$.

In the above expressions, $\mathrm{Ro}$ and $\mathrm{Ek}$ are the dimensionless Rossby and Ekman numbers, respectively, defined as

$$ (4) \qquad \mathrm{Ro} := \frac{\Omega}{(\pi \rho G)^{1/2}}, \quad \mathrm{Ek} := \frac{\nu}{a_1^2 \Omega}, $$

where $\rho$ is the constant fluid density, $G$ the gravitational constant, and $\nu$ the kinematic viscosity.

We introduce the Cauchy stress tensor of the fluid

$$ (5) \qquad \sigma_{ij} := -p\delta_{ij} + 2\mathrm{RoEk}\varepsilon_{ij} = -p\delta_{ij} + \mathrm{RoEk}(\partial_i u_j + \partial_j u_i), $$

where $\delta$ denotes the Kronecker delta function and $(i, j)$ label the components in the chosen coordinate system (the conversion of tensor components between two curvilinear bases is given in Appendix B). The surface of the ellipsoid being a free surface, we impose a dynamic no-stress boundary condition. Defining the outward unit normal and tangential vectors by $\mathbf{n}$
and $t$, this condition reads $n \cdot (\sigma \cdot n) = 0$ and $t \cdot (\sigma \cdot n) = 0$ (assuming vacuum outside the ellipsoid). In addition, the velocity field is required to satisfy a kinematic boundary condition enforcing continuity at the moving free surface.

We note that for $Ek \to 0$ (i.e. for the inviscid limit $\nu \to 0$), the no-shear stress condition becomes redundant and thus must be discarded. In the viscous case, however, consistency of the solution requires the presence of a thin boundary-layer adjacent to the inner surface of the ellipsoid, within which viscous effects remain significant.

The steady-state velocity field of Riemann ellipsoids is linear in the Cartesian coordinates. Explicitly, they are given by $u_0 = Sr$ (Riemann, 1860), where

$$ (6) \qquad S := \begin{pmatrix} 0 & Q_1 & 0 \\ Q_2 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \text{with } Q_1 := -\frac{\zeta}{1+\Gamma^2} \quad \text{and} \quad Q_2 := \frac{\zeta\Gamma^2}{1+\Gamma^2}. $$

We note that this solution is harmonic, i.e. it satisfies $\nabla^2 u_0 \equiv 0$. Therefore, the Riemann velocity field satisfies equally the Euler and the Navier-Stokes equations. Accordingly, the steady pressure and symmetric part of the velocity tensor are given by

$$ (7) \qquad p_0 = -P_0 F(\mathbf{r}), \quad \nabla u_0 + (\nabla u_0)^T = S + S^T = (Q_1 + Q_2) \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, $$

where $F(x, y, z)$ is the function defined in (1) and $P_0 > 0$ is recovered from the equations of hydrostatic balance. The transpose operator is denoted by the superscript $T$. We note that the strain-rate tensor vanishes identically for $\Gamma = 1$ (i.e. $a_1 = a_2$) or $\zeta = 0$, which corresponds to the Maclaurin and Jacobi families of rigidly rotating ellipsoids. The latter configuration and its stability have been notably studied in Chandrasekhar (1962).

To characterise admissible equilibrium configurations of Riemann ellipsoids, we follow the classical construction of Chandrasekhar (1965) and introduce the dimensionless ratio of internal vorticity to rigid rotation, $f := \zeta/R_0$. Fixing the value of $f$ defines a *Riemann sequence* of equilibria parameterised by the ellipsoidal geometry. In particular, the limiting cases $f = 0$ and $f \to \pm\infty$ correspond to the Jacobi and Dedekind sequences, respectively, representing purely rotational and purely internal-flow-supported equilibria. For such configurations to satisfy gravitational equilibrium, Chandrasekhar derived the coupled algebraic conditions

$$ (8) \qquad \Gamma^2 \left[ A_1 - A_2 + \frac{\Xi^2}{\Gamma^2} (1 - \Gamma^2) A_3 \right] f^2 + 2(1 + \Gamma^2)(A_1 - A_2\Gamma^2)f \\ \qquad \qquad + (1 + \Gamma^2)^2 \left[ A_1 - A_2 + \frac{\Xi^2}{\Gamma^2} (1 - \Gamma^2) A_3 \right] = 0, $$

and

$$ (9) \qquad \mathrm{Ro}^2 = \frac{2B_{12}}{1 + f^2 \Gamma^2 / (1 + \Gamma^2)^2}, $$

where the coefficients $A_i$ and $B_{ij}$ are the classical index symbols defined in (Chandrasekhar, 1969). For a prescribed value of $f$, solving (8) determines the admissible ellipsoidal geometry (through $\Gamma$ or $\Xi$), while (9) subsequently fixes the corresponding Rossby number. Fixing $f = 0$ and $\Gamma = 1$, the Rossby number of Maclaurin spheroids is expressed solely in terms of the eccentricity $e := d_\lambda$ (Chandrasekhar, 1969; Labarbe, 2021).
2.2. **Linear stability theory.** We restrict this work to slightly perturbed Riemann ellipsoids for the linear theory to remain valid. We denote by $\xi(\mathbf{r}, t)$ the Lagrangian displacement, which refers to a fluid parcel which would have been at time $t$ and position $\mathbf{r}$ in the unperturbed configuration. Accordingly, we introduce the Lagrangian difference operator as

$$ (10) \qquad \Delta f(\mathbf{r}, t) := f(\mathbf{r} + \xi(\mathbf{r}, t), t) - f_0(\mathbf{r}, t), $$

where $f$ refers to an arbitrary quantity evaluated in the perturbed configuration, whereas $f_0$ refers to the value of the same quantity but for the unperturbed flow — the *ghostly flow* solution in the terminology of Lynden-Bell and Ostriker (1967). The Eulerian change operator is recovered by setting $\xi \equiv 0$ in the definition (10), reading

$$ (11) \qquad \delta f(\mathbf{r}, t) := f(\mathbf{r}, t) - f_0(\mathbf{r}, t). $$

It is worth mentioning that Schutz and Friedman (1978a, 1978b) introduced a modified operator in their work on Lagrangian perturbation theory. The operator they define involves the Lie derivative with respect to the Lagrangian displacement and allows for a whole description in Lagrangian frame, without need to introduce a mixed Eulerian–Lagrangian formalism. Moreover, adopting the description of Friedman and Schutz eliminates the gauge freedom in the canonical Hamiltonian of the system and thus yields meaningful stability criteria. However, this Lagrangian picture is not suited to study viscous flows as the classical Laplacian would become a Laplace–Beltrami operator for the disturbed metric, rendering the problem highly nonlinear. Therefore, our linear theory is solely based on a mixed Lagrangian–Eulerian formulation to allow for viscosity to be taken into account, as done by Roberts and Stewartson (1963).

We apply operator (10) on (2)–(3) and expand the solution in Taylor series for arbitrarily small displacements. We recover the linearised equations of motion in a mixed Eulerian–Lagrangian description by retaining only first-order terms in disturbances

$$ (12a) \qquad \left( \frac{\partial}{\partial t} - \mathrm{RoEk} \nabla^2 + \mathbf{u}_0 \cdot \nabla + \mathbf{S} + 2\mathrm{Ro}e_z \times \right) \delta \mathbf{u} = \nabla (\delta \Phi - \delta p), $$

$$ (12b) \qquad \nabla \cdot \delta \mathbf{u} = 0, $$

$$ (12c) \qquad \delta \Phi(\mathbf{r}, t) = \frac{1}{\pi} \int_{\partial \mathcal{D}_0} \frac{\xi(\mathbf{r}', t) \cdot \mathbf{n}}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}', $$

where $\partial \mathcal{D}_0$ is the boundary of the equilibrium ellipsoid, i.e. the surface of equation $\lambda = 1$. As described in Appendix A, the outward unit normal vector is given by $\mathbf{n} := e_\lambda$.

Applying the same procedure as before, the linearised boundary conditions are written on the unperturbed free surface and in ellipsoidal coordinates (cf. Comins (1979) for details) as

$$ (13) \qquad \delta p + \xi \cdot \nabla p_0 - 2\mathrm{RoEk} \varepsilon_{\lambda\lambda} = 0, \quad \varepsilon_{\lambda\mu} = 0, \quad \varepsilon_{\lambda\varphi} = 0, $$

supplemented with kinematic condition. Despite the flow having internal motion with non-vanishing components of the stress tensor, their contributions to the linearised boundary conditions are not present since they remain uniform at all time and thus, are of the same order than the background state. Moreover, the disturbance fields are decomposed as $\delta \mathbf{u} := \delta \mathbf{u}^{(0)} + \delta \mathbf{u}^{(\nu)}$, $\delta p := \delta p^{(0)} + \delta p^{(\nu)}$ and $\xi := \xi^{(0)} + \xi^{(\nu)}$, where the superscripts denote whether these fields are solutions to the inviscid or the viscous equations.

In the following, we present an exhaustive analysis of the inviscid configuration, where we obtain a dispersion relation for arbitrary harmonic degree, extending the existing results of Chandrasekhar (1965, 1966) who used the tensor virial method. Subsequently, we present
a procedure to solve the viscous problem by means of a boundary-layer theory, as done by Roberts and Stewartson (1963) and Comins (1979) in the context of Maclaurin spheroids. Without loss of generality, we assume the disturbances to have harmonic motions of the form $\exp(-2i\omega t)$. The real and imaginary parts of the eigenfrequency $\omega$ correspond to, respectively, the growth rate (in magnitude) and the frequency of the disturbances.

### 3. INVISCID RIEMANN ELLIPSOIDS

The boundary condition for the inviscid problem is obtained by taking the limit $Ek \to 0$ in (13), which yields

$$ (14) \qquad [\delta p^{(0)}]_{\partial D_0} - 2P_0 [h_\lambda \xi_\lambda^{(0)}]_{\partial D_0} = 0, $$

where we have used the relation $\nabla p_0 = -|\nabla p_0|\mathbf{n} = -2P_0h_\lambda e_\lambda$. Our objective is to evaluate the different terms in (14) in order to derive a dispersion relation for the eigenfrequency $\omega$, and thereby assess the linear stability of the system for admissible configurations.

#### 3.1. Harmonic expansion of disturbance fields.
The product of the normal Lagrangian displacement and the corresponding metric scale factor, evaluated at the equilibrium boundary $\partial D_0$, is expanded in terms of surface ellipsoidal harmonics — whose definitions are given in Appendix C — as

$$ (15) \qquad [h_\lambda \xi_\lambda^{(0)}]_{\partial D_0} = \hat{\xi}_{nm} S_n^m(\mu, \varphi) e^{-2i\omega t}, $$

The corresponding expansion coefficients are given by

$$ (16) \qquad \hat{\xi}_{nm} = \frac{e^{2i\omega t}}{\varpi_n^m} \int [h_\lambda \xi_\lambda^{(0)}]_{\partial D_0} S_n^m(\mu, \varphi) dS, $$

where integration is performed over the surface of the equilibrium ellipsoid. Unless stated otherwise, all integrals for the surface element $dS$ appearing in the remainder of this article are taken over this same geometry.

Using the definition of the Lagrangian displacement,

$$ (17) \qquad \delta \mathbf{u}^{(0)} = \frac{\partial \xi^{(0)}}{\partial t} = -2i\omega \xi^{(0)}, $$

together with expression (63), the coefficients (16) can be expressed in terms of the Cartesian components of the inviscid velocity perturbation as

$$ (18) \qquad \hat{\xi}_{nm} = \frac{ie^{2i\omega t}}{2\omega\varpi_n^m} \int \left[ \frac{\partial x}{\partial\lambda} \delta u_x^{(0)} + \frac{\partial y}{\partial\lambda} \delta u_y^{(0)} + \frac{\partial z}{\partial\lambda} \delta u_z^{(0)} \right]_{\partial D_0} S_n^m(\mu, \varphi) dS. $$

Using the expansion (72) and exploiting the symmetry of the problem, the perturbation to the Newtonian gravitational potential can be expressed throughout space as (Chandrasekhar, 1969)

$$ (19) \qquad \delta\Phi^{(0)}(\lambda, \mu, \varphi, t) = \begin{cases} \hat{\Phi}_{nm} \frac{E_n^m(\lambda)}{E_n^m(1)} S_n^m(\mu, \varphi) e^{-2i\omega t}, & d_\mu \le \lambda \le 1, \\ \hat{\Phi}_{nm} \frac{F_n^m(\lambda)}{F_n^m(1)} S_n^m(\mu, \varphi) e^{-2i\omega t}, & 1 \le \lambda < \infty, \end{cases} $$

where the interior and exterior ellipsoidal harmonics are introduced in Appendix C.
The expansion coefficients are obtained by applying Gauss' theorem at the unperturbed boundary $\partial D_0$, namely

$$ (20) \qquad [\mathbf{n} \cdot \nabla \delta\Phi^{(0)}]_{\partial D_0^+} - [\mathbf{n} \cdot \nabla \delta\Phi^{(0)}]_{\partial D_0^-} = -4[\mathbf{n} \cdot \boldsymbol{\xi}^{(0)}]_{\partial D_0}, $$

where $\partial D_0^\pm$ denote the exterior and interior sides of the boundary, respectively. In ellipsoidal coordinates, this condition can be rewritten such that the right-hand side directly involves the quantity defined in (15). Explicitly,

$$ (21) \qquad \left[ \frac{\partial \delta \Phi^{(0)}}{\partial \lambda} \right]_{\partial D_0^+} - \left[ \frac{\partial \delta \Phi^{(0)}}{\partial \lambda} \right]_{\partial D_0^-} = -4 [h_\lambda \xi_\lambda^{(0)}]_{\partial D_0} = -4 \hat{\xi}_{nm} S_n^m e^{-2i\omega t}. $$

Substituting (19) into (21) and making use of the Wronskian relation (76), we obtain

$$ (22) \qquad [\delta\Phi^{(0)}]_{\partial D_0} = 4\hat{\xi}_{nm} W_n(1)^{-1} E_n^m(1) F_n^m(1) S_n^m(\mu, \varphi) e^{-2i\omega t} $$

$$ (23) \qquad = \frac{4\hat{\xi}_{nm}}{2n+1} \sqrt{(1-d_{\mu}^2)(1-d_{\varphi}^2)} E_n^m(1) F_n^m(1) S_n^m(\mu, \varphi) e^{-2i\omega t}. $$

To complete the inviscid boundary condition (14) and derive the dispersion relation for Riemann ellipsoids, it remains to obtain an explicit expression for the pressure perturbation. This is achieved in the following section by deriving a generalised Poincaré equation governing the normal modes of oscillation.

## 3.2. Normal modes of oscillation.
We first introduce an inviscid hydrodynamic potential, simultaneously taking into account the effects of pressure and gravity,

$$ (24) \qquad \delta W^{(0)} = \delta\Phi^{(0)} - \delta p^{(0)}. $$

Applying the divergence operator on (12a) and using (12b), we obtain a general Poincaré equation governing the small oscillations of the rotating fluid,

$$ (25) \qquad \mathcal{A}\{\delta\mathbf{u}^{(0)}\} = \mathcal{P}\{\delta W^{(0)}\}. $$

The Poincaré and advection operators are respectively given by

$$ (26) \qquad \mathcal{P}\{\delta W^{(0)}\} := \nabla \cdot (L^{-1}\nabla\delta W^{(0)}) := \text{tr}[(L^{-1})^T \text{Hess}(\delta W^{(0)})], $$

$$ (27) \qquad \mathcal{A}\{\delta\mathbf{u}^{(0)}\} := \nabla \cdot [L^{-1}(\mathbf{u}_0 \cdot \nabla\delta\mathbf{u}^{(0)})], $$

where we have used the identity $\nabla \cdot (A\mathbf{q}) = \text{tr}[A^T(\nabla \otimes \mathbf{q})]$ that is valid for any square matrix $A$ and vector $\mathbf{q}$. Hess denotes the Hessian matrix operator. The matrix $L$ is given by

$$ (28) \qquad L := \begin{pmatrix} -2i\omega & Q_1 - 2\text{Ro} & 0 \\ Q_2 + 2\text{Ro} & -2i\omega & 0 \\ 0 & 0 & -2i\omega \end{pmatrix}, $$

and is skew-Hermitian when $\omega \in \mathbb{R}$ and $Q_1 = Q_2 = 0$. In the absence of gravity, the spectrum of the Poincaré operator is known to be real and bounded in the interval $[-\omega, \omega]$, with polynomial eigenfunctions (Greenspan, 1965). When dynamical instabilities occur (usually through a Hamilton-Hopf bifurcation (Kirillov, 2021)), the spectrum becomes complex-valued and these properties no longer hold.

When $f = 0$, corresponding to a triaxial Jacobi ellipsoid, the problem reduces to $\mathcal{P}\{W\} = 0$, which is the classical Poincaré equation (Poincaré, 1885). Its solutions form a one-parameter family of homogeneous polynomials of the form $\alpha(x + iy)^n$, as established by Cartan (1922). Since the base flow $\mathbf{u}_0$ depends linearly on the Cartesian coordinates, the advection operator
A maps the space of Cartesian polynomials of degree $n-1$ onto itself. Accordingly, we seek solutions such that $\delta\mathbf{u}^{(0)} \in \mathbb{P}_{n-1}^3$ and $\delta W^{(0)} \in \mathbb{P}_n^1$, where $\mathbb{P}_n^d$ denotes the $d$-dimensional space of polynomials of degree $n$.

An extension of Cartan's polynomial family solving the generalised Poincaré equation (25) is found to be

$$ (29) \qquad \delta W^{(0)}(x, y, z) = -\alpha(x + iy)^n + \beta z^2 (x + \gamma y)^{n-2}, $$

where $(\alpha, \beta, \gamma)$ are real constants. It can be shown that $\beta$ depends linearly on $\alpha$, implying that (29) also constitutes a one-parameter family of solutions for Riemann ellipsoids.

Substituting this family into the general Poincaré equation (25), the inviscid velocity perturbation $\delta\mathbf{u}^{(0)} = (\delta u_x^{(0)}, \delta u_y^{(0)}, \delta u_z^{(0)})$ can be written in the form

$$ (30) \qquad \delta u_x^{(0)}(x, y, z) = \sum_{k=0}^{n-1} \hat{u}_k x^{n-1-k} y^k + \hat{u}_n z^2, \quad \delta u_y^{(0)}(x, y, z) = \sum_{k=0}^{n-1} \hat{v}_k x^{n-1-k} y^k + \hat{v}_n z^2, \\ \delta u_z^{(0)}(x, y, z) = z \sum_{k=0}^{n-2} \hat{w}_k x^{n-2-k} y^k, $$

where the polynomial coefficients $(\hat{u}_k, \hat{v}_k, \hat{w}_k)$ and the constants $(\beta, \gamma)$ are determined by term-by-term identification with the Euler and continuity equations. It can be shown that $(\beta, \gamma, \hat{w}_k)$ depend linearly on the background vorticity $\zeta$. Hence, in the limit $\zeta = 0$, the classical Cartan solutions parameterised by $\alpha$ are recovered. This family of solutions, for arbitrary harmonic degree and for axisymmetric configurations (i.e. Maclaurin spheroids with $\Gamma = 1$) is given by (Comins, 1979),

$$ (31) \qquad \delta u_x^{(0)} = -\frac{i n \alpha (x + iy)^{n-1}}{2 \mathrm{Ro}(\sigma + 1)}, \quad \delta u_y^{(0)} = \frac{n \alpha (x + iy)^{n-1}}{2 \mathrm{Ro}(\sigma + 1)}, \quad \delta u_z^{(0)} = 0. $$

For $n > 2$, the velocity perturbations of Riemann ellipsoids are not harmonic. This feature is a major obstacle to extend the classical viscous stability theory developed for Maclaurin spheroids, as historically done by Roberts and Stewartson (1963). Since the potential $\delta W^{(0)}$ is also non-harmonic in the triaxial case, the associated pressure perturbation

$$ (32) \qquad \delta p^{(0)} := \delta \Phi^{(0)} - \delta W^{(0)} = \delta \Phi^{(0)} + \alpha (x + iy)^n - \beta z^2 (x + \gamma y)^{n-2}, $$

is likewise non-harmonic for $n > 2$. Here, $\delta \Phi^{(0)}$ is expressed in terms of the velocity field (30) via relations (22) and (18).

### 3.3. Dispersion relation and stability results.
To evaluate the hydrodynamic potential at the boundary of the ellipsoid, we expand expression (29) in terms of surface harmonics,

$$ (33) \qquad [\delta W^{(0)}]_{\partial \mathcal{D}_0} = \hat{W}_{nm} \mathbb{S}_n^m e^{-2it\sigma\mathrm{Ro}}, $$

with expansion coefficients

$$ (34) \qquad \hat{W}_{nm} = \frac{e^{2it\sigma\mathrm{Ro}}}{\varpi_n^m} \int [-\alpha(x + iy)^n + \beta z^2 (x + \gamma y)^{n-2}]_{\partial \mathcal{D}_0} \mathbb{S}_n^m dS. $$

Substituting (15), (22), and (33) into the inviscid boundary condition (14) yields the dispersion relation for inviscid Riemann ellipsoids,

$$ (35) \qquad \mathcal{D}_n^m(\sigma) = 2\theta_n^m \hat{\xi}_{nm} + \hat{W}_{nm} = 0, $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Dispersion curves of the $m = n = 2$ inviscid Riemann ellipsoid over the planar aspect ratio $\Gamma$. Red (blue) lines denote the frequencies (growth rates) of the system, computed at different values of $f$ according to the legend.</figcaption></figure>

where

$$ (36) \qquad \theta_n^m := P_0 - 2W(1)^{-1}E_n^m(1)F_n^m(1). $$

Using expressions (18) and (34), the dispersion relation can be written explicitly as

$$ (37) \qquad \begin{aligned} \mathcal{D}^{(0)}(\sigma) = & \frac{i\theta_n^m}{\sigma \mathrm{Ro}} \int \left[ \frac{\partial x}{\partial \lambda} \delta u_x^{(0)} + \frac{\partial y}{\partial \lambda} \delta u_y^{(0)} + \frac{\partial z}{\partial \lambda} \delta u_z^{(0)} \right]_{\partial \mathcal{D}_0} \mathbb{S}_n^m dS \\ & - \int \left[ \alpha(x + iy)^n - \beta z^2 (x + \gamma y)^{n-2} \right]_{\partial \mathcal{D}_0} \mathbb{S}_n^m dS = 0. \end{aligned} $$

Once the degree and order of the ellipsoidal harmonic are specified, the associated Lamé functions are computed following the procedure described in Appendix C. The linear stability of the system then follows directly from solving (37) for the eigenfrequency $\sigma$.

Finally, setting $f = 0$ (or equivalently $\zeta = 0$) in (37) yields the classical dispersion relation for Jacobi ellipsoids (Chandrasekhar, 1962). In addition, setting $\Gamma = 1$ yields the characteristic equation for the stability of Maclaurin spheroids (Chandrasekhar, 1969; Comins, 1979)

$$ (38) \qquad \mathcal{D}_{\mathrm{Mc}}^{(0)}(\sigma) = \mathrm{Ro}^2\sigma(\sigma + 1) - \theta_n^m = 0. $$

We present in Figure 2 a representative set of dispersion curves for the case $m = n = 2$. This mode corresponds to the classical *bar mode*, which is known to be the most unstable configuration for Maclaurin spheroids (Chandrasekhar, 1969). In the limit $f = 0$, we recover the stability diagram of the Jacobi ellipsoid, including the bifurcation point located at $\Gamma_c \approx 0.4635$ and $\Xi_c \approx 0.3632$. The value of $\cos^{-1}(\Xi_c)$ at this critical point is tabulated in Chandrasekhar (1962), and our result agrees with this value to within less than $1^\circ$.

As $f$ departs from zero, similar bifurcations persist for $f > 0$; however, the ellipsoid remains linearly unstable beyond the critical threshold. For $f < 0$, additional instabilities arise through the crossings of distinct frequency branches. These unstable modes are attributed to the elliptic instability induced by the non-axisymmetric deformation of the internal streamlines, as discussed in (Kerswell, 2002).
## 4. VISCOUS RIEMANN ELLIPSOIDS

In this section, we investigate the effect of weak viscosity by assuming $Ek \ll 1$. We restrict our attention to the case $n = 2$, for which the inviscid velocity field (30) remains harmonic (see Section 3). For higher harmonic degrees ($n > 2$), the inviscid solutions no longer satisfy the viscous equations of motion, and $O(Ek)$ corrections to the velocity field are required for compatibility. By contrast, when $n = 2$, viscous effects enter the problem solely through the boundary conditions.

In the present framework, viscous corrections arise from the requirement that the stress-free boundary conditions be satisfied to first order in $Ek$. This is achieved by introducing a thin viscous boundary-layer at the inner edge of the Riemann ellipsoid and solving the associated Prandtl boundary-layer equations. This approach allows us to compute the leading-order viscous corrections to the normal stress condition and to derive an analytic dispersion relation, accounting for viscous effects.

### 4.1. Surface stress conditions.
We start by enforcing the normal stress boundary condition (13), which is decomposed as

$$ (39) \qquad [\delta p^{(0)} + \boldsymbol{\xi}^{(0)} \cdot \nabla p_0]_{\partial \mathcal{D}_0} + [\delta p^{(\nu)} + \boldsymbol{\xi}^{(\nu)} \cdot \nabla p_0 - 2\mathrm{RoEk} \varepsilon_{\lambda\lambda}^{(0)}] = 0, $$

where the first term in brackets corresponds to the inviscid contribution, which is determined by the dispersion relation obtained in Section 3. The second term in brackets is anticipated to be of order $O(Ek)$. As a consequence, the viscous contribution arising from the Cauchy stress tensor involves only the inviscid normal velocity field.

Expanding (39) in surface harmonics and integrating over the equilibrium boundary yields

$$ (40) \qquad \int [\delta p^{(0)} + \boldsymbol{\xi}^{(0)} \cdot \nabla p_0 - 2\mathrm{RoEk} \varepsilon_{\lambda\lambda}^{(0)}]_{\partial \mathcal{D}_0} \mathbb{S}_n^m dS + \int [\delta p^{(\nu)} + \boldsymbol{\xi}^{(\nu)} \cdot \nabla p_0]_{\partial \mathcal{D}_0} \mathbb{S}_n^m dS = 0, $$

where the first integral is evaluated using the inviscid analysis developed in Section 3.

Using the viscous counterpart of (17), together with the definition of the perturbation to the gravitational potential given in (22), the second integral may be rewritten as

$$ (41) \qquad \mathcal{J} = \int \left[ \delta p^{(\nu)} - \delta \Phi^{(\nu)} - \frac{i \theta_n^m}{\sigma \mathrm{Ro}} h_\lambda \delta u_\lambda^{(\nu)} \right]_{\partial \mathcal{D}_0} \mathbb{S}_n^m dS. $$

The objective of the present section is to evaluate this integral. To this end, we introduce a viscous boundary-layer formulation near the inner surface of the ellipsoid, following the approach previously employed for Maclaurin spheroids by Roberts and Stewartson (1963) and Comins (1979).

### 4.2. Boundary-layer theory.
We assume that advection in the direction normal to the boundary is negligible at leading-order in viscosity. In addition, we require the viscous correction terms to vanish at the inner surface of the ellipsoid and to asymptotically match the inviscid solution away from the boundary layer. Denoting the characteristic thickness of the layer by $d_\nu = O(Ek^{1/2})$, classical boundary-layer scaling arguments (Landau & Lifshitz, 1987) follow from enforcing consistency between the equations of motion and the boundary conditions at higher order in $Ek$. This yields

$$ \frac{\partial}{\partial \lambda} = O(Ek^{-1/2}), \quad \frac{\partial}{\partial \mu} = O(1), \quad \frac{\partial}{\partial \varphi} = O(1), $$
together with

$$
\begin{align*}
\delta u_{\lambda}^{(\nu)} &= O(\text{Ek}), \quad \delta u_{\mu}^{(\nu)} = O(\text{Ek}^{1/2}), \quad \delta u_{\varphi}^{(\nu)} = O(\text{Ek}^{1/2}), \quad \delta p^{(\nu)} = O(\text{Ek}), \quad \delta \Phi^{(\nu)} = O(\text{Ek}), \\
& u_{0\lambda} = O(\text{Ek}^{1/2}), \quad u_{0\mu} = O(1), \quad u_{0\varphi} = O(1),
\end{align*}
$$

where the scaling of the background velocity field $\mathbf{u}_0 := (u_{0\lambda}, u_{0\mu}, u_{0\varphi})$ follows from the
incompressibility constraint.

We now expand the incompressible Navier-Stokes equations (12a) together with the con-
tinuity equation (12b) in powers of Ek and retain the leading-order contributions in each
expression. Projecting the equations onto the normal direction and enforcing the solenoidal
condition yields, respectively,

$$
(42) \qquad h_{\mu}h_{\nu} \frac{\partial}{\partial \lambda} \left[ \delta p^{(\nu)} - \delta \Phi^{(\nu)} \right] = 2\operatorname{Roh}\lambda \left[ \frac{\partial z}{\partial \varphi} h_{\mu} \delta u_{\mu}^{(\nu)} - \frac{\partial z}{\partial \mu} h_{\varphi} \delta u_{\varphi}^{(\nu)} \right],
$$

$$
(43) \qquad h_{\mu} h_{\varphi} \frac{\partial}{\partial \lambda} \delta u_{\lambda}^{(\nu)} = - \left[ \frac{\partial}{\partial \mu} (h_{\lambda} h_{\varphi} \delta u_{\mu}^{(\nu)}) + \frac{\partial}{\partial \varphi} (h_{\lambda} h_{\mu} \delta u_{\varphi}^{(\nu)}) \right].
$$

Combining (42) and (43) so as to recover the integrand of $\mathcal{J}$, we obtain

$$
(44) \quad h_{\mu} h_{\varphi} \frac{\partial}{\partial \lambda} \left[ \delta p^{(\nu)} - \delta \Phi^{(\nu)} - \frac{i \theta_n^m}{\sigma \mathrm{Ro}} h_{\lambda} \delta u_{\lambda}^{(\nu)} \right] = 2 \mathrm{Roh} \lambda \left[ \frac{\partial z}{\partial \varphi} h_{\mu} \delta u_{\mu}^{(\nu)} - \frac{\partial z}{\partial \mu} h_{\varphi} \delta u_{\varphi}^{(\nu)} \right] \\
\qquad + \frac{i \theta_n^m}{\sigma \mathrm{Ro}} h_{\lambda} \left[ \frac{\partial}{\partial \mu} (h_{\lambda} h_{\varphi} \delta u_{\mu}^{(\nu)}) + \frac{\partial}{\partial \varphi} (h_{\lambda} h_{\mu} \delta u_{\varphi}^{(\nu)}) \right].
$$

It is important to notice that every geometric terms (e.g. the scale factors) are evaluated at
the boundary and thus, are independent of the normal coordinate $\lambda$.

We now expand the viscous corrections to the tangential velocity components in terms
of surface harmonics, allowing for an arbitrary dependence on the normal coordinate. This
yields

$$
(45) \qquad \delta u_{\mu}^{(\nu)} = \mathrm{Ek}^{1/2} \hat{u}_{\mu}(\lambda) \mathbb{S}_{n}^{m}(\mu, \varphi), \quad \delta u_{\varphi}^{(\nu)} = \mathrm{Ek}^{1/2} \hat{u}_{\varphi}(\lambda) \mathbb{S}_{n}^{m}(\mu, \varphi).
$$

Substituting this ansatz into (44), multiplying by the appropriate metric factors, integrating
over the surface element (73), and invoking Leibniz's rule (since all $\lambda$-dependent quantities
other than the velocity corrections are constant), we obtain

$$
(46) \qquad \frac{d\mathcal{J}}{d\lambda} = \mathrm{Ek}^{1/2} [\hat{u}_{\mu}(\lambda)\mathcal{J}_{\mu} + \hat{u}_{\varphi}(\lambda)\mathcal{J}_{\varphi}],
$$

where the coefficients are explicitly given by

$$
(47) \quad \mathcal{J}_{\mu} := 2\mathrm{Ro} \int \frac{h_{\lambda}}{h_{\varphi}} \frac{\partial z}{\partial \varphi} [\mathbb{S}_{n}^{m}]^2 dS + \frac{i \theta_{n}^{m}}{\sigma \mathrm{Ro}} \int \frac{h_{\lambda}}{h_{\mu} h_{\varphi}} \frac{\partial}{\partial \mu} (h_{\lambda} h_{\varphi} \mathbb{S}_{n}^{m}) \mathbb{S}_{n}^{m} dS,
$$

$$
(48) \quad \mathcal{J}_{\varphi} := 2\mathrm{Ro} \int \frac{h_{\lambda}}{h_{\mu}} \frac{\partial z}{\partial \mu} [\mathbb{S}_{n}^{m}]^2 dS + \frac{i \theta_{n}^{m}}{\sigma \mathrm{Ro}} \int \frac{h_{\lambda}}{h_{\mu} h_{\varphi}} \frac{\partial}{\partial \varphi} (h_{\lambda} h_{\mu} \mathbb{S}_{n}^{m}) \mathbb{S}_{n}^{m} dS.
$$

As anticipated, $\mathcal{J}$ is of order $O(\text{Ek})$.

We now seek explicit expressions for the expansion coefficients $(\hat{u}_{\mu}, \hat{u}_{\varphi})$. To this end, we
project the Navier-Stokes equations onto the tangent plane and retain the leading-order
viscous contributions. This procedure yields a coupled system of Prandtl-type equations,

$$ (49a) \quad \left[ -2i\sigma\mathrm{Ro} + \mathbf{u}_0 \cdot \nabla - \frac{\mathrm{Ro}\mathrm{Ek}}{h_{\lambda}^2} \frac{\partial^2}{\partial\lambda^2} + \frac{1}{h_{\mu}} \frac{\partial u_{0\mu}}{\partial\mu} \right] \delta u_{\mu}^{(\nu)} + \left[ \frac{1}{h_{\varphi}} \frac{\partial u_{0\mu}}{\partial\varphi} - \frac{2\mathrm{Ro}}{h_{\lambda}} \frac{\partial z}{\partial\lambda} \right] \delta u_{\varphi}^{(\nu)} = 0, $$

$$ (49b) \quad \left[ -2i\sigma\mathrm{Ro} + \mathbf{u}_0 \cdot \nabla - \frac{\mathrm{Ro}\mathrm{Ek}}{h_{\lambda}^2} \frac{\partial^2}{\partial\lambda^2} + \frac{1}{h_{\varphi}} \frac{\partial u_{0\varphi}}{\partial\varphi} \right] \delta u_{\varphi}^{(\nu)} + \left[ \frac{1}{h_{\mu}} \frac{\partial u_{0\varphi}}{\partial\mu} + \frac{2\mathrm{Ro}}{h_{\lambda}} \frac{\partial z}{\partial\lambda} \right] \delta u_{\mu}^{(\nu)} = 0. $$

In the absence of prevalent motion, the system (49) reduces to a pair of coupled non-homogeneous heat equations whose solutions are described analytically; see Roberts and Stewartson (1963) or Comins (1979).

### 4.3. Analytic solutions and secular instability.
We require the viscous correction fields to decay away from the boundary so that the velocity field asymptotically matches the inviscid solution in the bulk of the ellipsoid. This condition is expressed as

$$ (50) \qquad \lim_{\eta \to -\infty} [\hat{u}_{\mu}] = \lim_{\eta \to -\infty} [\hat{u}_{\varphi}] = 0, $$

where we introduced a stretched normal coordinate

$$ (51) \qquad \eta := d_{\nu}^{-1}(\lambda - 1) = \frac{\lambda - 1}{\mathrm{Ek}^{1/2}} \in (-\infty, 0], $$

which captures variations across the boundary layer.

Substituting (45) into the system (49) and integrating over the tangential coordinates yields a coupled system of linear ordinary differential equations with constant coefficients,

$$ (52a) \qquad a_0 \frac{d^2}{d\eta^2} \hat{u}_\mu(\eta) + a_1 \frac{d}{d\eta} \hat{u}_\mu(\eta) + a_2 \hat{u}_\mu(\eta) + a_3 \hat{u}_\varphi(\eta) = 0, $$

$$ (52b) \qquad b_0 \frac{d^2}{d\eta^2} \hat{u}_\varphi(\eta) + b_1 \frac{d}{d\eta} \hat{u}_\varphi(\eta) + b_2 \hat{u}_\varphi(\eta) + b_3 \hat{u}_\mu(\eta) = 0, $$

where the boundary-layer scaling introduced in Section 4.2 has been used to isolate the viscous contribution from the global dynamics.

Similarly, integrating the tangential stress boundary conditions (13) yields

$$ (53) \qquad \lim_{\eta=0} \left[ \varpi_n^m \frac{d}{d\eta} \hat{u}_\mu + \tau_\mu \right] = \lim_{\eta=0} \left[ \varpi_n^m \frac{d}{d\eta} \hat{u}_\varphi + \tau_\varphi \right] = 0. $$

Explicit expressions for the coefficients ($a_j, b_j, \tau_\mu, \tau_\varphi$) are provided in Appendix D.

Introducing the state vector $\xi := (\hat{u}_\mu, \hat{u}_\varphi)$, the system (52) can be written compactly as

$$ (54) \qquad \begin{pmatrix} a_0 & 0 \\ 0 & b_0 \end{pmatrix} \frac{d^2 \xi}{d\eta^2} + \begin{pmatrix} a_1 & 0 \\ 0 & b_1 \end{pmatrix} \frac{d \xi}{d\eta} + \begin{pmatrix} a_2 & a_3 \\ b_3 & b_2 \end{pmatrix} \xi := \mathbf{A}_0 \frac{d^2 \xi}{d\eta^2} + \mathbf{A}_1 \frac{d \xi}{d\eta} + \mathbf{A}_2 \xi = \mathbf{0}, $$

supplemented with boundary conditions (50) and (53). This linear system admits exponentially decaying solutions of the form $\xi \sim \exp(\kappa\eta)$. The corresponding characteristic equation is

$$ (55) \qquad P(\kappa) = \det (\kappa^2 \mathbf{A}_0 - \kappa \mathbf{A}_1 + \mathbf{A}_0) = 0, $$

which is a quartic polynomial in $\kappa$ whose roots generally possess both positive and negative real parts.
Since physical admissibility requires decay away from the boundary, only roots with positive real parts are retained. Denoting these roots by $\kappa_1$ and $\kappa_2$, the tangential velocity corrections may be written as

$$ (56) \qquad \hat{u}_{\mu}(\eta) = A_1 e^{\kappa_1 \eta} + A_2 e^{\kappa_2 \eta}, \quad \hat{u}_{\varphi}(\eta) = B_1 e^{\kappa_1 \eta} + B_2 e^{\kappa_2 \eta}. $$

The constants $K_j$ are obtained by enforcing the boundary conditions (50)–(53) and substituting (56) into (54), which yields the linear system

$$ (57) \qquad \begin{pmatrix} A_1 \\ A_2 \\ B_1 \\ B_2 \end{pmatrix} = \begin{pmatrix} P_a(\kappa_1) & P_a(\kappa_2) & a_3 & a_3 \\ b_3 & b_3 & P_b(\kappa_1) & P_b(\kappa_2) \\ -\bar{\omega}_n^m \kappa_1 & -\bar{\omega}_n^m \kappa_2 & 0 & 0 \\ 0 & 0 & -\bar{\omega}_n^m \kappa_1 & -\bar{\omega}_n^m \kappa_2 \end{pmatrix}^{-1} \begin{pmatrix} 0 \\ 0 \\ \tau_{\mu} \\ \tau_{\varphi} \end{pmatrix}, $$

where $P_a(\kappa) := a_0\kappa^2 + a_1\kappa + a_2$ and $P_b(\kappa) := b_0\kappa^2 + b_1\kappa + b_2$.

The integral $\mathcal{I}$ entering the normal stress condition is then obtained by integrating (46) over the stretched coordinate $\eta$ and substituting the expressions for $\kappa_j$ and $K_j$. This yields

$$ \mathcal{I} = \mathrm{Ek} \left[ \mathcal{I}_{\mu} \int_{-\infty}^{0} \hat{u}_{\mu} d\eta + \mathcal{I}_{\varphi} \int_{-\infty}^{0} \hat{u}_{\varphi} d\eta \right] = \frac{\mathrm{Ek}}{\kappa_1 \kappa_2} \left[ \mathcal{I}_{\mu}(\kappa_2 A_1 + \kappa_1 A_2) + \mathcal{I}_{\varphi}(\kappa_2 B_1 + \kappa_1 B_2) \right]. $$

The dispersion relation governing the stability of viscous Riemann ellipsoids is given by

$$ (58) \qquad \mathcal{D}(\sigma) = \mathcal{D}^{(0)}(\sigma) + \mathrm{Ek}\mathcal{D}^{(\nu)}(\sigma) = 0, $$

where $\mathcal{D}^{(0)}$ is obtained from (37) and

$$ (59) \qquad \mathcal{D}^{(\nu)}(\sigma) = -2\mathrm{Ro} \int \left[ \varepsilon_{\lambda\lambda}^{(0)} \right]_{\mathcal{D}_0} \mathcal{S}_n^m dS + \mathcal{I}_{\mu} \left( \frac{A_1}{\kappa_1} + \frac{A_2}{\kappa_2} \right) + \mathcal{I}_{\varphi} \left( \frac{B_1}{\kappa_1} + \frac{B_2}{\kappa_2} \right). $$

We summarise the procedure for performing the linear stability analysis as follows:

1. Fix the degree and order $(n, m)$ of the perturbation and compute the corresponding ellipsoidal harmonic $\mathcal{S}_n^m$, as described in Appendix C.

2. Determine the polynomial coefficients of the inviscid velocity field (30) following the procedure outlined in Section 3. The case $n = 2$ is already provided explicitly.

3. Compute the roots of the characteristic polynomial (55) with positive real parts and determine the associated coefficients $K_j$ by solving the linear system (57).

4. Substitute the resulting expressions into the viscous dispersion relation (58) and solve for the complex frequency $\sigma$ using a suitable root-finding algorithm. The stability of the system follows from the sign of the growth rate.

As mentioned at the beginning of this section, the above procedure applies directly only to the case $n = 2$, for which the inviscid velocity perturbations (30) satisfy the Navier-Stokes equations exactly. For $n > 2$, additional $O(\mathrm{Ek})$ corrections to the bulk velocity field would be required, since the inviscid solutions are no longer harmonic. Although the present analysis is not extended to that case, the methodology would remain formally unchanged.

Although the viscous formulation yields an explicit dispersion relation, the evaluation of the surface integrals entering (41) remains technically challenging. Both integral contributions in (40) exhibit endpoint singularities and must be reformulated in terms of complete elliptic integrals of the first, second, and third kinds. However, the second contribution involves integrals that do not admit closed analytical forms and therefore require numerical quadrature. By contrast, the first viscous contribution can be evaluated explicitly (using symbolic calculus softwares like Maple or Mathematica), without resorting to numerical discretisation. We
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Same panels as in Figure 2 but for the viscous case with Ek = $10^{-1}$, computed from the first term in expression (40).</figcaption></figure>

therefore restrict the present computations to the case where $\mathcal{J} \equiv 0$, which still captures a substantial viscous contribution. We emphasise that this restriction is motivated by analytical and numerical feasibility rather than asymptotic considerations, and that the neglected contribution is formally of the same order in Ek.

The same set of panels as in Figure 2 is shown in Figure 3 for the viscous case with Ek = $10^{-1}$. Viscosity primarily manifests itself through the appearance of *avoided crossings* in the neighbourhood of inviscid Hamilton–Hopf bifurcation points (Labarbe, 2021), while the global structure of the dispersion curves remains close to that of the inviscid case. Notably, we observe viscosity-driven secular instabilities for Riemann ellipsoids outside the inviscid interval of instability. This behaviour slightly contrasts with that of Maclaurin spheroids, for which arbitrarily small viscosity destabilises the equilibrium at which the axisymmetric sequence bifurcates toward a triaxial configuration (namely the Jacobi ellipsoids). Riemann ellipsoids being intrinsically non-axisymmetric, viscosity does induce an analogous loss of stability but not through symmetry breaking. In the viscous framework, Riemann ellipsoids of almost all admissible shapes (except for the irrotational sequence with $f = -2$) are found to be unstable to second-order harmonic perturbations.

## 5. DISCUSSION

In this work, we have developed a unified and tractable linear stability framework for S-type Riemann ellipsoids, addressing long-standing limitations of previous inviscid analyses (Chandrasekhar, 1965, 1966; Lebovitz & Lifschitz, 1996). Earlier studies based on the virial method led to global tensorial equations whose complexity grows rapidly with the harmonic order, thereby restricting practical computations to second- and third-order normal modes (Chandrasekhar, 1965). As noted by Lebovitz and Lifschitz (1996), higher-order disturbances are expected to be less strongly influenced by the gravitational field. However, a systematic treatment of this regime had remained unavailable, except through short-wavelength approximations. By contrast, the present formulation provides a consistent local approach that is valid for arbitrary harmonic perturbations, while remaining computationally efficient across the relevant parameter space.
In the inviscid limit, we derived a generalised Poincaré-type equation governing small-amplitude oscillations of a triaxial ellipsoid with internal strain. We showed that its solutions form a one-parameter family of polynomials, extending the classical results of Poincaré (1885) and Cartan (1922), which are recovered as special cases. This formulation offers a systematic alternative to the virial tensor method traditionally employed for Riemann ellipsoids (Chandrasekhar, 1965, 1966), while avoiding the wavelength restrictions inherent to short-wavelength (WKB) analyses (Lebovitz & Lifschitz, 1996). As a result, the full inviscid spectrum can be computed explicitly for arbitrary ellipsoidal harmonics, yielding dispersion relations in closed-form. This represents an extension of earlier investigations of differentially rotating self-gravitating bodies (Lynden-Bell & Ostriker, 1967). Notably, we show that almost all admissible S-type Riemann ellipsoids are unstable to second-order harmonic disturbances.

Building on these results, we examined the influence of weak viscosity by developing a small-Ekman-number asymptotic theory based on Prandtl's boundary-layer framework (Landau & Lifshitz, 1987). For the case $n = 2$, where the inviscid velocity field remains harmonic, viscous effects enter at leading order solely through the surface stress conditions. By resolving the associated boundary-layer dynamics analytically, we derived first-order viscous corrections to the inviscid dispersion relation and obtained an explicit criterion for secular instability. This analysis provides a systematic description of dissipation-induced instabilities in Riemann ellipsoids, closely following the classical treatments available for Maclaurin spheroids (Comins, 1979; Labarbe, 2021; Roberts & Stewartson, 1963). As in those systems, dissipation acts as a destabilising mechanism (Bloch et al., 1994), undermining the gyroscopic balance that sustains rotating self-gravitating configurations (Thompson & Tait, 1867). However, a definitive assessment of viscous stability requires the numerical evaluation of the second contribution in (40), which is formally of the same asymptotic order but involves surface integrals that do not admit closed-form expressions. The development of robust numerical methods to compute this term is left for future work.

A natural extension of this work concerns the Chandrasekhar–Friedman–Schutz instability (Schutz & Friedman, 1978a, 1978b). Addressing this problem requires incorporating relativistic corrections to the Newtonian gravitational potential, following the pioneering approach of Chandrasekhar (1970). The destabilisation of rapidly rotating self-gravitating stars by gravitational radiation reaction remains a compelling target for future gravitational-wave observations. Although direct detection from isolated stellar objects remains beyond the sensitivity of current interferometers (Abott et al., 2016), extending existing analyses (Labarbe, 2021) to include internal fluid motions may contribute to improved modelling and parameterisation for next-generation detectors.

Beyond their mathematical interest, the results we presented bear direct relevance for geo-physical and astrophysical fluid dynamics. Riemann ellipsoids serve as idealised yet physically meaningful models of rotating planets, stars, and fluid interiors with internal strain. Their stability properties play a central role in secular evolution, mode selection, and energy dissipation in such systems. The framework developed in this study opens the door for systematic investigations of linear instabilities in broader classes of rotating self-gravitating flows.

**Funding.** This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence Université Côte d'Azur under reference number ANR-15-IDEX-01.

**Declaration of interests.** The author reports no conflict of interest.
## APPENDIX A. DESCRIPTION OF THE ELLIPSOIDAL SYSTEM OF COORDINATES

Cartesian coordinates $(x, y, z)$ are expressed in terms of the confocal ellipsoidal coordinates $(\lambda, \mu, \varphi)$ through the mapping

$$ (60) \quad x^2 = \frac{\lambda^2 \mu^2 \varphi^2}{d_{\mu}^2 d_{\varphi}^2}, \quad y^2 = \frac{(\lambda^2 - d_{\varphi}^2)(\mu^2 - d_{\varphi}^2)(d_{\varphi}^2 - \varphi^2)}{d_{\lambda}^2 d_{\varphi}^2}, \quad z^2 = \frac{(\lambda^2 - d_{\mu}^2)(d_{\mu}^2 - \mu^2)(d_{\mu}^2 - \varphi^2)}{d_{\lambda}^2 d_{\mu}^2}, $$

where the positive semifocal distances are defined by

$$ (61) \quad d_{\lambda} := \sqrt{\Gamma^2 - \Xi^2}, \quad d_{\mu} := \sqrt{1 - \Xi^2}, \quad d_{\varphi} := \sqrt{1 - \Gamma^2}. $$

The geometric surfaces defined by constant values of $(\lambda, \mu, \varphi)$ correspond, respectively, to ellipsoids, one-sheet hyperboloids, and two-sheet hyperboloids. The surface of the unperturbed ellipsoid — defined by the equilibrium condition (1) — is obtained by setting $\lambda = 1$. Accordingly, the fluid domain occupied by a steady Riemann ellipsoid is given by

$$ \mathcal{D}_0 := \{(\lambda, \mu, \varphi) : d_{\mu} \le \lambda \le 1, d_{\varphi} \le \mu \le d_{\mu}, -d_{\varphi} \le \varphi \le d_{\varphi}\}. $$

The associated orthonormal basis $(e_{\lambda}, e_{\mu}, e_{\varphi})$ is particularly convenient for boundary-value problems, since the outward unit normal vector at the free surface is simply $n := e_{\lambda}$.

Because the coordinate system is orthogonal, the metric scale factors are given by the square roots of the diagonal metric components, namely

$$ (62) \quad h_{\lambda}^2 := \frac{(\lambda^2 - \mu^2)(\lambda^2 - \varphi^2)}{(\lambda^2 - d_{\mu}^2)(\lambda^2 - d_{\varphi}^2)}, \quad h_{\mu}^2 := \frac{(\mu^2 - \lambda^2)(\mu^2 - \varphi^2)}{(\mu^2 - d_{\mu}^2)(\mu^2 - d_{\varphi}^2)}, \quad h_{\varphi}^2 := \frac{(\varphi^2 - \lambda^2)(\varphi^2 - \mu^2)}{(\varphi^2 - d_{\mu}^2)(\varphi^2 - d_{\varphi}^2)}. $$

For axisymmetric configurations, it is often advantageous to work instead with oblate or prolate spheroidal coordinates. The correspondence between these coordinate systems and the general ellipsoidal geometry provides valuable simplifications in symmetry-reduced problems, as discussed in detail in the monograph by Dassios (2012) and references therein.

## APPENDIX B. GEOMETRIC RELATIONS IN CURVILINEAR COORDINATES

Let $f$ be an arbitrary vector in $\mathbb{R}^d$. Its components in two coordinate systems $\{u_i\}_{i=1,...,d}$ and $\{v_i\}_{i=1,...,d}$ are denoted by $\{f_u\}_i$ and $\{f_v\}_i$, respectively. The transformation of vector components between these two bases is given by (Feshbach & Morse, 1953)

$$ (63) \quad f_{v_i} = \frac{1}{h_{v_i}} \sum_{j=1}^{d} \frac{\partial u_j}{\partial v_i} f_{u_j}, \quad i = 1, \dots, d, $$

where $h_{v_i}$ denotes the metric scale factor associated with the coordinate $v_i$.

Let $T$ be a tensor of rank 2 whose components $T_{ij}$ are expressed in the first coordinate basis, as defined previously. The covariant transformation of this tensor on the second orthonormal curvilinear basis is given by (Feshbach & Morse, 1953)

$$ (64) \quad T_{v_i v_j} = \frac{1}{h_{v_i} h_{v_j}} \sum_{k=1}^{d} \sum_{l=1}^{d} \frac{\partial u_k}{\partial v_i} \frac{\partial u_l}{\partial v_j} T_{u_k u_l} \quad i = 1, \dots, d, \quad j = 1, \dots, d. $$

Accordingly, the gradient and Laplacian operators in an orthogonal curvilinear coordinate system $\{u_i\}_{i=1}^d$ take the standard form (Feshbach & Morse, 1953)

$$ (65) \quad \nabla := \sum_{j=1}^{d} \frac{e_{u_j}}{h_{u_j}} \frac{\partial}{\partial u_j}, \quad \nabla^2 := J^{-1} \sum_{j=1}^{d} \frac{\partial}{\partial u_j} \left( \frac{J}{h_{u_j}^2} \frac{\partial}{\partial u_j} \right), $$
where $J := \prod_{j=1}^{d} h_{u_j}$ is the Jacobian determinant.

## APPENDIX C. DEFINITION OF ELLIPSOIDAL HARMONICS

Ellipsoidal harmonics are the Lamé eigenfunctions of the Laplace operator expressed in ellipsoidal coordinates, see (65). They form an orthogonal basis and provide a natural framework for representing sufficiently regular scalar fields in ellipsoidal domains (Dassios, 2012).

We denote by $E_n^m(x)$ the ellipsoidal harmonic of degree $n$ and order $m$ of the first kind, associated with the coordinate $x$ (which may represent $\lambda$, $\mu$, or $\varphi$). These functions belong to one of four classical families of Lamé functions, given by

$$ (66) \qquad K_n^m(x) := \sum_{k=0}^{r} \hat{a}_k(p_m^K) x^{n-2k}, $$

$$ (67) \qquad L_n^m(x) := \sqrt{x^2 - d_\varphi^2} \sum_{k=0}^{r-1} \hat{b}_k(p_m^L) x^{n-1-2k}, $$

$$ (68) \qquad M_n^m(x) := \sqrt{x^2 - d_\mu^2} \sum_{k=0}^{r-1} \hat{b}_k(p_m^M) x^{n-1-2k}, $$

$$ (69) \qquad N_n^m(x) := \sqrt{x^2 - d_\varphi^2} \sqrt{x^2 - d_\mu^2} \sum_{k=0}^{r-1} \hat{c}_k(p_m^N) x^{n-2-2k}, $$

where $r := \lfloor n/2 \rfloor$ and $p_m$ denote the roots of the associated characteristic polynomials.

These characteristic polynomials $P_m(p)$ arise from tridiagonal matrix eigenvalue problems and satisfy the recurrence relation

$$ (70) \qquad P_m(p) := v^{(d)}(p)P_{m-1}(p) - v^{(u)}v^{(l)}P_{m-2}(p), $$

with $P_{-1} := 0$, $P_0 := 1$. The vectors $v^{(d)}$, $v^{(u)}$, and $v^{(l)}$ correspond to the diagonal, upper, and lower bands of the tridiagonal matrix, respectively.

For ellipsoidal harmonics, these coefficients take the explicit forms given in Dassios (2012), which we reproduce here for completeness. The roots $p_m^K$ are obtained from

$$ \begin{aligned} v_i^{(d)}(p) &:= -(d_\mu^2 + d_\varphi^2)[p - (n - 2i + 2)^2], \quad i = 1, \dots, r + 1, \\ v_i^{(u)}(p) &:= 2i(2n - 2i + 1), \quad v_i^{(l)}(p) := -d_\mu^2 d_\varphi^2(n - 2i + 1)(n - 2i + 2), \quad i = 1, \dots, r. \end{aligned} $$

The roots $p_m^L = p_m^M$ are determined from

$$ \begin{aligned} v_i^{(d)}(p) &:= -(d_\mu^2 + d_\varphi^2)[p - (n - 2i + 1)^2] + h_\mu^2(2n - 4i + 3), \quad i = 1, \dots, r + 1, \\ v_i^{(u)}(p) &:= 2i(2n - 2i + 1), \quad v_i^{(l)}(p) := -d_\mu^2 d_\varphi^2(n - 2i + 1)(n - 2i), \quad i = 1, \dots, r - 1, \end{aligned} $$

and the roots $p_m^N$ from

$$ \begin{aligned} v_i^{(d)}(p) &:= -(d_\mu^2 + d_\varphi^2)[p - (n - 2i + 1)^2], \quad i = 1, \dots, r, \\ v_i^{(u)}(p) &:= 2i(2n - 2i + 1), \quad v_i^{(l)}(p) := -d_\mu^2 d_\varphi^2(n - 2i - 1)(n - 2i), \quad i = 1, \dots, r - 1. \end{aligned} $$
The expansion coefficients $(\hat{a}, \hat{b}, \hat{c})$ appearing in (66)-(69) are obtained from recurrence relations. For $K_n^m$ and $k \ge 2$,

$$ \hat{a}_k := \frac{(d_\mu^2 + d_\varphi^2)[p - (n - 2k + 2)^2]\hat{a}_{k-1} + d_\mu^2 d_\varphi^2(n - 2k + 3)(n - 2k + 4)\hat{a}_{k-2}}{2k(2n - 2k + 1)}, $$

with $\hat{a}_0 := 1$ and $\hat{a}_1 := (d_\mu^2 + d_\varphi^2)(p - n^2)/[2(2n - 1)]$. For the Lamé functions $L_n^m$ and $M_n^m$ — since their expansion coefficients are the same — and $k \ge 3$, it yields

$$ \hat{b}_{k-1} := \frac{(d_{\mu}^{2} + d_{\varphi}^{2})[p - (n - 2k + 3)^{2}] - d_{\mu}^{2}(2n - 4k + 7)\hat{b}_{k-2}}{2(k-1)(2n - 2k + 3)} + \frac{d_{\mu}^{2}d_{\varphi}^{2}(n - 2k + 5)(n - 2k + 4)\hat{b}_{k-3}}{2(k-1)(2n - 2k + 3)}, $$

with $\hat{b}_0 := 1$ and $\hat{b}_1 := [(d_\mu^2 + d_\varphi^2)(p - (n-1)^2) - d_\mu^2(2n-1)]/[2(2n-1)]$. For the Lamé function $N_n^m$ and $k \ge 3$, it yields

$$ \hat{c}_{k-1} := \frac{(d_{\mu}^{2} + d_{\varphi}^{2})[p - (n - 2k + 3)^{2}]\hat{c}_{k-2} + d_{\mu}^{2}d_{\varphi}^{2}(n - 2k + 3)(n - 2k + 4)\hat{c}_{k-3}}{2(k-2)(2n - 2k + 3)}, $$

with $\hat{c}_0 := 1$ and $\hat{c}_1 := (d_\mu^2 + d_\varphi^2)[p - (n-1)^2]/[2(2n-1)]$. Efficient numerical strategies for computing the roots and coefficients are discussed extensively in Dassios (2012).

The ellipsoidal harmonics of the second kind, also known as exterior solutions, are defined by

$$ (71) \qquad F_n^m(x) := (2n + 1)E_n^m(x) \int_x^\infty \frac{ds}{[E_n^m(s)]^2(s^2 - d_\mu^2)^{1/2}(s^2 - d_\varphi^2)^{1/2}}. $$

Ellipsoidal harmonics provide a natural basis for expanding square-integrable functions on the ellipsoid. For any $f \in L^2(\mathcal{D}_0)$, we write

$$ (72) \qquad f(\lambda, \mu, \varphi, t) = \sum_{n=0}^{\infty} \sum_{m=0}^{2n+1} \left[ \hat{f}_{nm}^{(i)}(t) E_m^n(\lambda) + \hat{f}_{nm}^{(e)}(t) F_m^n(\lambda) \right] \mathbb{S}_n^m(\mu, \varphi), $$

where $\mathbb{S}_n^m(\mu, \varphi) := E_n^m(\mu)E_n^m(\varphi)$ denote the surface ellipsoidal harmonics.

These functions form a complete orthogonal set with respect to the surface measure

$$ (73) \qquad dS := wh_\mu h_\varphi d\mu d\varphi, \quad \text{with} \quad w := \frac{1}{\sqrt{(1 - \mu^2)(1 - \varphi^2)}}, $$

and satisfy the orthogonality relation

$$ (74) \qquad \int_{\mathcal{D}_0} \mathbb{S}_n^m(\mu, \varphi) \mathbb{S}_q^p(\mu, \varphi) dS = \varpi_n^m \delta_{nq} \delta_{mp}, $$

with

$$ (75) \qquad \varpi_n^m = \int_{\mathcal{D}_0} [\mathbb{S}_n^m(\mu, \varphi)]^2 dS. $$

Finally, the Wronskian relation between interior and exterior harmonics reads

$$ (76) \qquad W_n(\lambda) := F_n^m(\lambda) \frac{\partial}{\partial \lambda} E_n^m(\lambda) - E_n^m(\lambda) \frac{\partial}{\partial \lambda} F_n^m(\lambda) = \frac{2n+1}{(\lambda^2 - d_\mu^2)^{1/2}(\lambda^2 - d_\varphi^2)^{1/2}}. $$
## APPENDIX D. COEFFICIENTS OF THE INTEGRATED PRANDTL EQUATIONS

We give the coefficients of the integrated Prandtl equations as

$$a_0 = b_0 = -\bar{\omega}_n^m \text{Ro}, \quad a_1 = b_1 = \int [h_\lambda u_{0\lambda}]_{\partial D_0} [\mathbb{S}_n^m]^2 dS,$$

$$a_2 = \int h_\lambda^2 \left[ \left( \frac{1}{h_\mu} \frac{\partial u_{0\mu}}{\partial \mu} - 2i\sigma \text{Ro} \right) + \frac{u_{0\mu}}{2h_\mu} \frac{\partial}{\partial \mu} + \frac{u_{0\varphi}}{2h_\varphi} \frac{\partial}{\partial \varphi} \right]_{\partial D_0} [\mathbb{S}_n^m]^2 dS,$$

$$b_2 = \int h_\lambda^2 \left[ \left( \frac{1}{h_\varphi} \frac{\partial u_{0\varphi}}{\partial \varphi} - 2i\sigma \text{Ro} \right) + \frac{u_{0\mu}}{2h_\mu} \frac{\partial}{\partial \mu} + \frac{u_{0\varphi}}{2h_\varphi} \frac{\partial}{\partial \varphi} \right]_{\partial D_0} [\mathbb{S}_n^m]^2 dS,$$

$$a_3 = \int h_\lambda \left[ \frac{h_\lambda}{h_\varphi} \frac{\partial u_{0\mu}}{\partial \varphi} - 2\text{Ro} \frac{\partial z}{\partial \lambda} \right]_{\partial D_0} [\mathbb{S}_n^m]^2 dS, \quad b_3 = \int h_\lambda \left[ \frac{h_\lambda}{h_\mu} \frac{\partial u_{0\varphi}}{\partial \mu} + 2\text{Ro} \frac{\partial z}{\partial \lambda} \right]_{\partial D_0} [\mathbb{S}_n^m]^2 dS.$$

Similarly, the coefficients of the tangent stress conditions are

$$\tau_\mu = \int \left[ \frac{\partial}{\partial \lambda} \delta u_\mu^{(0)} + \frac{h_\lambda}{h_\mu} \frac{\partial}{\partial \mu} \delta u_\lambda^{(0)} \right]_{\partial D_0} \mathbb{S}_n^m dS, \quad \tau_\varphi = \int \left[ \frac{\partial}{\partial \lambda} \delta u_\varphi^{(0)} + \frac{h_\lambda}{h_\varphi} \frac{\partial}{\partial \varphi} \delta u_\lambda^{(0)} \right]_{\partial D_0} \mathbb{S}_n^m dS.$$

## REFERENCES

* Abott et al. (2016). Observation of gravitational waves from a binary black hole merger. *Phys. Rev. Lett.*, **116**(6), 061102.

* Backus, G., & Rieutord, M. (2017). Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid. *Phys. Rev. E*, **95**(5), 053116.

* Bayly, B. J. (1986). Three-dimensional instability of elliptical flow. *Phys. Rev. Lett.*, **57**(17), 2160-2163.

* Bayly, B. J., Orszag, S. A., & Herbert, T. (1988). Instability mechanisms in shear-flow transition. *Annu. Rev. Fluid Mech.*, **20**, 359-391.

* Bloch, A. M., Krishnaprasad, P. S., Marsden, J. E., & Ratiu, T. S. (1994). Dissipation induced instabilities. *Ann. Inst. H. Poincaré Anal. Non Linéaire*, **11**(1), 37-90.

* Bryan, G. (1889). The waves on a rotating liquid spheroid of finite ellipticity. *Philos. Trans. R. Soc. London*, **180**, 187-219.

* Cartan, E. (1922). Sur les petites oscillations d'une masse fluide. *Bull. Sci. Math.*, **46**, 317-352 and 356-369.

* Chandrasekhar, S. (1962). On the point of bifurcation along the sequence of the Jacobi ellipsoids. *Astrophys. J.*, **136**, 1048-1068.

* Chandrasekhar, S. (1965). The Equilibrum and the Stability of the Riemann Ellipsoids. I. *Astrophys. J.*, **142**, 890-961.

* Chandrasekhar, S. (1966). The Equilibrum and the Stability of the Riemann Ellipsoids. II. *Astrophys. J.*, **145**, 842-877.

* Chandrasekhar, S. (1969). *Ellipsoidal Figures of Equilibrium*. New Haven: Yale University Press.

* Chandrasekhar, S. (1970). Solutions of two problems in the theory of gravitational radiation. *Phys. Rev. Lett.*, **24**(11), 611-615.

* Comins, N. (1979). On secular instabilities of rigidly rotating stars in general relativity – I. Theoretical formalism. *Mon. Not. R. Astron. Soc.*, **189**(2), 233-253.

* Dassios, G. (2012). *Ellipsoidal Harmonics: Theory and Applications*. Cambridge: Cambridge University Press.
Dedekind, R. (1861). Zusatz zu der vorstehenden Abhandlung. *J.f. Reine und Angew Math.*, **58**, 217–228.

Feshbach, H., & Morse, P. M. (1953). *Methods of Theoretical Physics*. Cambridge: Cambridge University Press.

Greenspan, H. P. (1965). On the general theory of contained rotating fluid motions. *J. Fluid Mech.*, **22**(3), 449–462.

Jacobi, C. G. J. (1839). Note von der geodätischen Linie auf einem Ellipsoid und den verschiedenen Anwendungen einer merkwürdigen analytischen Substitution. *J. Reine Angew. Math.*, **19**, 309–313.

Kerswell, R. R. (2002). Elliptical instability. *Annu. Rev. Fluid Mech.*, **34**(1), 83–113.

Kirillov, O. N. (2021). *Nonconservative Stability Problems of Modern Physics* (2nd rev. and exten. edition ed., Vol. 14). Berlin, Boston: De Gruyter.

Labarbe, J. (2021). *Radiative and diffusive instabilities in moving fluids* (PhD thesis). Northumbria University.

Landau, L. D., & Lifshitz, E. M. (1987). *Fluid Mechanics*. Oxford: Pergamon Press.

Landman, M. J., & Saffman, P. G. (1987). The three-dimensional instability of strained vortices in a viscous fluid. *Phys. Fluids*, **30**(8), 2339–2342.

Lebovitz, N. R., & Lifschitz, A. (1996). Short-wavelength instabilities of Riemann ellipsoids. *Philos. Trans. R. Soc. A*, **354**(1709), 927–950.

Liouville, J. (1851). Mémoire sur les figures ellipsoïdales à trois axes inégaux, qui peuvent convenir à l'équilibre d'une masse liquide homogène, douée d'un mouvement de rotation. *J. Math. Pures Appl.*, **16**, 241–254.

Lynden-Bell, D., & Ostriker, J. P. (1967). On the stability of differentially rotating bodies. *Mon. Not. R. Astron. Soc.*, **136**(3), 293–310.

Maclaurin, C. (1742). *A Treatise of Fluxions: In Two Books. 1* (Vol. 72). Ruddimans.

Meyer, C. O. (1842). De aequilibrii formis ellipsoidicis. *J. Reine Angew. Math.*, **24**, 44–59.

Moore, D. W., & Saffman, P. G. (1975). The instability of a straight vortex filament in a strain field. *Proc. R. Soc. Lond. A*, **346**(1646), 413–425.

Morrison, P. J., Lebovitz, N. R., & Biello, J. A. (2009). The Hamiltonian description of incompressible fluid ellipsoids. *Ann. Phys.*, **324**(8), 1747–1762.

Pierrehumbert, R. T. (1986). Universal short-wave instability of two-dimensional eddies in an inviscid fluid. *Phys. Rev. Lett.*, **57**(17), 2157–2159.

Poincaré, H. (1885). Sur l'équilibre d'une masse fluide animée d'un mouvement de rotation. *Acta Math.*, **7**, 259–380.

Riemann, B. (1860). Ein Beitrag zu den Untersuchengen über die Bewegung eines flüssigen gleichartigen Ellipsoides. *Gött. Abh.*, IX, 3–36.

Roberts, P. H., & Stewartson, K. (1963). On the stability of a Maclaurin spheroid of small viscosity. *Astrophys. J.*, **137**, 777–790.

Schutz, B. F., & Friedman, J. L. (1978a). Langrangian perturbation theory of nonrelativistic fluids. *Astrophys. J.*, **221**, 937–957.

Schutz, B. F., & Friedman, J. L. (1978b). Secular instability of rotating Newtonian stars. *Astrophys. J.*, **222**, 281–296.

Thompson, W., & Tait, P. G. (1867). *A Treatise on Natural Philosophy – I*. Oxford: Clarendon Press.