Inertial Modes in Near-Spherical Geometries

TL;DR

This paper introduces a spectral numerical method to accurately compute inertial modes in near-spherical geometries, including oblate spheroids and triaxial ellipsoids, by solving a polynomial eigenvalue problem.

Contribution

It develops a spectral discretisation approach using spherical harmonics and Gegenbauer polynomials to solve the Poincare equation for inertial modes in near-spherical containers.

Findings

Accurately recovers inertial modes of oblate spheroids to machine precision.

Demonstrates convergence for triaxial ellipsoids with slightly deviating boundaries.

Provides a flexible method for near-spherical geometries in geophysical fluid dynamics.

Abstract

We propose a numerical method to compute the inertial modes of a container with near-spherical geometry based on the fully spectral discretisation of the angular and radial directions using spherical harmonics and Gegenbauer polynomial expansion respectively. This allows to solve simultaneously the Poincare equation and the no penetration condition as an algebraic polynomial eigenvalue problem. The inertial modes of an exact oblate spheroid are recovered to machine precision using an appropriate set of spheroidal coordinates. We show how other boundaries that deviate slightly from a sphere can be accommodated for with the technique of equivalent spherical boundary and we demonstrate the convergence properties of this approach for the triaxial ellipsoid.

Equations

$${\partial }_t\mathbf{u}+2\hat{\mathbf{z}}\times \mathbf{u}=-\nabla p,$$

$$\mathbf{u}$$

$$p$$

$$-{\lambda }^2{\nabla }^2{p}_0+(\hat{\mathbf{z}}\cdot \nabla {)}^2{p}_0=0,$$

$$-{\lambda }^2\hat{\mathbf{n}}\cdot \nabla p_0+i\lambda (\hat{\mathbf{z}}\times \hat{\mathbf{n}})\cdot \nabla p_0+(\hat{\mathbf{n}}\cdot \hat{\mathbf{z}})(\hat{\mathbf{z}}\cdot \nabla p_0){|}_{\partial \mathcal{V}}=0.$$

$${\mathbf{u}}_0\equiv \frac{1}{1-{\lambda }^2}\frac{1}{2}\left(\hat{\mathbf{z}}\times \nabla p_0-i\lambda \nabla p_0+\frac{i}{\lambda }(\hat{\mathbf{z}}\cdot \nabla p_0)\hat{\mathbf{z}}\right).$$

$$(A_0{\lambda }^0+A_1{\lambda }^1+A_2{\lambda }^2)\mathbf{x}=0.$$

$$\mathbf{u}=\nabla \times \nabla \times (P\mathbf{r})+\nabla \times (T\mathbf{r}),$$

$$p_0(r,\theta ,\varphi )=\sum _{\ell =0}^L\sum _{m=-\ell }^{\ell }{p}_{\ell ,m}(r)Y_{\ell }^m(\theta ,\varphi ),$$

$${\left\{p_{|m|,m},p_{|m+2|,m},p_{|m+4|,m},\dots \right\}}^{\text{T}}$$

$${\left\{p_{|m+1|,m},p_{|m+3|,m},p_{|m+5|,m},\dots \right\}}^{\text{T}}$$

$$Y_{\ell }^m(\theta +\pi ,\varphi )=(-1{)}^{\ell }{Y}_{\ell }^m(\theta ,\varphi ),$$

$$p_0(-r,\theta ,\varphi )=p_0(r,\pi -\theta ,\varphi +\pi ),$$

$$p_{\ell ,m}(-r)=(-1{)}^{\ell }{p}_{\ell ,m}(r).$$

$$R(\theta ,\varphi )$$

$$=R_0\left(1+\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{ϵ}_{\ell ,m}{Y}_{\ell }^m(\theta ,\varphi )\right).$$

$$\zeta {|}_R=\zeta {|}_{R_0}+\frac{d\zeta }{dr}{|}_{R_0}(R-R_0)+\mathcal{O}({ϵ}^2)=0.$$

$$\sum _{\ell =0}^L\sum _{m=-\ell }^{\ell }\left({\zeta }_{\ell ,m}{|}_{R_0}{Y}_{\ell }^m+\sum _{{\ell }^{\prime }=0}^L\sum _{m^{\prime }=-{\ell }^{\prime }}^{{\ell }^{\prime }}\frac{d{\zeta }_{\ell ,m}}{dr}{|}_{R_0}{ϵ}_{{\ell }^{\prime },m^{\prime }}{Y}_{\ell }^m{Y}_{{\ell }^{\prime }}^{m^{\prime }}\right)+\mathcal{O}({ϵ}^2)=0.$$

$$x^2+\frac{y^2}{(1-f^2)}+\frac{z^2}{(1-e^2)}=a^2,$$

$$x$$

$$y$$

$$z$$

$$\beta =\sqrt{\frac{1-{\lambda }^2}{{\lambda }^2}}\text{where}\lambda =\frac{\omega }{2\Omega }$$

$${\nabla }^2{p}_0-\frac{1}{{\lambda }^2}(\hat{\mathbf{z}}\cdot \nabla {)}^2{p}_0$$

$$={\left(\frac{{\partial }^2{p}_0}{\partial X^2}\right)}_{Y,Z}+{\left(\frac{{\partial }^2{p}_0}{\partial Y^2}\right)}_{X,Z}+{\left(\frac{{\partial }^2{p}_0}{\partial Z^2}\right)}_{X,Y}$$

$$=0$$

$$\Phi \equiv \frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}-1=0,$$

$$\mathbf{n}=\frac{x}{a^2}\hat{\mathbf{x}}+\frac{y}{b^2}\hat{\mathbf{y}}+\frac{z}{c^2}\hat{\mathbf{z}}.$$

$$-\left(\frac{iy\lambda }{b^2}+\frac{x{\lambda }^2}{a^2}\right){\left(\frac{\partial p_0}{\partial x}\right)}_{y,z}+\left(\frac{ix\lambda }{a^2}-\frac{y{\lambda }^2}{b^2}\right){\left(\frac{\partial p_0}{\partial y}\right)}_{x,z}+\left(\frac{z}{c^2}-\frac{z{\lambda }^2}{c^2}\right){\left(\frac{\partial p_0}{\partial z}\right)}_{x,y}=0.$$

$$[\mathbf{J}]=\left[\begin{array}{ccc}{\left(\frac{\partial X}{\partial u}\right)}_{v,w}& {\left(\frac{\partial X}{\partial v}\right)}_{u,w}& {\left(\frac{\partial X}{\partial w}\right)}_{u,v}\\ {\left(\frac{\partial Y}{\partial u}\right)}_{v,w}& {\left(\frac{\partial Y}{\partial v}\right)}_{u,w}& {\left(\frac{\partial Y}{\partial w}\right)}_{u,v}\\ \beta {\left(\frac{\partial Z}{\partial u}\right)}_{v,w}& \beta {\left(\frac{\partial Z}{\partial v}\right)}_{u,w}& \beta {\left(\frac{\partial Z}{\partial w}\right)}_{u,v}\end{array}\right]$$

$$\left[\begin{array}{ccc}{\left(\frac{\partial p_0}{\partial x}\right)}_{y,z}& {\left(\frac{\partial p_0}{\partial y}\right)}_{x,z}& {\left(\frac{\partial p_0}{\partial z}\right)}_{x,y}\end{array}\right]=\left[\begin{array}{ccc}{\left(\frac{\partial p_0}{\partial u}\right)}_{v,w}& {\left(\frac{\partial p_0}{\partial v}\right)}_{u,w}& {\left(\frac{\partial p_0}{\partial w}\right)}_{u,v}\end{array}\right]\cdot [\mathbf{J}{]}^{-1}$$

$$x$$

$$y$$

$$z$$

$$k=\sqrt{a^2+\frac{c^2}{{\beta }^2}}=a\sqrt{\frac{1-e^2{\lambda }^2}{1-{\lambda }^2}},$$

Full text

keywords : Planetary interiors; Core; Numerical solutions.

1. Introduction

Coriolis forces in a rotating fluid support oscillatory motions known as inertial waves. When the fluid is rotating inside a boundary, as is the case for many astrophysical objects including our planet, inertial modes can exist. Experimentally, these modes can be excited mechanically by means of libration of the bounding surface (Aldridge and Toomre, 1969), by precession (Malkus, 1968), by tidal forces (Morize et al., 2010), or by the differential rotation of a solid inner core (Kelley et al., 2007). Although the first theoretical studies of inertial modes date back towards the end of the 19th century, with the works of Thomson (1880); Poincaré (1885) and Bryan (1889), the role of inertial modes in the dynamics of rotating stars, planets and other astrophysical bodies remains largely unexplored.

With an eye towards the dynamics of bodies deformed by rotation and tides, it is relevant to obtain the solutions for inertial mode oscillations of rotating fluids within near-spherical boundaries such as spheroids or triaxial ellipsoids. In the analytical realm, implicit solutions for inviscid inertial eigenmodes in a spheroid were first obtained by Bryan (1889) using bi-spheroidal coordinates. In a triaxial ellipsoid the first relevant work is that of Hough (1895), who employed Lamé functions. Explicit solutions in the spheroid for few low degree modes were obtained much later by Kudlick (1966). Explicit solutions for all inertial eigenmodes in a spheroid were not available until the relatively recent work of Zhang et al. (2004), who also included the first order viscous corrections. In the case of a triaxial ellipsoid, Vantieghem (2014) provided an algorithm to construct analytical inviscid solutions that are linear and quadratic in the cartesian coordinates, with the possibility to include higher order solutions. The long standing question about the completeness of the inertial modes to represent any smooth fluid flow within a full rotating ellipsoid was answered affirmatively by the work of Backus and Rieutord (2017) and Ivers (2017a). It turns out that inertial modes are also complete when the rotation axis is arbitrary with respect to the principal axes of the ellipsoid (Ivers, 2017b).

In the numerical realm, studies of inertial eigenmodes in ellipsoids are based predominantly on finite element methods (e.g. Chan et al., 2010; Cébron et al., 2010; Vantieghem, 2014). The works of Schmitt (2006) and Schmitt and Jault (2004) are perhaps the only ones using a set of oblate spheroidal coordinates. All these numerical studies usually focus on the time evolution of the flow field inside the spinning ellipsoid, which is in turn subject to some kind of additional mechanical forcing (e.g. libration) to excite inertial modes.

The purpose of this paper is the numerical computation of the inertial modes of a fluid rotating inside an ellipsoidal boundary numerically using a fully spectral discretisation. The radial direction is discretised using a relatively recent technique which employs Chebyshev and Gegenbauer polynomials (Olver and Townsend, 2013). This has the advantage to result in a sparse matrix representation of differential operators. The discretisation in the angular directions is carried out using a spherical harmonics decomposition. For the special cases of an exactly spherical or oblate spheroidal container, the spectral decomposition can be used to recover the inertial modes with an accuracy that is limited only by numerical precision. In order to treat the triaxial case, we extend the equivalent spherical domain technique introduced by Smith (1974) and substantiate it with an analysis of numerical convergence. This makes it possible, at least in principle, to compute the inertial modes for any boundary that deviates only slightly from sphericity. We also demonstrate a technique to recover semi-analytical solutions in an oblate spheroid and in a triaxial ellipsoid using systems of bi-spheroidal and bi-ellipsoidal coordinates. This allows us to cross-validate both approaches.

Our numerical method paves the way for future extensions of the model to accommodate features that the analytical solutions are not capable of, such as the inclusion of viscosity, magnetic effects or stratification and, most prominently, the effect of boundary topography.

The inviscid problem that we present in this study is an interesting instance of a polynomial eigenvalue problem with a non-standard boundary condition. This problem also serves as a concrete example demonstrating the potential of a more general method that we have developed to treat similar problems consisting of general systems of Partial Differential Equations (PDEs) in near-spherical geometries.

The paper is structured as follows. In Sec. 2, we review the mathematical description of inertial waves inside a cavity before we explain the methods used to solve the resulting equations numerically. In Sec. 3, we demonstrate the validity of our approach by treating the case of a flow within an ellipsoidal boundary. Sec. 4 provides a discussion of how the method can easily be extended to more complex systems. The Appendix at the end of this paper provides details on how to compute the analytical solutions to the inertial modes as well as other useful technical information used throughout the present paper.

2. Method

2.1. Mathematical description of inertial modes

We consider an incompressible, homogeneous and inviscid fluid contained in a rigid cavity which rotates with angular velocity, $\Omega$ around the z-axis, ${\bf\hat{z}}$. In the frame attached rigidly to the rotating cavity and using $1/\Omega$ as the unit of time, the momentum balance reads

[TABLE]

where $\bf u$ represents the flow velocity and $p$ the reduced pressure (which includes the centrifugal and gravitational potential). We are assuming that the fluid velocity is small compared to the maximum tangential velocity of the cavity so that we can discard the non-linear $({\bf u}\cdot\nabla){\bf u}$ term. We further assume that the fluid motion is periodic in time:

[TABLE]

where $*$ denotes the complex conjugate. After some algebraic manipulation, Eq. (1) reduces to the so-called Poincaré equation for the pressure amplitude $p_{0}$ (Poincaré, 1885):

[TABLE]

where $\lambda$ is the half-frequency of the motion. We employ a no-penetration boundary condition for the velocity ${\bf u}\cdot{\bf\hat{n}}|_{\partial\mathcal{V}}=0$ at the bounding surface $\partial\mathcal{V}$, which translates into the following condition for the pressure field $p_{0}$ (Greenspan, 1968):

[TABLE]

The velocity field can be recovered from the solution for $p$ using the following expression:

[TABLE]

Finding analytical solutions to the problem just described proceeds in two steps. The first one is a change of coordinates that reduces Eq. (4) to a Laplace equation. The Laplace operator and the boundary condition must both be separable in these coordinates, something that is only possible in a limited number of coordinates. The second step is to inject the formal solution of the Laplace equation into the boundary condition and solve for $\lambda$. We give a detailed illustration for the cases of the sphere, the oblate spheroid and the triaxial ellipsoid in Appendix A. The solutions obtained will serve us to validate the numerical results of Sec. 3.

The system of Eqs. (4-5) can be cast numerically as an eigenvalue problem that is polynomial in the eigenvalue $\lambda$. There are methods to solve algebraic systems of this type. In order to use these, we first turn the above differential problem into an algebraic problem of the form

[TABLE]

where the $A_{i}$’s are square matrices, Fig. 1 gives a visual representation of these matrices. Sec. 2.2 explains how these can be obtained from Eqs. (4) and (5).

The direct numerical resolution of the Poincaré equation for the motion of a rotating fluid is not very common practice. The reason being that this formulation is not valid when viscosity is taken into account. When this is the case, one uses an approach based on the following decomposition of the (solenoidal) velocity field, $\boldsymbol{u}$ :

[TABLE]

where $P(\boldsymbol{r})$ and $T(\boldsymbol{r})$ are the poloidal and toroidal scalar fields respectively. After inserting the above in the momentum Eq. (1), one can write two independent scalar equations for $P$ and $T$ by separately considering the radial projections of the curl of the momentum equation and the curl of its curl (Rieutord and Valdettaro, 1997). This latter approach can be used as well when dealing with the inviscid problem so long as the domain of integration encloses the origin of coordinates, i.e., containers that have the topology of the spherical shell are prohibited (Rieutord et al., 2000).

Both methods have their own advantages and drawbacks. The formulation in $P$ and $T$, when discretised, leads to a linear eigenvalue problem and so allows the usage of state of the art codes designed for viscous computations with minor modifications. On the other hand the formulation in $p$ leads to a quadratic eigenvalue problem but the scalar nature of the Poincaré equation makes it straightforwardly adaptable to oblate spheroidal coordinates as discussed in Appendix C.

In the rest of this section, we give the details of the discretisation method assuming the formulation in terms of the Poincaré equation. Both formulations are used to obtain the results of Sec. 3.

2.2. Spectral discretisation

2.2.1. Angular discretisation

When working in spherical coordinates $\{r,\theta,\phi\}$, the Partial Differential Equation (PDE) Eq. (4) can be reduced to a set of coupled Ordinary Differential Equations (ODEs) using the method of spherical harmonics decomposition. The pressure field amplitude $p_{0}$ is written as the following truncated series

[TABLE]

where $L$ is an integer chosen as large as possible. The numerical task consists in finding an approximate expression for the radial functions $\{p_{\ell,m}(r)\}$. In the end, the resulting system will be turned into a fully algebraic (matrix) problem by discretisation of the radial direction.

We now review symmetry considerations which allow to decouple the problem further. The second term of Eq. (4) induces a coupling between each harmonic component ($\ell,m$) and its closest neighbours with a degree of the same parity ($\ell\pm 2,m$). The origin of the coupling traces back to the presence of the Coriolis force in the momentum equation. Appendix B gives the analytical expression of this term as well as the last two terms of Eq. (5) in spherical coordinates illustrating the case of a spherical container $({\bf\hat{n}}={\bf\hat{r}})$. The boundary condition induces no further coupling of the spherical harmonics components in that case. Moreover, components with different azimuthal $m$ numbers remain uncoupled, an important fact that is carried over to the spheroidal container case (or indeed any axisymmetric container) and allows to solve for modes with different $m$ independently. This classification of modes by their azimuthal $m$ number is no longer applicable for (non-axisymmetric) ellipsoidal containers.

The nature of the coupling of the $\ell$-numbers reflects the decoupling of modes with a pressure profile that is symmetric by reflection across the equatorial plane and those that are anti-symmetric. By using the property, $Y_{\ell}^{m}(\theta,\phi)=(-1)^{\ell+m}Y_{\ell}^{m}(\pi-\theta,\phi)$, one can show that modes with $(\ell+m)$ even are symmetric modes while those with $(\ell+m)$ odd are anti-symmetric modes.

Following the above considerations, the layout of the radial functions $p_{\ell,m}(r)$ (which will become the eigenvectors) prior to their discretisation in $r$ reads:

[TABLE]

2.2.2. Radial discretisation

For the radial discretisation we follow the method proposed by Olver and Townsend (2013) based on polynomial expansion on the (truncated) basis of Chebyshev polynomials for the radial functions and Gegenbauer polynomials for their radial derivatives. Its main advantage compared to the more common collocation methods is that the matrices that result from the discretisation of differential operators are small in size and sparse, as opposed to the collocation methods that lead to small but dense matrices for equal resolution. Finite difference methods lead to large sparse matrices and do not provide exponential convergence. The sparsity of the resulting matrices in the method of Olver and Townsend (2013) is particularly handy if very high numerical resolution is needed. The only drawback is that it is slightly less straightforward to implement. The use of Chebyshev polynomials is especially convenient due to the possibility to compute coefficients using a Fast Fourier Transform.

We deal with bounding surfaces that deviate slightly from sphericity, with $R$ representing the mean radius of the cavity. In its original form, the spectral method of Olver and Townsend (2013) is designed to deal with differential equations in a single spatial direction $x$, limited to the interval $x\in[-1,1]$. Special care is needed when we consider the radial domain $[0,R]$ of the fluid. One might naively map the radial fluid domain $[0,R]$ to the $[-1,1]$ interval and enforce the appropriate boundary condition (the regularity condition at the centre of coordinates) at $x=-1$. This is however a poor choice because the resulting radial functions are not necessarily compatible with the intrinsic symmetries of the spherical harmonics. A better choice is to extend the fluid domain to $[-R,R]$ and map it to the interval $[-1,1]$ to match the natural domain of the Chebyshev polynomials. Since the spherical harmonics satisfy

[TABLE]

the following identity holds for the pressure field $p_{0}({\bf r})$:

[TABLE]

wich implies

[TABLE]

In our approach each function $p_{\ell,m}(r)$ will be represented as a linear combination of Chebyshev polynomials, $T_{k}(x)$ : $p_{\ell,m}(r)=\sum_{k=0}^{N}p_{\ell,m}^{k}T_{k}(r/R)$. As the parity of each Chebyshev polynomial can be read from that of its integer index, $k$, the condition Eq. (13) can be enforced by only keeping the coefficients $p_{\ell,m}^{k}$ that have the correct parity.

Boundary conditions are enforced by row replacement (Olver and Townsend, 2013) after expansion on spherical harmonics.

One important detail that we need to point out is the fact that, when treating a given problem, one should always try to eliminate inverse powers of $r$ from the starting expression. This in order to keep the Chebyshev expansion of any prefactor in the equation as small as possible which greatly reduces the bandwidth of every block matrix (Olver and Townsend, 2013). For example, the first step in dealing with Eq. (4) is to multiply it by $r^{2}$.

The usage of spherical harmonics can be extended to equations written in oblate spheroidal coordinates. An explanation of how to do this is given in Appendix C.

Fig.1a shows an example of the final set of matrices making up the discretised version of Eqs. (4-5) for a spherical container using spherical coordinates. The rows dedicated to the enforcement of the boundary condition can be seen at the bottom of each matrix. Fig. 1b shows the same set of matrices for a spheroidal container using the oblate spheroidal coordinates. Notice that the number of blocks is larger compared to the spherical case of Fig.1a. This illustrates the more extensive coupling between the spherical harmonics components $p_{\ell,m}$. Each individual block also has a wider diagonal representing the increased number of coefficients needed in the polynomial expansion of each component.

2.3. Near-spherical boundaries

When dealing with non-spherical domains, the use of spherical harmonic decomposition must be adapted. We now present a general way to do so. In essence, it amounts to treat the position of the physical boundary as the result of a series expansion around the sphere. For this reason, this method is only efficient when dealing with near-spherical boundaries. The radius of the physical boundary becomes a function of $\theta$ and $\phi$ which we parametrise as

[TABLE]

The technique will work better if $\epsilon(\theta,\phi)\ll 1$ and the coefficients $\epsilon_{\ell,m}$ are small. The trick is to transform a boundary condition at $r=R(\theta,\phi)$, the physical boundary, into a condition on an equivalent spherical domain, the computational boundary (Smith, 1974). For example, suppose that one wishes to impose the Dirichlet boundary condition on a single scalar field $\zeta(r,\theta,\phi)$. At the physical boundary, the condition is simply $\zeta|_{R}=0$. Assuming a small deviation from the spherical boundary, one can write

[TABLE]

Using Eq. (14), one finds the following expansion in spherical harmonics

[TABLE]

The product of spherical harmonics can then be reduced to a sum using 3-j symbols. which will generally consist of $2\ell^{\prime}+1$ terms featuring spherical harmonics ranging from $Y_{|\ell-\ell^{\prime}|}^{m}$ to $Y_{\ell+\ell^{\prime}}^{m}$. Those terms effectively couple together these harmonics with $Y_{\ell}^{m}$ from the original expression. Imposing the Dirichlet boundary condition – in its series expansion form – is therefore equivalent to equating each spherical harmonics component of Eq. (16) to zero independently. It is possible to increase the precision of this scheme by keeping higher powers of $\epsilon$ in the Taylor expansion Eq. (15). This will bring out products of spherical harmonics that include more factors resulting in more extensive coupling in the spherical harmonics.

In general, the expression of the boundary condition may also contain vectors. It is actually the case for Eq. (5) which features the pressure gradient, ${\bf\nabla}p$, and the normal vector ${\bf\hat{n}}$. In such cases, the procedure is similar to what we have described, except that we now have to expand both ${\bf\hat{n}}$ and ${\bf\nabla}p$ onto a basis of vector spherical harmonics. The idea remains the same although the handling of symbolic expressions can become quite tedious. In practice, we use TenGSHui, a dedicated Mathematica package (Trinh, in prep.), for these manipulations.

The technique described above can be used to deal with ellipsoidal boundaries with small eccentricities. Such boundaries can be parametrised as

[TABLE]

where $a^{2}$ stands for the semi-major axis of the ellipsoid in the equatorial plane and $e^{2}$ and $f^{2}$ (both taken as $\ll 1$) respectively represent the polar and equatorial squared eccentricities in the principal axes of symmetry of the ellipsoid.

We have already argued that the spherical harmonics components with different azimuthal $m$ numbers decouple in the special case when $f^{2}=0$ (see Sec. 2.2.1). In the general (triaxial) case the spherical harmonics expansion will contain terms with different $m$ numbers. The extra amount of coupling greatly increases the computational cost when one increases the angular resolution, i.e the upper-bound on the $\ell$ number (see Eq. (9)). The considerations on the decoupling of the equatorially symmetric and antisymmetric modes, however remains valid with the effect that the shape function, $\epsilon_{\ell,m}$, (see Eq. (14)) contains only even $\ell$ and $m$ numbers.111the explicit expression of each coefficient $\epsilon_{\ell,m}$ is given in Appendix D. This ensures that modes with even $m$ in an axisymmetric container will have a triaxial counterpart whose expansion contains only even $m$ numbers. The same is true of axisymmetric modes with an odd $m$ number. Fig. 1c shows the matrices involved in the discretised versions of Eqs. (4-5) for a triaxial ellipsoid using the method of series expansion of the boundary condition. Notice the larger size of the matrices compared to Fig. 1a-1b due to the extra-coupling between the coefficients with different $m$ numbers in the boundary conditions.

2.4. Solver

We solve the algebraic problem resulting from the discretisation of Eqs. (4-5) using the SLEPc package (Hernandez et al., 2005). SLEPc is an open-source software built on top of PETSc, another open-source package dedicated to efficiently solve large matrix equations (Balay et al., 1997, 2017a, 2017b). SLEPc is the part that deals with eigenvalue problems. As the matrices involved in our problem can be quite large, it is impractical to solve for the whole eigenspectrum. Instead, we perform a shift-and-invert spectral transformation of the original problem which greatly improves the efficiency by limiting the computation of the spectrum to a small region around a given target eigenvalue (the pivot) provided as a guess by the user. Our own usage of SLEPc was greatly inspired by the one described by Vidal and Schaeffer (2015). The only difference being that we make use of the built-in PEP solver (Campos and Roman, 2016). The shift and invert method was also used by Rieutord and Valdettaro (1997) in a similar context.

3. Results

We apply the methods of the previous section to the computation of the inertial modes in an ellipsoid. We start with the spherical case before considering the oblate spheroid and finally the triaxial ellipsoid.

3.1. Sphere

Fig. 2 shows meridional cuts in the spatial velocity profiles of the first twelve inertial modes in order of increasing spatial complexity. Each mode is identified by its (single) azimuthal $m$ number, its maximum $\ell$ number, $\bar{\ell}$, and its physical frequency, $\omega=2\lambda$. The first of these modes which corresponds to $\bar{\ell}=2$ and $m=1$ is of particular importance and is called the spin-over mode. In the rotating reference frame, this mode corresponds to a solid body rotation around an axis, itself rotating within the equatorial plane with unit angular frequency. The computation was carried out using a set of Chebyshev polynomials with maximum degree $N=20$ and a maximum degree of spherical harmonics $L=10$. The values of the frequencies that we compute agree with the analytical predictions of Eq. (33) to numerical precision.

3.2. Oblate spheroid

Fig. 3 is analogous to Fig. 2 for a spheroidal container of squared eccentricity $e^{2}=\frac{9}{16}$. The resolution used is $N=20$, $L=10$. The values of the eigenfrequencies again agree with the analytical prediction of Eq. (33) to numerical precision.

These results were computed using the specially tailored set of oblate spheroidal coordinates (Appendix C).

The inertial modes of an oblate spheroid can also be computed using the technique of equivalent spherical boundary exposed in Sec. 2.3. Fig. 4 shows the difference between the numerical and analytical eigenvalues for the inertial modes of lowest maximum degree, $\bar{\ell}$. The integer, $n$, represents the maximum order of the Taylor expansion. The error made on the computed eigenfrequencies scales as some near-integer power of the squared eccentricity. We see that the precision on the numerical eigenvalues saturates to second order convergence for all modes except the spin-over ($\bar{\ell}=2,m=1$). We attribute this to the fact of working in finite arithmetic precision. The simple geometry of the spin-over mode somehow makes this issue less important. Interestingly, the error for $n=2$ for this mode already scales as $\sim(e^{2})^{4}$ so that nothing is gained by using $n=3$.

3.3. Triaxial ellipsoid

The Taylor expansion of the boundary condition is performed in the two eccentricity parameters $e^{2}$ and $f^{2}$ which are considered of the same order of magnitude. Fig. 5 shows the residual error on three eigenvalues in the special case where $e^{2}=f^{2}$. The error scales as expected for Taylor expansion to order $n=1$ and $n=2$ respectively.

Fig. 6a shows the error on the evaluation of the spin-over frequency as a function of $e^{2}$, this time setting $f^{2}=0.01$. The dip around $e^{2}\sim 0.01$ corresponds to the point where the error changes sign. This feature is also present on Fig. 6b which shows the error made by directly expanding the analytical formula of Eq. (40) in powers of $e^{2}$ and $f^{2}$. Decreasing the value of $e^{2}$ below the threshold $e^{2}=f^{2}$ has no effect on the accuracy of the solution which becomes limited by the value of $f^{2}$, hence the saturation observed on the left side of both graphs. Above the threshold, one recovers the expected power laws in $e^{2}$ indicating the good numerical convergence properties of the method.

4. Discussion

We have presented a new numerical method to compute the inviscid inertial modes of a rotating near-spherical container. This is based on the fully spectral discretisation of the angular and radial directions using spherical harmonics in the angular directions and Chebyshev as well as Gegenbauer polynomials in the radial direction. We have shown how the method can be used to solve the Poincaré equation and the no penetration condition simultaneously as an algebraic polynomial eigenvalue problem and we have employed it to recover the inertial modes of a sphere and an oblate spheroid numerically with an accuracy limited by machine precision only (Sec. 3.1 and 3.2). We have also shown how the method of equivalent spherical boundary introduced by Smith (1974) can be used to compute the inertial modes inside a boundary that deviates only slightly from a sphere. We substantiated this technique with an analysis of its numerical convergence for the inertial modes of lowest degree both in an oblate axisymmetric ellipsoid and a triaxial ellipsoid with small eccentricities (Sec. 3.2 and 3.3).

The demonstration of the well posedness and numerical convergence of the method exposed in this paper in the context of inviscid inertial modes is an important first step towards its application to other problems for which there exists no analytical solution. It is not limited to the study of the sole Poincaré equation and can, in principle, be applied to the resolution of any system of differential equations in near-spherical geometry. Such problems are ubiquitous to the fields of planetology and geophysical/astrophysical fluid dynamics. Future foreseen applications include the study of the impact of inertial modes on the global rotation of a planet via coupling with the Liouville equation of rotational dynamics. This will be the subject of a future paper by the authors (Triana et al., 2018). This study also investigates the effect of viscosity. The effects of stratification and magnetisation of the liquid core are also among possible applications. Such studies are natural extensions to the formalism of the present paper via inclusion of the equations for heat diffusion and magnetic induction. With minor modifications, the spectral discretisation described here can also accommodate multi-layered physical configurations. This would allow to couple the dynamics of the liquid core to that of a visco-elastic and self gravitating mantle with topography. These studies are long-terms endeavours that are currently being investigated by the authors.

Acknowledgement

The authors would like to thank J. Vidal and N. Schaeffer for the useful discussions in the first stage of development of the method which lead to the present paper. We would also like to thank S. Olver and the PETSc and SLEPc support teams for their help with technical details in the making of this work. The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Advanced Grant agreement No 670874).

Appendix A Analytical solutions

The idea is to introduce a new (frequency-dependent) set of coordinates $\{u,v,w\}$

[TABLE]

with the rescaling factor

[TABLE]

so as to turn the Poincaré equation (Eq. (4)) into a Laplace equation:

[TABLE]

The set of coordinates $\{u,v,w\}$ should be chosen so that the Laplace equation is separable and the no-penetration boundary condition (Eq. (5)) is easy to implement.

From the implicit expression for the general ellipsoid in $\mathbb{R}^{3}$

[TABLE]

the (unnormalised) normal vector ${\bf n}\sim{\bf\nabla}\Phi$ at a surface point $x{\bf\hat{x}}+y{\bf\hat{y}}+z{\bf\hat{z}}$ writes

[TABLE]

Therefore, the boundary condition Eq. (5) in Cartesian coordinates reads

[TABLE]

evaluated at a surface point $x{\bf\hat{x}}+y{\bf\hat{y}}+z{\bf\hat{z}}$. To rewrite the boundary condition Eq. (25) in terms of the new set of coordinates $\{u,v,w\}$, we compute the Jacobian of the transformation of coordinates

[TABLE]

and deduce

[TABLE]

We now show that such a convenient set of coordinates can be found in spherical, spheroidal, and ellipsoidal geometries. Note that in the sequel we implicitly assume $0<\lvert\lambda\rvert<1$. Sec. 2.7 of Greenspan [1968] shows that all inertial modes in general bounded geometries have a half-frequency $-1\leq\lambda\leq 1$, that $\lambda=\pm 1$ is not an eigenfrequency, and that $\lambda=0$ is an eigenfrequency with infinite multiplicity.

A.1. Sphere and oblate spheroid

The appropriate set of coordinates in this case is the (frequency-dependent) set of bi-spheroidal coordinates $\{\xi,\mu,\phi\}$222Note that these are different from the orthogonal oblate spheroidal coordinates presented in Appendix C., which are defined, in terms of the Cartesian coordinates, as

[TABLE]

with $\beta$ as previously (Eq. (19)) and

[TABLE]

where $0<\xi<1$, $-\xi<\mu<\xi$, and $0<\phi<2\pi$. Note that these coordinates only map a bounded domain shaped as a pair of cuberdons joined along their base (corresponding to $\xi=\mu$). For convenience, we allow a double covering of this domain by extending $\mu$ to $-1<\mu<1$. Then, the level surfaces of $\xi$ are (prolate or oblate) spheroids, those of $\mu$ are half-spheroids, and those of $\phi$ are vertical half-planes. Fig. 7 gives a graphical representation of these coordinates in the $xz$-plane. Notice how this system of coordinates is non-orthogonal.

Under the above transformation, the Poincaré equation takes the form of a Laplace equation, of which the general solutions are solid bi-spheroidal harmonics,

[TABLE]

where the $p_{\ell m}$ are constants and $\mathrm{P}_{\ell}^{m}$ is the degree-$\ell$ order-$m$ associated Legendre function of the first kind. Regularity conditions prohibit the appearance of associated Legendre functions of the second kind $\mathrm{Q}_{\ell}^{m}$.

One interest of bi-spheroidal coordinates is that the surface of the oblate spheroid is a level surface of $\xi$,

[TABLE]

which makes it easier to impose the no-penetration boundary condition. The latter requires the normal velocity to vanish, and can be expressed in terms of the pressure as Eq. (25), written in bi-spheroidal coordinates. Plug Eq. (30) into Eq. (25):

[TABLE]

The functions $\mathrm{P}_{\ell}^{m}(\mu)e^{im\phi}$ form a basis of surface spheroidal harmonics over the surface of the oblate spheroid. The boundary condition therefore leaves the coefficients $p_{\ell m}$ uncoupled, another advantage of bi-spheroidal coordinates. By orthogonality of the surface spheroidal harmonics, each harmonic coefficient has to vanish. The frequencies of the inertial modes in an oblate spheroid correspond to those values of $\lambda$ or $\omega$ that allow non-zero values for $p_{\ell m}$, i.e. one of the multiplicative factors vanishes. After some rearrangement,

[TABLE]

The sphere corresponds to the special case where $e^{2}=0$, in which case one recovers Eq. (2.12.8) of Greenspan [1968].

The so-called planetary modes [Rieutord, 2014], or r-modes, are a special family of solutions that satisfy $\ell=\lvert m\rvert+1$. Inserting the latter into Eq. (33) reduces it to an equation of the first degree in $\lambda$ with solution [Zhang et al., 2004],

[TABLE]

The spin-over mode is the simplest member of this family, corresponding to $m=\pm 1$ and $\ell=2$.

A.2. Triaxial ellipsoid

The appropriate set of coordinates in this case is the (frequency-dependent) set of bi-ellipsoidal coordinates $\{\rho,\mu,\nu\}$, which are defined, in terms of the Cartesian coordinates, as

[TABLE]

with $\beta$ and $k$ as previously (Eqs. (19-29)) and

[TABLE]

where $h<\rho<k$, $h<\mu<\rho$, and $0<\nu<h$. Note that these coordinates only map one octant of a bounded domain shaped as a pair of cuberdons joined along their base (corresponding to $\rho=\mu$), hence the $\pm$ signs. For convenience, we allow a double covering of this domain by extending $\mu$ to $h<\mu<k$. Then, the level surfaces of $\rho$ and $\mu$ are (sections of) ellipsoids, those of $\nu$ are (sections of) 1-sheeted hyperboloids. Fig. 8 gives a graphical representation of these coordinates in the $xz$ and $xy$ planes. Notice how this system of coordinates is non-orthogonal.

Under the above transformation, the Poincaré equation takes the form of a Laplace equation, of which the general solutions are solid bi-ellipsoidal harmonics,

[TABLE]

where the $c_{np}$ are constants and $\mathrm{E}_{n,p}$ is the degree-$n$ index-$p$ Lamé function of the first kind333see e.g. pp.1304-1309 of Morse and Feshbach [1953] for reference. (slightly modified to ensure that it always evaluates to a real number). Regularity conditions prohibit the appearance of Lamé functions of the second kind $\mathrm{F}_{n,p}$.

One interest of bi-ellipsoidal coordinates is that the surface of the triaxial ellipsoid is a level surface of $\rho$,

[TABLE]

which makes it easier to impose the no-penetration boundary condition. The latter requires the normal velocity to vanish, and can be reexpressed in terms of the pressure as Eq. (25), written in bi-ellipsoidal coordinates. Plug Eq. (37) into Eq. (25):

[TABLE]

The functions $\mathrm{E}_{n,p}(\mu)\mathrm{E}_{n,p}(\nu)$ form a basis of surface ellipsoidal harmonics over the surface of the triaxial ellipsoid. The boundary condition therefore only slightly couples the coefficients $c_{np}$, another advantage of bi-ellipsoidal coordinates. The second term of Eq. (39) only mixes terms of the same degree and symmetry. Indeed, the second-line terms can be rewritten as a sum of surface bi-ellipsoidal harmonics [Hough, 1895]. There is, however, no general expression and each degree and symmetry has to be considered individually. By orthogonality of the surface ellipsoidal harmonics, each harmonic coefficient has to vanish. This can be rewritten as a number of homogeneous matrix equations (one for each degree and symmetry) for the coefficients $c_{np}$. The frequencies of the inertial modes in a triaxial ellipsoid correspond to those values of $\lambda$ or $\omega$ that allow non-zero values for $c_{np}$, i.e. that cancel the determinant of the matrix.

We now list the half-frequency of the first few inertial modes in a triaxial ellipsoid.

For the spin-over mode444reference? (an equatorially antisymmetric mode of degree $n=2$)

[TABLE]

The equatorially symmetric modes of degree 3 are

[TABLE]

The equatorially antisymmetric modes of degree 3 are

[TABLE]

Appendix B Spherical Harmonics decompositions

We give the spherical harmonics coefficients of the expressions that appear in Eq. (4-5) in spherical coordinates.555Those expressions depend on the choice of normalisation. Here, we use the semi-normalised spherical harmonics for which one has $Y_{0}^{0}=1$.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix C Oblate spheroidal coordinates

It is possible to extend the usage of spherical harmonics to the treatment of equations written in oblate spheroidal coordinates:

[TABLE]

The level surfaces labelled with different constant values of $\xi$ correspond to a family of confocal spheroids with foci separated by a distance $2c>0$ in the $xy$-plane.

A scalar function given in oblate spheroidal coordinates can be expanded on a series of spherical harmonics just as easily as if it were given in spherical coordinates. The only difference being that, the coordinate $\theta$ no longer represents the polar angle (geocentric colatitude) but rather the angle between the basis vector ${\bf\hat{\theta}}$ and the horizontal as illustrated on Fig. 9 (geodetic colatitude). The $\xi$ coordinate now plays a role similar to that of the radial spherical coordinate $r$. In oblate spheroidal coordinates, the (truncated) expansion of the reduced pressure reads

[TABLE]

The trick to using spherical harmonics expansion effectively for the resolution of differential equations in these coordinates consists in realising that the result of any differential operator acting on a single harmonic function, $Y_{\ell}^{m}$, can be written as a finite sum of harmonic functions multiplied by some power of a common prefactor:

[TABLE]

When the equation to solve consists of the sum of differential operators, each with its own power of the prefactor (49), the idea is to reduce the whole equation to the same denominator. This generates products of the spherical harmonics in each numerator with some power of $\cos^{2}\theta$ (which has a simple harmonic expansion in terms of $Y_{2}^{0}$ and $Y_{0}^{0}$ only). These products can be reduced to a finite sum over spherical harmonics using the usual rules. This technique critically reduces the coupling between each resulting ODE which would otherwise be infinite. In practice, we use the Mathematica package TenGSHui to expand the equations in oblate spheroidal harmonics.

Appendix D series expansion of an ellipsoid

The coefficients of the expansion of the triaxial ellipsoid surface in spherical harmonics which appear in Eq. (14) are given below to second order in the squared eccentricities $e^{2}$ and $f^{2}$.

[TABLE]