# Spatiospectral Concentration on a Sphere

Frederik Simons, F. Dahlen, Mark Wieczorek

► To cite this version:

Frederik Simons, F. Dahlen, Mark Wieczorek. Spatiospectral Concentration on a Sphere. SIAM Review, 2006, 48 (3), pp.504-536. 10.1137/S0036144504445765. hal-02458541

HAL Id: hal-02458541

https://hal.science/hal-02458541

Submitted on 29 Jan 2023

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

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

<figure><img src="image_1.png" /></figure>
# Spatiospectral Concentration on a Sphere*

Frederik J. Simons†

F. A. Dahlen‡

Mark A. Wieczorek§

**Abstract.** We pose and solve the analogue of Slepian's time-frequency concentration problem on the surface of the unit sphere to determine an orthogonal family of strictly bandlimited functions that are optimally concentrated within a closed region of the sphere or, alternatively, of strictly spacelimited functions that are optimally concentrated in the spherical harmonic domain. Such a basis of simultaneously spatially and spectrally concentrated functions should be a useful data analysis and representation tool in a variety of geophysical and planetary applications, as well as in medical imaging, computer science, cosmology, and numerical analysis. The spherical Slepian functions can be found by solving either an algebraic eigenvalue problem in the spectral domain or a Fredholm integral equation in the spatial domain. The associated eigenvalues are a measure of the spatiospectral concentration. When the concentration region is an axisymmetric polar cap, the spatiospectral projection operator commutes with a Sturm–Liouville operator; this enables the eigenfunctions to be computed extremely accurately and efficiently, even when their area-bandwidth product, or Shannon number, is large. In the asymptotic limit of a small spatial region and a large spherical harmonic bandwidth, the spherical concentration problem reduces to its planar equivalent, which exhibits self-similarity when the Shannon number is kept invariant.

**Key words.** bandlimited function, commuting differential operator, concentration problem, eigenvalue problem, multitaper spectral analysis, spherical harmonics

**AMS subject classifications.** 42B05, 42B35, 45B05, 47B32

**DOI.** 10.1137/S0036144504445765

**I. Introduction.** In a classic series of papers published in the 1960s, Slepian, Landau, and Pollak solved a fundamental problem in information theory, namely, that of optimally concentrating a given signal in both the time and frequency domains [35, 36, 58, 59]. The orthogonal family of data windows, or tapers, that arise in this context, and their discrete and multidimensional extensions [5, 22, 38, 56, 57], form the basis of the multitaper method of spectral analysis [69, 70, 71], which has enjoyed application in a wide range of physical, computational, and biomedical disciplines (e.g., climatology, communications engineering, geodesy, neurology, optics,

<footnote>*Received by the editors August 23, 2004; accepted for publication (in revised form) August 23, 2005; published electronically August 1, 2006. Financial support for this work was provided by the U.S. National Science Foundation under grant EAR-0105387. This is IPGP contribution 2084. http://www.siam.org/journals/sirev/48-3/44576.html</footnote>

<footnote>†Department of Geosciences, Princeton University, Princeton, NJ 08544. Current address: Department of Earth Sciences, University College London, London WC1E 6BT, UK (fjsimons@alum.mit.edu).</footnote>

<footnote>‡Department of Geosciences, Princeton University, Princeton, NJ 08544.</footnote>

<footnote>§Département de Géophysique Spatiale et Planétaire, Institut de Physique du Globe de Paris, 94701 St. Maur-des-Fossés, France.</footnote>
seismology, and speech recognition). Time-frequency and time-scale concentration in more general settings and a variety of geometries has subsequently been studied by several authors, e.g., [7, 11, 13, 14, 16, 17, 18, 21, 37, 40, 43, 44, 51].

Adhering to Slepian's original quadratic concentration criterion, we consider the simultaneous spatial and spectral localization of a real-valued function of geographical position on the surface of the unit sphere. Apart from Gilbert and Slepian's work on concentrated Legendre polynomials [20] and the extremely insightful analysis by Grünbaum and his colleagues [21], Slepian's concentration problem on a sphere has been revisited relatively rarely, by workers in geodesy [1], medical imaging [50], and planetary science [76]. These studies all restrict attention to the special case of an axisymmetric concentration region; here, we pose and solve the spherical spatiospectral concentration problem for a geographical region of arbitrary shape.

The optimal orthogonal family of spherical multitapers that we derive here will be useful in a variety of geophysical, planetary, cosmological, and other applications; the single tapers that have recently been developed for some of the purposes listed above, e.g., [19, 32, 54, 63, 68], are inferior in the extraction of localized statistical information [49, 74, 76]. In geophysics, the local thickness or elastic strength of the terrestrial lithosphere can be estimated from the cross-spectral properties of Earth's surface topography and associated gravitational field [73]; in planetary science, the lithospheric properties of other bodies in the solar system may be investigated in the same way. The data used in these studies are most commonly available as a bandlimited set of spherical harmonic coefficients measured by artificial satellites. In many if not most cases, planetary curvature prohibits the use of locally flat approximations [75]. Thus, the determination of spatially localized estimates of terrestrial and planetary lithospheric properties requires spatiospectral localization methods that go beyond those available in the plane, e.g., [53]. The construction of bandlimited spherical basis functions that are optimally concentrated within a specific geographical region also has applications in satellite hydrology [64], which aims to constrain fluxes of surface and ground water by measurements of the associated temporal variations in Earth's gravity field, and in geodesy, where one is frequently required to determine the spherical harmonic coefficients of a planet's gravity or magnetic field using data from an incompletely sampled sphere [61]. In astronomy and cosmology, spherical multitapers will be useful in estimating the power spectra of extragalactic objects cataloged in sky surveys [23, 48] as well as the spectrum of the cosmic microwave background radiation, either from ground-based surveys of a limited region of the sky or from space-based measurements that are contaminated by emission from the galactic plane [66, 67].

**2. Slepian's Concentration Problem.** We begin with a brief review of the one-dimensional, continuous-continuous, time-frequency concentration problem in order to provide a template for the spherical spatiospectral concentration problem, which we shall consider in the remainder of the paper. We use $t$ and $\omega$ to denote time and angular frequency, respectively, and adopt a normalization convention in which a real-valued time-domain signal $f(t)$ and its Fourier transform $F(\omega)$ are related by

$$ (2.1) \qquad f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} d\omega, \qquad F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} dt. $$

The problem considered by Slepian and Pollak [59] is that of optimally concentrating a strictly bandlimited signal $g(t)$, with a spectrum $G(\omega)$ that vanishes for frequencies $|\omega| > W$, into a time interval $|t| \le T$. No such bandlimited signal $g(t)$ can be
completely concentrated within a finite interval by virtue of the Paley–Wiener theorem [12, 39]. Although alternative criteria have been analyzed, e.g., [19, 51], the optimally concentrated signal is considered to be the one with the least energy outside the interval:

$$ (2.2) \qquad \lambda = \frac{\int_{-T}^{T} g^2(t) \, dt}{\int_{-\infty}^{\infty} g^2(t) \, dt} = \text{maximum.} $$

Bandlimited signals $g(t)$ satisfying the variational problem (2.2) have spectra $G(\omega)$ that satisfy the frequency-domain convolutional integral eigenvalue equation

$$ (2.3) \qquad \int_{-W}^{W} \frac{\sin T(\omega - \omega')}{\pi(\omega - \omega')} G(\omega') \, d\omega' = \lambda G(\omega), \quad |\omega| \le W. $$

A closely related problem is that of concentrating the spectrum $H(\omega)$ of a strictly timelimited function $h(t)$, which vanishes for times $|t| > T$, into a spectral interval $|\omega| \le W$. Slepian's measure of concentration in this case is

$$ (2.4) \qquad \lambda = \frac{\int_{-W}^{W} |H(\omega)|^2 \, d\omega}{\int_{-\infty}^{\infty} |H(\omega)|^2 \, d\omega} = \text{maximum.} $$

Timelimited signals $h(t)$ whose spectra satisfy the variational problem (2.4) themselves satisfy the time-domain eigenvalue equation

$$ (2.5) \qquad \int_{-T}^{T} \frac{\sin W(t-t')}{\pi(t-t')} h(t') \, dt' = \lambda h(t), \quad |t| \le T. $$

Both (2.3) and (2.5) have the same eigenvalues $1 > \lambda_1 > \lambda_2 > \dots > 0$, with associated time-domain eigentapers $g_1(t), g_2(t), \dots$ and $h_1(t), h_2(t), \dots$, which coincide, to within a multiplicative constant, within the interval $|t| \le T$, and eigenspectra $G_1(\omega), G_2(\omega), \dots$ and $H_1(\omega), H_2(\omega), \dots$, which coincide within the interval $|\omega| \le W$.

A change of both the independent and dependent variables transforms both (2.3) and (2.5) into the same dimensionless eigenvalue equation:

$$ (2.6) \qquad \int_{-1}^{1} \frac{\sin TW(x-x')}{\pi(x-x')} \psi(x') \, dx' = \lambda \psi(x), \quad |x| \le 1. $$

Equation (2.6) shows that the eigenvalues $\lambda_1, \lambda_2, \dots$ and suitably scaled eigenfunctions $\psi_1(x), \psi_2(x), \dots$ depend only upon the time-bandwidth product $TW$. The sum of the eigenvalues is related to this product by

$$ (2.7) \qquad N = \sum_{\alpha=1}^{\infty} \lambda_{\alpha} = \frac{2TW}{\pi}. $$

Because of the characteristic step shape of the eigenvalue spectrum [33, 60], this so-called Shannon number [49] is a good estimate of the number of significant eigenvalues
or, roughly speaking, the number of signals $f(t)$ that can be simultaneously well concentrated into a finite time interval $|t| \le T$ and a finite frequency interval $|\omega| \le W$.

Serendipitously, the integral operator acting upon $\psi$ in (2.6) commutes with a second-order differential operator, $\mathcal{P} = \frac{d}{dx}(1-x^2)\frac{d}{dx} - T^2W^2x^2$, which arises in the separation of the three-dimensional scalar wave equation in prolate spheroidal coordinates [55]. Thanks to this “lucky accident” [58], it is also possible to find the scaled eigenfunctions $\psi_1(x), \psi_2(x), \dots$ by solving the Sturm–Liouville equation

$$ (2.8) \qquad \frac{d}{dx}(1-x^2)\frac{d\psi}{dx} + (\chi - T^2W^2x^2)\psi = 0, \quad |x| \le 1, $$

where $\chi \neq \lambda$ is the associated eigenvalue.

The bandlimited prolate spheroidal eigentapers may be chosen to be orthonormal over the infinite time interval $|t| \le \infty$ and orthogonal over the finite interval $|t| \le T$:

$$ (2.9) \qquad \int_{-\infty}^{\infty} g_{\alpha}g_{\beta} dt = \delta_{\alpha\beta}, \qquad \int_{-T}^{T} g_{\alpha}g_{\beta} dt = \lambda_{\alpha}\delta_{\alpha\beta}. $$

Almost all of the above results have analogues in the spati spectral concentration problem for functions on the unit sphere. As we shall see, this two-dimensional problem is enriched by the arbitrary shape of the region of spatial concentration.

**3. Preliminaries.** The geometry of the unit sphere $\Omega = \{\hat{\mathbf{r}} : \|\hat{\mathbf{r}}\| = 1\}$ is depicted in Figure 3.1. We denote the colatitude of a geographical point $\hat{\mathbf{r}}$ by $0 \le \theta \le \pi$ and

<figure><img src="image_2.png" /><figcaption>**Fig. 3.1** Sketch illustrating the geometry of the spherical concentration problem. Lower right shows an axisymmetric polar cap of colatitudinal radius $\Theta$, treated in Section 5. The area of the region of concentration, $R = R_1 \cup R_2 \cup \dots$, is denoted by A in the text.</figcaption></figure>
the longitude by $0 \le \phi < 2\pi$; the geodesic angular distance between two points $\hat{\mathbf{r}}$ and $\hat{\mathbf{r}}'$ will be denoted by $\Delta$, where $\cos \Delta = \hat{\mathbf{r}} \cdot \hat{\mathbf{r}}' = \cos \theta \cos \theta' + \sin \theta \sin \theta' \cos(\phi - \phi')$.
We use $R$ to denote a region of $\Omega$, of area $A = \int_R d\Omega$, within which we seek to concentrate a bandlimited function of position $\hat{\mathbf{r}}$. The region may consist of a number of unconnected subregions, $R = R_1 \cup R_2 \cup \dots$, and it may have an irregularly shaped boundary, as shown. The region complementary to $R$ will be denoted by $\Omega - R$.

**3.1. Spherical Harmonics.** Since we restrict our attention to real-valued functions, we use real surface spherical harmonics $Y_{lm}(\hat{\mathbf{r}}) = Y_{lm}(\theta, \phi)$ defined by [10, 15]

$$ (3.1) \qquad Y_{lm}(\theta, \phi) = \begin{cases} \sqrt{2}X_{lm}(\theta) \cos m\phi & \text{if } -l \le m < 0, \\ X_{l0}(\theta) & \text{if } m = 0, \\ \sqrt{2}X_{lm}(\theta) \sin m\phi & \text{if } 0 < m \le l, \end{cases} $$

$$ (3.2) \qquad X_{lm}(\theta) = (-1)^m \left(\frac{2l+1}{4\pi}\right)^{1/2} \left[\frac{(l-m)!}{(l+m)!}\right]^{1/2} P_{lm}(\cos \theta), $$

$$ (3.3) \qquad P_{lm}(\mu) = \frac{1}{2^l l!} (1 - \mu^2)^{m/2} \left(\frac{d}{d\mu}\right)^{l+m} (\mu^2 - 1)^l. $$

The quantity $0 \le l \le \infty$ is known as the angular degree of the spherical harmonic, and $-l \le m \le l$ is its angular order. The $l \to \infty$ asymptotic wavenumber associated with a harmonic of degree $l$ is $\sqrt{l(l+1)} \approx l + 1/2$ [29]. The function $P_{lm}(\mu)$ defined in (3.3) is the associated Legendre function of integer degree $l$ and order $m$. The spherical harmonics $Y_{lm}(\hat{\mathbf{r}})$ are eigenfunctions of the Laplace-Beltrami operator, $\nabla^2 Y_{lm} = -l(l+1)Y_{lm}$, where $\nabla^2 = \partial_\theta^2 + \cot \theta \partial_\theta + (\sin \theta)^{-2} \partial_\phi^2$. Our choice of the constants in (3.1)-(3.2) orthonormalizes the harmonics on the unit sphere:

$$ (3.4) \qquad \int_{\Omega} Y_{lm} Y_{l'm'} d\Omega = \delta_{ll'} \delta_{mm'}. $$

The corresponding fixed-order orthogonality relations for $X_{lm}(\theta)$ and $P_{lm}(\mu)$ are

$$ (3.5) \qquad \int_{0}^{\pi} X_{lm} X_{l'm} \sin \theta \, d\theta = \frac{1}{2\pi} \delta_{ll'}, \qquad \int_{-1}^{1} P_{lm} P_{l'm} \, d\mu = \frac{2}{2l+1} \frac{(l+m)!}{(l-m)!} \delta_{ll'}. $$

The integral of a Legendre polynomial $P_l(\mu) = P_{l0}(\mu)$ over a cap $\cos \Theta \le \mu \le 1$ is [6]

$$ (3.6) \qquad \int_{\cos \Theta}^{1} P_l \, d\mu = \frac{1}{2l+1} [P_{l-1}(\cos \Theta) - P_{l+1}(\cos \Theta)], $$

where $P_{-1}(\mu) \equiv 1$, and the product of two fixed-order Legendre functions is

$$ (3.7) \qquad X_{lm}(\theta) X_{l'm}(\theta) = (-1)^m \sum_{n=|l-l'|}^{l+l'} \sqrt{\frac{(2n+1)(2l+1)(2l'+1)}{4\pi}} \\ \qquad \times \begin{pmatrix} l & n & l' \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l & n & l' \\ m & 0 & -m \end{pmatrix} X_{n0}(\theta), $$

where the arrays of indices are Wigner 3-j symbols [15, 41]. We shall use two recursion relations involving the associated Legendre functions and their derivatives, namely,

$$ (3.8a) \qquad (2l+1)\mu P_{lm} = (l-m+1)P_{l+1m} + (l+m)P_{l-1m}, $$

$$ (3.8b) \qquad (1 - \mu^2) \frac{dP_{lm}}{d\mu} = (l+1)\mu P_{lm} - (l-m+1)P_{l+1m}. $$
Finally, two relations involving sums of products of Legendre functions evaluated at different arguments are useful in the discussion that follows. The first is the well-known spherical harmonic addition theorem [15]

$$ (3.9) \qquad \sum_{m=-l}^{l} Y_{lm}(\hat{\mathbf{r}})Y_{lm}(\hat{\mathbf{r}}') = \left(\frac{2l+1}{4\pi}\right) P_l(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}') = \left(\frac{2l+1}{4\pi}\right) P_l(\cos \Delta), $$

and the second is the Legendre version of the Christoffel–Darboux identity [28, 62, 65]

$$ (3.10) \qquad (\mu - \mu') \sum_{l=m}^{L} (2l+1) \frac{(l-m)!}{(l+m)!} P_{lm}(\mu)P_{lm}(\mu') \\ = \frac{(L-m+1)!}{(L+m)!} [P_{L+1m}(\mu)P_{Lm}(\mu') - P_{Lm}(\mu)P_{L+1m}(\mu')], $$

which is a straightforward consequence of the three-term relation (3.8a).

**3.2. Functions on the Sphere.** Let $f(\hat{\mathbf{r}})$ be a real-valued, square-integrable function on the unit sphere $\Omega$. Any such function can be expanded as

$$ (3.11) \qquad f = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} f_{lm} Y_{lm}, \qquad f_{lm} = \int_{\Omega} f Y_{lm} d\Omega. $$

Equation (3.11) is the spherical analogue of the one-dimensional Fourier transform pair (2.1). The finite character of the unit sphere quantizes the colatitudinal and longitudinal “frequencies” $0 \le l \le \infty$ and $-l \le m \le l$. We use a sans serif $\mathbf{f}$ to denote the ordered column vector of spherical harmonic coefficients: $\mathbf{f} = (\cdots f_{lm} \cdots)^{\mathrm{T}}$. The norms of a spatial-domain function $f(\hat{\mathbf{r}})$ and its spectral-domain equivalent $\mathbf{f}$ will be denoted by

$$ (3.12) \qquad \|f\|_{\Omega}^{2} = \int_{\Omega} f^{2} d\Omega, \qquad \|f\|_{\infty}^{2} = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} f_{lm}^{2}. $$

Using this notation, we can write Parseval’s relation in the form $\|f\|_{\Omega}^{2} = \|f\|_{\infty}^{2}$. The power spectral density or variance per spherical harmonic degree $l$ and per unit area of a function $f(\hat{\mathbf{r}})$ is defined by $\langle f_l^2 \rangle = (2l+1)^{-1} \sum_{m=-l}^{l} f_{lm}^2$. We use $\delta(\hat{\mathbf{r}}, \hat{\mathbf{r}}')$ for the Dirac delta function on the sphere, with the sifting property

$$ (3.13) \qquad \int_{\Omega} \delta(\hat{\mathbf{r}}, \hat{\mathbf{r}}') f(\hat{\mathbf{r}}) d\Omega = f(\hat{\mathbf{r}}'). $$

The spherical harmonic representation of $\delta(\hat{\mathbf{r}}, \hat{\mathbf{r}}') = (\sin \theta)^{-1} \delta(\theta - \theta') \delta(\phi - \phi')$ is

$$ (3.14) \qquad \delta(\hat{\mathbf{r}}, \hat{\mathbf{r}}') = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} Y_{lm}(\hat{\mathbf{r}})Y_{lm}(\hat{\mathbf{r}}') = \sum_{l=0}^{\infty} \left(\frac{2l+1}{4\pi}\right) P_l(\cos \Delta). $$

The power spectral density of $\delta(\hat{\mathbf{r}}, \hat{\mathbf{r}}')$ is white: $\langle \delta_l^2 \rangle = 1/(4\pi)$ for all $0 \le l \le \infty$.
**3.3. Bandlimited and Spacelimited Functions.** We are interested in two subspaces of the space of all square-integrable functions on the unit sphere $\Omega$. We use $S_L = \{g: \langle g_l^2 \rangle = 0 \text{ for } L < l \le \infty\}$ to denote the space of all bandlimited functions

$$ (3.15) \qquad g = \sum_{l=0}^{L} \sum_{m=-l}^{l} g_{lm} Y_{lm}, $$

with no power above the bandwidth $L$, and $S_R = \{h: h = 0 \text{ in } \Omega - R\}$ to denote the space of all spacelimited functions $h(\hat{\mathbf{r}})$ that are strictly contained within a region $R$. The space $S_R$ is infinite-dimensional but $\dim S_L = (L+1)^2$, since the vector of spherical harmonic coefficients $\mathbf{g} = (g_{00} \cdots g_{LL})^\top$ associated with a function $g(\hat{\mathbf{r}})$ of the form (3.15) has $\sum_{l=0}^L (2l+1) = (L+1)^2$ entries. Spatial and spectral (semi)norms analogous to (3.12) are defined as

$$ (3.16) \qquad \|f\|_R^2 = \int_R f^2 d\Omega, \qquad \|f\|_L^2 = \sum_{l=0}^L \sum_{m=-l}^{l} f_{lm}^2. $$

**4. Concentration within an Arbitrarily Shaped Region.** No function can be strictly spacelimited as well as strictly bandlimited, i.e., no $f(\hat{\mathbf{r}})$ can be in both subspaces $S_R$ and $S_L$ simultaneously. The objective of this paper is to determine those bandlimited functions $g(\hat{\mathbf{r}}) \in S_L$ that are optimally concentrated within a spatial region $R$ and those spacelimited functions $h(\hat{\mathbf{r}}) \in S_R$ whose spectrum is optimally concentrated within an interval $0 \le l \le L$. As in the time-frequency case, these two spatiospectral concentration problems will be shown to be each other's duals.

**4.1. Spatial Concentration of a Bandlimited Function.** To maximize the spatial concentration of a bandlimited function $g(\hat{\mathbf{r}}) \in S_L$ within a region $R$, we maximize the ratio of the (semi)norms:

$$ (4.1) \qquad \lambda = \frac{\|g\|_R^2}{\|g\|_{\Omega}^2} = \frac{\int_R g^2 d\Omega}{\int_{\Omega} g^2 d\Omega} = \text{maximum.} $$

The two-dimensional variational problem (4.1) is analogous to the one-dimensional problem (2.2). Here, as there, the ratio $0 < \lambda < 1$ is a measure of the spatial concentration. Upon inserting the representation (3.15) of $g(\hat{\mathbf{r}})$ into (4.1) and interchanging the order of summation and integration, we can express $\lambda$ in the form

$$ (4.2) \qquad \lambda = \frac{\sum_{l=0}^{L} \sum_{m=-l}^{l} g_{lm} \sum_{l'=0}^{L} \sum_{m'=-l'}^{l'} D_{lm, l'm'} g_{l'm'}}{\sum_{l=0}^{L} \sum_{m=-l}^{l} g_{lm}^2}, $$

where we have used the orthonormality relation (3.4) and defined the quantities

$$ (4.3) \qquad D_{lm,l'm'} = \int_R Y_{lm} Y_{l'm'} d\Omega. $$
Upon introducing the $(L+1)^2 \times (L+1)^2$ matrix

$$ (4.4) \qquad \mathbf{D} = \begin{pmatrix} D_{00,00} & \cdots & D_{00,LL} \\ \vdots & \ddots & \vdots \\ D_{LL,00} & \cdots & D_{LL,LL} \end{pmatrix}, $$

with elements $D_{lm,l'm'}$, $0 \le l \le L$ and $-l \le m \le l$, we can rewrite (4.1) as a classical matrix variational problem in the space of $(L+1)^2$-tuples [25]:

$$ (4.5) \qquad \lambda = \frac{\mathbf{g}^\mathrm{T} \mathbf{D} \mathbf{g}}{\mathbf{g}^\mathrm{T} \mathbf{g}} = \text{maximum.} $$

Column vectors $\mathbf{g}$ that render the Rayleigh quotient $\lambda$ in (4.5) stationary are solutions of the $(L+1)^2 \times (L+1)^2$ algebraic eigenvalue problem

$$ (4.6) \qquad \mathbf{D}\mathbf{g} = \lambda\mathbf{g}. $$

Equation (4.6) is the discrete spherical analogue of the one-dimensional frequency-domain equation (2.3). The matrix (4.4) is real, symmetric ($\mathbf{D}^\mathrm{T} = \mathbf{D}$), and positive definite ($\mathbf{g}^\mathrm{T}\mathbf{D}\mathbf{g} > 0$ for all $\mathbf{g} \neq \mathbf{0}$), so the $(L+1)^2$ eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_{(L+1)^2}$ and eigenvectors $\mathbf{g}$ are always real [25]. We order the eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_{(L+1)^2}$ and eigenvectors $\mathbf{g}_1, \mathbf{g}_2, \dots, \mathbf{g}_{(L+1)^2}$ so that $1 > \lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_{(L+1)^2} > 0$. Every spectral-domain eigenvector $\mathbf{g}_\alpha$ gives rise to an associated bandlimited spatial eigenfunction $g_\alpha(\hat{\mathbf{r}})$ defined by (3.15). The largest eigenvalue $\lambda_1$ is strictly smaller than one because no bandlimited function can be completely confined within the region $R$, and the smallest eigenvalue $\lambda_{(L+1)^2}$ is strictly greater than zero because of the positive definiteness of the matrix $\mathbf{D}$.

The symmetry $\mathbf{D}^\mathrm{T} = \mathbf{D}$ also guarantees that the eigenvectors $\mathbf{g}_1, \mathbf{g}_2, \dots, \mathbf{g}_{(L+1)^2}$ are mutually orthogonal [25]. We choose them to be orthonormal:

$$ (4.7) \qquad \mathbf{g}_\alpha^\mathrm{T} \mathbf{g}_\beta = \delta_{\alpha\beta}, \qquad \mathbf{g}_\alpha^\mathrm{T} \mathbf{D} \mathbf{g}_\beta = \lambda_\alpha \delta_{\alpha\beta}. $$

The associated spatial eigenfunctions $g_1(\hat{\mathbf{r}}), g_2(\hat{\mathbf{r}}), \dots, g_{(L+1)^2}(\hat{\mathbf{r}})$ form a basis for $S_L$. They are orthonormal over the whole sphere $\Omega$ and orthogonal over the region $R$:

$$ (4.8) \qquad \int_{\Omega} g_{\alpha} g_{\beta} \, d\Omega = \delta_{\alpha\beta}, \qquad \int_{R} g_{\alpha} g_{\beta} \, d\Omega = \lambda_{\alpha} \delta_{\alpha\beta}. $$

The two spatial-domain relations in (4.8) are equivalent to the corresponding matrix spectral relations (4.7), and they are analogous to the one-dimensional orthogonality relations (2.9). The eigenfunction $g_1(\hat{\mathbf{r}})$ associated with the largest eigenvalue $\lambda_1$ is the member of the space $S_L$ of bandlimited functions that is most spatially concentrated within region $R$; the eigenfunction $g_2(\hat{\mathbf{r}})$ is the next best concentrated function in $S_L$ orthogonal to $g_1(\hat{\mathbf{r}})$ over both $\Omega$ and $R$; and so on.

Written out in full using index notation, the matrix eigenvalue equation (4.6) is

$$ (4.9) \qquad \sum_{l'=0}^{L} \sum_{m'=-l'}^{l'} D_{lm,l'm'} g_{l'm'} = \lambda g_{lm}. $$

Upon multiplying (4.9) by $Y_{lm}(\hat{\mathbf{r}})$ and summing over all $0 \le l \le L$ and $-l \le m \le l$, we obtain a spatial-domain eigenvalue equation:

$$ (4.10) \qquad \int_{R} D(\hat{\mathbf{r}}, \hat{\mathbf{r}}') g(\hat{\mathbf{r}}') \, d\Omega' = \lambda g(\hat{\mathbf{r}}), \quad \hat{\mathbf{r}} \in \Omega, $$
where we have defined the bandlimited Dirac delta function

$$ (4.11) \qquad D(\hat{\mathbf{r}}, \hat{\mathbf{r}}') = \sum_{l=0}^{L} \sum_{m=-l}^{l} Y_{lm}(\hat{\mathbf{r}}) Y_{lm}(\hat{\mathbf{r}}') = \sum_{l=0}^{L} \left( \frac{2l+1}{4\pi} \right) P_l(\cos \Delta). $$

Equation (4.10) is a homogeneous Fredholm integral equation of the second kind, with a finite-rank, symmetric, separable kernel [31, 72]. Upon inserting the representations (3.15) and (4.11) into (4.10), we recover the matrix equation (4.6), so that the spectral-domain eigenvalue problem for $\mathbf{g}$ and the spatial-domain eigenvalue problem for a bandlimited $g(\hat{\mathbf{r}}) \in S_L$ are completely equivalent.

In summary, we can find an orthogonal family of bandlimited eigenfunctions that are optimally concentrated within a region $R$ on the unit sphere $\Omega$ by solving either the $(L+1)^2 \times (L+1)^2$ matrix eigenvalue problem (4.6) for the spectral-domain eigenvectors $\mathbf{g}_1, \mathbf{g}_2, \dots, \mathbf{g}_{(L+1)^2}$ or the Fredholm integral equation (4.10) for the associated spatial-domain eigenfunctions $g_1, g_2, \dots, g_{(L+1)^2}$. Either method determines the optimally concentrated eigenfunctions at all points $\hat{\mathbf{r}} \in \Omega$, i.e., both in the region $R$, where they are concentrated, and in the complementary region $\Omega - R$, where they exhibit inevitable leakage.

**4.2. Spectral Concentration of a Spacelimited Function.** Instead of seeking to concentrate a bandlimited function $g(\hat{\mathbf{r}}) \in S_L$ within a spatial region $R$, we may seek to concentrate a spacelimited function $h(\hat{\mathbf{r}}) \in S_R$ within a spectral interval $0 \le l \le L$. A suitable measure of concentration is then the spectral (semi)norm ratio, analogous to the one-dimensional ratio (2.4):

$$ (4.12) \qquad \lambda = \frac{\|\mathbf{h}\|_L^2}{\|\mathbf{h}\|_\infty^2} = \frac{\sum_{l=0}^{L} \sum_{m=-l}^{l} h_{lm}^2}{\sum_{l=0}^{\infty} \sum_{m=-l}^{l} h_{lm}^2} = \text{maximum.} $$

Upon inserting the representation of the spherical harmonic expansion coefficients,

$$ (4.13) \qquad h_{lm} = \int_R h Y_{lm} d\Omega, $$

and interchanging the order of summation and integration, we can rewrite the variational problem (4.12) in the form

$$ (4.14) \qquad \lambda = \frac{\int_{R} \int_{R} h(\hat{\mathbf{r}}) D(\hat{\mathbf{r}}, \hat{\mathbf{r}}') h(\hat{\mathbf{r}}') d\Omega d\Omega'}{\int_{R} h^2(\hat{\mathbf{r}}) d\Omega} = \text{maximum,} $$

where we have made use of the replication property (3.13) of the delta function (3.14) and the definition (4.11) of the kernel $D(\hat{\mathbf{r}}, \hat{\mathbf{r}}')$. Functions $h(\hat{\mathbf{r}}) \in S_R$ that render the Rayleigh quotient (4.14) stationary are solutions of the Fredholm integral equation

$$ (4.15) \qquad \int_{R} D(\hat{\mathbf{r}}, \hat{\mathbf{r}}') h(\hat{\mathbf{r}}') d\Omega' = \lambda h(\hat{\mathbf{r}}), \quad \hat{\mathbf{r}} \in R. $$

Equation (4.15) is the spherical analogue of the one-dimensional time-domain eigenvalue equation (2.5). In fact, this equation for $h(\hat{\mathbf{r}}) \in S_R$ is identical to (4.10) for
$g(\hat{\mathbf{r}}) \in S_L$. The only difference is that (4.10) is applicable on the entire sphere $\Omega$, whereas the domain of (4.15) is limited to the region $R$, within which $h(\hat{\mathbf{r}}) \neq 0$. Evidently, the eigenfunctions $h(\hat{\mathbf{r}})$ that maximize the spectral norm ratio (4.12) are identical, within the region $R$, to the eigenfunctions $g(\hat{\mathbf{r}})$ that maximize the spatial norm ratio (4.1). We normalize such that

$$ (4.16) \qquad h(\hat{\mathbf{r}}) = \begin{cases} g(\hat{\mathbf{r}}) & \text{if } \hat{\mathbf{r}} \in R, \\ 0 & \text{otherwise.} \end{cases} $$

Every bandlimited eigenfunction $g_\alpha \in S_L$ gives rise to a spacelimited eigenfunction $h_\alpha \in S_R$ defined by the restriction (4.16). The associated eigenvalues $\lambda_\alpha$ measure the spatiospectral concentration; the fractional spatial energy $1 - \lambda_\alpha$ leaked by $g_\alpha$ to the region $\Omega - R$ is identical to the fractional spectral energy leaked by $h_\alpha$ into the degrees $L < l \le \infty$. Had we started with the variational prescription (4.12) rather than (4.1), we could have obtained the integral equation (4.10) governing a bandlimited $g(\hat{\mathbf{r}}) \in S_L$ by extending the domain of (4.15) to the whole sphere $\Omega$.

The spacelimited eigenfunctions $h_1(\hat{\mathbf{r}})$, $h_2(\hat{\mathbf{r}})$, ..., $h_{(L+1)^2}(\hat{\mathbf{r}})$ defined by (4.16) are orthogonal (but not orthonormal) over both the whole sphere $\Omega$ and the region $R$:

$$ (4.17) \qquad \int_{\Omega} h_{\alpha} h_{\beta} d\Omega = \int_{R} h_{\alpha} h_{\beta} d\Omega = \lambda_{\alpha} \delta_{\alpha\beta}. $$

The relation $h_{lm} = \sum_{l'=0}^{L} \sum_{m'=-l'}^{l'} D_{lm,l'm'} g_{l'm'}$ expresses the coefficients $h_{lm}$, where $0 \le l \le \infty$, in terms of the coefficients $g_{lm}$, with $0 \le l \le L$. By (4.9), this amounts to $h_{lm} = \lambda g_{lm}$ when $0 \le l \le L$. In addition to the $(L+1)^2$ eigenfunctions with nonzero eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_{(L+1)^2}$, (4.15) has an infinite-dimensional null space of eigenfunctions with associated eigenvalue $\lambda = 0$. Every function $h(\hat{\mathbf{r}})$ that vanishes in $\Omega - R$ and has no power in the interval $0 \le l \le L$ is a member of this null space.

**4.3. Significant and Insignificant Eigenvalues.** The sum of the eigenvalues of the matrix $D$ defined in (4.4) is given by

$$ (4.18) \qquad N = \sum_{\alpha=1}^{(L+1)^2} \lambda_{\alpha} = \text{tr} D = \sum_{l=0}^{L} \sum_{m=-l}^{l} D_{lm,lm} = \int_{R} D(\hat{\mathbf{r}}, \hat{\mathbf{r}}) d\Omega = (L+1)^2 \frac{A}{4\pi}, $$

where we have substituted for the diagonal matrix elements $D_{lm,lm}$ from (4.3) and used (4.11) and the identity $P_l(1) = 1$. The quantity $N$ in (4.18) is the spherical analogue of the Shannon number (2.7) in Slepian's one-dimensional concentration problem. Eigenfunctions $g_\alpha(\hat{\mathbf{r}})$ that are well concentrated within the region $R$ will have eigenvalues $\lambda_\alpha$ near unity, whereas those that are poorly concentrated will have eigenvalues $\lambda_\alpha$ near zero. If, as in the one-dimensional problem, the spectrum of eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_{(L+1)^2}$ has a relatively narrow transition band from values near unity to values near zero [60], then the total number of significant ($\lambda_\alpha \approx 1$) eigenvalues will be well approximated by the rounded sum (4.18). For this reason, we expect $N$ to be a good estimate of the number of significant eigenvalues. Roughly speaking, the spherical Shannon number (4.18) is the dimension of the space of two-dimensional functions $f(\hat{\mathbf{r}})$ that are both approximately limited in the spectral domain to spherical harmonic degrees $0 \le l \le L$ and approximately limited in the spatial domain to an arbitrarily shaped region $R$ of area $A$ [33, 34].

Rather than seeking a bandlimited function $g(\hat{\mathbf{r}}) \in S_L$ that is optimally concentrated within a spatial region $R$, we could have decided to seek one that is optimally
excluded from $R$, i.e., one that is optimally concentrated within the complementary region $\Omega - R$. In that case, we would have sought to minimize rather than maximize the Rayleigh quotient (4.1). In fact, all that we have found are the bandlimited functions $g(\hat{\mathbf{r}}) \in S_L$ that render the eigenvalue $\lambda$ stationary, so we have actually solved the concentration and exclusion problems simultaneously. The optimally concentrated eigenfunctions and the optimally excluded eigenfunctions are identical, but with the ordering indices reversed, i.e., the bandlimited function that is most excluded from $R$ is $g_{(L+1)^2}$, the next most excluded is $g_{(L+1)^2-1}$, and so on. Since $\lambda_\alpha$ is the fractional power of the stationary eigenfunction $g_\alpha$ within $R$, its fractional power within $\Omega - R$ is $1 - \lambda_\alpha$. Whenever the area $A$ of the region $R$ is a small fraction of the area of the sphere, $A \ll 4\pi$, there will be many more well-excluded eigenfunctions with insignificant ($\lambda_\alpha \approx 0$) eigenvalues than well-concentrated eigenfunctions with significant ($\lambda_\alpha \approx 1$) eigenvalues, i.e., $N \ll (L+1)^2$.

The kernel $D(\hat{\mathbf{r}}, \hat{\mathbf{r}}')$ in the integral eigenvalue equation (4.10) can be expressed in terms of the spatial-domain eigenfunctions $g_1, g_2, \dots, g_{(L+1)^2}$ in the form

$$ (4.19) \qquad D(\hat{\mathbf{r}}, \hat{\mathbf{r}}') = \sum_{\alpha=1}^{(L+1)^2} g_{\alpha}(\hat{\mathbf{r}})g_{\alpha}(\hat{\mathbf{r}}'). $$

To verify that (4.19) is equivalent to the original representation (4.11), it suffices to note that both $Y_{lm}, 0 \le l \le L, -l \le m \le l$, and $g_\alpha, \alpha = 1, 2, \dots, (L+1)^2$, are $(L+1)^2$-dimensional orthonormal bases for $S_L$. The transformed representation (4.19) is the spherical version of Mercer's theorem [17, 31, 72]. Upon setting $\hat{\mathbf{r}}' = \hat{\mathbf{r}}$ in (4.19), we deduce that the sum of the squares of the $(L+1)^2$ bandlimited eigenfunctions $g_\alpha(\hat{\mathbf{r}})$ is a constant, independent of position $\hat{\mathbf{r}}$ on the sphere $\Omega$; in fact,

$$ (4.20) \qquad \sum_{\alpha=1}^{(L+1)^2} g_{\alpha}^2(\hat{\mathbf{r}}) = D(\hat{\mathbf{r}}, \hat{\mathbf{r}}) = \frac{(L+1)^2}{4\pi} = \frac{N}{A}. $$

If the first $N$ eigenfunctions $g_1, g_2, \dots, g_N$ have eigenvalues near unity and lie mostly within $R$, and the remainder $g_{N+1}, g_{N+2}, \dots, g_{(L+1)^2}$ have eigenvalues near zero and lie mostly in $\Omega - R$, then we expect the eigenvalue-weighted sum of squares to be

$$ (4.21) \qquad \sum_{\alpha=1}^{(L+1)^2} \lambda_{\alpha} g_{\alpha}^2(\hat{\mathbf{r}}) \approx \sum_{\alpha=1}^{N} \lambda_{\alpha} g_{\alpha}^2(\hat{\mathbf{r}}) \approx \begin{cases} N/A & \text{if } \hat{\mathbf{r}} \in R, \\ 0 & \text{otherwise.} \end{cases} $$

The terms with $N+1 \le \alpha \le (L+1)^2$ should be negligible, so it is immaterial whether they are included in the sum (4.21) or not. Taken together, the first $N$ orthogonal eigenfunctions $g_\alpha, \alpha = 1, 2, \dots, N$, with significant eigenvalues $\lambda_\alpha \approx 1$, provide an essentially uniform coverage of the region $R$. This is really the essence of the spatiospectral concentration problem: the number of degrees of freedom is reduced from $\dim S_L = (L+1)^2$ to the Shannon number $N = (L+1)^2 A/(4\pi)$.

**4.4. Abstract Operator Formulation.** We conclude this section on the concentration problem for an arbitrarily shaped region by reiterating the above results using an abstract operator notation. We use $\mathcal{H}$ to denote the operator that acts upon square-integrable functions $f(\hat{\mathbf{r}})$ in the spatial domain to produce the associated infinite-dimensional column vectors $\mathbf{f}$ of spherical harmonic coefficients $f_{lm}$ in the spectral domain, and we use $\mathcal{H}^{-1}$ to denote its inverse, so that $\mathcal{H}\mathbf{f} = \mathbf{f}$ and
$\mathcal{H}^{-1}\mathbf{f} = \mathbf{f}$. We introduce two operators, $\mathcal{R}$ and $\mathcal{L}$, which project onto the space $S_R$ of spacelimited functions and the space $S_L$ of bandlimited functions, respectively. The first of these acts to spatially restrict functions $f(\hat{\mathbf{r}})$ in the spatial domain,

$$ (4.22) \qquad \mathcal{R}f(\hat{\mathbf{r}}) = \begin{cases} f(\hat{\mathbf{r}}) & \text{if } \hat{\mathbf{r}} \in \mathcal{R}, \\ 0 & \text{otherwise,} \end{cases} $$

whereas the second acts to bandlimit column vectors $\mathbf{f}$ in the spectral domain,

$$ (4.23) \qquad \mathcal{L}\mathbf{f} = \mathcal{L} \begin{pmatrix} f_{00} \\ \vdots \\ f_{\infty\infty} \end{pmatrix} = \begin{pmatrix} f_{00} \\ \vdots \\ f_{LL} \end{pmatrix}. $$

The operator product $\mathcal{H}^{-1}\mathcal{L}\mathcal{H}$ acts to bandlimit an arbitrary function $f(\hat{\mathbf{r}})$. In numerical analysis [4, 28, 46, 62], this operation is referred to as spherical or uniform-resolution filtering or as triangular truncation. We denote the inner products in the two domains by $\langle f, f' \rangle_{\Omega} = \int_{\Omega} f f' d\Omega$ and $\langle \mathbf{f}, \mathbf{f}' \rangle_{\infty} = \mathbf{f}^{\top}\mathbf{f}'$. The spatial and spectral norms introduced in (3.12) are given by $\|f\|_{\Omega}^2 = \langle f, f \rangle_{\Omega}$ and $\|f\|_{\infty}^2 = \langle \mathbf{f}, \mathbf{f} \rangle_{\infty}$.

The spatial concentration variational problem (4.1) and the spectral concentration variational problem (4.12) can be written using this operator notation in the form

$$ (4.24) \qquad \lambda = \frac{\langle \mathcal{R}\mathcal{H}^{-1}\mathcal{L}\mathbf{f}, \mathcal{R}\mathcal{H}^{-1}\mathcal{L}\mathbf{f} \rangle_{\Omega}}{\langle \mathcal{H}^{-1}\mathcal{L}\mathbf{f}, \mathcal{H}^{-1}\mathcal{L}\mathbf{f} \rangle_{\Omega}} = \frac{\langle \mathcal{L}\mathcal{H}\mathcal{R}\mathbf{f}, \mathcal{L}\mathcal{H}\mathcal{R}\mathbf{f} \rangle_{\infty}}{\langle \mathcal{H}\mathcal{R}\mathbf{f}, \mathcal{H}\mathcal{R}\mathbf{f} \rangle_{\infty}} = \text{maximum.} $$

The associated spectral-domain and spatial-domain eigenvalue equations are

$$ (4.25) \qquad (\mathcal{L}\mathcal{H}\mathcal{R}\mathcal{H}^{-1}\mathcal{L})(\mathcal{L}\mathbf{f}) = \lambda(\mathcal{L}\mathbf{f}), \qquad (\mathcal{R}\mathcal{H}^{-1}\mathcal{L}\mathcal{H}\mathcal{R})(\mathcal{R}\mathbf{f}) = \lambda(\mathcal{R}\mathbf{f}), $$

where we have made use of the facts that $\mathcal{H}$ and $\mathcal{H}^{-1}$ are each others' adjoints, that both $\mathcal{R}$ and $\mathcal{L}$ are self-adjoint, and that $\mathcal{R}^2 = \mathcal{R}$ and $\mathcal{L}^2 = \mathcal{L}$. The equations (4.25) are the operator equivalents of the algebraic eigenvalue equation (4.6) and the integral eigenvalue equation (4.15). Any solution of (4.24) is a bandlimited column vector of the form $\mathbf{g} = \mathcal{L}\mathbf{f}$, whereas any solution of (4.25) is a spacelimited function of the form $\mathbf{h} = \mathcal{R}\mathbf{f}$. Both the spectral-domain operator $\mathcal{L}\mathcal{H}\mathcal{R}\mathcal{H}^{-1}\mathcal{L}$ and the spatial-domain operator $\mathcal{R}\mathcal{H}^{-1}\mathcal{L}\mathcal{H}\mathcal{R}$ are symmetric by inspection.

**5. Concentration within an Axisymmetric Polar Cap.** We turn our attention next to the special but important case in which the region of concentration is a circularly symmetric cap of colatitudinal radius $\Theta$, centered on the north pole, as illustrated in the lower right of Figure 3.1. In practical applications, the eigenfunctions that are optimally concentrated within such a polar cap $R = \{\theta: 0 \le \theta \le \Theta\}$ can be rotated to an arbitrarily positioned circular cap on the unit sphere using standard spherical harmonic rotation formulae [3, 10, 15].

**5.1. Decomposition of the Spectral-Domain Eigenvalue Problem.** The matrix elements (4.3) in this axisymmetric case reduce to

$$ (5.1) \qquad D_{lm,l'm'} = 2\pi \delta_{mm'} \int_{0}^{\Theta} X_{lm}X_{l'm} \sin\theta \, d\theta. $$

The Kronecker delta $\delta_{mm'}$ renders the $(L+1)^2 \times (L+1)^2$ matrix $\mathbf{D}$ of (4.4) block diagonal: $\mathbf{D} = \text{diag}(\mathbf{D}_0, \mathbf{D}_1, \mathbf{D}_1, \dots, \mathbf{D}_L, \mathbf{D}_L)$, where every submatrix $\mathbf{D}_m \neq \mathbf{D}_0$ occurs
twice as a result of the doublet degeneracy associated with $\pm m$. Rather than solving the full eigenvalue equation (4.6), we solve a series of $(L-m+1) \times (L-m+1)$ spectral-domain algebraic eigenvalue problems, $\mathbf{D}_m\mathbf{g}_m = \lambda_m\mathbf{g}_m$, one for each non-negative order $m$. In the remainder of this section, we drop the identifying subscript $m$ and write each fixed-order eigenvalue problem as, simply,

$$ (5.2) \qquad \mathbf{D}\mathbf{g} = \lambda\mathbf{g}. $$

The matrix $\mathbf{D}$ and the column vector $\mathbf{g}$ in (5.2) are of the form

$$ (5.3) \qquad \mathbf{D} = \begin{pmatrix} D_{mm} & \cdots & D_{mL} \\ \vdots & \ddots & \vdots \\ D_{Lm} & \cdots & D_{LL} \end{pmatrix}, \qquad \mathbf{g} = \begin{pmatrix} g_m \\ \vdots \\ g_L \end{pmatrix}, $$

where, for a particular order $0 \le m \le L$ and $m \le l \le L$,

$$ (5.4) \qquad D_{ll'} = 2\pi \int_0^\Theta X_{lm} X_{l'm} \sin \theta \, d\theta. $$

A simple analytical expression for this integral exists when $m=0$ and $l \neq l'$ [76]; more generally, (5.4) can be evaluated with the aid of (3.6)-(3.7):

$$ (5.5) \qquad D_{ll'} = (-1)^m \frac{\sqrt{(2l+1)(2l'+1)}}{2} \sum_{n=|l-l'|}^{l+l'} \begin{pmatrix} l & n & l' \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} l & n & l' \\ m & 0 & -m \end{pmatrix} \times [P_{n-1}(\cos \Theta) - P_{n+1}(\cos \Theta)]. $$

We rank order the distinct $L-m+1$ eigenvalues obtained by solving each of the fixed-order eigenvalue problems (5.2) so that $1 > \lambda_1 > \lambda_2 > \cdots > \lambda_{L-m+1} > 0$, and we orthonormalize the associated eigenvectors $\mathbf{g}_1, \mathbf{g}_2, \dots, \mathbf{g}_{L-m+1}$ as in (4.7) so that

$$ (5.6) \qquad \mathbf{g}_{\alpha}^{\top} \mathbf{g}_{\beta} = \delta_{\alpha\beta}, \qquad \mathbf{g}_{\alpha}^{\top} \mathbf{D} \mathbf{g}_{\beta} = \lambda_{\alpha} \delta_{\alpha\beta}. $$

The associated bandlimited eigenfunctions $g_1(\theta), g_2(\theta), \dots, g_{L-m+1}(\theta)$, defined by

$$ (5.7) \qquad g = \sum_{l=m}^{L} g_l X_{lm}, $$

then satisfy the colatitudinal orthogonality relations

$$ (5.8) \qquad 2\pi \int_{0}^{\pi} g_{\alpha} g_{\beta} \sin \theta \, d\theta = \delta_{\alpha\beta}, \qquad 2\pi \int_{0}^{\Theta} g_{\alpha} g_{\beta} \sin \theta \, d\theta = \lambda_{\alpha} \delta_{\alpha\beta}. $$

The optimally concentrated spatial eigenfunctions $g(\hat{\mathbf{r}})$ for a given $-L \le m \le L$ are expressed in terms of the fixed-order colatitudinal eigenfunctions (5.7) by

$$ (5.9) \qquad g(\theta, \phi) = \begin{cases} \sqrt{2} g(\theta) \cos m\phi & \text{if } -L \le m < 0, \\ g(\theta) & \text{if } m = 0, \\ \sqrt{2} g(\theta) \sin m\phi & \text{if } 0 < m \le L. \end{cases} $$
<figure><img src="image_3.png" /><figcaption>**Fig. 5.1** Colatitudinal dependence of the first four bandlimited eigenfunctions $g_\alpha(\theta)$, $\alpha = 1, 2, 3, 4$, of fixed order $m = 0$ (top) to $m = 4$ (bottom). The radius of the polar cap is $\Theta = 40^\circ$, and the bandwidth is $L = 18$. Black curves show the concentration within the cap $0^\circ \le \theta \le 40^\circ$; gray curves show the leakage into the rest of the sphere, $40^\circ < \theta \le 180^\circ$. Labels show the eigenvalues $\lambda_\alpha$, which express the quality of the spatial concentration. The corresponding spacelimited spectral-domain eigenfunctions are shown in Figure 5.2.</figcaption></figure>

The four most optimally concentrated eigenfunctions $g_1(\theta)$, $g_2(\theta)$, $g_3(\theta)$, $g_4(\theta)$ for orders $0 \le m \le 4$ are plotted in Figure 5.1. The associated eigenvalues $\lambda_1, \lambda_2, \lambda_3, \lambda_4$ are listed to six-figure accuracy. The radius of the polar cap in this example is $\Theta = 40^\circ$, the bandwidth is $L = 18$, and the rounded Shannon number is $N = 42$. The first zonal ($m=0$) eigenfunction, $g_1(\theta)$, has no nodes within the cap $0^\circ \le \theta \le 40^\circ$; the second, $g_2(\theta)$, has one node; and so on. The nonzonal ($m > 0$) eigenfunctions all vanish at the north pole, $\theta = 0^\circ$. The first four zonal eigenfunctions, the first three $m = 1$ and $m = 2$ eigenfunctions, and the first two $m = 3$ and $m = 4$ eigenfunctions are all very well concentrated ($\lambda > 0.9$), whereas the fourth $m = 3$ and $m = 4$ eigenfunctions exhibit significant leakage ($\lambda < 0.1$).

**5.2. Decomposition of the Spatial-Domain Eigenvalue Problem.** The integral eigenvalue problem (4.10) in the spatial domain likewise decomposes into a series of
fixed-order, one-dimensional Fredholm eigenvalue equations,

$$ (5.10) \qquad \int_0^\Theta D(\theta, \theta') g(\theta') \sin \theta' d\theta' = \lambda g(\theta), \quad 0 \le \theta \le \pi, $$

each with an m-dependent, separable, symmetric kernel

$$ (5.11) \qquad D(\theta, \theta') = 2\pi \sum_{l=m}^{L} X_{lm}(\theta) X_{lm}(\theta'). $$

The results (5.10)–(5.11) can be obtained either by multiplying $\sum_{l'=m}^L D_{ll'}g_{l'} = \lambda g_l$, the index form of (5.2), by $X_{lm}(\theta)$ and summing over $m \le l \le L$, or by substituting the representation (5.9) into (4.10) and using the orthogonality of the longitudinal functions ..., $\sqrt{2}\cos m\phi$, ..., $1$, ..., $\sqrt{2}\sin m\phi$, ... over $0 \le \phi < 2\pi$. The matrix eigenvalue problem (5.2) can in turn be derived starting from the separable Fredholm equation (5.10), so that the fixed-order spectral and spatial eigenvalue problems are completely equivalent. The fixed-order spacelimited eigenfunctions

$$ (5.12) \qquad h(\theta) = \begin{cases} g(\theta) & \text{if } 0 \le \theta \le \Theta, \\ 0 & \text{otherwise} \end{cases} $$

satisfy an equation identical to (5.10), but only within the polar cap itself:

$$ (5.13) \qquad \int_0^\Theta D(\theta, \theta') h(\theta') \sin \theta' d\theta' = \lambda h(\theta), \quad 0 \le \theta \le \Theta. $$

The eigenvalue $\lambda$ is a measure of both the spatial concentration of $g(\theta) \in S_L$ within the cap $0 \le \theta \le \Theta$ and the spectral concentration of $h(\theta) \in S_R$ within the interval $0 \le l \le L$. The substitution $\mu = \cos \theta$ converts (5.10) and (5.13) into

$$ (5.14a) \qquad \int_{\cos \Theta}^1 D(\mu, \mu') g(\mu') d\mu' = \lambda g(\mu), \quad -1 \le \mu \le 1, $$

$$ (5.14b) \qquad \int_{\cos \Theta}^1 D(\mu, \mu') h(\mu') d\mu' = \lambda h(\mu), \quad \cos \Theta \le \mu \le 1. $$

The kernel $D(\mu, \mu')$ can be simplified using the Christoffel–Darboux identity (3.10):

$$ \begin{align} D(\mu, \mu') &= \frac{1}{2} \sum_{l=m}^{L} (2l+1) \frac{(l-m)!}{(l+m)!} P_{lm}(\mu) P_{lm}(\mu') \nonumber \\ (5.15) \qquad \qquad \qquad \qquad &= \frac{(L-m+1)!}{2(L+m)!} \left[ \frac{P_{L+1,m}(\mu) P_{L,m}(\mu') - P_{L,m}(\mu) P_{L+1,m}(\mu')}{\mu - \mu'} \right], \end{align} $$

with L'Hôpital's rule for $\mu = \mu'$. Equation (5.15) generalizes the $m = 0$ result of [20].

The squared spherical harmonic coefficients $h_l^2$ of the four best concentrated space-limited eigenfunctions $h_1(\theta), h_2(\theta), h_3(\theta), h_4(\theta)$ for $0 \le m \le 4$ are plotted versus spherical harmonic degree $l$ in Figure 5.2. The cap radius $\Theta = 40^\circ$, spectral concentration interval $L = 18$, and layout are the same as in Figure 5.1. The maximum contribution to the $\alpha$th zonal ($m = 0$) eigenfunction comes from the degree $l$ satisfying $2\pi\alpha/\sqrt{l(l+1)} \approx 2\Theta$; physically, this corresponds to fitting an integral number of asymptotic wavelengths within the cap of diameter $2\Theta$.
<figure><img src="image_4.png" /><figcaption>**Fig. 5.2** Squared spherical harmonic coefficients $h_l^2$ of the first four spacelimited eigenfunctions $h_\alpha(\theta)$, $\alpha = 1, 2, 3, 4$, of fixed order $m = 0$ (top) to $m = 4$ (bottom). The radius of the polar cap is $\Theta = 40^\circ$, and the spectral concentration interval is $L = 18$. Black curves show the power within the interval $0 \le l \le 18$; gray curves show the power leaked to $19 \le l \le 127$. Values of $h_l^2$ are in decibels, normalized to zero at the individual maxima. Labels show the eigenvalues $\lambda_\alpha$, which express the quality of the spectral concentration. The corresponding bandlimited spatial-domain eigenfunctions are shown in Figure 5.1.</figcaption></figure>

**5.3. Eigenvalue Spectrum and Eigenfunctions.** For each of the fixed-order eigenvalue problems (5.2), (5.10), (5.13), and (5.14), the number of significant eigenvalues, or partial Shannon number, can be computed using any of the formulae

$$ (5.16) \qquad N_m = \sum_{\alpha=1}^{L-m+1} \lambda_{\alpha} = \sum_{l=m}^{L} D_{ll} = \int_{0}^{\Theta} D(\theta, \theta) d\theta. $$

We can write the final relation in (5.16) using (5.15) in the form

$$ (5.17) \qquad N_m = \frac{(L-m+1)!}{2(L+m)!} \int_{\cos \Theta}^{1} [P'_{L+1, m} P_{L, m} - P'_{L, m} P_{L+1, m}] d\mu, $$

where the prime denotes differentiation with respect to $\mu$.

Once the $L+1$ sequences of fixed-order eigenvalues have been found, they can be re-sorted to exhibit an overall mixed-order ranking. The total number of significant
<figure><img src="image_5.png" /><figcaption>**Fig. 5.3** Reordered eigenvalue spectra ($\lambda_\alpha$ versus rank $\alpha$) for axisymmetric polar caps of colatitudinal radii $\Theta = 10^\circ, 20^\circ, 30^\circ, 40^\circ$ and a common maximal spherical harmonic degree $L = 18$. The total number of eigenvalues is $(L+1)^2 = 381$; only $\lambda_1$ through $\lambda_{60}$ are shown. Different symbols are used to plot the various orders $-11 \le m \le 11$; juxtaposed identical symbols are $\pm m$ doublets. Vertical gridlines and top labels specify the rounded Shannon numbers $N = 3, 11, 24, 42$.</figcaption></figure>

eigenvalues (4.18) is then $N = N_0 + 2 \sum_{m=1}^{L} N_m$, where the factor of two accounts for the $\pm m$ degeneracy. In Figure 5.3, we show the reordered, mixed-$m$ eigenvalue spectra for four different polar caps, with colatitudinal radii $\Theta = 10^\circ, 20^\circ, 30^\circ, 40^\circ$; the maximal spherical harmonic degree is $L = 18$. The spectra have a characteristic step shape [33, 49, 60], showing significant ($\lambda \approx 1$) and insignificant ($\lambda \approx 0$) eigenvalues separated by a narrow transition band. The rounded Shannon numbers $N = 3, 11, 24, 42$ roughly separate the reasonably well concentrated eigensolutions ($\lambda \ge 0.5$) from the more poorly concentrated ones ($\lambda < 0.5$) in all four cases.

Figure 5.4 shows a polar plot of the first 32 eigenfunctions $g(\theta, \phi)$ concentrated within a cap of radius $\Theta = 40^\circ$, defined by (5.9). The maximal spherical harmonic degree is $L = 18$ and the Shannon number is $N = 42$, as in Figure 5.1. The eigenvalue ranking is mixed order, as in Figure 5.3, and all degenerate $\sqrt{2} \cos m\phi$, $\sqrt{2} \sin m\phi$ doublets are shown. The concentration factors $1 < \lambda \le 0.849$ and orders $m$ of each eigenfunction are indicated. Blue and red colors represent positive and negative values, respectively; however, all signs could be reversed without violating the quadratic concentration criteria (4.1) and (4.12).

Finally, in Figure 5.5, we illustrate the eigenvalue-weighted pointwise sums of squares $\sum_{\alpha} \lambda_{\alpha} g_{\alpha}^2(\theta, \phi)$ for polar caps of radii $\Theta = 10^\circ, 20^\circ, 30^\circ, 40^\circ$ and bandwidth $L = 18$. The sums are concentrated within the polar cap $0 \le \theta \le \Theta$; solid lines of differing shades of gray distinguish the sums carried up to the first $N$ or all $(L+1)^2$
<figure><img src="image_6.png" /><figcaption>**Fig. 5.4** Bandlimited eigenfunctions $g(\theta, \phi)$ that are optimally concentrated within a circular cap of colatitudinal radius $\Theta = 40^{\circ}$. Dashed circles denote the cap boundary. The bandwidth is $L = 18$, and the rounded Shannon number is $N = 42$. Subscripts on the eigenvalues $\lambda_\alpha$ specify the fixed-order rank. The eigenvalues have been re-sorted into a mixed-order ranking, with the best concentrated eigenfunction plotted on the top left and the 32nd best on the lower right. Blue is positive and red is negative; regions in which the absolute value is less than one hundredth of the maximum value on the sphere are left white.</figcaption></figure>

possible terms. Our heuristic expectation (4.21) is confirmed: inside the cap, the weighted sums rapidly approach $N/A = (L+1)^2/(4\pi)$. In contrast, the full unweighted sum $\sum_\alpha g_\alpha^2(\theta, \phi)$ of $(L+1)^2$ terms (dashed lines) is exactly $N/A$ over the entire sphere in accordance with (4.20).
<figure><img src="image_7.png" /><figcaption>**Fig. 5.5** Cumulative energy of the eigenfunctions concentrated within circularly symmetric polar caps of colatitudinal radii $\Theta = 10^\circ, 20^\circ, 30^\circ, 40^\circ$. The maximal spherical harmonic degree is $L = 18$; the rounded Shannon numbers are $N = 3, 11, 24, 42$. The sums of squares $g_1^2(\theta, \phi) + g_2^2(\theta, \phi) + \cdots$ and $\lambda_1 g_1^2(\theta, \phi) + \lambda_2 g_2^2(\theta, \phi) + \cdots$ are plotted versus colatitude $\theta$ along a fixed arbitrary meridian $\phi$. Dashed lines show the full unweighted sums of $(L+1)^2$ terms, which attain the constant value $N/A$ over the entire sphere $0 \le \theta \le \pi$. Solid lines show the eigenvalue-weighted partial sums of $N$ terms and the full sums of $(L+1)^2$ terms, which are very nearly equal and concentrated within $0 \le \theta \le \Theta$.</figcaption></figure>

**5.4. Commuting Differential Operator.** The analysis of the time-frequency concentration problem was advanced considerably by Slepian's discovery of commuting prolate spheroidal differential operator. Remarkably, there is also a differential operator that commutes with the integral operator in (5.13) and (5.14b) [20]. Grünbaum, Longhi, and Perlstadt [21] show that, for any $0 \le m \le L$, it is of the form

$$ (5.18) \qquad \mathcal{G} = (\cos \Theta - \cos \theta) \nabla_m^2 + \sin \theta \frac{d}{d\theta} - L(L+2) \cos \theta, $$

where $\nabla_m^2 = d^2/d\theta^2 + \cot \theta (d/d\theta) - m^2(\sin \theta)^{-2}$ is the fixed-order Laplace-Beltrami operator. Rewritten in terms of $\mu = \cos \theta$, the Grünbaum operator (5.18) is

$$ (5.19) \qquad \mathcal{G} = \frac{d}{d\mu} \left[ (\cos \Theta - \mu)(1 - \mu^2) \frac{d}{d\mu} \right] - L(L+2)\mu - \frac{m^2(\cos \Theta - \mu)}{1 - \mu^2}. $$

Since commuting operators have the same eigenfunctions, we can find the spacelimited, fixed-order eigenfunctions $h(\theta)$ or $h(\mu)$ by solving the differential eigenvalue equation $\mathcal{G}h = \chi h$, where $\chi \neq \lambda$ is the associated Grünbaum eigenvalue.
To confirm that the differential operator $\mathcal{G}$ in (5.19) commutes with the integral operator acting on $h(\mu')$ in (5.14b), we are required to show that

$$ (5.20) \qquad \int_{\cos \Theta}^{1} D(\mu, \mu') \mathcal{G}_{\mu'} h(\mu') d\mu' = \int_{\cos \Theta}^{1} \mathcal{G}_{\mu} D(\mu, \mu') h(\mu') d\mu' $$

for an arbitrary spacelimited function $h(\mu')$. Following Grünbaum, Longhi, and Perl-stadt [21], we first show that the left side of (5.20) can be rewritten as

$$ (5.21) \qquad \int_{\cos \Theta}^{1} D(\mu, \mu') \mathcal{G}_{\mu'} h(\mu') d\mu' = \int_{\cos \Theta}^{1} \mathcal{G}_{\mu'} D(\mu, \mu') h(\mu') d\mu', $$

and then we verify that

$$ (5.22) \qquad \mathcal{G}_{\mu} D(\mu, \mu') = \mathcal{G}_{\mu'} D(\mu, \mu'). $$

The first result (5.21) is a straightforward exercise in integration by parts: for any two functions $\zeta(\mu)$ and $\eta(\mu)$, it may be easily shown that

$$ (5.23) \qquad \int_{\cos \Theta}^{1} \zeta(\mathcal{G}\eta) d\mu = - \int_{\cos \Theta}^{1} \left[ (\cos \Theta - \mu)(1 - \mu^2) \frac{d\zeta}{d\mu} \frac{d\eta}{d\mu} + L(L+2)\mu\zeta\eta \right. \\ \qquad \left. + m^2(\cos \Theta - \mu)(1 - \mu^2)^{-1}\zeta\eta \right] d\mu = \int_{\cos \Theta}^{1} (\mathcal{G}\zeta) \eta \, du. $$

Although it is not needed to prove (5.20), we note for future reference that (5.23) is also valid if the integrations are carried out over the full interval $-1 \le \mu \le 1$. To verify the second result (5.22) we use the definition (5.15) of $D(\mu, \mu')$, the identity $\nabla_m^2 P_{lm} = -l(l+1)P_{lm}$, and the notation $A_{lm} = (2l+1)(l-m)/(l+m)!$ to write

$$ (5.24) \qquad (\mathcal{G}_{\mu} - \mathcal{G}_{\mu'}) D(\mu, \mu') = \frac{1}{2}(\mu - \mu') \sum_{l=m}^{L} A_{lm} P_{lm}(\mu) P_{lm}(\mu') [l(l+1) - L(L+2)] \\ \qquad \qquad - \frac{1}{2}(1 - \mu^2) \sum_{l=m}^{L} A_{lm} \frac{d}{d\mu} P_{lm}(\mu) P_{lm}(\mu') \\ \qquad \qquad + \frac{1}{2}(1 - \mu'^2) \sum_{l=m}^{L} A_{lm} P_{lm}(\mu) \frac{d}{d\mu'} P_{lm}(\mu'). $$

An application of the Legendre derivative identity (3.8b) transforms (5.24) into

$$ (5.25) \qquad (\mathcal{G}_{\mu} - \mathcal{G}_{\mu'}) D(\mu, \mu') = \frac{1}{2}(\mu - \mu') \sum_{l=m}^{L} A_{lm} P_{lm}(\mu) P_{lm}(\mu') [l^2 - (L+1)^2] \\ \qquad \qquad + \frac{1}{2} \sum_{l=m}^{L} A_{lm}(l-m+1) [P_{l+1, m}(\mu) P_{lm}(\mu') - P_{lm}(\mu) P_{l+1, m}(\mu')], $$

and the Christoffel–Darboux identity (3.10) transforms (5.25) into
$$ (G_{\mu} - G_{\mu'})D(\mu, \mu') = \frac{1}{2}(\mu - \mu') \sum_{l=m}^{L} A_{lm}P_{lm}(\mu)P_{lm}(\mu') [l^2 - (L+1)^2] \\ + \frac{1}{2}(\mu - \mu') \sum_{l=m}^{L} (2l+1) \sum_{n=m}^{l} A_{nm}P_{nm}(\mu)P_{nm}(\mu') \\ (5.26) \qquad = \frac{1}{2}(\mu - \mu') \sum_{l=m}^{L} A_{lm}P_{lm}(\mu)P_{lm}(\mu') \left[ l^2 - (L+1)^2 + \sum_{n=l}^{L} (2n+1) \right], $$

where we obtained the last equality by interchanging the order of summation and relabeling the sums. The final sum over $n$ is equal to $-l^2 + (L+1)^2$; therefore, $(G_{\mu} - G_{\mu'})D(\mu, \mu') = 0$, and the commutation relation (5.20) is confirmed.

**5.5. Grünbaum's Equation.** The above argument shows that we can compute the fixed-order, spacelimited, colatitudinal eigenfunctions $h_1(\theta), h_2(\theta), \dots, h_{L-m+1}(\theta)$ by solving either the integral equation (5.13) or the differential equation

$$ (5.27) \qquad (\cos \Theta - \cos \theta) \nabla_m^2 h + \sin \theta \frac{dh}{d\theta} - L(L+2) \cos \theta h = \chi h, \quad 0 \le \theta \le \Theta. $$

The equivalent equation in terms of $\mu = \cos \theta$ is in standard Sturm-Liouville form [8]:

$$ (5.28) \qquad (ph')' - qh + \chi\rho h = 0, \quad \cos \Theta \le \mu \le 1, $$

where $p(\mu) = (\mu - \cos \Theta)(1 - \mu^2)$, $q(\mu) = m^2(1 - \mu^2)^{-1}(\mu - \cos \Theta) - L(L+2)\mu$, $\rho(\mu) = 1$, and the prime denotes differentiation with respect to $\mu$. Equation (5.28) must be solved subject to the requirement that $h(\mu)$ remain finite at the endpoints $\mu = \cos \Theta$ and $\mu = 1$. The associated variational problem is [8]

$$ (5.29) \qquad \chi = \frac{\int_{\cos \Theta}^{1} (ph'^2 + qh^2) d\mu}{\int_{\cos \Theta}^{1} \rho h^2 d\mu} = \text{minimum.} $$

All of the familiar Sturm-Liouville theorems apply. In particular, we know that (5.28) has a simple spectrum, with an infinite number of distinct eigenvalues $\chi_1 < \chi_2 < \dots$, having an accumulation point at infinity. The rank orderings of the eigenvalues $\chi_1, \chi_2, \dots$ and the spatiowspectral concentration factors $\lambda_1, \lambda_2, \dots, \lambda_{L-m+1}$ are reversed, so that the eigenfunction $h_1(\theta)$ associated with the numerically smallest eigenvalue $\chi_1$, which has no nodes in the polar cap $0 \le \theta \le \Theta$, is the best concentrated fixed-order eigenfunction; $h_2(\theta)$, which has exactly one node, is the next best concentrated; and so on. Only the first $L-m+1$ eigenfunctions $h_1(\theta), h_2(\theta), \dots, h_{L-m+1}(\theta)$ with nonzero eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_{L-m+1}$ are of interest in most applications. The remaining eigenfunctions $h_{L-m+2}(\theta), h_{L-m+3}(\theta), \dots$ are in the null space of the integral equation (5.13) [42].

**5.6. Commuting Tridiagonal Matrix.** As in the case of (5.10) and (5.13), we are free to extend the domain of (5.27) to the entire domain $0 \le \theta \le \pi$; in that case, the unknown function must be bandlimited rather than spacelimited:

$$ (5.30) \qquad (\cos \Theta - \cos \theta) \nabla_m^2 g + \sin \theta \frac{dg}{d\theta} - L(L+2) \cos \theta g = \chi g, \quad 0 \le \theta \le \pi. $$
Upon substituting the harmonic representation (5.7) of $g(\theta)$ into (5.30), multiplying both sides by $2\pi \sin\theta X_{lm}(\theta)$, integrating over $0 \le \theta \le \pi$, and invoking the orthogonality relation (3.5), we obtain the algebraic eigenvalue equation

$$ (5.31) \qquad \mathbf{Gg} = \chi \mathbf{g}, \quad \text{where} \quad \mathbf{G} = \begin{pmatrix} G_{mm} & \cdots & G_{mL} \\ \vdots & \ddots & \vdots \\ G_{Lm} & \cdots & G_{LL} \end{pmatrix} $$

is the $(L-m+1) \times (L-m+1)$ matrix with elements

$$ (5.32) \qquad G_{ll'} = 2\pi \int_0^\pi X_{lm}(\mathcal{G}X_{l'm}) \sin\theta \, d\theta. $$

Equation (5.31) is the spectral-domain version of the differential eigenvalue equation (5.30) just as (5.2) is the spectral equivalent of the integral equation (5.10).

The symmetry relation (5.23) is valid even if the interval of integration is extended to $0 \le \theta \le \pi$, as in (5.32). This shows that the Grünbaum matrix $\mathbf{G}$ is symmetric: $\mathbf{G} = \mathbf{G}^\mathrm{T}$. In addition, the matrices $\mathbf{D}$ and $\mathbf{G}$ commute, $\mathbf{DG} = \mathbf{GD}$, so they have identical eigenvectors. The index version of the commutation relation is

$$ (5.33) \qquad \sum_{n=m}^{L} D_{ln} G_{nl'} = 2\pi \int_{0}^{\Theta} X_{lm}(\mathcal{G} X_{l'm}) \sin\theta \, d\theta = \sum_{n=m}^{L} G_{ln} D_{nl'}. $$

The interior expression involving both the operator $\mathcal{G}$ and integration over the region of concentration $0 \le \theta \le \Theta$ is the $ll'$ or $l'l$ element of the symmetric matrix product $\mathbf{DG} = (\mathbf{DG})^\mathrm{T}$. Verification of the intermediate steps requires the use of the orthogonality relation (3.5), the operator identity (5.22), and both the symmetry relation (5.23) and its extension to the interval $0 \le \theta \le \pi$.

There are a number of ways to evaluate the matrix elements (5.32), perhaps the most straightforward of which is to make use of the relation $\nabla_m^2 X_{lm} = -l(l+1)X_{lm}$ and the Legendre identities (3.8). In fact, the Grünbaum matrix $\mathbf{G}$ is tridiagonal:

$$ (5.34a) \qquad G_{ll'} = -l(l+1) \cos \Theta, $$

$$ (5.34b) \qquad G_{l l+1} = G_{l+1 l} = [l(l+2) - L(L+2)] \sqrt{\frac{(l+1)^2 - m^2}{(2l+1)(2l+3)}}, $$

$$ (5.34c) \qquad G_{ll'} = 0 \quad \text{otherwise.} $$

The solution of equation (5.31) offers a particularly attractive means of computing the $(L-m+1)$-dimensional eigenvectors $\mathbf{g}$ and thus the optimally concentrated polar cap eigenfunctions $g(\theta) \in S_L$, because it only requires the numerical diagonalization of a tridiagonal matrix $\mathbf{G}$ with analytically prescribed elements (5.34) and a spectrum of eigenvalues $\chi$ that is guaranteed to be simple.

Among other things, Grünbaum's procedure enables the stable computation of bandlimited functions that are optimally concentrated in a large rather than a small region of $\Omega$. To illustrate this, we show in Figure 5.6 the first four $(g_1, g_2, g_3, g_4)$ and the last four $(g_{16}, g_{17}, g_{18}, g_{19})$ zonal ($m=0$) eigenfunctions for a polar cap of radius $\Theta = 40^\circ$ and a maximal harmonic degree $L=18$. As noted in section 4.3, the eigenfunctions $g_{16}, g_{17}, g_{18}, g_{19}$ that are optimally excluded from the polar cap $\Theta = 40^\circ$ are optimally concentrated within the much larger antipodal cap $\Theta = 140^\circ$. The actual eigenvalues $\lambda_{16}, \lambda_{17}, \lambda_{18}, \lambda_{19}$ are many orders of magnitude smaller than the listed values, which simply represent the noise floor of our double-precision computations.
<figure><img src="image_8.png" /><figcaption>**Fig. 5.6** Optimally concentrated (top row) and optimally excluded (bottom row) zonal ($m = 0$) eigenfunctions for a circular polar cap of colatitudinal radius $\Theta = 40^{\circ}$ and a maximal spherical harmonic degree $L = 18$. The optimally excluded eigenfunctions cannot be accurately computed by double-precision diagonalization of the matrix D. Solution of the concentration problem for a polar cap of radius $\Theta = 140^{\circ}$ gives rise to the same eigenfunctions $g_1(\theta), g_2(\theta), g_3(\theta), g_4(\theta), \dots, g_{16}(\theta), g_{17}(\theta), g_{18}(\theta), g_{19}(\theta)$, but in reverse order.</figcaption></figure>

**5.7. Abstract Operator Formulation.** The spatial-domain commutation relation (5.20) can be expressed using the operator notation of section 4.4 as

$$ (5.35) \qquad (\mathcal{R}\mathcal{H}^{-1}\mathcal{L}\mathcal{H}\mathcal{R})\mathcal{G} = \mathcal{G}(\mathcal{R}\mathcal{H}^{-1}\mathcal{L}\mathcal{H}\mathcal{R}). $$

Upon premultiplying (5.35) by $\mathcal{L}\mathcal{H}$, postmultiplying it by $\mathcal{H}^{-1}\mathcal{L}$, and making use of the relations $\mathcal{G} = \mathcal{G}\mathcal{R} = \mathcal{R}\mathcal{G}$ and the fact that $\mathcal{L}^2 = \mathcal{L}$, we obtain

$$ (5.36) \qquad (\mathcal{L}\mathcal{H}\mathcal{R}\mathcal{H}^{-1}\mathcal{L})(\mathcal{L}\mathcal{H}\mathcal{G}\mathcal{H}^{-1}\mathcal{L}) = (\mathcal{L}\mathcal{H}\mathcal{G}\mathcal{H}^{-1}\mathcal{L})(\mathcal{L}\mathcal{H}\mathcal{R}\mathcal{H}^{-1}\mathcal{L}). $$

Equation (5.36) is the abstract operator formulation of the spectral-domain matrix commutation relation $\mathcal{D}\mathcal{G} = \mathcal{G}\mathcal{D}$. The operator equivalents of the spectral-domain eigenvalue equation (5.31) and the spatial-domain eigenvalue equation (5.27) are

$$ (5.37) \qquad (\mathcal{L}\mathcal{H}\mathcal{G}\mathcal{H}^{-1}\mathcal{L})(\mathcal{L}\mathcal{f}) = \chi(\mathcal{L}\mathcal{f}), \qquad \mathcal{G}(\mathcal{R}\mathcal{f}) = \chi(\mathcal{R}\mathcal{f}), $$

where $\mathcal{f}(\theta)$ is an arbitrary colatitudinal function that is neither bandlimited nor space-limited. Because of the commutation relations (5.35) and (5.36), we are free to solve (5.37) rather than the fixed-order version of (4.25) to find the bandlimited eigenvectors $\mathcal{g} = \mathcal{L}\mathcal{f}$ and the spacelimited eigenfunctions $\mathcal{h}(\theta) = \mathcal{R}\mathcal{f}(\theta)$.

**6. Continental Concentration.** To illustrate the theory for irregularly shaped regions, we consider the spatiospectral concentration in six of Earth's continental regions, listed in Table 6.1 with their rounded Shannon numbers $N = (L+1)^2A/(4\pi)$ for different bandwidths. The spherical Slepian functions should be ideally suited to the spectral analysis of data within either Earth's continents or oceans, as required in geodesy, geophysics, and oceanography, e.g., [26, 27, 47, 52].
<table><caption>**Table 6.1** Fractional areas, Shannon numbers, and bandwidths for the continental concentration problem.</caption><thead><tr><th rowspan="2">Continental region</th><th rowspan="2">Fractional area<br>A/(4π) in %</th><th colspan="4">Shannon number N</th></tr><tr><th>L = 6</th><th>L = 12</th><th>L = 18</th><th>L = 24</th></tr></thead><tbody><tr><td>Greenland</td><td>0.44</td><td>0</td><td>1</td><td>2</td><td>3</td></tr><tr><td>Australia</td><td>1.50</td><td>1</td><td>3</td><td>5</td><td>9</td></tr><tr><td>South America</td><td>3.50</td><td>2</td><td>6</td><td>13</td><td>22</td></tr><tr><td>North America</td><td>3.98</td><td>2</td><td>7</td><td>14</td><td>25</td></tr><tr><td>Africa</td><td>5.78</td><td>3</td><td>10</td><td>21</td><td>36</td></tr><tr><td>Eurasia</td><td>9.98</td><td>5</td><td>17</td><td>36</td><td>62</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.1** Eigenvalue spectra for Greenland, Australia, North America, Africa, and Eurasia. From upper left to lower right, four different bandwidths, L = 6, 12, 18, 24, are considered. The horizontal axis in each panel is truncated; the total number of eigenvalues (L + 1)² = 49, 169, 361, 625 appears to the right of the arrow. Vertical gridlines and the five leftmost ordinate labels specify the rounded Shannon numbers N.</figcaption></figure>

Figure 6.1 shows the eigenvalue spectra for five of the six regions (Greenland, Australia, North America, Africa, and Eurasia) and four maximal spherical harmonic degrees, L = 6, 12, 18, 24, corresponding to (L + 1)² = 49, 169, 391, 625 eigenfunctions each. The minimum wavelength associated with a bandwidth limit L is 2π/√L(L + 1)/ ≈ 2π/(L + 1/2) multiplied by Earth's radius [29]. The cutoff wavelengths corresponding to the choices L = 6, 12, 18, and 24 are 6200, 3200, 2200, and 1600 km, respectively. Only the largest continent, Eurasia, is sizable enough to exhibit at least one nearly perfectly concentrated eigenfunction for the smallest bandwidth, L = 6, and the smallest region considered, Greenland, is too tiny to exhibit even a single eigenfunction with a concentration factor λ near unity for the largest bandwidth,
<figure><img src="image_10.png" /><figcaption>**Fig. 6.2** Bandlimited $L = 18$ eigenfunctions $g_1, g_2, \dots, g_{12}$ that are optimally concentrated within Australia. The concentration factors $\lambda_1, \lambda_2, \dots, \lambda_{12}$ are indicated; the rounded Shannon number is $N = 5$. Order of concentration is left to right, top to bottom, as with Flemish text. Positive values are blue and negative values are red (though, as we have noted, these could be reversed, since the sign of an eigenfunction is arbitrary). Regions in which the absolute value is less than one hundredth of the maximum value on the sphere are left white.</figcaption></figure>

$L = 24$. As in the case of a polar cap (Figure 5.3), the rounded Shannon numbers $N$ roughly separate the well-concentrated eigenfunctions with associated eigenvalues $\lambda \ge 0.5$ from the poorly concentrated ones with eigenvalues $\lambda < 0.5$.

In Figures 6.2 and 6.3 we show map views of the first twelve $L = 18$ eigenfunctions $g_1(\hat{\mathbf{r}}), g_2(\hat{\mathbf{r}}), \dots, g_{12}(\hat{\mathbf{r}})$ that are optimally concentrated within the continents Australia and Africa. In the case of Australia (Figure 6.2) the first five eigenfunctions are reasonably well concentrated within the continental boundaries ($\lambda_5 = 0.607$); however, the concentration factors $\lambda$ diminish rapidly thereafter, so that $g_{12}$ is far more excluded than concentrated ($\lambda_{12} = 0.049$). With a limiting bandwidth $L = 18$, and thus a cutoff wavelength of 2200 km, it is only possible to concentrate $N = 5$ orthogonal bandlimited eigenfunctions $g_1, g_2, g_3, g_4, g_5$ into a continent which, across its north-south waist, is only about 1500 km wide. This situation is improved in the case of Africa (Figure 6.3), which has an area that is 3.9 times larger than that of Australia. In fact, Africa has $N = 21$ reasonably well concentrated $L = 18$ eigenfunctions, of which only the first twelve ($\lambda_{12} = 0.887$) are shown. The first eigenfunction, $g_1$, is a roughly circular dome centered in the middle of the continent, as in the case of Australia. Subsequent orthogonal eigenfunctions $g_2, g_3, \dots$ exhibit lobes in previously uncovered
<figure><img src="image_11.png" /><figcaption>**Fig. 6.3** Bandlimited $L = 18$ eigenfunctions $g_1, g_2, \dots, g_{12}$ that are optimally concentrated within Africa. The concentration factors $\lambda_1, \lambda_2, \dots, \lambda_{12}$ are indicated; the rounded Shannon number is $N = 21$. Format is identical to that in Figure 6.2.</figcaption></figure>

regions. West Africa is uncovered by $g_1$ and $g_2$, but becomes reasonably well covered by $g_3$ and $g_5$; likewise, Southern Africa is uncovered until $g_4$ and $g_5$. Other, smaller, geographical features successively become well covered by the increasingly oscillatory orthogonal eigenfunctions (e.g., Egypt by $g_7$ and $g_{12}$).

Finally, Figure 6.4 shows the eigenvalue-weighted sum of squares $\sum_{\alpha} \lambda_{\alpha} g_{\alpha}^2(\hat{\mathbf{r}})$ of the bandlimited $L = 18$ eigenfunctions of all six of Earth's landmasses (excluding Antarctica). We find the eigenfunctions $g_1(\hat{\mathbf{r}}), g_2(\hat{\mathbf{r}}), \dots, g_{(L+1)^2}(\hat{\mathbf{r}})$ by diagonalizing the $(L+1)^2 \times (L+1)^2$ matrix (4.4) formed by summing the corresponding matrices $\mathbf{D}_{\text{Eurasia}} + \mathbf{D}_{\text{Africa}} + \dots$ of each of the regions. The combined fractional area of all six regions is $A/(4\pi) = 25.2\%$, and the rounded Shannon number is $N = 91$; the partial sums of the first $N/4$, $N/2$, and $N$ terms, as well as the full sum of all $(L+1)^2 = 361$ terms, are shown. The ability of the first $N$ eigenfunctions to uniformly cover the target region is evident; as in Figure 5.5, the coverage is only marginally improved by adding the remaining, poorly concentrated, $(L+1)^2 - N = 250$ terms. Due to their small size, Australia and Greenland do not appear until the $1 \to N/2$ and $1 \to N$ partial sums, respectively. Even then, the coverage of Greenland is imperfect, an expected consequence of its small Shannon number ($N = 2$ for a bandwidth $L = 18$).
<figure><img src="image_12.png" /><figcaption>**Fig. 6.4** Cumulative eigenvalue-weighted energy of the first N/4, N/2, N and all (L + 1)² eigenfunctions that are optimally concentrated within the ensemble of Eurasia, Africa, North America, South America, Australia, and Greenland. The bandwidth is L = 18; the cumulative fractional area is A/(4π) = 25.2%; the rounded Shannon number is N = 91. The darkest blue on the color bar corresponds to the expected value (4.21) of the sum, as shown. Regions where the value is smaller than one hundredth of the maximum value on the sphere are left white.</figcaption></figure>

**7. Asymptotic Scaling.** As we have noted, the eigenvalues $\lambda_1, \lambda_2, \dots$ and suitably scaled eigenfunctions $\psi_1(x), \psi(2)_x, \dots$ of the original Slepian concentration problem (2.6) depend only upon the Shannon number $N = 2TW/\pi$. This scaling is the only important feature of the one-dimensional problem that does not carry over to the spatospectral concentration problem on a sphere. Fundamentally, this lack of scaling is a consequence of the fact that it is not possible to shrink or magnify a region, such as Africa, on a sphere $\Omega$ of fixed radius $\|\hat{\mathbf{r}}\| = 1$, while keeping the angular relationships among all of the interior points the same. Shannon-number scaling on a sphere is exhibited only asymptotically in the limit

$$ (7.1) \qquad A \to 0, \quad L \to \infty, \quad \text{with } N = (L+1)^2 \frac{A}{4\pi} \quad \text{held fixed.} $$

In that limit of a small concentration area A and a large bandwidth L, the curvature of the sphere becomes negligible and the spherical concentration problem becomes identical to the two-dimensional concentration problem in the plane [56].

**7.1. Hilb Approximation and Poisson Sum Formula.** Two results underlie the consideration of the “flat-Earth” limit (7.1), which we undertake in this section. The first is Hilb’s asymptotic approximation for the Legendre functions [2, 9, 24, 65],

$$ (7.2) \qquad X_{lm}(\theta) \approx (-1)^m \sqrt{\frac{l+1/2}{2\pi}} \sqrt{\frac{\theta}{\sin\theta}} J_m\left[(l+1/2)\theta\right], \quad 0 \le \theta \ll \pi, $$
where $J_m(x)$ is the Bessel function of the first kind, and the second is the truncated Poisson sum formula, valid for an arbitrary $2\pi$-periodic function $f(x)$:

$$ (7.3) \qquad \sum_{l=0}^{L} f(l + 1/2) = \sum_{s=-\infty}^{\infty} (-1)^s \int_{0}^{L+1} f(k)e^{-2\pi i sk} dk. $$

**7.2. Scaled Integral Equation for an Arbitrary Region.** Making use of (7.2) and (7.3), we can approximate the kernel $D(\hat{\mathbf{r}}, \hat{\mathbf{r}}')$ in (4.11) by

$$ (7.4) \qquad D(\Delta) \approx \frac{1}{2\pi} \sqrt{\frac{\Delta}{\sin \Delta}} \sum_{s=-\infty}^{\infty} (-1)^s \int_{0}^{L+1} J_0(k\Delta)e^{-2\pi i sk} k \, dk. $$

Upon substituting $k = (L+1)p$ and taking the limit $L \to \infty, \Delta \to 0$, with the product $L\Delta$ held fixed, (7.4) reduces to

$$ (7.5) \qquad D(\Delta) \approx \frac{(L+1)^2}{2\pi} \int_0^1 J_0[(L+1)p\Delta] p \, dp = \frac{(L+1)J_1[(L+1)\Delta]}{2\pi\Delta}, $$

where we have made the approximation $\Delta/\sin\Delta \approx 1$ and used the Riemann-Lebesgue lemma [45] to eliminate the $s \neq 0$ terms involving the highly oscillatory factors $e^{-2\pi i s(L+1)p}$. In the limit $x \to 0$, the ratio $J_1(x)/x \to 1/2$, so the $\Delta \to 0$ limit of the kernel (7.5) is $D(0) = (L+1)^2/(4\pi)$, guaranteeing that the Shannon number, or sum of the eigenvalues (4.18), is still given in this asymptotic approximation by

$$ (7.6) \qquad N = \int_R D(0) \, d\Omega = (L+1)^2 \frac{A}{4\pi}. $$

To obtain a scaled version of (4.10) dependent only upon the Shannon number $N$, we introduce the independent and dependent variable transformations

$$ (7.7) \qquad \mathbf{x} = \sqrt{\frac{4\pi}{A}} \hat{\mathbf{r}}, \quad \mathbf{x}' = \sqrt{\frac{4\pi}{A}} \hat{\mathbf{r}}', \quad \psi(\mathbf{x}) = g(\hat{\mathbf{r}}), \quad \psi(\mathbf{x}') = g(\hat{\mathbf{r}}'). $$

The scaled coordinates $\mathbf{x}, \mathbf{x}'$ are the projections of the points $\hat{\mathbf{r}}, \hat{\mathbf{r}}' \in \Omega$ onto a large sphere $\Omega_*$ of squared radius $||\mathbf{x}||^2 = 4\pi/A$. The geodesic distance between the scaled points $\mathbf{x}, \mathbf{x}' \in \Omega_*$ and the differential surface area on $\Omega_*$ are

$$ (7.8) \qquad \|\mathbf{x} - \mathbf{x}'\| = \sqrt{\frac{4\pi}{A}} \Delta, \quad d\Omega_* = \frac{4\pi}{A} d\Omega. $$

Upon making the substitutions (7.7)–(7.8), equations (4.10) and (7.5) reduce to

$$ (7.9) \qquad \int_{R_*} D_*(\mathbf{x}, \mathbf{x}') \psi(\mathbf{x}') \, d\Omega_*' = \lambda \psi(\mathbf{x}), $$

where $R_*$, of area $4\pi$, is the projection of the region $R$ onto the scaled sphere $\Omega_*$, and

$$ (7.10) \qquad D_*(\mathbf{x}, \mathbf{x}') = \frac{\sqrt{N} J_1(\sqrt{N} \|\mathbf{x} - \mathbf{x}'\|)}{2\pi \|\mathbf{x} - \mathbf{x}'\|} $$

is the symmetric, $N$-dependent Fredholm kernel. Equations (7.9)–(7.10) are the two-dimensional planar analogue of the one-dimensional scaled eigenvalue equation (2.6). The flat-Earth eigenvalues $\lambda_1, \lambda_2, \dots$ and scaled eigenfunctions $\psi_1(\mathbf{x}), \psi_2(\mathbf{x}), \dots$ depend upon the maximal degree $L$ and the area $A$ only through the Shannon number $N = (L+1)^2 A/(4\pi)$. As in the case of (4.10) and (4.15), we are free to solve equations (7.9)–(7.10) either on all of $\Omega_*$, in which case the eigenfunctions are bandlimited, or only in the region of concentration $R_*$, in which case they are spacelimited.
**7.3. Scaled Eigenvalue Equation for an Axisymmetric Polar Cap.** The flat-Earth asymptotic version of the fixed-order colatitudinal eigenvalue problem (5.10) can be obtained in two different ways: either by an application of the Hilb approximation (7.2) and the Poisson sum formula (7.3) to the kernel $D(\theta, \theta')$ given in (5.11), or by using the addition theorem for Bessel functions [30],

$$ (7.11) \qquad J_0(k\Delta) = J_0(k\theta)J_0(k\theta') + 2 \sum_{m=1}^{\infty} J_m(k\theta)J_m(k\theta') \cos m(\phi - \phi'), $$

together with the representation (5.9); the orthonormality of the longitudinal factors ..., $\sqrt{2} \cos m\phi$, ..., $1$, ..., $\sqrt{2} \sin m\phi$, ... over the interval $0 \le \phi < 2\pi$; and the first expression for the asymptotic kernel of (7.5) to decompose (4.10) into a series of individual eigenvalue problems, one for each order $0 \le m \le L$. Using either method, we find that (5.10) can be approximated in the limit (7.1) by

$$ (7.12) \qquad \int_0^\Theta D(\theta, \theta') g(\theta') \theta' d\theta' = \lambda g(\theta), $$

where we have approximated $\sin \theta' d\theta' \approx \theta' d\theta'$, and where

$$ (7.13) \qquad D(\theta, \theta') = (L+1)^2 \int_0^1 J_m[(L+1)p\theta] J_m[(L+1)p\theta'] p \, dp. $$

It is convenient in the present instance to approximate the area of the small polar cap by $A \approx \pi\Theta^2$, such that $N \approx (L+1)^2\Theta^2/4$, and to introduce scaled coordinates that are slightly different from those in (7.7), namely, $x = \theta/\Theta$, $x' = \theta'/\Theta$, and $\psi(x) = g(\theta)$, $\psi(x') = g(\theta')$. This leads to a scaled, fixed-order eigenvalue problem,

$$ (7.14) \qquad \int_0^1 D_*(x, x') \psi(x') x' dx' = \lambda \psi(x), $$

with an associated kernel

$$ (7.15) \qquad D_*(x, x') = 4N \int_0^1 J_m(2\sqrt{N}px) J_m(2\sqrt{N}px') p \, dp, $$

whose eigenvalues $\lambda_1, \lambda_2, \dots$ and scaled eigenfunctions $\psi_1(x), \psi_2(x), \dots$ depend upon the maximal spherical harmonic degree $L$ and the cap radius $\Theta$ only through the small polar-cap Shannon number $N$. Slepian [56] has noted that (7.14)–(7.15) are an iterated version of the equivalent “square root” equation

$$ (7.16) \qquad 2\sqrt{N} \int_0^1 J_m(2\sqrt{N}xx') \psi(x') x' dx' = \sqrt{\lambda} \psi(x). $$

In principle, these asymptotic results would enable the determination of approximate polar cap eigenfunctions $g(\theta)$ for varying values of $L$ and $\Theta$ by scaling a precomputed catalog of fixed-$N$ eigenfunctions. In practice, the construction and diagonalization of the tridiagonal Grünbaum matrix (5.34) is so straightforward and efficient that it is preferable to simply compute the optimally concentrated eigenfunctions $g(\theta)$ exactly.
**7.4. Asymptotic Fixed-Order Shannon Number.** The asymptotic approximation to the number of significant eigenvalues associated with a given order *m* is

$$ (7.17) \qquad \begin{aligned} N_m &= \int_0^1 D_*(x, x) x \, dx = 4N \int_0^1 \int_0^1 J_m^2(2\sqrt{N}px) \, p \, dp \, x \, dx \\ &= 2N \left[ J_m^2(2\sqrt{N}) + J_{m+1}^2(2\sqrt{N}) \right] - (2m+1)\sqrt{N}J_m(2\sqrt{N})J_{m+1}(2\sqrt{N}) \\ & \quad - \frac{m}{2} \left[ 1 - J_0^2(2\sqrt{N}) - 2 \sum_{n=1}^m J_n^2(2\sqrt{N}) \right]. \end{aligned} $$

The relationship $N = N_0 + 2 \sum_{m=1}^{\infty} N_m$ between the total number of significant eigenvalues and the number associated with each order *m* is preserved in this asymptotic approximation by virtue of the identity $J_0^2(x) + 2 \sum_{m=1}^{\infty} J_m^2(x) = 1$. The number of significant *m* = 0 eigenvalues can be even more simply approximated by the relation $N_0 \approx 2\sqrt{N}/\pi \approx (L+1)\Theta/\pi$ (see [76]), which can be derived from (7.17) using the large-argument asymptotic expansion of the Bessel function [30, 45]. The lengthier result (7.17) is exact in the case of concentration within a two-dimensional plane.

**8. Conclusion.** An orthogonal family of bandlimited spherical harmonic expansions that are optimally concentrated within a finite region *R* of the unit sphere can be computed by solving either a symmetric matrix eigenvalue problem in the spectral domain or an equivalent Fredholm integral eigenvalue problem in the spatial domain. Every eigenvalue $0 < \lambda < 1$ is a measure of both the spatial concentration of the bandlimited eigenfunction $g(\hat{\mathbf{r}})$ and the spectral concentration of the spacelimited eigenfunction $h(\hat{\mathbf{r}})$ that coincides with $g(\hat{\mathbf{r}})$ inside the region of concentration. The number of well-concentrated eigenfunctions is $N = (L+1)^2 A/(4\pi)$, where *L* is the bandwidth and *A* is the area of the region of concentration. Roughly speaking, this Shannon number *N* is the dimension of the space of functions $f(\hat{\mathbf{r}})$ that can be concentrated both within a finite region *R* of the sphere and within a spectral interval $0 \le l \le L$. For a small region $A \ll 4\pi$ and a moderate maximal spherical harmonic degree *L*, the optimally concentrated bandlimited eigenfunctions $g(\hat{\mathbf{r}})$ and associated spacelimited eigenfunctions $h(\hat{\mathbf{r}})$ can be computed accurately, even for an irregularly shaped region *R*. In the special, but important, case of a circular polar cap, every eigenfunction can be computed accurately by numerical diagonalization of a commuting tridiagonal matrix, which has a simple Sturm-Liouville spectrum. Just as Slepian's prolate spheroidal eigentapers have proven to be extremely useful for spectral analysis in Cartesian geometry, we expect the spherical eigenfunctions developed here to have a wide variety of spatospectral data analysis applications in fields such as geophysics, planetary science, and cosmology.

**Acknowledgments.** F.J.S. thanks Sofia Akber-Knutson, Ingrid Daubechies, Peter Harris, Partha Mitra, Jean Steiner, Bill Symes, and Carl Wunsch for comments and discussions and the Département de Géophysique Spatiale et Planétaire at the Institut de Physique du Globe de Paris for its hospitality and financial support. We appreciate the reviews by David Thomson, John Wahr, and an anonymous third referee and thank Margaret Wright for her expert editorial handling of the manuscript, as well as Julia Cochrane and Louis Primus for their skillful copyediting.
REFERENCES

[1] A. ALBERTELLA, F. SANSÒ, AND N. SNEEUW, *Band-limited functions on a bounded spherical domain: The Slepian problem on the sphere*, J. Geodesy, 73 (1999), pp. 436–447.

[2] R. D. AMADO, K. STRICKER-BAUER, AND D. A. SPARROW, *Semiclassical methods and the summation of the scattering partial wave series*, Phys. Rev. C, 32 (1985), pp. 329–332.

[3] M. A. BLANCO, M. FLÓREZ, AND M. BERMEJO, *Evaluation of the rotation matrices in the basis of real spherical harmonics*, J. Mol. Struct. (Theochem.), 419 (1997), pp. 19–27.

[4] M. BÖHME AND D. POTTS, *A fast algorithm for filtering and wavelet decomposition on the sphere*, Electron. Trans. Numer. Anal., 16 (2003), pp. 70–93.

[5] T. P. BRONEZ, *Spectral estimation of irregularly sampled multidimensional processes by generalized prolate spheroidal sequences*, IEEE Trans. Acoust. Speech Signal Process., 36 (1988), pp. 1862–1873.

[6] W. E. BYERLY, *An Elementary Treatise on Fourier's Series and Spherical, Cylindrical, and Ellipsoidal Harmonics*, Ginn & Co., Boston, 1893.

[7] L. COLZANI, *Regularity of spherical means and localization of spherical harmonic expansions*, J. Aust. Math. Soc. Ser. A, 41 (1986), pp. 287–297.

[8] R. COURANT AND D. HILBERT, *Methods of Mathematical Physics*, Interscience, New York, 1953.

[9] F. A. DAHLEN, *A uniformly valid asymptotic representation of normal mode multiplet spectra on a laterally heterogeneous Earth*, Geophys. J. R. Astron. Soc., 62 (1980), pp. 225–247.

[10] F. A. DAHLEN AND J. TROMP, *Theoretical Global Seismology*, Princeton University Press, Princeton, NJ, 1998.

[11] I. DAUBECHIES, *Time-frequency localization operators: A geometric phase space approach*, IEEE Trans. Inform. Theory, 34 (1988), pp. 605–612.

[12] I. DAUBECHIES, *Ten Lectures on Wavelets*, SIAM, Philadelphia, 1992.

[13] I. DAUBECHIES AND T. PAUL, *Time-frequency localisation operators—A geometric phase space approach: II. The use of dilations*, Inverse Problems, 4 (1988), pp. 661–680.

[14] D. L. DONOHO AND P. B. STARK, *Uncertainty principles and signal recovery*, SIAM J. Appl. Math., 49 (1989), pp. 906–931.

[15] A. R. EDMONDS, *Angular Momentum in Quantum Mechanics*, Princeton University Press, Princeton, NJ, 1996.

[16] N. L. FERNÁNDEZ, *Polynomial bases on the sphere*, in Advanced Problems in Constructive Approximation, M. D. Buhmann and D. H. Mache, eds., Int. Ser. Numer. Math., 142, Birkhäuser, Basel, 2002, pp. 39–52.

[17] P. FLANDRIN, *Time-Frequency/Time-Scale Analysis*, Academic Press, San Diego, CA, 1999.

[18] W. FREEDEN AND V. MICHEL, *Orthogonal zonal, tesseral and sectorial wavelets on the sphere for the analysis of satellite data*, Adv. Comput. Math., 21 (2004), pp. 181–217.

[19] W. FREEDEN AND U. WINDHEUSER, *Combined spherical harmonic and wavelet expansion—A future concept in Earth’s gravitational determination*, Appl. Comput. Harmon. Anal., 4 (1997), pp. 1–37.

[20] E. N. GILBERT AND D. SLEPIAN, *Doubly orthogonal concentrated polynomials*, SIAM J. Math. Anal., 8 (1977), pp. 290–319.

[21] F. A. GRÜNBAUM, L. LONGHI, AND M. PERLSTADT, *Differential operators commuting with finite convolution integral operators: Some non-Abelian examples*, SIAM J. Appl. Math., 42 (1982), pp. 941–955.

[22] A. HANSSEN, *Multidimensional multitaper spectral estimation*, Signal Process., 58 (1997), pp. 327–332.

[23] M. G. HAUSER AND P. J. E. PEEBLES, *Statistical analysis of catalogs of extragalactic objects. II. The Abell catalog of rich clusters*, Astrophys. J., 185 (1973), pp. 757–785.

[24] E. HILB, *Über die Laplacesche Reihe*, Math. Z., 5 (1919), p. 17.

[25] R. A. HORN AND C. R. JOHNSON, *Matrix Analysis*, Cambridge University Press, Cambridge, UK, 1990.

[26] C. HWANG, *Spectral analysis using orthonormal functions with a case study on sea surface topography*, Geophys. J. Int., 115 (1993), pp. 1148–1160.

[27] C. HWANG AND S.-K. CHEN, *Fully normalized spherical cap harmonics: Application to the analysis of sea-level data from TOPEX/POSEIDON and ERS-1*, Geophys. J. Int., 129 (1997), pp. 450–460.

[28] R. JAKOB-CHIEN AND B. K. ALPERT, *A fast spherical filter with uniform resolution*, J. Comput. Phys., 136 (1997), pp. 580–584.

[29] J. JEANS, *The propagation of earthquake waves*, Philos. Trans. R. Soc. Lond. Ser. A, 102 (1923), pp. 554–574.

[30] H. JEFFREYS AND B. S. JEFFREYS, *Methods of Mathematical Physics*, 3rd ed., Cambridge University Press, Cambridge, UK, 1988.
[31] R. P. KANWAL, *Linear Integral Equations: Theory and Technique*, Academic Press, New York, 1971.

[32] M. KIDO, D. A. YUEN, AND A. P. VINCENT, *Continuous wavelet-like filter for a spherical surface and its application to localized admittance function on Mars*, Phys. Earth Planet. Int., 135 (2003), pp. 1-14.

[33] H. J. LANDAU, *On the eigenvalue behavior of certain convolution equations*, Trans. Amer. Math. Soc., 115 (1965), pp. 242-256.

[34] H. J. LANDAU, *Necessary density conditions for sampling and interpolation of certain entire functions*, Acta Math. Uppsala, 117 (1967), pp. 37-52.

[35] H. J. LANDAU AND H. O. POLLAK, *Prolate spheroidal wave functions, Fourier analysis and uncertainty*-II, Bell Syst. Tech. J., 40 (1960), pp. 65-84.

[36] H. J. LANDAU AND H. O. POLLAK, *Prolate spheroidal wave functions, Fourier analysis and uncertainty*—III: *The dimension of the space of essentially time- and band-limited signals*, Bell Syst. Tech. J., 41 (1962), pp. 1295-1336.

[37] J. M. LILLY AND J. PARK, *Multiwavelet spectral and polarization analyses of seismic records*, Geophys. J. Int., 122 (1995), pp. 1001-1021.

[38] T.-C. LIU AND B. D. VAN VEEN, *Multiple window based minimum variance spectrum estimation for multidimensional random fields*, IEEE Trans. Signal Process., 40 (1992), pp. 578-589.

[39] S. MALLAT, *A Wavelet Tour of Signal Processing*, Academic Press, San Diego, CA, 1998.

[40] C. MEANEY, *Localization of spherical harmonic expansions*, Monatsh. Math., 98 (1984), pp. 65-74.

[41] A. MESSIAH, *Quantum Mechanics*, Dover, New York, 2000.

[42] L. MIRANIAN, *Slepian functions on the sphere, generalized Gaussian quadrature rule*, Inverse Problems, 20 (2004), pp. 877-892.

[43] F. J. NARCOWICH AND J. D. WARD, *Nonstationary wavelets on the m-sphere for scattered data*, Appl. Comput. Harmon. Anal., 3 (1996), pp. 324-336.

[44] S. OLHEDE AND A. T. WALDEN, *Generalized Morse wavelets*, IEEE Trans. Signal Process., 50 (2002), pp. 2661-2670.

[45] F. W. J. OLVER, *Asymptotics and Special Functions*, A. K. Peters, Wellesley, MA, 1997.

[46] S. M. OULD KABER, *A Legendre pseudospectral viscosity method*, J. Comput. Phys., 128 (1996), pp. 165-180.

[47] R. PAIL, G. PLANK, AND W.-D. SCHUH, *Spatially restricted data distributions on the sphere: The method of orthonormalized functions and applications*, J. Geodesy, 75 (2001), pp. 44-56.

[48] P. J. E. PEEBLES, *Statistical analysis of catalogs of extragalactic objects*. I. *Theory*, Astrophys. J., 185 (1973), pp. 413-440.

[49] D. B. PERCIVAL AND A. T. WALDEN, *Spectral Analysis for Physical Applications, Multitaper and Conventional Univariate Techniques*, Cambridge University Press, New York, 1993.

[50] A. S. POLYAKOV, *Local Basis Expansions for Linear Inverse Problem*, Ph.D. thesis, New York University, 2002.

[51] K. S. RIEDEL AND A. SIDORENKO, *Minimum bias multiple taper spectral estimation*, IEEE Trans. Signal Process., 43 (1995), pp. 188-195.

[52] S. S. SHAPIRO, B. H. HAGER, AND T. H. JORDAN, *The continental tectosphere and Earth's long-wavelength gravity field*, Lithos, 48 (1999), pp. 135-152.

[53] F. J. SIMONS, R. D. VAN DER HILST, AND M. T. ZUBER, *Spatio-spectral localization of isostatic coherence anisotropy in Australia and its relation to seismic anisotropy: Implications for lithospheric deformation*, J. Geophys. Res., 108 (2003), article 2250.

[54] M. SIMONS, S. C. SOLOMON, AND B. H. HAGER, *Localization of gravity and topography: Constraints on the tectonics and mantle dynamics of Venus*, Geophys. J. Int., 131 (1997), pp. 24-44.

[55] D. SLEPIAN, *Estimation of signal parameters in the presence of noise*, Trans. IRE Prof. Gr. Inform. Theory, 3 (1954), pp. 68-89.

[56] D. SLEPIAN, *Prolate spheroidal wave functions, Fourier analysis and uncertainty*—IV: *Extensions to many dimensions; generalized prolate spheroidal functions*, Bell Syst. Tech. J., 43 (1964), pp. 3009-3057.

[57] D. SLEPIAN, *Prolate spheroidal wave functions, Fourier analysis and uncertainty*—V: *The discrete case*, Bell Syst. Tech. J., 57 (1978), pp. 1371-1429.

[58] D. SLEPIAN, *Some comments on Fourier analysis, uncertainty and modeling*, SIAM Rev., 25 (1983), pp. 379-393.

[59] D. SLEPIAN AND H. O. POLLAK, *Prolate spheroidal wave functions, Fourier analysis and uncertainty*—I, Bell Syst. Tech. J., 40 (1960), pp. 43-63.

[60] D. SLEPIAN AND E. SONNENBLICK, *Eigenvalues associated with prolate spheroidal wave functions of zero order*, Bell Syst. Tech. J., 44 (1965), pp. 1745-1759.
[61] N. SNEEUW AND M. V. GELDEREN, *The polar gap*, in *Geodetic Boundary Value Problems in View of the One Centimeter Geoid*, F. Sansò and R. Rummel, eds., Lecture Notes in Earth Sci. 65, Springer, Berlin, 1997, pp. 559–568.

[62] P. N. SWARZTRAUBER AND W. F. SPOTZ, *Generalized discrete spherical harmonic transforms*, J. Comput. Phys., 159 (2000), pp. 213–230.

[63] S. SWENSON AND J. WAHR, *Methods for inferring regional surface-mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time-variable gravity*, J. Geophys. Res., 107 (2002), article 2193.

[64] S. SWENSON, J. WAHR, AND P. C. D. MILLY, *Estimated accuracies of regional water storage variations inferred from the Gravity Recovery and Climate Experiment (GRACE)*, Water Resources Res., 39 (2003), article 1223.

[65] G. SZEGÖ, *Orthogonal Polynomials*, 4th ed., AMS, Providence, RI, 1975.

[66] M. TEGMARK, *A method for extracting maximum resolution power spectra from galaxy surveys*, Astrophys. J., 455 (1995), pp. 429–438.

[67] M. TEGMARK, *The angular power spectrum of the four-year COBE data*, Astrophys. J., 464 (1996), pp. L35–L38.

[68] M. TEGMARK, *A method for extracting maximum resolution power spectra from microwave sky maps*, Mon. Not. R. Astron. Soc, 280 (1996), pp. 299–308.

[69] D. J. THOMSON, *Spectrum estimation and harmonic analysis*, Proc. IEEE, 70 (1982), pp. 1055–1096.

[70] D. J. THOMSON, *Quadratic-inverse spectrum estimates: Applications to paleoclimatology*, Philos. Trans. R. Soc. Lond. Ser. A, 332 (1990), pp. 539–597.

[71] D. J. THOMSON, *Multitaper analysis of nonstationary and nonlinear time series data*, in *Nonlinear and Nonstationary Signal Processing*, R. L. Smith, W. J. Fitzgerald, A. T. Walden, and P. C. Young, eds., Cambridge University Press, Cambridge, UK, 2001, Chap. 11, pp. 317–394.

[72] F. G. TRICOMI, *Integral Equations*, 5th ed., Interscience, New York, 1970.

[73] D. L. TURCOTTE, R. J. WILLEMANN, W. F. HAXBY, AND J. NORBERRY, *Role of membrane stresses in the support of planetary topography*, J. Geophys. Res., 86 (1981), pp. 3951–3959.

[74] A. T. WALDEN, E. J. MCCOY, AND D. B. PERCIVAL, *The variance of multitaper spectrum estimates for real Gaussian processes*, IEEE Trans. Signal Process., 2 (1994), pp. 479–482.

[75] M. A. WIECZOREK AND R. J. PHILLIPS, *Potential anomalies on a sphere: Applications to the thickness of the lunar crust*, J. Geophys. Res., 103 (1998), pp. 1715–1724.

[76] M. A. WIECZOREK AND F. J. SIMONS, *Localized spectral analysis on the sphere*, Geophys. J. Int., 162 (2005), pp. 655–675.