Inertial Waves in a Rotating Spherical Shell with Homogeneous Boundary Conditions

Abstract

We find the analytical form of inertial waves in an incompressible, rotating fluid constrained by concentric inner and outer spherical surfaces with homogeneous boundary conditions on the normal components of velocity and vorticity. These fields are represented by Galerkin expansions whose basis consists of toroidal and poloidal vector functions, i.e., products and curls of products of spherical Bessel functions and vector spherical harmonics. These vector basis functions also satisfy the Helmholtz equation and this has the benefit of providing each basis function with a well-defined wavenumber. Eigenmodes and associated eigenfrequencies are determined for both the ideal and dissipative cases. These eigenmodes are formed from linear combinations of the Galerkin expansion basis functions. The system is truncated to numerically study inertial wave structure, varying the number of eigenmodes. The largest system considered in detail is a 25 eigenmode system and a graphical depiction is presented of the five lowest dissipation eigenmodes, all of which are non-oscillatory. These results may be useful in understanding data produced by numerical simulations of fluid and magnetofluid turbulence in a spherical shell that use a Galerkin, toroidal–poloidal basis as well as qualitative features of liquids confined by a spherical shell.

1. Introduction

The desire to understand fluid motion in the Earth’s outer core, and ultimately the geodynamo, provides the physical motivation for our study. The magnetofluid in the outer core is primarily liquid iron and although its electrical conductivity and convective turbulence give rise to magnetic fields, we will ignore these here, but will include dissipation. Furthermore, the fluid will be treated as incompressible, as is usually done in geodynamo modeling (e.g., [1,2]). Velocity and vorticity will be expanded in terms of ‘toroidal’ and ‘poloidal’ functions, as they were first called and defined by [3]; these will be detailed presently. Here, we use Galerkin toroidal–poloidal (T–P) expansions to focus on velocity and vorticity, where ‘Galerkin’ means each term satisfies designated boundary conditions (e.g., see [4]). (We leave magnetic fields for future study. Understanding magnetohydrodynamic inertial modes through Galerkin T–P expansions is a further goal and the present work is viewed as a necessary step in that direction).

Our T–P model is based on that of [5] and is similar to that used by [6,7], although these authors work with positive and negative helicity Chandrasekhar–Kendall (C–K) functions [8], which are linear combinations of T–P functions, as will be seen explicitly in Equation (10) below; also, our model system has an inner boundary, while theirs did not. When T–P functions are put into C–K form, again as seen in Equation (10), we have ‘eigenfunctions of the curl’ [9,10] that also satisfy the Helmholtz equation and are Galerkin. The focus of either the T–P or C–K approach is on internal, free-stream turbulence and neither pretends to accurately address boundary layer motion. However, the work of [11] indicates that interior solutions for spherical shell flows are independent of b.c.s in the limit of vanishing Ekman number. In addition, little is known about the nature of the core–mantle boundary or the inner core boundary, so it is uncertain as to what kind of boundary conditions actually apply. Thus, we have the latitude to choose b.c.s that are appropriate for a T–P expansion of an incompressible fluid confined within a spherical shell.

A critical feature of the velocity and vorticity expansions is that each term will be required to individually satisfy the Helmholtz equation as well as the b.c.s. Mathematically, this means that our T–P basis functions are a product of spherical Bessel functions and spherical harmonics. The benefit of this is that each term in the expansions has a unique wavenumber associated with it. However, this affects the nature of the boundary conditions that can be imposed. The first boundary condition will be, of course, that the normal component of velocity vanishes at the boundaries. Next, the choice of T–P basis functions that satisfy the Helmholtz equation in a spherical shell leads to the second boundary condition: the normal component of vorticity also vanishes at the boundaries. What is left on the bounding surfaces is a two-dimensional potential flow, as will be seen presently in (29). If the fluid in the spherical shell was compressible, then velocity on the boundaries could be completely zeroed out by adding another term (related to compressibility) to the T–P expansion and the flow could be made to obey a no-slip condition (more details will appear in Section 4) However, we wish to treat the flow as incompressible here, to make our work commensurate with previous work on the geodynamo.

Thus, we do not use the no-slip or stress-free or inviscid boundary conditions that were compared by [11], but something else, which we will call homogeneous boundary conditions, i.e., the normal components of velocity and vorticity vanish at the boundaries. In the case of no-slip or stress-free, there are conditions on all three velocity components, while for the inviscid case, there is only one condition, that on the normal component of velocity. Here, we have the two boundary conditions stated above plus a third boundary condition intrinsically linked to the condition of incompressibility, i.e., that there is a 2-D potential flow on the boundaries. (Again, more details are given in Section 4).

At this point, let us set the context for our mathematical model and T–P vector basis functions (or equivalent C–K functions) that we use, followed by our results.

2. Historical Overview

Historically, it has long been recognized that stars and planets contain turbulent fluids. Turbulence theory and simulation began by using Fourier transformations by which the velocity and magnetic fields are represented by their Fourier coefficients. Fourier analysis, in using expansions in terms of sines and cosines, implies that physical motion is taking place in periodic box, i.e., a 3-torus which has no boundaries and hence no need for boundary conditions. Although usually called ‘periodic boundary conditions,’ this is a misnomer because, again, a 3-torus has no boundaries. The motivation for using Fourier analysis is well-described by [12], though the topological ramifications are not. A 3-torus is a self-contained compact, closed space without boundary and there is no ‘outside’ to it, topologically speaking [13]. A 3-torus is a mathematically useful space for studying fluid motion, though unphysical (except perhaps for our universe as a whole, which some conjecture may, in fact, be a 3-torus [14]). However, a 3-torus is problematic when a magnetofluid is the object of study, as there is no outside into which the model magnetic field can emerge to be compared with geophysical or astrophysical observations. Nevertheless, the majority of research on fluid and magnetofluid turbulence has relied on ‘periodic boundary conditions.’ Although inertial waves in a periodic box have been examined using Fourier methods [15], the vector basis functions are all uncoupled in the Fourier case, but, as we will see, in the rotating spherical shell, T–P functions are strongly coupled.

In regard to Galerkin methods, an important move away from Fourier analysis had been taken by [6,7], who numerically examined incompressible MHD inside a sphere (no inner boundary) using expansions in terms of positive and negative helicity C–K functions [8]. (However, [1,2] and others had used non-Galerkin T–P methods previously, using Chebyshev polynomials instead of Bessel functions for radial dependence). C–K functions have also been used to treat astrophysical problems in a spherical geometry, rather a periodic box [16]. Again, C–K and T–P functions satisfy ‘homogeneous’ boundary conditions (b.c.s): the velocity and magnetic field, as well as their curls (and the curls of those, and of those, ad infinitum) all have zero normal components on the boundary. This is analogous to Fourier analysis, where the velocity and all of its derivatives of every order are still periodic. A further analogy is that both the Fourier and C–K or T–P basis functions satisfy the Helmholtz equation, which is essential as it allows a specific wavenumber to be associated with each basis function. This, in turn, allows us to define precisely what is meant by the energy or power spectrum in fluid mechanics.

Extending the Galerkin approach of [6,7] was done in [5] by providing the sphere with an inner boundary upon which homogeneous b.c.s are also applied, but using a T–P form rather than a C–K form. This system can be set rotating mathematically, giving us opportunity here to study waves within a spherical shell with homogeneous b.c.s. Using Galerkin T–P expansions to analyze the linearized Navier–Stokes equation for flows within a spherical shell with homogeneous b.c.s, leads us to explicit formulas for inertial waves, which are the new results that we present.

Rotation enters the Navier–Stokes equation through the Coriolis term and linearization gives the basic equation for inertial waves. Traditionally, for the case of a sphere with no-slip conditions, these waves have been studied and solutions provided by [17,18], building on the work of [19], while [20] have used a perturbation expansion to examine pressure oscillations, finding solutions that are ‘pathological.’ Moreover, in the case of flow in a rotating spherical shell, Ref. [17] has stated that, “Explicit formulas for all the eigenmodes in the spherical annulus do not exist. Bryan’s transformation [p. 64], the essential step in the solution procedure for the spherical container, does not apply and for this reason very little is known about the inertial waves.” This sentiment concerning analytical solutions is echoed by [21,22], as well as by [23], who, in their examination of inertial waves in a rotating spherical shell, used non-Galerkin T–P methods, as radial cofactors have primarily been Chebyshev polynomials, rather than Bessel functions. This may have its benefits numerically, but such basis functions do not satisfy the Helmholtz equation and wavenumbers cannot be precisely associated with them.

Our approach here varies in obvious and significant details from the traditional approaches that use, for example, no-slip [17] or stress-free conditions [24]. First, we use ‘homogeneous b.c.s’ and choose to avoid conjecturing about the nature of the boundary between, say, the Earth’s outer core and mantle or between the outer core and inner core, boundaries which are presently unknown, and possibly unknowable. Second, velocity and vorticity are represented by Galerkin T–P functions (or the equivalent C–K form), i.e., vector basis functions that have as radial cofactors a linear combination of spherical Bessel functions and vector spherical harmonics.

The results presented here are extensions of our previous work [5,25] in which T–P expansions using a Galerkin basis as defined by [3] were used, while others may also have used a general T–P formulation, these were designed for computational efficacy and to satisfy no-slip b.c.s. For example, radial functions can be Chebyshev polynomials [1,2,23] or be represented numerically by finite differences [26], but again, these are not Galerkin methods. Here, we use a Galerkin expansion with spherical Bessel of the first and second kind (the second kind are also called Neumann functions) that describe radial variation.

Using our model, with velocity and vorticity expanded in terms of Galerkin T–P vector functions, we are able to find explicit formulas for all the eigenmodes of inertial waves in an incompressible, rotating fluid contained between concentric spherical surfaces with homogeneous b.c.s, rather than traditional no-slip or stress-free boundary conditions. Our method of solution appears robust as we find no pathologies. We determine eigenfrequencies which we compute using Earth-like parameters. We present these results following a brief discussion of the basic equations and the Galerkin method of solution. As an example, we consider a dissipative 25 eigenmode system in detail; the method can, of course, be expanded to any desired number of eigenmodes given sufficient computational resources.

3. Basic Equations

In non-dimensional form, the linearized incompressible Navier–Stokes equation in a rotating frame of reference with constant angular velocity ${\Omega }_o$ is

(1)$\frac{\partial \mathit{u}}{\partial t}=\mathbf{-}\nabla p+2\mathit{u}\times \widehat{\mathit{z}}+E{\nabla }^2\mathit{u},\nabla ·\mathit{u}=0.$

For the Earth, ${\Omega }_o\cong 7.27\times {10}^{\mathbf{-}5}$ rad/s and $E\sim {10}^{\mathbf{-}14}$ is the Ekman number for the outer core [27]. Although E is very small, the dissipation term ${\nabla }^2\mathit{u}$ will get large at a small enough length-scale, i.e., large enough wavenumber, which needs to be estimated.

We approximate the outer core to be a spherical shell with inner radius $R_1=1.22\times {10}^6$ m and outer radius $R_2=3.482\times {10}^6$ m [27]; $R_1$ is the characteristic distance and the ratio $r_o=R_2/R_1=2.85$. The non-dimensional Equation (1) will be examined in a spherical polar coordinate system $(r,\theta ,\phi )$, with $1\le r\le r_o$, $0\le \theta \le \pi$ and $0\le \phi <2\pi$. The boundary conditions at $r=1$ and $r_o$ are that the normal components of $\mathit{u}$ and $\mathit{\omega }=\nabla \times \mathit{u}$ are zero:

(2)$\widehat{\mathit{r}}·\mathit{u}{{|}_{1,r_o}=0,\widehat{\mathit{r}}·\mathit{\omega }|}_{1,r_o}=0.$

Thus, we appear to have two conditions: (a) a Dirichlet condition on $u_r$ and (b) a Neumann condition linking angular derivatives of $u_{\theta }$ and $u_{\phi }$ on the model boundaries (we use the appellations ‘Dirichlet condition’ and ‘Neumann condition’ in a general sense); conditions (a) and (b) are implicit in the ‘homogeneous b.c.s’ given by (2). It should be remembered, though, that the curl of the vorticity, the curl of that, etc., also satisfy homogeneous b.c.s; this is simply a property of T–P (or C–K) expansions. As stated earlier, a third condition could be imposed if the flow was not considered incompressible.

Again, many different boundary conditions have been employed for geodynamo studies, although the particular set used often does not appear to have significant effects on final results, as discussed by [11,28]. Mathematically, the choice (2) allows us to employ the Galerkin expansion functions that are used for determining the structure of inertial waves in a rotating spherical shell. As mentioned before, the individual basis functions of T–P (or C–K) expansions each satisfy a different Helmholtz equation, allowing each to have a well-defined wavenumber.

4. Galerkin Expansion

It is well-known (e.g., see [29]) that a divergence-less vector field in a spherical geometry can be represented as a ‘toroidal–poloidal decomposition’:

(3)$\begin{array}{ccc}\hfill \mathit{u}(\mathit{x},t)& =& \nabla \times \left[\mathit{r}U\right(r,\theta ,\phi ,t\left)\right]+\nabla \times (\nabla \times [\mathit{r}W(r,\theta ,\phi ,t)\left]\right).\hfill \end{array}$

In order to satisfy the Helmholtz equation in a spherical shell, the functions U and W in (3) must take the form

(4)$\begin{array}{ccc}\hfill U(r,\theta ,\varphi ,t)=\sum _{l,m,n}{u}_{lmn}\left(t\right)g_l\left(k_{ln}r\right){\mathrm{Y}}_{lm}(\theta ,\phi ),& & W(r,\theta ,\varphi ,t)=\sum _{l,m,n}{w}_{lmn}\left(t\right)g_l\left(k_{ln}r\right){\mathrm{Y}}_{lm}(\theta ,\phi ).\hfill \end{array}$

Here, ${\mathrm{Y}}_{lm}$ are spherical harmonics and the $g_l\left(k_{ln}r\right)$ are linear combinations of the spherical Bessel functions of the first kind, $j_l\left(k_{ln}r\right)$, and of the second kind, $n_l\left(k_{ln}r\right)$ (also called spherical Neumann functions). The exact form of the $g_l\left(k_{ln}r\right)$ will be given below, while the wavenumbers $k_{ln}$ will be determined by the b.c.s. Note that if the flow were compressible, we would need to add a term $\nabla \Phi$ to (3), where $\Phi$ would have an expansion similar to those in (4), with radial derivative zero on the boundaries. Here, we only wish to outline a novel approach to studying compressible flow, while our main aim is to focus on incompressible flow.

Now, putting (4) into (3), we can write the toroidal–poloidal decomposition in the form:

(5)$\begin{array}{ccc}\hfill \mathit{u}(\mathit{x},t)& =& \sum _{l,m,n}\left[u_{lmn}\left(t\right){\mathit{T}}_{lmn}\left(\mathit{x}\right)+w_{lmn}\left(t\right){\mathit{P}}_{lmn}\left(\mathit{x}\right)\right],\hfill \end{array}$

(6)${\mathit{T}}_{lmn}\left(\mathit{x}\right)=c_{lmn}{g}_l\left(k_{ln}r\right)\mathit{r}\times \nabla {\mathrm{Y}}_{lm}(\theta ,\phi ),$

(7)$\nabla \times {\mathit{T}}_{lmn}\left(\mathit{x}\right)={\mathit{P}}_{lmn}\left(\mathit{x}\right),\nabla \times {\mathit{P}}_{lmn}\left(\mathit{x}\right)=k_{ln}^2{\mathit{T}}_{lmn}\left(\mathit{x}\right),$

(8)$u_{l,\mathbf{-}m,n}\left(t\right)={(\mathbf{-}1)}^m{u}_{lmn}^{*}\left(t\right),w_{l,\mathbf{-}m,n}\left(t\right)={(\mathbf{-}1)}^m{w}_{lmn}^{*}\left(t\right).$

Here, $c_{lmn}$ is a normalizing constant, to be defined presently; also, ${}^{*}$ denotes complex conjugation. Furthermore, in (5), the overall minus sign that occurred because $\nabla \times \mathit{r}{\mathrm{Y}}_{lm}=\mathbf{-}\mathit{r}\times \nabla {\mathrm{Y}}_{lm}$ has been dropped. Furthermore, if one wanted to approach the problem with no-slip b.c.s. and compressibility, then one would still use the vector expansion functions ${\mathit{T}}_{lmn}\left(\mathit{x}\right)$ and ${\mathit{P}}_{lmn}\left(\mathit{x}\right)$ in addition to introducing a new set of expansion functions for the compressible term $\nabla \Phi$.

Obviously, $\nabla ·{\mathit{T}}_{lmn}\left(\mathit{x}\right)=\nabla ·{\mathit{P}}_{lmn}\left(\mathit{x}\right)=0$, while (8) is required by the realness of $\mathit{u}$ and the definition of the toroidal and poloidal vector functions, ${\mathit{T}}_{lmn}\left(\mathit{x}\right)$ and ${\mathit{P}}_{lmn}\left(\mathit{x}\right)$ given below, where it will be seen that the ${\mathit{T}}_{lmn}\left(\mathit{x}\right)$ and ${\mathit{P}}_{lmn}\left(\mathit{x}\right)$ form an orthogonal set for $1\le r\le r_o$, each member of which will be seen to satisfy the boundary conditions (so that we have a Galerkin expansion), while the $k_{ln}^2$ are numerical constants and the summation indices $l,m,n$ range over prescribed, finite range of integer values. It is evident that we have two sets of essential, generally complex, coefficients ($u_{lmn}$ and $w_{lmn}$) to represent the velocity. The flow we are studying has vorticity $\mathit{\omega }=\nabla \times \mathit{u}$ and using the relations given above, this is

(9)$\begin{array}{ccc}\hfill \mathit{\omega }(\mathit{x},t)& =& \sum _{l,m,n}\left[k_{ln}^2{w}_{lmn}\left(t\right){\mathit{T}}_{lmn}\left(\mathit{x}\right)+u_{lmn}\left(t\right){\mathit{P}}_{lmn}\left(\mathit{x}\right)\right].\hfill \end{array}$

We mention the vorticity because it is used to define the kinetic helicity which, along with energy, is an ideal invariant for (1) with homogeneous boundary conditions [30]. In what follows, we will use $\omega$ to represent angular frequency, but there should be no confusion as the magnitude of $\mathit{\omega }$ will not appear here, although its component ${\omega }_{\phi }$ will.

Note that we could have written (5) and (9) in terms of an explicitly helical set of vector basis functions:

(10)$\begin{array}{c}\hfill {\mathit{J}}_{lmn}^{\pm }\left(\mathit{x}\right)=\pm k_{ln}{\mathit{T}}_{lmn}\left(\mathit{x}\right)+{\mathit{P}}_{lmn}\left(\mathit{x}\right),\nabla \times {\mathit{J}}_{lmn}^{\pm }\left(\mathit{x}\right)=\pm k_{ln}{\mathit{J}}_{lmn}^{\pm }\left(\mathit{x}\right).\end{array}$

The ${\mathit{J}}_{lmn}^{\pm }\left(\mathit{x}\right)$ are Chandrasekhar-Kendall functions [8], and are eigenfunctions of the curl operator [9,10]. Although we choose to use ${\mathit{T}}_{lmn}\left(\mathit{x}\right)$ and ${\mathit{P}}_{lmn}\left(\mathit{x}\right)$ instead, our toroidal-poloidal vector basis functions are essentially the same as those in (10), differing only by a linear transformation.

In general, the various quantities used above can be defined as follows. First,

(11)$\begin{array}{ccc}\hfill {\mathit{T}}_{lmn}\left(\mathit{x}\right)& =& {\widehat{g}}_l\left(k_{ln}r\right){\widehat{\Phi }}_{lm}(\theta ,\phi ),1\le l\le L,\mathbf{-}l\le m\le l,1\le n\le N,\hfill \\ \hfill \end{array}$

(12)$\begin{array}{ccc}\hfill {\mathit{P}}_{lmn}\left(\mathit{x}\right)& =& \mathbf{-}\frac{1}{r}\left[\sqrt{l(l+1)}{\widehat{g}}_l\left(k_{ln}r\right){\widehat{\mathit{Y}}}_{lm}(\theta ,\phi )+\frac{\mathrm{d}}{\mathrm{d}r}\left[r{\widehat{g}}_l\left(k_{ln}r\right)\right]{\widehat{\Psi }}_{lm}(\theta ,\phi )\right].\hfill \end{array}$

The orthonormal vector functions ${\widehat{\mathit{Y}}}_{lm}(\theta ,\phi )$, ${\widehat{\Psi }}_{lm}(\theta ,\phi )$ and ${\widehat{\Phi }}_{lm}(\theta ,\phi )$ are the (normalized) vector spherical harmonics of [31] and are based on the spherical harmonics ${\mathrm{Y}}_{l,m}(\theta ,\phi )$ described in [32]:

(13)$\begin{array}{ccc}\hfill {\mathrm{Y}}_{l,\mathbf{-}m}(\theta ,\phi )& =& {(\mathbf{-}1)}^m{\mathrm{Y}}_{lm}^{*}(\theta ,\phi ),\hfill \\ \hfill \end{array}$

(14)$\begin{array}{ccc}\hfill {\widehat{\mathit{Y}}}_{lm}(\theta ,\phi )& =& \widehat{\mathit{r}}{\mathrm{Y}}_{lm}(\theta ,\phi ),\hfill \\ \hfill \end{array}$

(15)$\begin{array}{ccc}\hfill {\widehat{\Psi }}_{lm}(\theta ,\phi )& =& \frac{r}{\sqrt{l(l+1)}}\nabla {\mathrm{Y}}_{lm}(\theta ,\phi ),\hfill \\ \hfill \end{array}$

(16)$\begin{array}{ccc}\hfill {\widehat{\Phi }}_{lm}(\theta ,\phi )& =& \widehat{\mathit{r}}\times {\widehat{\Psi }}_{lm}(\theta ,\phi ).\hfill \end{array}$

The radial functions ${\widehat{g}}_l\left(k_{ln}r\right)$ in (11) and (12) form an orthonormal set for $n=1,2,\dots ,N$ at each value of l over the range $1\le r\le r_o$, and are linear combinations of the spherical Bessel and Neumann functions, $j_l\left(z\right)$ and $n_l\left(z\right)={(\mathbf{-}1)}^{l+1}{j}_{\mathbf{-}l\mathbf{-}1}\left(z\right)$ (see, e.g., [33]):

(17)$\begin{array}{ccc}\hfill {\widehat{g}}_l\left(k_{ln}r\right)& =& N_{ln}^{\mathbf{-}1}{g}_l\left(k_{ln}r\right),\hfill \\ \hfill \end{array}$

(18)$\begin{array}{ccc}\hfill g_l\left(k_{ln}r\right)& =& n_l\left(k_{ln}{r}_o\right)j_l\left(k_{ln}r\right)\mathbf{-}{j}_l\left(k_{ln}{r}_o\right)n_l\left(k_{ln}r\right),\hfill \\ \hfill \end{array}$

(19)$\begin{array}{ccc}\hfill N_{ln}^2& =& \frac{1}{2k_{ln}^4}\left(\frac{1}{r_o}\mathbf{-}{\left[\frac{n_l\left(k_{ln}{r}_o\right)}{n_l\left(k_{ln}\right)}\right]}^2\right),\hfill \\ \hfill \end{array}$

(20)$\begin{array}{ccc}\hfill g_l\left(k_{ln}\right)& =& 0,n=1,2,\dots ,N;l=1,2,\dots ,L.\hfill \end{array}$

Therefore, ${\widehat{g}}_l\left(k_{ln}\right)={\widehat{g}}_l\left(k_{ln}{r}_o\right)=0$ and the zeroes $k_{ln}$ of $g_l\left(z\right)$ satisfy

(21)$j_l\left(k_{ln}\right)n_l\left(k_{ln}{r}_o\right)=n_l\left(k_{ln}\right)j_l\left(k_{ln}{r}_o\right).$

(An equivalent form of this equation has been studied earlier by [34]). Using the results above, along with (15) and (16), the normalizing constant in (6) is seen to be $c_{lmn}={\left(\sqrt{l(l+1)}{N}_{ln}\right)}^{\mathbf{-}1}$. The $k_{ln}$ in (21) clearly depend on the ratio $r_o$; again, we use $r_o\equiv 2.85$, which is appropriate for the Earth’s outer core. The $k_{ln}$ are easily found numerically, and for $l,n=1,\dots ,10$, they appear in Table 1.

Let us define integrals over radius r as ${⟨F⟩}_r$ and over solid angle $\sigma$ as ${⟨G⟩}_{\sigma }$, where

(22)${⟨F⟩}_r={\int }_1^{r_o}F\left(r\right)r^2\mathrm{d}r,{⟨G⟩}_{\sigma }={\int }_{\theta =0}^{\pi }{\int }_{\phi =0}^{2\pi }G(\theta ,\phi )sin\theta \mathrm{d}\theta \mathrm{d}\phi .$

Thus, the integral over the volume of the spherical shell interior is

(23)$⟨F\left(r\right)G(\theta ,\phi )⟩={\int }_1^{r_o}{\int }_{\theta =0}^{\pi }{\int }_{\phi =0}^{2\pi }F\left(r\right)G(\theta ,\phi )r^{2}sin\theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\phi ={⟨F⟩}_r{⟨G⟩}_{\sigma }.$

The orthonormality properties of the Galerkin basis functions are

(24)$\begin{array}{ccc}\hfill {⟨{\widehat{g}}_l\left(k_{ln}r\right){\widehat{g}}_l\left(k_{lk}r\right)⟩}_r& =& {\delta }_{nk},{⟨{\widehat{\mathit{V}}}_{lm}^{\left(i\right)}·{\widehat{\mathit{V}}}_{pq}^{\left(j\right)*}⟩}_{\sigma }={\delta }_{ij}{\delta }_{lp}{\delta }_{mq},\hfill \end{array}$

(25)$\begin{array}{ccc}\hfill {\widehat{\mathit{V}}}_{lm}^{\left(1\right)}& =& {\widehat{\mathit{Y}}}_{lm}(\theta ,\phi ),{\widehat{\mathit{V}}}_{lm}^{\left(2\right)}={\widehat{\Psi }}_{lm}(\theta ,\phi ),{\widehat{\mathit{V}}}_{lm}^{\left(3\right)}={\widehat{\Phi }}_{lm}(\theta ,\phi ).\hfill \end{array}$

Using these and the properties given in (11) and (12), we see that the integrals of the inner products of ${\mathit{T}}_{lmn}\left(\mathit{x}\right)$ and ${\mathit{P}}_{lmn}\left(\mathit{x}\right)$ over the volume enclosed by the homogeneous boundaries are

(26)$\begin{array}{ccc}\hfill ⟨{\widehat{\mathit{T}}}_{pqs}·{\widehat{\mathit{T}}}_{lmn}^{*}⟩& =& {⟨{\widehat{g}}_p\left(k_{ps}r\right){\widehat{g}}_l\left(k_{ln}r\right)⟩}_r{⟨{\widehat{\Phi }}_{pqs}·{\widehat{\Phi }}_{lmn}^{*}⟩}_{\sigma }={\delta }_{pl}{\delta }_{qm}{\delta }_{sn},\hfill \\ \hfill \end{array}$

(27)$\begin{array}{ccc}\hfill ⟨{\widehat{\mathit{P}}}_{pqs}·{\widehat{\mathit{P}}}_{lmn}^{*}⟩& =& k_{ln}^2⟨{\widehat{\mathit{T}}}_{pqs}·{\widehat{\mathit{T}}}_{lmn}^{*}⟩,\hfill \\ \hfill \end{array}$

(28)$\begin{array}{ccc}\hfill ⟨{\widehat{\mathit{T}}}_{pqs}·{\widehat{\mathit{P}}}_{lmn}^{*}⟩& =& 0.\hfill \end{array}$

These will be useful in determining the linear inertial waves of the system.

While we use a model with an inner core, refs. [6,7] perform numerical simulations on a model with no inner core. The results we present do not depend on whether there is an inner core or not. The only adjustment necessary is that the radial expansion functions (18) become $g_l\left(k_{ln}r\right)=j_l\left(k_{ln}r\right)$ when there is no inner core. In this case, the numerical values of the $k_{ln}$ in Table 1 change slightly, but this does not affect our analysis, which is independent of the exact values of the $k_{ln}$, as these are always distinct and well-ordered solutions of $g_l\left(k_{ln}\right)=0$, as discussed by [35].

As to the flow on the boundary, consider (11) and (12) at $r=r_o,1$:

(29)${\mathit{T}}_{lmn}{\left(\mathit{x}\right)|}_{r_o,1}=0,{\mathit{P}}_{lmn}\left(\mathit{x}\right){{|}_{r_o,1}=\mathbf{-}\frac{k_{ln}}{\sqrt{l(l+1)}}r{\widehat{g}}_l^{\prime }\left(k_{ln}r\right)\nabla {\mathrm{Y}}_{lm}(\theta ,\phi )|}_{r_o,1}.$

This indicates a 2-D potential flow on the bounding surfaces that could be negated if compressibility was assumed for the fluid, as mentioned earlier in the discussion following (3) and (4). The development of this possibility is beyond the scope of the present work and will be deferred.

5. Inertial Waves in a Spherical Shell

If we place Equation (5) into the Navier–Stokes Equation (1), take the inner product of both sides with ${\widehat{\mathit{T}}}_{lmn}^{*}$ and ${\widehat{\mathit{P}}}_{lmn}^{*}$, then integrate over the spherical shell volume, using the orthogonality relations (24) and (25), we get

(30)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{u}_{lmn}}{\mathrm{d}t}& =& \mathbf{-}2\widehat{\mathit{z}}·\sum _{p,q,s}\left[u_{pqs}⟨{\widehat{\mathit{T}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩+w_{pqs}⟨{\widehat{\mathit{P}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩\right]\mathbf{-}E{k}_{ln}^2{u}_{lmn},\hfill \\ \hfill \end{array}$

(31)$\begin{array}{ccc}\hfill k_{ln}^2\frac{\mathrm{d}{w}_{lmn}}{\mathrm{d}t}& =& \mathbf{-}2\widehat{\mathit{z}}·\sum _{p,q,s}\left[u_{pqs}⟨{\widehat{\mathit{T}}}_{pqs}\times {\widehat{\mathit{P}}}_{lmn}^{*}⟩+w_{pqs}⟨{\widehat{\mathit{P}}}_{pqs}\times {\widehat{\mathit{P}}}_{lmn}^{*}⟩\right]\mathbf{-}E{k}_{ln}^4{w}_{lmn}.\hfill \end{array}$

Note that pressure p does not appear above because $⟨(\nabla p)·{\mathit{Q}}^{*}⟩=⟨\nabla ·\left(p{\mathit{Q}}^{*}\right)⟩=0$, where $\mathit{Q}$ is either ${\widehat{\mathit{T}}}_{lmn}$ or ${\widehat{\mathit{P}}}_{lmn}$. Here, we see that the dissipative term may be neglected in a model system with given L and N if $E{k}_{ln}^2\ll 1$ for $l\le L$ and $n\le N$. For Earth-like parameters, the dissipative term is extremely small and dissipation wavenumber $k_D\sim k_{LN}\sim E^{\mathbf{-}1/2}\sim {10}^7$. If the system of Equations (30) and (31) is truncated at moderate sizes of N and L, where $k_{LN}\ll {10}^7$ the dissipative term may be ignored for earth-like E, though we will study model systems here by making E large enough to have significant effect.

Now, consider the integral $⟨{\widehat{\mathit{P}}}_{pqs}\times {\widehat{\mathit{P}}}_{lmn}^{*}⟩$. We can use (11) and (12), as well as standard vector identities, along with Gauss’s divergence theorem and the boundary conditions to determine that

(32)$⟨{\widehat{\mathit{P}}}_{pqs}\times {\widehat{\mathit{P}}}_{lmn}^{*}⟩=k_{ln}^2⟨{\widehat{\mathit{T}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩=k_{qs}^2⟨{\widehat{\mathit{T}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩.$

These equations lead immediately to

(33)$\left(k_{ps}^2\mathbf{-}{k}_{ln}^2\right)⟨{\widehat{\mathit{T}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩=0.$

Since the wavenumbers $k_{ln}$ are all distinct, and using the fact that ${\widehat{\mathit{T}}}_{pqs}$ only has $\phi$-dependence through a factor $exp\left(\mathrm{i}q\phi \right)$, where $\mathrm{i}=\sqrt{\mathbf{-}1}$, we arrive at

(34)$\begin{array}{ccc}\hfill ⟨{\widehat{\mathit{T}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩& =& ⟨{\widehat{\mathit{T}}}_{lmn}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩{\delta }_{pl}{\delta }_{qm}{\delta }_{sn},\hfill \\ \hfill \end{array}$

(35)$\begin{array}{ccc}\hfill ⟨{\widehat{\mathit{P}}}_{pqs}\times {\widehat{\mathit{P}}}_{lmn}^{*}⟩& =& k_{ln}^2⟨{\widehat{\mathit{T}}}_{lmn}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩{\delta }_{pl}{\delta }_{qm}{\delta }_{sn}.\hfill \end{array}$

Next, we must determine $\widehat{\mathit{z}}·⟨{\widehat{\mathit{T}}}_{lmn}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩$ and then $\widehat{\mathit{z}}·⟨{\widehat{\mathit{P}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩$.

The quantum mechanical angular momentum operator $\mathit{L}$, whose properties are well known, is (with $\hslash =1$)

(36)$\mathit{L}=\mathbf{-}\mathrm{i}\mathit{r}\times \nabla .$

Using this, along with (11) and (15), allows ${\widehat{\mathit{T}}}_{lmn}$ to be expressed as

(37)${\widehat{\mathit{T}}}_{lmn}=\mathrm{i}\frac{{\widehat{g}}_l\left(k_{ln}\right)}{\sqrt{l(l+1)}}\mathit{L}{Y}_{lm}.$

The differential operator $\mathit{L}$ has the properties

(38)${\mathit{L}}^{*}=\mathbf{-}\mathit{L},\mathit{L}\times {\mathit{L}}^{*}=\mathrm{i}{\mathit{L}}^{*},\widehat{\mathit{z}}·\mathit{L}{Y}_{lm}=L_z{Y}_{lm}=m{Y}_{lm}.$

Thus, we have

(39)$\begin{array}{ccc}\hfill \widehat{\mathit{z}}·⟨{\widehat{\mathit{T}}}_{lmn}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩& =& \frac{1}{l(l+1)}\widehat{\mathit{z}}·{⟨\mathit{L}{Y}_{lm}\times {\mathit{L}}^{*}{Y}_{lm}^{*}⟩}_{\sigma }\hfill \\ \hfill \end{array}$

(40)$\begin{array}{ccc}& =& \frac{1}{l(l+1)}\widehat{\mathit{z}}·{⟨\mathit{L}\times \left(Y_{lm}{\mathit{L}}^{*}{Y}_{lm}^{*}\right)\mathbf{-}{Y}_{lm}\mathit{L}\times {\mathit{L}}^{*}{Y}_{lm}^{*}⟩}_{\sigma }\hfill \\ \hfill \end{array}$

(41)$\begin{array}{ccc}& =& \frac{\mathbf{-}\mathrm{i}}{l(l+1)}\widehat{\mathit{z}}·{⟨Y_{lm}{\mathit{L}}^{*}{Y}_{lm}^{*}⟩}_{\sigma }\hfill \\ \hfill \end{array}$

(42)$\begin{array}{ccc}& =& \frac{\mathbf{-}\mathrm{i}m}{l(l+1)}.\hfill \end{array}$

The first term in (40) vanishes because of the identity ${\oint }_S(d\mathit{S}\times \nabla )\times \mathit{A}=0$.

Now, to determine the cross terms $\widehat{\mathit{z}}·⟨{\widehat{\mathit{P}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩$. Using (11), (12), (15) and (16), along with the recurrence relations of the spherical Bessel, Neumann and associated Legendre functions [33], and after much algebra, we find

(43)$\begin{array}{ccc}\hfill \widehat{\mathit{z}}·⟨{\widehat{\mathit{P}}}_{pqs}\times {\widehat{\mathit{T}}}_{lmn}^{*}⟩& =& {\delta }_{mq}\left(G_{psln}{C}_{lm}{\delta }_{l,p+1}\mathbf{-}{G}_{lnps}{C}_{pm}{\delta }_{p,l+1}\right)\hfill \\ \hfill \end{array}$

(44)$\begin{array}{ccc}\hfill G_{psln}& =& {⟨\begin{array}{c}{\widehat{g}}_p\left(k_{ps}r\right)\left\{l{\widehat{g}}_l\left(k_{ln}r\right)+\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left[r{\widehat{g}}_l\left(k_{ln}r\right)\right]\right\}\end{array}⟩}_r,\hfill \\ \hfill \end{array}$

(45)$\begin{array}{ccc}\hfill C_{lm}& =& \frac{\sqrt{l^2\mathbf{-}1}}{l}\sqrt{\frac{l^2\mathbf{-}{m}^2}{4l^2\mathbf{-}1}}.\hfill \end{array}$

Note that the vector product is antisymmetric when $pqs$ and $lmn$ are interchanged, as it must be, and that $0\le C_{lm}<1/2$ for $l\ge 1$ and $\left|m\right|<l$; for $\left|m\right|\ge l$, $C_{lm}=0$. Furthermore, note that these cross terms do not vanish when $m=0$. Furthermore,

(46)$\begin{array}{ccc}\hfill \frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left[r{\widehat{g}}_l\left(k_{ln}r\right)\right]& =& \left(\frac{l+1}{2l+1}{\delta }_{a,l\mathbf{-}1}\mathbf{-}\frac{l}{2l+1}{\delta }_{a,l+1}\right)D_{lna}\left(r\right),\hfill \\ \hfill \end{array}$

(47)$\begin{array}{ccc}\hfill D_{lna}\left(r\right)& =& \frac{k_{ln}}{N_{ln}}\left[n_l\left(k_{ln}{r}_o\right)j_a\left(k_{ln}r\right)\mathbf{-}{j}_l\left(k_{ln}{r}_o\right)n_a\left(k_{ln}r\right)\right].\hfill \end{array}$

When needed, this will be of use in computing (44).

The cross terms in (30) and (31), which have just been evaluated, pertain to a rotating spherical shell. Ref. [7] consider a rotating sphere and in that case called these cross terms ‘coupling coefficients’ and although they stated that these were computed numerically, no analytical determination similar to ours was presented. In our previous work [5], we did not consider these cross terms at all and in the work presented here, we make up for this deficiency.

Now, if we put all of these results into the set of ordinary differential Equations (30) and (31), we correct and expand the Coriolis terms given in [5]:

(48)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{u}_{lmn}}{\mathrm{d}t}& =& \mathrm{i}\frac{2m}{l(l+1)}{u}_{lmn}\mathbf{-}2\sum _{s=1}^N\left(G_{l\mathbf{-}1,sln}{C}_{lm}{w}_{l\mathbf{-}1,ms}\mathbf{-}{G}_{ln,l+1,s}{C}_{l+1,m}{w}_{l+1,ms}\right),\hfill \\ \hfill \end{array}$

(49)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{w}_{lmn}}{\mathrm{d}t}& =& \mathrm{i}\frac{2m}{l(l+1)}{w}_{lmn}\mathbf{-}\frac{2}{k_{ln}^2}\sum _{s=1}^N\left(G_{l\mathbf{-}1,sln}{C}_{lm}{u}_{l\mathbf{-}1,ms}\mathbf{-}{G}_{ln,l+1,s}{C}_{l+1,m}{u}_{l+1,ms}\right).\hfill \end{array}$

Here, we have set $E=0$ to eliminate the dissipative term, which can be placed back later, if desired. The equations above are qualitatively different from the inertial wave equations that are obtained by Fourier analysis [17,24]. In particular, we see that moving from a periodic box to a more realistic spherical shell model reveals the intrinsic coupling amongst the toroidal and poloidal linear modes. For each m, there are two separate, disjoint sets. In one set, toroidal $u_{l_{e}mn}$ with $l_e$ even couple to poloidal $w_{l_{o}ms}$ with $l_o$ odd; in the other set, toroidal $u_{l_{o}mn}$ couple to poloidal $w_{l_{e}ms}$. This result is fundamentally different to the one derived from Fourier analysis where there is no coupling.

Again, in (48) and (49), modes $u_{lmn}$ and $w_{l^{\prime }ms}$ are only coupled to modes with the same m and opposite parity for l and $l^{\prime }$. Toroidal modes $l,m,n$ are coupled to poloidal modes $l\pm 1,m,s$, and all modes with same m are coupled to one another in their respective set; for a given N and L, each set will have $NL$ eigenmodes. We can examine these coupled systems by restricting (48) and (49) to $l=1,2$ and $n=s=1$. For $m=0$, we have, using the notation ${\overline{w}}_{lmn}\equiv k_{ln}{w}_{lmn}$,

(50)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{u}_{101}}{\mathrm{d}t}& =& \frac{2}{\sqrt{5}}\frac{G_{1121}}{k_{21}}{\overline{w}}_{201},\frac{\mathrm{d}{\overline{w}}_{201}}{\mathrm{d}t}=\mathbf{-}\frac{2}{\sqrt{5}}\frac{G_{1121}}{k_{21}}{u}_{101},\hfill \\ \hfill \end{array}$

(51)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{\overline{w}}_{101}}{\mathrm{d}t}& =& \frac{2}{\sqrt{5}}\frac{G_{1121}}{k_{11}}{u}_{201},\frac{\mathrm{d}{u}_{201}}{\mathrm{d}t}=\mathbf{-}\frac{2}{\sqrt{5}}\frac{G_{1121}}{k_{11}}{\overline{w}}_{101}.\hfill \end{array}$

For the wave numbers given in Table 1, $G_{11,21}=2.0671$. Using this, the set of coupled Equations (50) and (51) gives the variables $u_{101}$ and ${\overline{w}}_{201}$ a dimensional oscillation frequency of ${\omega }_{21}^{\left(0\right)}\approx \pm 0.860{\Omega }_o$ and the set (50) and (51) gives $u_{201}$ and ${\overline{w}}_{101}$ an dimensional oscillation frequency of ${\overline{\omega }}_{21}^{\left(0\right)}\approx \pm 0.992{\Omega }_o$. (Again, ${\Omega }_o\cong 7.27\times {10}^{\mathbf{-}5}$ rad/s). Note that $u_{101}$ and ${\overline{w}}_{101}$ are not coupled and that $u_{201}$ and ${\overline{w}}_{201}$ are not coupled either. Instead, the Coriolis force links dipole and quadrupole components here.

Now, let us consider the variables with $m=1$ in the reduced system. Their equations are

(52)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{u}_{111}}{\mathrm{d}t}& =& \mathrm{i}{u}_{111}+\sqrt{\frac{3}{5}}\frac{G_{1121}}{k_{21}}{\overline{w}}_{211},\frac{\mathrm{d}{\overline{w}}_{211}}{\mathrm{d}t}=\frac{\mathrm{i}}{3}{\overline{w}}_{211}\mathbf{-}\sqrt{\frac{3}{5}}\frac{G_{1121}}{k_{21}}{u}_{111},\hfill \\ \hfill \end{array}$

(53)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{\overline{w}}_{111}}{\mathrm{d}t}& =& \mathrm{i}{\overline{w}}_{111}+\sqrt{\frac{3}{5}}\frac{G_{1121}}{k_{11}}{u}_{211},\frac{\mathrm{d}{u}_{211}}{\mathrm{d}t}=\frac{\mathrm{i}}{3}{u}_{211}\mathbf{-}\sqrt{\frac{3}{5}}\frac{G_{1121}}{k_{11}}{\overline{w}}_{111}.\hfill \end{array}$

Compared to (50) and (51), these equations have an extra term on the right; without the coupling terms, the 111 coefficients would have a dimensional angular frequency of ${\Omega }_o$ and the 211 coefficients have a frequency of ${\Omega }_o/3$. However, they are coupled and the eigenvalues of the coupling matrices provide the actual frequencies. Again using $G_{11,21}=2.0671$, the set of coupled Equation (52) gives the variables $u_{111}$ and ${\overline{w}}_{211}$ frequencies of ${\omega }_{21}^{\left(1\right)}\approx 1.59{\Omega }_o$ and $\mathbf{-}0.255{\Omega }_o$, while the set (50) and (51) gives $u_{211}$ and ${\overline{w}}_{111}$ the frequencies are ${\overline{\omega }}_{21}^{\left(1\right)}\approx 1.48{\Omega }_o$ and $\mathbf{-}0.149{\Omega }_o$. (For $m=\mathbf{-}1$, these eigenfrequencies have their signs reversed). Again, the Coriolis force links dipole and quadrupole components and allows them to oscillate at frequencies determined by the set of linked multipole equations.

The $m=0$ Equations (50) and (51) and $m=1$ Equations (52) and (53) define four independent coupled linear systems of two variables each. However, these systems were created by ignoring all variables with $l>2$. For example, the full equation for $u_{201}$, rather than (51), couples it to ${\overline{w}}_{301}$, so that we have

(54)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{u}_{201}}{\mathrm{d}t}& =& \mathbf{-}2\left(\frac{G_{1121}{C}_{20}}{k_{11}}{\overline{w}}_{101}\mathbf{-}\frac{G_{2131}{C}_{30}}{k_{31}}{\overline{w}}_{301}\right),\hfill \\ \hfill \end{array}$

(55)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{\overline{w}}_{301}}{\mathrm{d}t}& =& \mathbf{-}2\left(\frac{G_{2131}{C}_{30}}{k_{31}}{u}_{201}\mathbf{-}\frac{G_{3141}{C}_{40}}{k_{31}}{u}_{401}\right).\hfill \end{array}$

Thus, even though we can create a closed system for $u_{lmn}$ and ${\overline{w}}_{lmn}$ with all $l>2$ variables set to zero, a system can be set up with larger L and N, where $1\le l\le L$ and $1\le n\le N$. The values of L and N, which set the spatial resolution, can be made as large as computationally practicable, while the closure problem for an ideal system can be addressed by using a suitable $E>0$ to determine the dissipation wave number $k_D$, which (as in fluid turbulence) should be commensurate with the largest wavenumber of the system, $k_{LN}$. We now consider these two topics: resolution and dissipation.

6. General Solution

A uniform approximation of fluid velocity in a spherical shell would ostensibly require $L=N$; However, it may be beneficial to produce formulas that allow L and N to be independent. In this regard, let us define the following matrix elements (reintroducing dissipation):

(56)$\begin{array}{ccc}\hfill D_{ns}^{\left(lm\right)}& =& \left(\frac{m}{l(l+1)}+\mathrm{i}\frac{E}{2}{k}_{ln}^2\right){\delta }_{ns},A_{ns}^{\left(lm\right)}=\frac{G_{ls,l+1,n}{C}_{l+1,m}}{k_{ln}},B_{ns}^{\left(lm\right)}=\frac{G_{ln,l+1,s}{C}_{l+1,m}}{k_{l+1,s}}.\hfill \end{array}$

Using these, (48) and (49) can be written as

(57)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{u}_{lmn}}{\mathrm{d}t}& =& 2\mathrm{i}\sum _{s=1}^N\left(D_{ns}^{\left(lm\right)}{u}_{lmn}+\mathrm{i}{A}_{sn}^{(l\mathbf{-}1,m)}{\overline{w}}_{l\mathbf{-}1,ms}\mathbf{-}\mathrm{i}{B}_{ns}^{\left(lm\right)}{\overline{w}}_{l+1,ms}\right),\hfill \\ \hfill \end{array}$

(58)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{\overline{w}}_{l+1,mn}}{\mathrm{d}t}& =& 2\mathrm{i}\sum _{s=1}^N\left(D_{ns}^{(l+1,m)}{\overline{w}}_{l+1,mn}+\mathrm{i}{B}_{sn}^{\left(lm\right)}{u}_{lms}\mathbf{-}\mathrm{i}{A}_{ns}^{(l+1,m)}{u}_{l+2,ms}\right).\hfill \end{array}$

(Remember that ${\overline{w}}_{lmn}\equiv k_{ln}{w}_{lmn}$). Let us further define the $N\times 1$ vectors (‘${}^T$’ is transpose),

(59)$\begin{array}{ccc}\hfill {\mathsf{U}}_{lm}^{\left(N\right)}& =& {\left[u_{lm1}{u}_{lm2}\cdots u_{lmN}\right]}^T,{\mathsf{W}}_{lm}^{\left(N\right)}={\left[{\overline{w}}_{lm1}{\overline{w}}_{lm2}\cdots {\overline{w}}_{lmN}\right]}^T,\hfill \end{array}$

and $N\times N$ matrices,

(60)$\begin{array}{ccc}\hfill {\mathsf{D}}_N^{\left(lm\right)}& =& \left[D_{ns}^{\left(lm\right)}\right],{\mathsf{A}}_N^{\left(lm\right)}=\left[A_{ns}^{\left(lm\right)}\right],{\mathsf{B}}_N^{\left(lm\right)}=\left[B_{ns}^{\left(lm\right)}\right],1\le n,s\le N.\hfill \end{array}$

Note that $u_{lmn}={\overline{w}}_{lmn}=0$ and $A_{ns}^{\left(lm\right)}=B_{ns}^{\left(lm\right)}=0$ for $m>l$.

At this point, we can write (57) and (58) as

(61)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{\mathsf{U}}_{lm}^{\left(N\right)}}{\mathrm{d}t}& =& 2\mathrm{i}\left({\mathsf{D}}_N^{\left(lm\right)}{\mathsf{U}}_{lm}^{\left(N\right)}+\mathrm{i}{\mathsf{A}}_N^{(l\mathbf{-}1,m)T}{\mathsf{W}}_{l\mathbf{-}1,m}^{\left(N\right)}\mathbf{-}\mathrm{i}{\mathsf{B}}_N^{\left(lm\right)}{\mathsf{W}}_{l+1,m}^{\left(N\right)}\right),\hfill \\ \hfill \end{array}$

(62)$\begin{array}{ccc}\hfill \frac{\mathrm{d}{\mathsf{W}}_{l+1,m}^{\left(N\right)}}{\mathrm{d}t}& =& 2\mathrm{i}\left({\mathsf{D}}_N^{(l+1,m)}{\mathsf{W}}_{l+1,m}^{\left(N\right)}+\mathrm{i}{\mathsf{B}}_N^{\left(lm\right)T}{\mathsf{U}}_{lm}^{\left(N\right)}\mathbf{-}\mathrm{i}{\mathsf{A}}_N^{(l+1,m)}{\mathsf{U}}_{l+2,m}^{\left(N\right)}\right).\hfill \end{array}$

Continuing along this line, we construct the following ‘super vectors’:

(63)$\begin{array}{ccc}\hfill {\mathbb{U}}_{LN}^{\left(m\right)}& =& \left[\begin{array}{c}{\mathsf{U}}_{1m}^{\left(N\right)}\hfill \\ {\mathsf{W}}_{2m}^{\left(N\right)}\hfill \\ ⋮\hfill \\ {\mathsf{U}}_{L\mathbf{-}1,m}^{\left(N\right)}\hfill \\ {\mathsf{W}}_{Lm}^{\left(N\right)}\hfill \end{array}\right],{\overline{\mathbb{U}}}_{LN}^{\left(m\right)}=\left[\begin{array}{c}{\mathsf{W}}_{1m}^{\left(N\right)}\hfill \\ {\mathsf{U}}_{2m}^{\left(N\right)}\hfill \\ ⋮\hfill \\ {\mathsf{W}}_{L\mathbf{-}1,m}^{\left(N\right)}\hfill \\ {\mathsf{U}}_{Lm}^{\left(N\right)}\hfill \end{array}\right],\left(L\mathrm{even}\right);\hfill \\ \hfill \\ \hfill \end{array}$

(64)$\begin{array}{ccc}\hfill {\mathbb{U}}_{LN}^{\left(m\right)}& =& \left[\begin{array}{c}{\mathsf{U}}_{1m}^{\left(N\right)}\hfill \\ {\mathsf{W}}_{2m}^{\left(N\right)}\hfill \\ ⋮\hfill \\ {\mathsf{W}}_{L\mathbf{-}1,m}^{\left(N\right)}\hfill \\ {\mathsf{U}}_{Lm}^{\left(N\right)}\hfill \end{array}\right],{\overline{\mathbb{U}}}_{LN}^{\left(m\right)}=\left[\begin{array}{c}{\mathsf{W}}_{1m}^{\left(N\right)}\hfill \\ {\mathsf{U}}_{2m}^{\left(N\right)}\hfill \\ ⋮\hfill \\ {\mathsf{U}}_{L\mathbf{-}1,m}^{\left(N\right)}\hfill \\ {\mathsf{W}}_{Lm}^{\left(N\right)}\hfill \end{array}\right],\left(L\mathrm{odd}\right).\hfill \end{array}$

We also define the following ‘super matrices’, which are Hermitian if $E=0$:

(65)$\begin{array}{ccc}\hfill {\mathbb{M}}_{LN}^{\left(m\right)}& =& \left[\begin{array}{ccccccc}{\mathsf{D}}_N^{\left(1m\right)}& \mathbf{-}\mathrm{i}{\mathsf{B}}_N^{\left(1m\right)}& 0& 0& 0& \cdots & \cdots \\ \mathrm{i}{\mathsf{B}}_N^{\left(1m\right)T}& {\mathsf{D}}_N^{\left(2m\right)}& \mathbf{-}\mathrm{i}{\mathsf{A}}_N^{\left(2m\right)}& 0& 0& \cdots & \cdots \\ 0& \mathrm{i}{\mathsf{A}}_N^{\left(2m\right)T}& {\mathsf{D}}_N^{\left(3m\right)}& \mathbf{-}\mathrm{i}{\mathsf{B}}_N^{\left(3m\right)}& 0& \cdots & \cdots \\ 0& 0& \mathrm{i}{\mathsf{B}}_N^{\left(3m\right)T}& {\mathsf{D}}_N^{\left(4m\right)}& \mathbf{-}\mathrm{i}{\mathsf{A}}_N^{\left(4m\right)}& \cdots & \cdots \\ 0& 0& 0& \mathrm{i}{\mathsf{A}}_N^{\left(4m\right)T}& {\mathsf{D}}_N^{\left(5m\right)}& \cdots & \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& {\mathsf{D}}_N^{\left(Lm\right)}\end{array}\right],\hfill \\ \hfill \\ \hfill \end{array}$

(66)$\begin{array}{ccc}\hfill {\overline{\mathbb{M}}}_{LN}^{\left(m\right)}& =& \left[\begin{array}{ccccccc}{\mathsf{D}}_N^{\left(1m\right)}& \mathbf{-}\mathrm{i}{\mathsf{A}}_N^{\left(1m\right)}& 0& 0& 0& \cdots & \cdots \\ \mathrm{i}{\mathsf{A}}_N^{\left(1m\right)T}& {\mathsf{D}}_N^{\left(2m\right)}& \mathbf{-}\mathrm{i}{\mathsf{B}}_N^{\left(2m\right)}& 0& 0& \cdots & \cdots \\ 0& \mathrm{i}{\mathsf{B}}_N^{\left(2m\right)T}& {\mathsf{D}}_N^{\left(3m\right)}& \mathbf{-}\mathrm{i}{\mathsf{A}}_N^{\left(3m\right)}& 0& \cdots & \cdots \\ 0& 0& \mathrm{i}{\mathsf{A}}_N^{\left(3m\right)T}& {\mathsf{D}}_N^{\left(4m\right)}& \mathbf{-}\mathrm{i}{\mathsf{B}}_N^{\left(4m\right)}& \cdots & \cdots \\ 0& 0& 0& \mathrm{i}{\mathsf{B}}_N^{\left(4m\right)T}& {\mathsf{D}}_N^{\left(5m\right)}& \cdots & \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& {\mathsf{D}}_N^{\left(Lm\right)}\end{array}\right].\hfill \end{array}$

Finally, we get to the canonical forms:

(67)$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbb{U}}_{LN}^{\left(m\right)}}{\mathrm{d}t}=2\mathrm{i}{\mathbb{M}}_{LN}^{\left(m\right)}{\mathbb{U}}_{LN}^{\left(m\right)},\frac{\mathrm{d}{\overline{\mathbb{U}}}_{LN}^{\left(m\right)}}{\mathrm{d}t}=2\mathrm{i}{\overline{\mathbb{M}}}_{LN}^{\left(m\right)}{\overline{\mathbb{U}}}_{LN}^{\left(m\right)}.\end{array}$

In the ideal case $E=0$, the matrices ${\mathbb{M}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{M}}}_{LN}^{\left(m\right)}$ are Hermitian and thus have real eigenvalues. Furthermore, for each m there are three ideal invariants, the energies $E_{LN}^{\left(m\right)}$ and ${\overline{E}}_{LN}^{\left(m\right)}$, as well as the kinetic helicity $H_{LN}^{\left(m\right)}$ because of (56):

(68)$\begin{array}{ccc}\hfill E_{LN}^{\left(m\right)}& =& \frac{1}{2}{\mathbb{U}}_{LN}^{\left(m\right)†}{\mathbb{U}}_{LN}^{\left(m\right)},{\overline{E}}_{LN}^{\left(m\right)}=\frac{1}{2}{\overline{\mathbb{U}}}_{LN}^{\left(m\right)†}{\overline{\mathbb{U}}}_{LN}^{\left(m\right)},H_{LN}^{\left(m\right)}=\frac{1}{2}\left({\overline{\mathbb{U}}}_{LN}^{\left(m\right)†}{\mathbb{K}}_{LN}{\mathbb{U}}_{LN}^{\left(m\right)}+{\mathbb{U}}_{LN}^{\left(m\right)†}{\mathbb{K}}_{LN}{\overline{\mathbb{U}}}_{LN}^{\left(m\right)}\right).\hfill \end{array}$

Here, $H_{LN}^{\left(m\right)}$ is ‘kinetic helicity’ and ${\mathbb{K}}_{LN}$ is the diagonal matrix

(69)$\begin{array}{ccc}\hfill {\mathbb{K}}_{LN}& =& \mathrm{diag}\left[k_{11}{k}_{12}\cdots k_{1N}\cdots k_{L1}{k}_{L2}\cdots k_{LN}\right].\hfill \end{array}$

Each inertial wave vector ${\mathbb{U}}_{LN}^{\left(m\right)}$ or ${\overline{\mathbb{U}}}_{LN}^{\left(m\right)}$ as a whole conserves energy when $E=0$, but their individual modes do not. In the theory of ideal fluid turbulence [36], the ideal invariance of total energy and kinetic helicity allow for a statistical description of a model system; here, the lack of an inertial term $\mathit{u}·\nabla \mathit{u}$ in (1) precludes turbulence. However, it turns out that there are $2LN+1$ ideal invariants for each m, rather than the three in (68) and (69).

Let ${\mathbb{S}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{S}}}_{LN}^{\left(m\right)}$ be the unitary matrices that diagonalize ${\mathbb{M}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{M}}}_{LN}^{\left(m\right)}$:

(70)$\begin{array}{c}\hfill {\Lambda }_{LN}^{\left(m\right)}={\mathbb{S}}_{LN}^{\left(m\right)}{\mathbb{M}}_{LN}^{\left(m\right)}{\mathbb{S}}_{LN}^{\left(m\right)†},{\overline{\Lambda }}_{LN}^{\left(m\right)}={\overline{\mathbb{S}}}_{LN}^{\left(m\right)}{\overline{\mathbb{M}}}_{LN}^{\left(m\right)}{\overline{\mathbb{S}}}_{LN}^{\left(m\right)†}.\end{array}$

The superscript ‘${}^{†}$’ denotes, as usual, a Hermitian adjoint. The diagonal matrices ${\Lambda }_{LN}^{\left(m\right)}$ and ${\overline{\Lambda }}_{LN}^{\left(m\right)}$ contain the eigenvalues, which are complex when $E>0$. Multiplying these eigenvalues by 2 yields the non-dimensional eigenfrequencies; when these, in turn are multiplied by ${\Omega }_o$, they are the dimensional eigenfrequencies of the eigenmodes of the system.

The unitary matrices ${\mathbb{S}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{S}}}_{LN}^{\left(m\right)}$ in (70) also produce the eigenmodes of the system by transforming the vectors ${\mathbb{U}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{U}}}_{LN}^{\left(m\right)}$ defined in (63) and (64):

(71)$\begin{array}{c}\hfill {\mathbb{V}}_{LN}^{\left(m\right)}={\mathbb{S}}_{LN}^{\left(m\right)}{\mathbb{U}}_{LN}^{\left(m\right)},{\overline{\mathbb{V}}}_{LN}^{\left(m\right)}={\overline{\mathbb{S}}}_{LN}^{\left(m\right)}{\overline{\mathbb{U}}}_{LN}^{\left(m\right)}.\end{array}$

The eigenmodes are the components of ${\mathbb{V}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{V}}}_{LN}^{\left(m\right)}$ and they are linear combinations of the components of ${\mathbb{U}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{U}}}_{LN}^{\left(m\right)}$, respectively. Thus, the eigenmodes are combinations of the various multipoles (dipole, quadrupole, etc.) which comprise the inertial waves of the system. The extent of this mixing will be considered in the next section for the case $L=N=5$.

In terms of ${\mathbb{V}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{V}}}_{LN}^{\left(m\right)}$, Equation (67) become

(72)$\begin{array}{c}\hfill \frac{\mathrm{d}{\mathbb{V}}_{LN}^{\left(m\right)}}{\mathrm{d}t}=2\mathrm{i}{\Lambda }_{LN}^{\left(m\right)}{\mathbb{V}}_{LN}^{\left(m\right)},\frac{\mathrm{d}{\overline{\mathbb{V}}}_{LN}^{\left(m\right)}}{\mathrm{d}t}=2\mathrm{i}{\overline{\Lambda }}_{LN}^{\left(m\right)}{\overline{\mathbb{V}}}_{LN}^{\left(m\right)}.\end{array}$

These have immediate solution:

(73)$\begin{array}{c}\hfill {\mathbb{V}}_{LN}^{\left(m\right)}\left(t\right)=exp\left(2\mathrm{i}{\Lambda }_{LN}^{\left(m\right)}\right){\mathbb{V}}_{LN}^{\left(m\right)}\left(0\right),{\overline{\mathbb{V}}}_{LN}^{\left(m\right)}\left(t\right)=exp\left(2\mathrm{i}{\overline{\Lambda }}_{LN}^{\left(m\right)}\right){\overline{\mathbb{V}}}_{LN}^{\left(m\right)}\left(0\right).\end{array}$

The eigenvectors ${\mathbb{V}}_{LN}^{\left(m\right)}\left(t\right)$ and ${\overline{\mathbb{V}}}_{LN}^{\left(m\right)}\left(t\right)$ have components (eigenmodes) with constant magnitude and periodically varying phase, whereas including an inertial term in (1) would connect them and subject both to statistical fluctuations. These $2LN$ magnitudes plus the kinetic helicity (69) are the ideal invariants of the model system, i.e., when $E=0$.

The case where $L=2$ and $N=1$ have been expressed in (50)–(53), where it was seen that the eigenfrequencies for the former are different from those for the latter. In general. For large L and N, the matrices (65) and (66) present a formidable eigenvalue problem, particularly with respect to accuracy. However, for small L and N, the eigenvalues are easily found numerically and we now turn to a discussion of these.

7. Eigenfrequencies and Eigenmodes

We have seen some examples already, in 2-dimensional systems (50) and (52) for modes $u_{lmn}$${\overline{w}}_{lmn}$ with $l\le L=2$ and $n=N=1$. Again, a more uniform approximation of fluid velocity in a spherical shell would ostensibly require $L=N$. So, let us first truncate (63)–(66) to $l,n\le L=N=3$ and look at the systems

(74)$\begin{array}{ccc}\hfill {\mathbb{U}}_3^{\left(m\right)}& =& {\left[u_{1m1}{u}_{1m2}{u}_{1m3}{\overline{w}}_{2m1}{\overline{w}}_{2m2}{\overline{w}}_{2m3}{u}_{3m1}{u}_{3m2}{u}_{3m3}\right]}^T,\hfill \\ \hfill \end{array}$

(75)$\begin{array}{ccc}\hfill {\overline{\mathbb{U}}}_3^{\left(m\right)}& =& {\left[{\overline{w}}_{1m1}{\overline{w}}_{1m2}{\overline{w}}_{1m3}{u}_{2m1}{u}_{2m2}{u}_{2m3}{\overline{w}}_{3m1}{\overline{w}}_{3m2}{\overline{w}}_{3m3}\right]}^T,\hfill \end{array}$

(76)$\begin{array}{ccc}\hfill {\mathbb{M}}_{33}^{\left(m\right)}& =& \left[\begin{array}{ccc}{\mathsf{D}}_3^{\left(1m\right)}& \mathbf{-}\mathrm{i}{\mathsf{B}}_3^{\left(1m\right)}& 0\\ \mathrm{i}{\mathsf{B}}_3^{\left(1m\right)T}& {\mathsf{D}}_3^{\left(2m\right)}& \mathbf{-}\mathrm{i}{\mathsf{A}}_3^{\left(2m\right)}\\ 0& \mathrm{i}{\mathsf{A}}_3^{\left(2m\right)T}& {\mathsf{D}}_3^{\left(3m\right)}\end{array}\right],{\overline{\mathbb{M}}}_{33}^{\left(m\right)}=\left[\begin{array}{ccc}{\mathsf{D}}_3^{\left(1m\right)}& \mathbf{-}\mathrm{i}{\mathsf{A}}_3^{\left(1m\right)}& 0\\ \mathrm{i}{\mathsf{A}}_3^{\left(1m\right)T}& {\mathsf{D}}_3^{\left(2m\right)}& \mathbf{-}\mathrm{i}{\mathsf{B}}_3^{\left(2m\right)}\\ 0& \mathrm{i}{\mathsf{B}}_3^{\left(2m\right)T}& {\mathsf{D}}_3^{\left(3m\right)}\end{array}\right].\hfill \end{array}$

Let us set $E=0$ and $m=0$. Using the matrix elements defined in (56), along with Table 1 and formulas (44)–(47), we compute (using MatLab) the following, where ${\mathsf{O}}_3$ is the zero $3\times 3$ matrix,

(77)$\begin{array}{ccc}\hfill {\mathsf{D}}_3^{\left(10\right)}& =& {\mathsf{O}}_3,{\mathsf{D}}_3^{\left(20\right)}={\mathsf{O}}_3,{\mathsf{D}}_3^{\left(30\right)}={\mathsf{O}}_3,\hfill \end{array}$

(78)$\begin{array}{ccc}\hfill {\mathsf{A}}_3^{\left(10\right)}& \cong & \left[\begin{array}{ccc}\hfill 0.4960& \hfill \mathbf{-}0.3273& \hfill \mathbf{-}0.0408\\ \hfill 0.1728& \hfill 0.2764& \hfill \mathbf{-}0.3249\\ \hfill 0.0018& \hfill 0.2177& \hfill 0.1898\end{array}\right],{\mathsf{A}}_3^{\left(20\right)}\cong \left[\begin{array}{ccc}\hfill 0.6872& \hfill \mathbf{-}0.2839& \hfill \mathbf{-}0.0676\\ \hfill 0.1643& \hfill 0.3899& \hfill \mathbf{-}0.3185\\ \hfill 0.0087& \hfill 0.2190& \hfill 0.2963\end{array}\right],\hfill \end{array}$

(79)$\begin{array}{ccc}\hfill {\mathsf{B}}_3^{\left(10\right)}& \cong & \left[\begin{array}{ccc}\hfill 0.4300& \hfill 0.1641& \hfill 0.0017\\ \hfill \mathbf{-}0.2837& \hfill 0.2624& \hfill 0.2123\\ \hfill \mathbf{-}0.0354& \hfill \mathbf{-}0.3085& \hfill 0.1851\end{array}\right],{\mathsf{B}}_3^{\left(20\right)}\cong \left[\begin{array}{ccc}\hfill 0.5899& \hfill 0.1534& \hfill 0.0084\\ \hfill \mathbf{-}0.2437& \hfill 0.3639& \hfill 0.2113\\ \hfill \mathbf{-}0.0580& \hfill \mathbf{-}0.2973& \hfill 0.2859\end{array}\right].\hfill \end{array}$

Now, let us keep $E=0$ but set $m=1$. Using the matrix elements defined in (56), along with Table 1 and Equations (44)–(47), we compute (using MatLab) the following, where ${\mathsf{I}}_3$ is the unit $3\times 3$ matrix,

(80)$\begin{array}{ccc}\hfill {\mathsf{D}}_3^{\left(11\right)}& =& \frac{1}{2}{\mathsf{I}}_3,{\mathsf{D}}_3^{\left(21\right)}=\frac{1}{6}{\mathsf{I}}_3,{\mathsf{D}}_3^{\left(31\right)}=\frac{1}{12}{\mathsf{I}}_3,\hfill \end{array}$

(81)$\begin{array}{ccc}\hfill {\mathsf{A}}_3^{\left(11\right)}& \cong & \left[\begin{array}{ccc}\hfill 0.4296& \hfill \mathbf{-}0.2834& \hfill \mathbf{-}0.0353\\ \hfill 0.1497& \hfill 0.2393& \hfill \mathbf{-}0.2814\\ \hfill 0.0015& \hfill 0.1885& \hfill 0.1643\end{array}\right],{\mathsf{A}}_3^{\left(21\right)}\cong \left[\begin{array}{ccc}\hfill 0.6479& \hfill \mathbf{-}0.2677& \hfill \mathbf{-}0.0637\\ \hfill 0.1549& \hfill 0.3676& \hfill \mathbf{-}0.3002\\ \hfill 0.0082& \hfill 0.2065& \hfill 0.2793\end{array}\right].\hfill \end{array}$

(82)$\begin{array}{ccc}\hfill {\mathsf{B}}_3^{\left(11\right)}& \cong & \left[\begin{array}{ccc}\hfill 0.3724& \hfill 0.1421& \hfill 0.0015\\ \hfill \mathbf{-}0.2457& \hfill 0.2272& \hfill 0.1838\\ \hfill \mathbf{-}0.0306& \hfill \mathbf{-}0.2672& \hfill 0.1603\end{array}\right],{\mathsf{B}}_3^{\left(21\right)}\cong \left[\begin{array}{ccc}\hfill 0.5562& \hfill 0.1446& \hfill 0.0080\\ \hfill \mathbf{-}0.2298& \hfill 0.3431& \hfill 0.1993\\ \hfill \mathbf{-}0.0547& \hfill \mathbf{-}0.2803& \hfill 0.2695\end{array}\right].\hfill \end{array}$

The numerical matrix elements given above are, of course, approximate values and presented for purposes of illustration and an example for verification. One feature to notice in the matrices (78), (79), (81) and (82) is that the first element $n=s=1$ in each has the largest magnitude. This property appears to hold for all ${\mathsf{A}}_N^{\left(lm\right)}$ and ${\mathsf{B}}_N^{\left(lm\right)}$, with $N\ge 1$.

From the above results, we can form the matrices given in (75) and (76), i.e., ${\mathbb{M}}_{33}^{\left(0\right)}$, ${\overline{\mathbb{M}}}_{33}^{\left(0\right)}$, ${\mathbb{M}}_{33}^{\left(1\right)}$ and ${\overline{\mathbb{M}}}_{33}^{\left(1\right)}$. (Again, $E=0$ here). The eigenvalues of these, multiplied by 2, are the eigenfrequencies of associated truncated systems; these eigenfrequencies, written as vectors, are

(83)$\begin{array}{ccc}\hfill {\omega }_{33}^{\left(0\right)}& =& \left[1.853,1.369,0.857,0.000,0.000,0.000,\mathbf{-}0.857,\mathbf{-}1.369,\mathbf{-}1.853\right],\hfill \\ \hfill \end{array}$

(84)$\begin{array}{ccc}\hfill {\overline{\omega }}_{33}^{\left(0\right)}& =& \left[1.785,1.352,0.844,0.000,0.000,0.000,\mathbf{-}0.844,\mathbf{-}1.352,\mathbf{-}1.785\right],\hfill \\ \hfill \end{array}$

(85)$\begin{array}{ccc}\hfill {\omega }_{33}^{\left(1\right)}& =& \left[2.130,1.738,1.305,0.749,0.663,0.647,\mathbf{-}0.473,\mathbf{-}0.898,\mathbf{-}1.362\right],\hfill \\ \hfill \end{array}$

(86)$\begin{array}{ccc}\hfill {\overline{\omega }}_{33}^{\left(1\right)}& =& \left[2.113,1.743,1.308,0.681,0.612,0.592,\mathbf{-}0.452,\mathbf{-}0.857,\mathbf{-}1.239\right].\hfill \end{array}$

Above, ${\omega }_{33}^{\left(m\right)}$ are the eigenfrequencies associated with ${\mathbb{M}}_{33}^{\left(m\right)}$ and ${\overline{\omega }}_{33}^{\left(m\right)}$ are the eigenfrequencies associated with ${\overline{\mathbb{M}}}_{33}^{\left(m\right)}$ for $m=0$ and 1; those associated with $m=\mathbf{-}1$ are the negative of those in (86). Here, these are all real because $E=0$; they become complex, with positive imaginary parts when $E>0$. (The ordering of the eigenfrequencies and their corresponding eigenvectors is arbitrary; above, the choice is most positive to most negative eigenfrequency, but this is not essential).

The eigenmodes can be found using (71), but this will not be done here; however, the case $L=N=5$ will be considered presently and the structure of the eigenmodes will be examined graphically. Note the appearance of zero eigenfrequencies in (83) and (84); N zero values occur when L is odd, but none when L is even or when $\left|m\right|>0$ (or when $E>0$, though the real part will be zero for any odd L and become zero if E is large enough, as will be seen presently).

Earlier, in the discussion following (50) and (51), the eigenfrequencies associated with ${\mathbb{M}}_{21}^{\left(0\right)}$, ${\overline{\mathbb{M}}}_{21}^{\left(0\right)}$ were given, and similarly, following (52) and (53), the eigenfrequencies associated with ${\mathbb{M}}_{21}^{\left(1\right)}$, ${\overline{\mathbb{M}}}_{21}^{\left(1\right)}$ were given, while above $L=N=3$ was considered. Of course, we can find the eigenvalues of ${\mathbb{M}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{M}}}_{LN}^{\left(m\right)}$ for any values of L and N. As L and N increase, the number of eigenfrequencies (and eigenmodes) is $L\times N$, and the presentation of either as an array of numbers would perhaps be overwhelming. Instead, results will be presented graphically (Figure 1, 5 and 6) in what follows and will be limited to $L=N$ for $L\le 5$ and $m=0,\pm 1$.

8. Undamped and Damped Inertial Waves

Consider the cases with $E=0$ and $L=N\le 5$. The case $L=1$ is simple: for $m=0$, the eigenfrequencies are zero and for $m=\pm 1$, they are $\pm 1$ (in units of ${\Omega }_o$). For $2\le L\le 5$, we use MatLab to find the eigenfrequencies from ${\mathbb{M}}_{LL}^{\left(m\right)}$ and ${\overline{\mathbb{M}}}_{LL}^{\left(m\right)}$. In the case where $m=0$, the eigenfrequencies are presented in Figure 2. There it can be seen that the undamped eigenfrequencies appear to increase without bound for $E=0$. For $L=N$ = 1, 3, 5 there are 1, 3 and 5 zero eigenfrequencies, respectively, although they are on top of each other for $E=0$.

However, in physical body such as the Earth, E, though small, is not zero, but using $E={10}^{\mathbf{-}14}$ would require $L\sim {10}^7$, which from a standpoint of accuracy and computation time is not very feasible. So, for example, we look at the cases where $L=5$ and $0.01\le E\le 1.25$; the results of doing so gives us the sequence of plots shown in Figure 3. There, each of the eigenfrequencies has gained a positive imaginary part, leading to damping of the eigenmode. As E increases, more and more eigenfrequencies get a real part of zero, indicating overdamping. Thus, as L increases beyond what is shown in Figure 2, the bounding curve will reach a limit and eventually return to $\Re {\omega }_{LN}^{\left(m\right)}=0$.

In Figure 4, the damped eigenfrequencies for $E=0.01$ and the cases $m=0,\pm 1$ are presented. Notice that in both Figure 3 and Figure 4, the five zeroes that were superimposed for the $L=5$ case in Figure 2 have now separated due to each gaining a different damping factor $\Im {\omega }_{55}^{\left(0\right)}$. The same will occur for the multiple zeros of any odd L when damping is introduced.

In (71), the matrices ${\mathbb{S}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{S}}}_{LN}^{\left(m\right)}$ used to transform the primary vectors ${\mathbb{U}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{U}}}_{LN}^{\left(m\right)}$ in eigenvectors ${\mathbb{V}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{V}}}_{LN}^{\left(m\right)}$, each of the components of which is an eigenmode. These eigenmodes are combinations of the components (modes) of ${\mathbb{U}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{U}}}_{LN}^{\left(m\right)}$, as defined by (59), (60), (63) and (64); the specific manner in which each mode goes in each eigenmode is determined by the $LN$ rows ($I=1,\dots ,LN$) $S^I$ of ${\mathbb{S}}_{LN}^{\left(m\right)}$ and ${\overline{S}}^I$ of ${\overline{\mathbb{S}}}_{LN}^{\left(m\right)}$. The $LN$ elements ($J=1,\dots ,LN$) $S_J^I$ of $S^I$ and ${\overline{S}}_J^I$ of ${\overline{S}}^I$ are generally complex and for $L=N=5$ would give an array of 25 complex numbers, too many perhaps to be very illuminating.

Instead, for $E=0.01$ and $m=0$, we give the eigenfrequencies in Figure 1 and Table 2, and plot the magnitudes $|S_J^I|$ and $|{\overline{S}}_J^I|$ for the $L=N=5$ case in Figure 5, which contains the 25 magnitudes for each $I=1,\dots ,15$; the corresponding eigenfrequencies are given in Table 2. In the $m=0$ case, there is no need to give the magnitudes and eigenfrequencies for $I=16,\dots ,25$ because these are the same as those found in the range $I=1,\dots ,10$ (although the complex values for these $S_J^I$ and ${\overline{S}}_J^I$ are different). Similar figures can be made for the $m\ne 0$ eigenmodes, but this will not be done here, where the $m=0$ case is chosen as an example.

One interesting feature seen in plots $I=11,\dots ,15$, of Figure 5 is that the $|S_J^I|$ and $|{\overline{S}}_J^I|$ have only three significant magnitudes, while the plots, $I=1,\dots ,10$, have many. Table 3 explores this feature in more detail by showing the actual values, rather than just the magnitudes, of these large magnitude $S_J^I$ and $|S_J^I|$. In Table 3, we see that actual values are all real and positive. Remember that for $I=11,\dots ,15$, $\Re {\omega }_{LN}^{\left(0\right)}\left(I\right)=0$ and $\Re {\overline{\omega }}_{LN}^{\left(0\right)}\left(I\right)=0$; thus, these are non-oscillatory eigenmodes.

If we consider the structure of ${\mathbb{U}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{U}}}_{LN}^{\left(m\right)}$ given in (63) and (64), we see that the three significant magnitudes for the eigenmodes $I=11,\dots ,15$, in Figure 5 which have $\Re {\omega }_{LN}^{\left(0\right)}\left(I\right)=0$, are the toroidal parts of ${\mathbb{U}}_{LN}^{\left(0\right)}$ and those which have $\Re {\overline{\omega }}_{LN}^{\left(0\right)}\left(I\right)=0$ are the poloidal parts of ${\overline{\mathbb{U}}}_{LN}^{\left(0\right)}$, while the eigenmodes whose frequencies have zero real part are almost purely either toroidal or poloidal, the eigenmodes $I=1,\dots ,10$, have a mixture of toroidal and poloidal components; according to Table 2, these also have $\Re {\omega }_{LN}^{\left(0\right)}\left(I\right)\ne 0$ and $\Re {\overline{\omega }}_{LN}^{\left(0\right)}\left(I\right)\ne 0$, and are thus time-dependent eigenmodes.

9. Physical Space Eigenfunctions

In order to visualize the non-oscillatory eigenmodes, we find their functional form in physical space. First, we define arrays containing the expansion functions ${\mathit{T}}_{lmn}\left(\mathit{x}\right)$ and ${\mathit{P}}_{lmn}\left(\mathit{x}\right)$:

(87)$\begin{array}{ccc}\hfill {\mathsf{T}}_{lm}^{\left(N\right)}& =& {\left[{\mathit{T}}_{lm1}{\mathit{T}}_{lm2}\cdots {\mathit{T}}_{lmN}\right]}^T,\hfill \\ \hfill \end{array}$

(88)$\begin{array}{ccc}\hfill {\mathsf{P}}_{lm}^{\left(N\right)}& =& {\left[k_{lm1}^{\mathbf{-}1}{\mathit{P}}_{lm1}{k}_{lm2}^{\mathbf{-}1}{\mathit{P}}_{lm2}\cdots k_{lmN}^{\mathbf{-}1}{\mathit{P}}_{lmN}\right]}^T.\hfill \end{array}$

We use these to define larger arrays analogous to ${\mathbb{U}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{U}}}_{LN}^{\left(m\right)}$:

(89)$\begin{array}{ccc}\hfill {\mathbb{Q}}_{LN}^{\left(m\right)}& =& \left[\begin{array}{c}{\mathsf{T}}_{1m}^{\left(N\right)}\hfill \\ {\mathsf{P}}_{2m}^{\left(N\right)}\hfill \\ ⋮\hfill \\ {\mathsf{T}}_{L\mathbf{-}1,m}^{\left(N\right)}\hfill \\ {\mathsf{P}}_{Lm}^{\left(N\right)}\hfill \end{array}\right],{\overline{\mathbb{Q}}}_{LN}^{\left(m\right)}=\left[\begin{array}{c}{\mathsf{P}}_{1m}^{\left(N\right)}\hfill \\ {\mathsf{T}}_{2m}^{\left(N\right)}\hfill \\ ⋮\hfill \\ {\mathsf{P}}_{L\mathbf{-}1,m}^{\left(N\right)}\hfill \\ {\mathsf{T}}_{Lm}^{\left(N\right)}\hfill \end{array}\right],\left(L\mathrm{even}\right);\hfill \\ \hfill \end{array}$

(90)$\begin{array}{ccc}\hfill {\mathbb{Q}}_{LN}^{\left(m\right)}& =& \left[\begin{array}{c}{\mathsf{T}}_{1m}^{\left(N\right)}\hfill \\ {\mathsf{P}}_{2m}^{\left(N\right)}\hfill \\ ⋮\hfill \\ {\mathsf{P}}_{L\mathbf{-}1,m}^{\left(N\right)}\hfill \\ {\mathsf{T}}_{Lm}^{\left(N\right)}\hfill \end{array}\right],{\overline{\mathbb{Q}}}_{LN}^{\left(m\right)}=\left[\begin{array}{c}{\mathsf{P}}_{1m}^{\left(N\right)}\hfill \\ {\mathsf{T}}_{2m}^{\left(N\right)}\hfill \\ ⋮\hfill \\ {\mathsf{T}}_{L\mathbf{-}1,m}^{\left(N\right)}\hfill \\ {\mathsf{P}}_{Lm}^{\left(N\right)}\hfill \end{array}\right],\left(L\mathrm{odd}\right).\hfill \end{array}$

Now, Equations (1) and (9) for $\mathit{u}$ and $\mathit{\omega }=\nabla \times \mathit{u}$ can be written, using (69), as

(91)$\begin{array}{ccc}\hfill \mathit{u}& =& \sum _m{(\mathbf{-}1)}^m\left({\mathbb{U}}_{LN}^{(\mathbf{-}m)†}{\mathbb{Q}}_{LN}^{\left(m\right)}+{\overline{\mathbb{U}}}_{LN}^{(\mathbf{-}m)†}{\overline{\mathbb{Q}}}_{LN}^{\left(m\right)}\right),\hfill \\ \hfill \end{array}$

(92)$\begin{array}{ccc}\hfill \mathit{\omega }& =& \sum _m{(\mathbf{-}1)}^m\left({\overline{\mathbb{U}}}_{LN}^{(\mathbf{-}m)†}{\mathbb{K}}_{LN}{\mathbb{Q}}_{LN}^{\left(m\right)}+{\mathbb{U}}_{LN}^{(\mathbf{-}m)†}{\mathbb{K}}_{LN}{\overline{\mathbb{Q}}}_{LN}^{\left(m\right)}\right).\hfill \end{array}$

Here, the factor of ${(\mathbf{-}1)}^m$ appears because of (8) and (13).

The arrays ${\mathbb{Q}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{Q}}}_{LN}^{\left(m\right)}$ transform differently from the rule (71) that applies to ${\mathbb{U}}_{LN}^{\left(m\right)}$ and ${\overline{\mathbb{U}}}_{LN}^{\left(m\right)}$:

(93)$\begin{array}{ccc}\hfill {\mathbb{R}}_{LN}^{\left(m\right)}={\mathbb{S}}_{LN}^{\left(m\right)†}{\mathbb{Q}}_{LN}^{\left(m\right)},& & {\overline{\mathbb{R}}}_{LN}^{\left(m\right)}={\overline{\mathbb{S}}}_{LN}^{\left(m\right)†}{\overline{\mathbb{Q}}}_{LN}^{\left(m\right)},\hfill \\ \hfill \end{array}$

(94)$\begin{array}{ccc}\hfill {\mathbb{Z}}_{LN}^{\left(m\right)}={\overline{\mathbb{S}}}_{LN}^{\left(m\right)†}{\mathbb{K}}_{LN}{\mathbb{Q}}_{LN}^{\left(m\right)},& & {\overline{\mathbb{Z}}}_{LN}^{\left(m\right)}={\mathbb{S}}_{LN}^{\left(m\right)†}{\mathbb{K}}_{LN}{\overline{\mathbb{Q}}}_{LN}^{\left(m\right)}.\hfill \end{array}$

In the case $m=0$, the elements of ${\mathbb{U}}_{LN}^{\left(0\right)}$ and ${\overline{\mathbb{U}}}_{LN}^{\left(0\right)}$ are real, while those of ${\mathbb{S}}_{LN}^{\left(0\right)}$ and ${\overline{\mathbb{S}}}_{LN}^{\left(0\right)}$ are generally complex. Under these transformations, (91) and (92) become

(95)$\begin{array}{ccc}\hfill \mathit{u}& =& \sum _m{(\mathbf{-}1)}^m\left({\mathbb{V}}_{LN}^{(\mathbf{-}m)†}{\mathbb{R}}_{LN}^{\left(m\right)}+{\overline{\mathbb{V}}}_{LN}^{(\mathbf{-}m)†}{\overline{\mathbb{R}}}_{LN}^{\left(m\right)}\right),\hfill \\ \hfill \end{array}$

(96)$\begin{array}{ccc}\hfill \mathit{\omega }& =& \sum _m{(\mathbf{-}1)}^m\left({\overline{\mathbb{V}}}_{LN}^{(\mathbf{-}m)†}{\mathbb{Z}}_{LN}^{\left(m\right)}+{\mathbb{V}}_{LN}^{(\mathbf{-}m)†}{\overline{\mathbb{Z}}}_{LN}^{\left(m\right)}\right).\hfill \end{array}$

However, in this example with $L=N=5$ and $m=0$, the important elements of ${\mathbb{S}}_{55}^{\left(0\right)}$ and ${\overline{\mathbb{S}}}_{55}^{\left(0\right)}$ associated with $\Re {\omega }_{55}^{\left(0\right)}=0$ and $\Re {\overline{\omega }}_{55}^{\left(0\right)}=0$ are real (and positive) as Table 3 shows. We can use these values and (93) to find the approximate form of the physical space eigenfunctions ${\mathbb{R}}_{55}^{\left(0\right)}\left(I\right)$ and ${\overline{\mathbb{R}}}_{55}^{\left(0\right)}\left(I\right)$ associated with each ${\omega }_{55}^{\left(0\right)}\left(I\right)$, and ${\mathbb{Z}}_{55}^{\left(0\right)}\left(I\right)$ and ${\overline{\mathbb{Z}}}_{55}^{\left(0\right)}\left(I\right)$ associated with each ${\overline{\omega }}_{55}^{\left(0\right)}\left(I\right)$, $I=11,\dots ,15$. Since $m=0$, these have no $\phi$-dependence, plotting their shape on any $\phi$ = constant slice of the spherical shell will reveal their structure. For this case, we make use of (11)–(15) to see that $\widehat{\mathit{\phi }}·{\mathit{P}}_{l0n}=0$; thus, only the toroidal ${\mathit{T}}_{l0n}$ in (95) and (96) will contribute as $\widehat{\mathit{\phi }}·{\overline{\mathbb{R}}}_{55}^{\left(0\right)}=0$ and $\widehat{\mathit{\phi }}·{\overline{\mathbb{Z}}}_{55}^{\left(0\right)}=0$. Let us define

(97)$F_I(r,\theta )=\widehat{\mathit{\phi }}·{\mathbb{R}}_{55}^{\left(0\right)}\left(I\right),G_I(r,\theta )=\widehat{\mathit{\phi }}·{\mathbb{Z}}_{55}^{\left(0\right)}\left(I\right),$

(98)$h_{ln}(r,\theta )=\widehat{\mathit{\phi }}·{\mathit{T}}_{l0n}=\sqrt{\frac{2l+1}{l(l+1)}}\frac{{\widehat{g}}_l\left(k_{ln}r\right)}{\sqrt{4\pi }}\frac{\mathrm{d}}{\mathrm{d}\theta }{P}_l(cos\theta ),{h^{\prime }}_{ln}=k_{ln}{h}_{ln}.$

Here, $P_l\left(x\right)$ is, of course, a Legendre polynomial.

Now, we can examine the physical space eigenfunctions associated with the $\Re {\omega }_{55}^{\left(0\right)}=0$ and $\Re {\overline{\omega }}_{55}^{\left(0\right)}=0$ eigenmodes contained in the expansions for velocity (91) and vorticity (92) using the data in Table 3. In any $(r,\theta )$ slice, the velocity eigenfunctions will be denoted $f_n(r,\theta )$, $n=1,\dots ,5$, where

(99)$f_1(r,\theta )=F_{15}(r,\theta )=S_1^{15}{h}_{11}(r,\theta )+S_{11}^{15}{h}_{31}(r,\theta )+S_{21}^{15}{h}_{51}(r,\theta ),$

(100)$f_2(r,\theta )=F_{14}(r,\theta )=S_2^{14}{h}_{12}(r,\theta )+S_{12}^{14}{h}_{32}(r,\theta )+S_{22}^{14}{h}_{52}(r,\theta ),$

(101)$f_3(r,\theta )=F_{13}(r,\theta )=S_3^{13}{h}_{13}(r,\theta )+S_{13}^{13}{h}_{33}(r,\theta )+S_{23}^{13}{h}_{53}(r,\theta ),$

(102)$f_2(r,\theta )=F_{12}(r,\theta )=S_4^{12}{h}_{14}(r,\theta )+S_{14}^{12}{h}_{34}(r,\theta )+S_{24}^{12}{h}_{54}(r,\theta ),$

(103)$f_5(r,\theta )=F_{11}(r,\theta )=S_5^{11}{h}_{11}(r,\theta )+S_{15}^{11}{h}_{35}(r,\theta )+S_{25}^{11}{h}_{55}(r,\theta ).$

The vorticity eigenfunctions associated with these are ${f^{\prime }}_n(r,\theta )$, $n=1,\dots ,5$,

(104)$\begin{array}{ccc}\hfill {f^{\prime }}_1(r,\theta )=G_{15}(r,\theta )& =& {\overline{S}}_1^{15}{h^{\prime }}_{11}(r,\theta )+{\overline{S}}_{11}^{15}{h^{\prime }}_{31}(r,\theta )+{\overline{S}}_{21}^{15}{h^{\prime }}_{51}(r,\theta ),\hfill \\ \hfill \end{array}$

(105)$\begin{array}{ccc}\hfill {f^{\prime }}_2(r,\theta )=G_{14}(r,\theta )& =& {\overline{S}}_2^{14}{h^{\prime }}_{12}(r,\theta )+{\overline{S}}_{12}^{14}{h^{\prime }}_{32}(r,\theta )+{\overline{S}}_{22}^{14}{h^{\prime }}_{52}(r,\theta ),\hfill \\ \hfill \end{array}$

(106)$\begin{array}{ccc}\hfill {f^{\prime }}_3(r,\theta )=G_{13}(r,\theta )& =& {\overline{S}}_3^{13}{h^{\prime }}_{13}(r,\theta )+{\overline{S}}_{13}^{13}{h^{\prime }}_{33}(r,\theta )+{\overline{S}}_{23}^{13}{h^{\prime }}_{53}(r,\theta ),\hfill \\ \hfill \end{array}$

(107)$\begin{array}{ccc}\hfill {f^{\prime }}_2(r,\theta )=G_{12}(r,\theta )& =& {\overline{S}}_4^{12}{h^{\prime }}_{14}(r,\theta )+{\overline{S}}_{14}^{12}{h^{\prime }}_{34}(r,\theta )+{\overline{S}}_{24}^{12}{h^{\prime }}_{54}(r,\theta ),\hfill \\ \hfill \end{array}$

(108)$\begin{array}{ccc}\hfill {f^{\prime }}_5(r,\theta )=G_{11}(r,\theta )& =& {\overline{S}}_5^{11}{h^{\prime }}_{15}(r,\theta )+{\overline{S}}_{14}^{11}{h^{\prime }}_{34}(r,\theta )+{\overline{S}}_{24}^{11}{h^{\prime }}_{54}(r,\theta ).\hfill \end{array}$

Again, the numerical values of the $S_J^I$ and ${\overline{S}}_J^I$ used above can be found in Table 3.

The eigenfunctions $f_n(r,\theta )$ in (99)–(103) and ${f^{\prime }}_n(r,\theta )$ in (104)–(108) that represent the $n=1,\dots ,5$ physical space eigenmodes of velocity and vorticity are pictured in Figure 6. These are non-oscillatory flow patterns specific to the $L=N=5$ system that was used as an example; they are ordered by the smallest ($n=1$) to the largest ($n=5$) dissipation factors of the non-oscillatory eigenmodes, which are also the five smallest dissipation factors of the 25 eigenmodes of the whole system. What they indicate is that the non-oscillatory inertial wave eigenmodes take the form of circular vortex structures that surround the rotation axis of the spherical shell. This is in contrast to the columnar, axially aligned vortices that numerical geodynamo simulations of a turbulent magnetofluid show occurring outside the so-called tangent cylinder of the rotating shell [27,28,37]. However, it must be remembered that here we are dealing with linear inertial hydrodynamic waves, while the turbulent geodynamo simulations study a fully nonlinear magnetohydrodynamic (MHD) system. A connection to nonlinear MHD systems must await an analysis of the MHD analog of the hydrodynamic inertial wave Equation (1).

10. Conclusions

In this paper, we determined the inertial waves found in an incompressible fluid in a rotating spherical shell with ‘homogeneous’ boundary conditions (b.c.s). A dynamical system has been defined and it has been shown that this system, with definite values of L and N, has exact solutions, e.g., definite eigenfrequencies and eigenmodes. If L and N are changed, then we have a new system which may have different values for the eigenfrequencies and shapes of eigenmodes. In this context, the question can be asked as to whether eigenfrequencies converge, at least for the largest scales, to definite values for a given E, as L and N get larger and larger. The two largest-scale modes for L = N = 5 are identified by $I=15$ as seen in Figure 5; the eigenfrequencies of these two modes are $0+i0.056$ and $0+i0.085$ as given in Table 2; these have the smallest magnitude of all eigenfrequencies in the case L = N = 5. Furthermore, in Figure 2, it is evident that the smallest magnitude $\omega$ have Re $\omega$ = 0 for all odd L, while for increasing even values of L, Figure 2 suggests that smallest magnitude $\omega$ also have Re $\omega \to 0$. This indicates the convergence of the eigenfrequencies of the largest scale modes: $\omega \to 0$. This is interesting, as it appears that the largest-scale inertial waves are or become stationary as L and N increase. A fuller examination of this topic is beyond the scope of the present work and must be deferred.

In summary, we have found the analytical form of inertial waves in an incompressible, rotating fluid constrained by concentric inner and outer spherical surfaces with homogeneous boundary conditions by using a toroidal-poloidal basis equivalent to Chandrasekhar-Kendall functions. We have discovered formulas and procedures for determining eigenmodes and associated eigenfrequencies for both the ideal and dissipative cases. These eigenmodes are formed from linear combinations of the multipole components of the Bessel function-spherical harmonic series that represents the fluid velocity. We illustrated the use of our results with specific low-dimensional examples and discussed what may occur as one moves to larger systems. In particular, we have found in a specific example that the lowest dissipation eigenmodes, all of which are non-oscillatory, are vortex ring structures that encircle the rotation axis. Application of the approach presented here to the problem of determining magnetohydrodynamic inertial wave structure in a spherical shell with homogeneous b.c.s is expected to further increase our understanding of the geodynamo.

Finally, let us mention a few other points. First, we have outlined an approach to the case with no-slip b.c.s. and compressibility, where one would still use the vector expansion functions ${\mathit{T}}_{lmn}\left(\mathit{x}\right)$ and ${\mathit{P}}_{lmn}\left(\mathit{x}\right)$ in addition to introducing a new set of expansion functions suitable for including compressibility effects and enforcing the no-slip boundary condition. Second, the Galerkin basis used here is unique in that each term satisfies a Helmholtz equation, for the given b.c.s and incompressibility. One could use Chebyshev polynomials as is done by [1,2], or perhaps Legendre polynomials, but the individual terms no longer satisfy the Helmholtz equation, i.e., there are no longer unique wavenumbers for the modes. (Please see [38] for formulas related to Chebyshev and Legendre methods). Third, the results developed here could be used to analyze data from numerical simulations, most directly from simulations based on the Galerkin vector basis functions defined in the paper, but also, data from finite-difference or finite-volume methods could also examined, particularly, but not necessarily, if these other methods employ the same b.c.s.

Funding

This research received no external funding.

Data Availability Statement

The data presented here can be generated from the formulas given.

Conflicts of Interest

The author declares no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (colour online)The eigenfrequencies (a) [Forumla omitted. See PDF.] and (b) [Forumla omitted. See PDF.]. Here, the numbering with respect to I corresponds to that in Table 2 and Figure 5. The eigenfrequencies with negative values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], are found by reflecting [Forumla omitted. See PDF.] in (a) and (b) across the vertical axes [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 2. (colour online) Eigenfrequencies [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] of undamped eigenmodes associated with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] to 5. The smooth curve [Forumla omitted. See PDF.] fits the outer data but is only used for purpose of illustrating that the eigenfrequencies appear unbounded as [Forumla omitted. See PDF.]. However, dissipation ([Forumla omitted. See PDF.]) is eventually important and provides a bounding mechanism, as shown in Figure 3.

Figure 3. (colour online) Above, [Forumla omitted. See PDF.] signifies either set of eigenfrequencies (∘) [Forumla omitted. See PDF.] or (+) [Forumla omitted. See PDF.] of the damped eigenmodes associated with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. Notice that as E increases from (a–f), more and more eigenmodes become overdamped, i.e., have [Forumla omitted. See PDF.] become purely imaginary. When E is slightly larger than 1.25, all eigenmodes are overdamped; as E becomes smaller, overdamping starts to occur when [Forumla omitted. See PDF.].

Figure 4. (colour online) Above, [Forumla omitted. See PDF.] signifies either set of eigenfrequencies (∘) [Forumla omitted. See PDF.] or (+) [Forumla omitted. See PDF.] of the damped eigenmodes associated with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.], comparing the cases (a) [Forumla omitted. See PDF.] with (b) [Forumla omitted. See PDF.] and (c) [Forumla omitted. See PDF.].

Figure 5. (colour online) Magnitudes of components [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], of rows [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], of transformation matrices [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], respectively. These are for the case [Forumla omitted. See PDF.] and the index I = 1 to 15 is ordered to correspond to the largest values of [Forumla omitted. See PDF.] or [Forumla omitted. See PDF.] to their smallest values.

Figure 6. Non-oscillatory eigenmodes: Axisymmetric inertial wave structures in a rotating spherical shell. (a–e) Azimuthal components of velocity [Forumla omitted. See PDF.] and (f–j) corresponding azimuthal vorticity component [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], associated with the [Forumla omitted. See PDF.] eigenmodes I = 15–11 of Table 2 and Figure 5. These are helical flows encircling the rotation axis, which is the left vertical edge in each plot.

Table 1

Values for $k_{ln}$, $g_l\left(k_{ln}\right)=0$, for $l,n=1,\dots ,10$.

n =	1	2	3	4	5	6	B	8	9	10
$l=1$	1.8638	3.4929	5.1612	6.8434	8.5316	10.223	11.916	13.611	15.306	17.002
2	2.1497	3.6788	5.2927	6.9440	8.6129	10.291	11.975	13.662	15.352	17.043
3	2.5042	3.9411	5.4851	7.0931	8.7339	10.393	12.062	13.739	15.420	17.105
4	2.8910	4.2636	5.7330	7.2884	8.8935	10.527	12.178	13.841	15.511	17.187
5	3.2898	4.6294	6.0300	7.5271	9.0903	10.694	12.323	13.968	15.624	17.289
6	3.6911	5.0227	6.3684	7.8062	9.3227	10.892	12.494	14.119	15.759	17.411
7	4.0909	5.4310	6.7396	8.1218	9.5890	11.119	12.692	14.294	15.916	17.553
8	4.4882	5.8457	7.1347	8.4695	9.8870	11.376	12.917	14.493	16.094	17.714
9	4.8828	6.2618	7.5452	8.8441	10.214	11.661	13.166	14.715	16.293	17.894
10	5.2749	6.6767	7.9642	9.2399	10.569	11.972	13.441	14.959	16.513	18.093

Table 2

For the case $L=N=5$, $m=0$ and $E=0.01$: Values for the eigenfrequencies ${\omega }_{55}^{\left(0\right)}\left(I\right)$ and ${\overline{\omega }}_{55}^{\left(0\right)}\left(I\right)$ for $I=1,\dots ,15$. The index I = 1 to 15 is ordered to correspond to the largest values of $|{\omega }_{55}^{\left(0\right)}|$ or $|{\overline{\omega }}_{55}^{\left(0\right)}|$ to their smallest values. We define ${\omega }_{55}^{\left(0\right)}\left(I\right)=\mathbf{-}{\omega }_{55}^{\left(0\right)*}(j\mathbf{-}15)$ and ${\overline{\omega }}_{55}^{\left(0\right)}\left(I\right)=\mathbf{-}{\overline{\omega }}_{55}^{\left(0\right)*}(I\mathbf{-}15)$, $I=16,\dots ,25$.

I =	1	2	3	4	5
${\omega }_{55}^{\left(0\right)}\left(I\right)=$	2.63 + i0.09	2.03 + i0.20	1.84 + i0.39	1.62 + i0.54	1.42 + i0.093
${\overline{\omega }}_{55}^{\left(0\right)}\left(I\right)=$	2.61 + i0.09	2.02 + i0.21	1.82 + i0.39	1.67 + i0.55	1.35 + i0.087
$\mathit{I}=$	6	7	8	9	10
${\omega }_{55}^{\left(0\right)}\left(I\right)=$	1.03 + i0.61	1.11 + i0.20	1.01 + i0.40	0.88 + i0.55	0.56 + i0.55
${\overline{\omega }}_{55}^{\left(0\right)}\left(I\right)=$	1.06 + i0.59	0.99 + i0.63	1.08 + i0.20	0.96 + i0.37	0.35 + i0.59
$\mathit{I}=$	11	12	13	14	15
${\omega }_{55}^{\left(0\right)}\left(I\right)=$	0 + i0.76	0 + i0.51	0 + i0.30	0 + i0.15	0 + i0.056
${\overline{\omega }}_{55}^{\left(0\right)}\left(I\right)=$	0 + i0.63	0 + i0.51	0 + i0.31	0 + i0.17	0 + i0.085

Table 3

First, to minimize confusion, remember that $L=5$ and $N=5$; second, I is just an arbitrary numbering of the $L\times N=25$ eigenmodes; third, $J=N(l\mathbf{-}1)+n$, with $1\le l\le L$ and $1\le n\le N$. For $I=11,\cdots ,15$ of Figure 1 and Figure 5, this table has the values for the largest three $S_J^I$ and ${\overline{S}}_J^I$ for each I. These are the major matrix elements that define the $\Re {\omega }_{55}^{\left(0\right)}=0$ and $\Re {\overline{\omega }}_{55}^{\left(0\right)}=0$ eigenmodes, along with their associated wavenumbers $k_{ln}$, which are drawn from Table 1.

I	${\mathit{k}}_{\mathit{ln}}$	${\mathit{S}}_{\mathit{J}}^{\mathit{I}}$	${\mathit{k}}_{\mathit{ln}}$	${\mathit{S}}_{\mathit{J}}^{\mathit{I}}$	${\mathit{k}}_{\mathit{ln}}$	${\mathit{S}}_{\mathit{J}}^{\mathit{I}}$
15	$k_{11}=1.864$,	$S_1^{15}=0.7612$	$k_{31}=2.504$,	$S_{11}^{15}=0.5007$	$k_{51}=3.290$,	$S_{21}^{15}=0.4009$
14	$k_{12}=3.493$,	$S_2^{14}=0.6889$	$k_{32}=3.941$,	$S_{12}^{14}=0.5501$	$k_{52}=4.629$,	$S_{22}^{14}=0.4547$
13	$k_{13}=5.161$,	$S_3^{13}=0.6797$	$k_{33}=5.485$,	$S_{13}^{13}=0.5463$	$k_{53}=6.030$,	$S_{23}^{13}=0.4648$
12	$k_{14}=6.843$,	$S_4^{12}=0.6494$	$k_{34}=7.093$,	$S_{14}^{12}=0.5597$	$k_{54}=7.527$,	$S_{24}^{12}=0.4977$
11	$k_{15}=8.532$,	$S_5^{11}=0.6583$	$k_{35}=8.734$,	$S_{15}^{11}=0.5528$	$k_{55}=9.090$,	$S_{25}^{11}=0.4833$
I	$k_{ln}$	${\overline{S}}_J^I$	$k_{ln}$	${\overline{S}}_J^I$	$k_{ln}$	${\overline{S}}_J^I$
15	$k_{11}=1.864$,	${\overline{S}}_1^{15}=0.5851$	$k_{31}=2.504$,	${\overline{S}}_{11}^{15}=0.5302$	$k_{51}=3.290$,	${\overline{S}}_{21}^{15}=0.5608$
14	$k_{12}=3.493$,	${\overline{S}}_2^{14}=0.6136$	$k_{32}=3.941$,	${\overline{S}}_{12}^{14}=0.5466$	$k_{52}=4.629$,	${\overline{S}}_{22}^{14}=0.5322$
13	$k_{13}=5.161$,	${\overline{S}}_3^{13}=0.6491$	$k_{33}=5.485$,	${\overline{S}}_{13}^{13}=0.5359$	$k_{53}=6.030$,	${\overline{S}}_{23}^{13}=0.5063$
12	$k_{14}=6.843$,	${\overline{S}}_4^{12}=0.6859$	$k_{34}=7.093$,	${\overline{S}}_{14}^{12}=0.4970$	$k_{54}=7.527$,	${\overline{S}}_{24}^{12}=0.4551$
11	$k_{15}=8.532$,	${\overline{S}}_5^{11}=0.7996$	$k_{34}=7.093$,	${\overline{S}}_{14}^{11}=0.3103$	$k_{54}=7.527$,	${\overline{S}}_{24}^{11}=0.3080$

References

1. Glatzmaier, G.A.; Roberts, P.H. A three-dimensional self-consistent computer simulation of a geomagnetic field reversal. Nature; 1995; 377, pp. 203-209. [DOI: https://dx.doi.org/10.1038/377203a0]

2. Glatzmaier, G.A.; Roberts, P.H. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle. Phys. Earth Planet. Int.; 1995; 91, pp. 63-75. [DOI: https://dx.doi.org/10.1016/0031-9201(95)03049-3]

3. Elsässer, W.M. Induction Effects in Terrestrial Magnetism Part I. Theory. Phys. Rev.; 1946; 69, pp. 106-112. [DOI: https://dx.doi.org/10.1103/PhysRev.69.106]

4. Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. Spectral Methods in Fluid Dynamics; Springer: New York, NY, USA, 1987.

5. Shebalin, J.V. Broken ergodicity, magnetic helicity, and the MHD dynamo. Geophys. Astrophys. Fluid Dyn.; 2013; 107, pp. 353-375. [DOI: https://dx.doi.org/10.1080/03091929.2012.689299]

6. Mininni, P.D.; Montgomery, D. Magnetohydrodynamic activity inside a sphere. Phys. Fluids; 2006; 18, 116602. [DOI: https://dx.doi.org/10.1063/1.2393438]

7. Mininni, P.D.; Montgomery, D.; Turner, L. Hydrodynamic and magnetohydrodynamic computations inside a rotating sphere. New J. Phys.; 2007; 9, pp. 303-330. [DOI: https://dx.doi.org/10.1088/1367-2630/9/8/303]

8. Chandrasekhar, S.; Kendall, P.C. On Force-Free Magnetic Fields. Astrophys. J.; 1957; 12, pp. 457-460. [DOI: https://dx.doi.org/10.1086/146413]

9. Morse, E.C. Eigenfunctions of the curl in cylindrical geometry. J. Math. Phys.; 2005; 46, 113511. [DOI: https://dx.doi.org/10.1063/1.2118447]

10. Morse, E.C. Eigenfunctions of the curl in annular cylindrical and rectangular geometry. J. Math. Phys.; 2007; 48, 083504. [DOI: https://dx.doi.org/10.1063/1.2760391]

11. Livermore, P.W.; Bailey, L.M.; Hollerbach, R. A comparison of no-slip, stress-free and inviscid models of rapidly rotating fluid in a spherical shell. Sci. Rep.; 2016; 6, 22812. [DOI: https://dx.doi.org/10.1038/srep22812]

12. Batchelor, G.K. The Theory of Homogeneous Turbulence; Cambridge University Press: Cambridge, UK, 1953; pp. 28-33.

13. Ciufolini, I.; Wheeler, J.A. Gravitation and Inertia; Princeton University Press: Princeton, NJ, USA, 1995; pp. 224-226.

14. Luminet, J.-P. The Wraparound Universe; Peters, A.K. CRC Press: Wellsley, MA, USA, 2008.

15. di Leoni, P.; Cobelli, P.J.; Mininni, P.D.; Dmitruk, P.; Matthaeus, W.H. Quantification of the strength of inertial waves in a rotating turbulent flow. Phys. Fluids; 2014; 26, 035106. [DOI: https://dx.doi.org/10.1063/1.4868280]

16. Bucciantini, N. Polarization properties of turbulent synchrotron bubbles: An approach based on Chandrasekhar–Kendall functions. Mon. Not. Roy. Astron. Soc.; 2017; 471, pp. 4885-4893. [DOI: https://dx.doi.org/10.1093/mnras/stx1927]

17. Greenspan, H.G. The Theory of Rotating Fluids; Cambridge University Press: Cambridge, UK, 1968; pp. 63-68.

18. Zhang, K.; Earnshaw, P.; Xinhao, L.; Busse, F.H. On inertial waves in a rotating fluid sphere. J. Fluid Mech.; 2001; 437, pp. 103-119. [DOI: https://dx.doi.org/10.1017/S0022112001004049]

19. Bryan, G.H. The waves on a rotating liquid spheroid of finite ellipticity. Phil. Trans. R. Soc. Lond. A; 1889; 180, pp. 187-219.

20. Stewartson, K.; Rickard, J.A. Inertial waves in a rotating spherical shell: Attractors and asymptotic spectrum. J. Fluid Mech.; 2001; 435, pp. 103-144.

21. Aldridge, K.D. Axisymmetric Inertial Oscillations of a Fluid in a Rotating Spherical Shell. Mathematika; 1972; 19, pp. 163-168. [DOI: https://dx.doi.org/10.1112/S002557930000560X]

22. Aldridge, K.D. Inertial Waves and the Earth’s Outer Core. Geophys. J. R. Astro. Soc.; 1975; 42, pp. 337-345. [DOI: https://dx.doi.org/10.1111/j.1365-246X.1975.tb05865.x]

23. Rieutord, M.; Valdettaro, L. Inertial waves in a rotating spherical shell. J. Fluid Mech.; 1997; 341, pp. 77-99. [DOI: https://dx.doi.org/10.1017/S0022112097005491]

24. Finlay, C.C. Waves in the Presence of Magnetic Fields, Rotation and Convection. Dynamos: Lecture Notes of the Les Houches Summer School 2007; Cardin, P.; Cugliandolo, L.F. Elsevier: Amsterdam, The Netherlands, 2008; pp. 403-450.

25. Shebalin, J.V. Magnetic Helicity and the Geodynamo. Fluids; 2021; 6, 99. [DOI: https://dx.doi.org/10.3390/fluids6030099]

26. Kuang, W.; Bloxham, J. An Earth-like numerical dynamo model. Nature; 1997; 389, pp. 371-374. [DOI: https://dx.doi.org/10.1038/38712]

27. Olson, P. Overview. Core Dynamics; Schubert, G.; Olson, P. Elsevier: Amsterdam, The Netherlands, 2009; 6.

28. Christensen, U.R.; Wicht, J. Numerical Dynamo Simulations. Core Dynamics, Schubert, G., Olson, P., Eds.; Elsevier: Amsterdam, The Netherlands, 2009; pp. 245-282.

29. Dormy, E.; Soward, A.M. Mathematical Aspects of Natural Dynamos; CRC Press: London, UK, 2007.

30. Betchov, R. Semi-isotropic turbulence and helicoidal flows. Phys. Fluids; 1961; 4, pp. 925-926. [DOI: https://dx.doi.org/10.1063/1.1706426]

31. Barrera, R.G. Estévez, G.A.; Giraldo, J. Vector spherical harmonics and their application to magnetostatics. Eur. J. Phys.; 1995; 6, pp. 287-294. [DOI: https://dx.doi.org/10.1088/0143-0807/6/4/014]

32. Jackson, J.D. Classical Electrodynamics; 2nd ed. Wiley: New York, NY, USA, 1975; pp. 99-100.

33. Arfken, G.B.; Weber, H.J. Mathematical Methods for Physicists; 5th ed. Academic Press: New York, NY, USA, 2005; pp. 725-739.

34. Chandrasekhar, S.; Elbert, D. The Roots of J_−(l+½)(λη)J_l+½(λ)−J_l+½(λη)J_−(l+½)(λ) = 0. Mon. Proc. Camb. Phil. Soc.; 1953; 49, pp. 446-448. [DOI: https://dx.doi.org/10.1017/S0305004100028619]

35. Liu, H.; Zou, J. Zeros of the Bessel and spherical Bessel functions and their applications for uniqueness in inverse acoustic obstacle scattering. IMA J. Appl. Math.; 2007; 72, pp. 817-831. [DOI: https://dx.doi.org/10.1093/imamat/hxm013]

36. Kraichnan, R.H. Helical turbulence and absolute equilibrium. J. Fluid Mech.; 1973; 59, pp. 745-752. [DOI: https://dx.doi.org/10.1017/S0022112073001837]

37. Tilgner, A. Rotational Dynamics of the Core. Core Dynamics; Schubert, G.; Olson, P. Elsevier: Amsterdam, The Netherlands, 2009; pp. 207-243.

38. Gottlieb, D.; Orszag, S.A. Numerical Analysis of Spectral Methods: Theory and Applications; SIAM: Philadelphia, PA, USA, 1977.