Sparse spectral methods for partial differential equations on spherical caps

In recent years, sparse spectral methods for solving partial differential equations have been derived using hierarchies of classical orthogonal polynomials on intervals, disks, disk-slices and triangles. In this work we extend the methodology to a hierarchy of non-classical multivariate orthogonal polynomials on spherical caps. The entries of discretisations of partial differential operators can be effectively computed using formulae in terms of (non-classical) univariate orthogonal polynomials. We demonstrate the results on partial differential equations involving the spherical Laplacian and biharmonic operators, showing spectral convergence.


Introduction
This paper develops sparse spectral methods for solving linear partial differential equations on certain subsets of the sphere-specifically spherical caps.More precisely, we consider the solution of partial differential equations on the spherical cap Ω defined by where α ∈ (−1, 1) and β := 1.
Remark: For simplicity we focus on the case of a spherical cap, though there is an extension to a spherical band by taking β ∈ (α, 1).The methods presented here translate to the spherical band case by including the necessary adjustments to the weights and recurrence relations we present in this paper.These adjustments make the mathematics more involved, which is why they are omitted here, but the approach is the same.
We advocate using a basis that is polynomial in cartesian coordinates, that is, polynomial in x, y, and z, and orthogonal with respect to a prescribed weight: that is, multivariate orthogonal polynomials, whose construction was considered in [12].Equivalently, we can think of these as polynomials modulo the vanishing ideal {x 2 + y 2 + z 2 = 1}, or simply as a linear recombination of spherical harmonics that are orthogonalised on a subset of the sphere.This is in contrast to more standard approaches based on mapping the geometry to a simpler one (e.g., a rectangle or disk) and using orthogonal polynomials in the mapped coordinates (e.g., a basis that is polynomial in the spherical coordinates ϕ and θ).The benefit of the new approach is that we do not need to resolve Jacobians, and thereby we can achieve sparse discretisations for partial differential operators, including those with polynomial variable coefficients.Further, we avoid the singular nature at the poles or as α approaches 0 that such a projection may give, since our new approach yields a smooth polynomial basis for all α ∈ [−1, 1).
On the spherical cap, the family of weights we consider are of the form W (a) (x, y, z) := (z − α) a , for (x, y, z) ∈ Ω, noting that W (a) (x, y, z) = 0 for (x, y, z) ∈ ∂Ω when a > 0. The corresponding OPs denoted Q (a) n,k,i (x, y, z), where n denotes the polynomial degree, 0 ≤ k ≤ n and i ∈ {0, min(1, k)}.We define these to be orthogonalised lexicographically, that is, n,k,i (x, y, z) = C n,k,i x k−i y i z n−k + (lower order terms) where C n,k,i = 0 and "lower order terms" includes degree n polynomials of the form x j−i y i z n−j where j < k.The precise normalization arises from their definition in terms of one-dimensional OPs in Definition 2.
where ρ(z) := √ 1 − z 2 , we have i.e. ∆ S f (x) = ∆f ( x x ) for x := (x, y, z) ∈ R 3 .We do so by considering the component operators ρ ∂ ∂ϕ and ∂ ∂θ applied to OPs with a specific choices of weight so that their discretisation is sparse, see Theorem 1. Sparsity comes from expanding the domain and range of an operator using different choices of the parameter a, a la the ultraspherical spectral method for intervals [9], triangles [10] and disk-slices and trapeziums [14], and the related work on sparse discretisations on disks [16] and spheres [17,6].As in the disk-slice case in 2D [14], we use an integration-by-parts argument to deduce the sparsity structure.
The three-dimensional orthogonal polynomials defined here involve the same non-classical (in fact, semi-classical [7, §5]) 1D OPs as those outlined for the disk-slice, and so methods for calculating these 1D OP recurrence coefficients and integrals have already been outlined [14].In particular, by exploiting the connection with these 1D OPs we can construct discretizations of general partial differential operators of size (p + 1) 2 × (p + 1) 2 in O(p 3 ) operations, where p is the total polynomial degree.This clearly compares favourably to proceeding in a naïve approach where one would require O(p 6 ) operations.
Note that we consider partial differential operators that are not necessarily rotational invariant: for example, one can use these techniques for Schrödinger operators ∆ S + v(x, y, z) where v is first approximated by a polynomial.A nice feature though is that if the partial differential operator is invariant with respect to rotation around the z axis (e.g., a Schrödinger operator with potential v(z)) the discretisation decouples, and can be reordered as a block-diagonal matrix.This improves the complexity further to an optimal O(p 2 ), which is demonstrated in Figure 4 with v(x, y, z) = cos z.

An overview of the paper is as follows:
Section 2: We present our definition of a (one-parameter) family of 3D orthogonal polynomials (OPs) on the spherical cap domain Ω, by combining 1D OPs on the interval (α, 1) with Chebyshev polynomials, to form 3D OPs on the spherical cap surface.We show that these families will lead to sparse Jacobi operators for multiplication by x, y, z and demonstrate how to obtain the 3D OPs.
Section 3: We define several partial differential operators such as spherical Laplacians and show that these will be sparse when applied to a suitable choice of expansions in bases built from OPs on the spherical cap.We can exactly calculate the non-zero entries of these sparse operators using the quadrature rule associated with the non-classical 1D OPs.Section 4: We derive a quadrature rule on the spherical cap that can be used to expand a function in the OP basis up to a given order N , and demonstrate how to evaluate a function using the Clenshaw algorithm using the coefficients of its expansion.

Orthogonal polynomials on spherical caps
In this section we outline the construction and some basic properties of Q (a) n,k,i (x, y, z).

Explicit construction
We can construct the 3D orthogonal polynomials on Ω from 1D orthogonal polynomials on the interval [α, β], and from Chebyshev polynomials.We do so in terms of Fourier series, which, following [12], we write here as orthogonal polynomials in x and y: , and define the parameter θ for each (x, y) ∈ ω by x = cos θ, y = sin θ.Define the polynomials {Y k,i } for k = 0, 1, . . ., i = 0, 1 on (x, y) ∈ ω by p(x(θ)) q(x(θ)) dθ Note that we have defined Y 0 so as to ensure orthonormality.
For the spherical cap, we can use Proposition 1 to create our one-parameter family of OPs.We first introduce notation for our family of non-classical univariate OPs that will be used as the r n polynomials above.

Definition 2 ([14]
).Let w (a,b) R (x) be a weight function on the interval (α, 1) given by: and define the associated inner product by: where }, orthonormal with respect to the inner product defined in (1).
We can now define the 3D OPs for the spherical cap.
Definition 3. Define the one-parameter 3D orthogonal polynomials via: By construction, {Q n,k,i } are orthogonal with respect to the inner product We note that the weight w (a,b) R (z) has been used in the construction of 2D orthogonal polynomials on disk-slices and trapeziums [14], where a method for obtaining recurrence coefficients and evaluating integrals was established (the weight is in fact semi-classical, and is equivalent to a generalized Jacobi weight [7, §5]).

Jacobi matrices
We can express the three-term recurrences associated with R where the coefficients are calculatable (see [14]).We can use (5) to determine the 3D recurrences for n,k,i (x, y, z).Importantly, we can deduce sparsity in the recurrence relationships.We first require the following lemma.
Lemma 1.The following identities hold for k = 2, 3, . . ., j = 0, 1, . . .and i, h ∈ {0, 1}: Proof.Each follows from the definitions of Y k,i and Y 0 , as well as the relationships: Lemma 2. Define n,k,i (x, y, z) satisfy the following recurrences: n,k,3 n,k,4 n,k,5 if (i = 0 and j is odd) or (i = 1 and j is even) Remark: For z multiplication, note that different Fourier modes do not interact.This is because z is rotationally invariant.

These coefficients are given by
where we show the non-zero coefficients that result are the α . Then for m = 0, . . ., n + 1, j = 0, . . ., m, using a change of variables (cos θ sin ϕ, sin θ sin ϕ, cos ϕ) = (x, y, z): where δ k,j is the standard Kronecker delta function, using Lemma 1.Similarly, for the recurrence for multiplication by y, we can expand m,j,h (x, y, z).
These coefficients are given by where we show the non-zero coefficients that result are the β n,k,6 in the lemma: where again δ k,j is the standard Kronecker delta function, and we have used Lemma 1.
The recurrences in Lemma 2 lead to Jacobi operators that correspond to multiplication by x, y and z.In later sections we will use an ordering of the OPs so that they are grouped by Fourier mode k, which is convenient for the application of differential and other operators to the vector of coefficients of a given function's expansion (some operators will exploit this ordering for operators where Fourier modes do not interact, and thus will be blockdiagonal).Before that though, the ordering we will use in the remainder of this section is convenient for establishing Jacobi operators for multiplication by x, y and z, and hence building the OPs and importantly obtaining the associated recurrence coefficient matrices necessary for efficient function evaluation using the Clenshaw algorithm.In practice, it is simply a matter of converting coefficients between the two orderings.To this end, we define our OP-building ordering as follows.For n = 0, 1, 2, . . .: x , J y , J (a) z as the Jacobi matrices corresponding to where x/y/z,0 A (a) x/y/z,0 x/y/z,1 B (a) x/y/z,1 A (a) x/y/z,1 x/y/z,2 B (a) x/y/z,2 A (a) x/y/z,2 x , J y , J (a) z are banded-block-banded matrices: , and sub-block-bandwidths (λ, µ) if all blocks A i,j are banded with bandwidths (λ, µ).A matrix where the block-bandwidths and sub-block-bandwidths are small compared to the dimensions is referred to as a banded-block-banded matrix.
Note that the sparsity of the Jacobi matrices (in particular the sparsity of the sub-blocks) comes from the natural sparsity of the three-term recurrences of the 1D OPs and the circular harmonics, meaning that the sparsity is not limited to the specific spherical cap, and would extend to the spherical band.

Building the OPs
Following the triangle case [10], we use a multivariate analogue of Clenshaw's algorithm for evaluation, where we combine each system in ( 7) into a block-tridiagonal system, for any (x, y, z) ∈ Ω: , where we note Y 0 , and for each n = 0, 1, 2 . . ., For each n = 0, 1, 2 . . .let D n be any matrix that is a left inverse of A n , i.e. such that D n A n = I 2n+3 .Multiplying our system by the preconditioner matrix that is given by the block diagonal matrix of the D n 's, we obtain a lower triangular system [5, p78], which can be expanded to obtain the recurrence: n−1 (x, y, z), n = 0, 1, 2, . . . .
Note that we can define an explicit D n as follows: for n = 1, 2, . . .where again Γ (a) n,k,3 are defined in equations (8,9) for k = 0, . . ., n, and where η 0 , η 1 ∈ R 3(2n+1) with entries given by For n = 0, we can simply take It follows that we can apply D n in O(n) complexity, and thereby calculate Q

Sparse partial differential operators
In this section we will derive the entries of spherical partial differential operators applied to our basis, demonstrating their sparsity in the process.To this end, as alluded to in Section 2.2, we introduce new notation for a different ordering of the OP vector, in order to exploit the orthogonality the polynomials Y k,i will bring and thus ensure the operators will be block-diagonal.Let N ∈ N and define: 0,0,0 (x, y, z) . . .
N,0,0 (x, y, z) We further denote the weighted set of OPs on Ω by The operator matrices we derive here act on coefficient vectors, that represent a function f (x, y, z) defined on Ω in spectral space -such a function is approximated by its expansion up to degree N : n,k,i (x, y, z), where f = (f n,k,i ) is the coefficients vector for the function f .Definition 6.Let a be a nonnegative parameter, and ã ≥ 2 be a positive integer.Define the operator matrices W according to: W f , (for a = 1 only) The incrementing and decrementing of parameters as seen here is analogous to other well known orthogonal polynomial families' derivatives, for example the Jacobi polynomials on the interval, as seen in the DLMF [8, (18.9. 3)], on the triangle [11], and on the disk-slice [14].The operators we define here are for partial derivatives with respect to the spherical coordinates (ϕ, θ), so that we can more easily apply the operators to PDEs on the surface of a sphere (for example, surface Laplacian operator in the Poisson equation).With the OP ordering by Fourier mode k defined in equations (10,11,12) these rotationally invariant operators are block-diagonal, meaning simple and parallelisable practical application.
Theorem 1.The operator matrices W from Definition 6 are sparse, with banded-block-banded structure.More specifically: ϕ is block-diagonal with sub-block-bandwidths (2, 4) ϕ is block-diagonal with sub-block-bandwidths (4, 2) • L (a)→(a+ã) is block-diagonal with sub-block-bandwidths (0, 4) W is block-diagonal with sub-block-bandwidths (2, 2) In order to show the last part of Theorem 1, we require the following short lemma.Lemma 3.For any general parameter a and any n = 0, 1, . . ., k = 0, . . ., n we have that for some coefficients c{n,k},m .These coefficients are given by c{n,k},m = 1 We show that these are zero for m < n − k − 1 by integrating twice by parts: Proof of Theorem 1.For the operator D θ for partial differentiation by θ, we simply have that We now proceed with the case for the operator D ϕ for partial differentiation by ϕ.The entries of the operator are given by the coefficients in the expansion where the coefficients are . Now, note that: Then, dz which is zero for j = k, h = i, and m < n − 2 by orthogonality.
Similarly for the operator W (a) ϕ for partial differentiation by ϕ on the weighted space, the entries of the operator are given by the coefficients in the expansion ρ ∂ ∂ϕ (w m,j,h , where the coefficients are dz which is zero for j = k, h = i, and m < n − 1 by orthogonality.
We move on to the spherical Laplacian operators.Note that the Laplacian acting on the weighted and non-weighted spherical cap OP Q ∆ S w For the operator L (a)→(a+ã) for the surface Laplacian on a non-weighted space, the entries of the operator are given by the coefficients in the expansion where the coefficients are m,j,h Q (a+ã) .
Using equation ( 13), and integrating by parts twice, we then have that where r m−k+ã is a degree m − k + ã polynomial in z, and so the above is zero for n For the operator L for the surface Laplacian on a weighted space, the entries of the operator are given by the coefficients in the expansion where the coefficients are Using equation ( 14), and integrating by parts thrice, we then have that where r m−k is a degree m − k polynomial in z, and so the above is zero for n Finally, fix a = 1.For the operator ∆ W for the Laplacian on the weighted space, the entries of the operator are given by the coefficients in the expansion ∆ S w m,j,h , where the coefficients are given by m,j,h Q (1) .
Using equation ( 14) with a = 1, and Lemma 3, we then have that By applying these differential operators, we are (in some cases) incrementing or decrementing the parameter value a.It is therefore necessary to also be able to raise or lower the parameter by way of an independent operator.There exist conversion matrix operators that do exactly this, transforming the OPs from one (weighted or non-weighted) parameter space to another.Definition 7. Define the operator matrices T (a)→(a+ã) , T (a)→(a−ã) W for conversion between non-weighted spaces and weighted spaces respectively according to Lemma 4. The operator matrices in Definition 7 are sparse, with banded-block-banded structure.More specifically: • T (a)→(a+ã) is block-diagonal with sub-block bandwidths (0, 2ã) is block-diagonal with sub-block bandwidths (2ã, 0) Proof.We proceed with the case for the non-weighted operators T (a)→(a+ã) .Since {Q (a+ã) m,j,h } for m = 0, . . ., n, j = 0, . . ., m, h = 0, 1 is an orthogonal basis for any degree n polynomial, we can expand m,j,h .The coefficients of the expansion are then the entries of the operator matrix.We will show that the only non-zero coefficients are for k = j, i = h and m ≥ n − ã.Note that where The sparsity argument for the weighted parameter transformation operator follows similarly.

Further partial differential operators
General linear partial differential operators with polynomial variable coefficients can be constructed by composing the sparse representations for partial derivatives, conversion between bases, and Jacobi operators.As a canonical example, we can obtain the matrix operator for the ρ 2 -factored spherical Laplacian ρ(z) 2 ∆ S , that will take us from coefficients for expansion in the weighted space W (1) N (x, y, z) to coefficients in the non-weighted space Q(1) N (x, y, z).Note that this construction will ensure the imposition of the Dirichlet zero boundary conditions on Ω, similar to how the Dirichlet zero boundary conditions would be imposed for the operator ∆ W in Definition 6.The matrix operator for this ρ 2 -factored spherical Laplacian acting on the coefficients vector is then given by Importantly, this operator will have banded-block-banded structure, and hence will be sparse, as seen in Figure 1.
Another desirable operator is the Biharmonic operator ∆ 2 S , for which we assume zero Dirichlet and Neumann conditions.That is, The Laplace-Beltrami operator ∆ ) T where ∂Ω is the z = α boundary, and n(x, y, z) is the outward unit normal vector at the point (x, y, z) on the boundary, i.e. n(x, y, z) = n( x) := x x = x.The matrix operator for the Biharmonic operator will take us from coefficients in the space W (2) (x, y, z) to coefficients in the space Q(2) N (x, y, z).To construct this, we can simply multiply together two of the spherical Laplacian operators defined in Definition 6, namely L (0)→ (2)  and L (2)→(0) W :

Since the operator L
(2)→(0) W acts on coefficients in the W (2) (x, y, z) space, we ensure that we satisfy the zero Dirichlet and Neumann boundary conditions -such a function could be written u(x, y, z) = w (2,0) R (z) ũ(x, y, z) and thus its spherical gradient would be zero on the boundary z = α.This allows us to apply the L (0)→ (2) operator after, safe in the knowledge that boundary conditions have been accounted for.The sparsity and structure of this biharmonic operator are seen in Figure 1.

Computational aspects
In this section we discuss how to expand and evaluate functions in our proposed basis, and take advantage of the sparsity structure in partial differential operators in practical computational applications.} OPs in (5), see [14], by careful application of the Christoffel-Darboux formula [8, 18.2.12].

Quadrature rule on the spherical cap
In this section we construct a quadrature rule exact for polynomials on the spherical cap Ω that can be used to expand functions in the OPs Q (a) n,k,i (x, y, z) for a given parameter a.

Remark:
Note that the Gauss quadrature nodes and weights (t j , w (t) j ) will have to be calculated, however the Gauss quadrature nodes and weights (s j , w (s) j ) are simply the Chebyshev-Gauss quadrature nodes and weights given explicitly [8, 3.5.23]as s j := Proof.Let f : Ω → R. Define the functions f e , f o : Ω → R by so that x → f e (x, z) for fixed z is an even function, and x → f o (x, z) for fixed z is an odd function.Note that if f is a polynomial, then f e (ρ(t)x, ρ(t)y, t) is a polynomial in t ∈ [α, 1] for fixed (x, y) ∈ R 2 .
Firstly, we note that for some function g.Then, integrating the even function f e we have Ω f e (x, y, z) w l f e ρ(z)s l , ρ(z) 1 − s 2 l , z dz ( ) l f e ρ(t j )s l , ρ(t j ) 1 − s 2 l , t j ( ) w j f e (x j , y j , z j ).
Suppose f is a polynomial in x, y, z of degree N , and hence that f e is a degree ≤ N polynomial.It follows that s → f e ρ(z)s, ρ(z) √ 1 − s 2 , z for fixed z is then a polynomial of degree ≤ N .We therefore achieve equality at ( ) if 2M 2 −1 ≥ N and we achieve equality at ( ) if also 2M 1 − 1 ≥ N .
Recall from equation ( 4) that Using the quadrature rule detailed in Section 4.2 for the inner product, we can calculate the coefficients f n,k,i for each n = 0, . . ., N , k = 0, . . ., n, i = 0, 1: n,k,i (x j , y j , z j ) + f (−x j , −y j , z j )Q ).

Function evaluation
For a function f , with coefficients vector f for expansion in the {Q n,k,i } basis as determined via the method in Section 4.3 up to order N , we can use the Clenshaw algorithm to evaluate the function at a point (x, y, z) ∈ Ω as follows.Let A n , B n , D n , C n be the Clenshaw matrices from Definition 5, and define the rearranged coefficients vector f via The trivariate Clenshaw algorithm works similar to the bivariate Clenshaw algorithm introduced in [10] for expansions in the triangle: 2) For n = N : −1 : 0

Calculating non-zero entries of the operator matrices
The proofs of Theorem 1 and Lemma 4 provide a way to calculate the non-zero entries of the operator matrices given in Definition 6 and Definition 7. We can simply use quadrature to calculate the 1D inner products, which has a complexity of O(N 2 ).This proves much cheaper computationally than using the 3D quadrature rule to calculate the surface inner products, which has a complexity of O(N 3 ).

Obtaining operator matrices for variable coefficients
The Clenshaw algorithm outlined in Section 4.4 can also be used with Jacobi matrices J (a) x , J y , J (a) z replacing the point (x, y, z).Let v : Ω → R be the function that we wish to obtain an operator matrix V for v, so that i.e.V f is the coefficients vector for the function v(x, y, z) f (x, y, z).
To this end, let ṽ be the coefficients for expansion up to order N in the {Q n,k,i } basis of v (rearranged as in Section 4.4 so that v(x, y, z) = v Q (a) (x, y, z)).Denote X := (J (a) x ) , Y := (J (a) y ) , Z := (J (a) z ) .The operator V is then the result of the following: where at each iteration, ξ n is a vector of matrices.

Examples on spherical caps with zero Dirichlet conditions
We now demonstrate how the sparse linear systems constructed as above can be used to efficiently solve PDEs with zero Dirichlet conditions on the spherical cap defined by Ω.We  right hand sides we choose here are given by for differing choices of -this parameter serves to alter the distance from which we would have a singularity.In the plot, a "block" is simply the group of coefficients corresponding to OPs of the same degree, and so the plot shows how the norms of these blocks decay as the degree of the expansion increases.Thus, the rate of decay in the coefficients is a proxy for the rate of convergence of the computed solution: as typical of spectral methods, we expect the numerical scheme to converge at the same rate as the coefficients decay.We see that we achieve spectral convergence for these examples.

Inhomogeneous variable-coefficient Helmholtz
Find u(x, y) given functions v, f : Ω → R such that: where k ∈ R, noting the imposition of zero Dirichlet boundary conditions on u.
We can tackle the problem as follows.Denote the coefficient vector for expansion of u in the W N OP basis up to degree N by u, and the coefficient vector for expansion of f in the Q(1) N OP basis up to degree N by f .Since f is known, we can obtain the coefficients f using the quadrature rule in Section 4.3.
Figure 4: Time in seconds to build and solve the system ∆S + v(x, y, z) u(x, y, z) = f (x, y, z), for a rotationally invariant v(x, y, z) = v(z).This demonstrates that the approach is roughly of order O(N 2 ), where N is the degree to which we approximate the solution.Here, we used f = −2e x yz(2 + x) + (z − α)e x (y 3 + z 2 y − 4xy − 2y) and v(x, y, z) = v(z) = cos(z).
Define X := (J (0) x ) , Y := (J (0) y ) , Z := (J (0) z ) .We can obtain the matrix operator for the variable-coefficient function v(x, y, z) by using the Clenshaw algorithm with matrix inputs as the Jacobi matrices X, Y, Z, yielding an operator matrix of the same dimension as the input Jacobi matrices a la the procedure introduced in [10].We can denote the resulting operator acting on coefficients in the Q(0) N space by v(X, Y, Z).In matrix-vector notation, our system hence becomes: (∆ W + k 2 T (0)→(1) V T (1)→(0) W ) u = f which can be solved to find u.We can see the sparsity and structure of this matrix system in Figure 1 with v(x, y, z) = zxy 2 as an example.In Figure 3 we see the solution to the inhomogeneous variable-coefficient Helmholtz equation with zero boundary conditions given in (16) in the spherical cap Ω, with f (x, y, z) = ye x w (1,0) R (z), v(x, y, z) = 1 − (3(x − x 0 ) 2 + 5(y − y 0 ) 2 + 2(z − z 0 ) 2 ) where (x 0 , z 0 ) := (0.7, 0.2), y 0 := 1 − x 2 0 − z 2 0 and k = 100.In Figure 3 we also show the norms of each block of calculated coefficients for the approximation of the solution to the inhomogeneous variable-coefficient Helmholtz equation with various k values.Here, we use N = 200, that is, (N + 1) 2 = 40, 401 unknowns.Once again, the rate of decay in the coefficients is a proxy for the rate of convergence of the computed solution, and we see that we achieve spectral convergence.
In Figure 4 we plot the time taken 1 to construct the operator for ∆ S + v(x, y, z), with a rotationally invariant v(x, y, z) = v(z) = cos z, and solve a zero boundary condition Helmholtz problem.The plot demonstrates that as we increase the degree of approximation N , we achieve a complexity of an optimal O(N 2 ).

Biharmonic equation
Our last erxample is the biharmonic equation: find u(x, y, z) given a function f (x, y, z) such that: ∆ 2 S u(x, y, z) = f (x, y, z) in Ω u(x, y, z) = 0, ∂u ∂n (x, y, z) = ∇ S u(x, y, z) • n(x, y, z) = 0 on ∂Ω where ∆ 2 S is the Biharmonic operator, noting the imposition of zero Dirichlet and Neumann boundary conditions on u.For clarity, we reiterate that the unit normal vector in this sense is simply n(x, y, z) = n( x) := x x = x (see Section 3.1).In Figure 5 we see the solution to the Biharmonic equation (17) in the spherical cap Ω.In Figure 5 we also show the norms of each block of calculated coefficients of the approximation for four more complex right-hand sides of the biharmonic equation with N = 200, that is, (N + 1) 2 = 40, 401 unknowns.Once again, the rate of decay in the coefficients is a proxy for the rate of convergence of the computed solution, and we see that we achieve exponential convergence for these more complex functions.

Conclusions
We have shown that trivariate orthogonal polynomials can lead to sparse discretizations of general linear PDEs on spherical cap domains, with Dirichlet boundary conditions on the z = α ∈ (0, 1) boundary.We have provided a detailed practical framework for the application of the methods described for quadratic surfaces of revolution [12], by utilising the non-classical 1D OPs on the interval [α, 1] with the weight (z −α) a (1−z 2 ) b/2 defined for the disk-slice case [14].Generalisation to spherical bands (α ≤ z ≤ β) is straightforward.This work thus forms a building block in developing an hp−finite element method to solve PDEs on the sphere by using spherical band and spherical cap shaped elements.This work also serves as a stepping stone to constructing similar methods to solve partial differential equations on other 3D sub-domains of the sphere-it is clear from the construction in this paper that discretizations of spherical gradients and Laplacian's are sparse on other suitable sub-components of the sphere.The resulting sparsity in high-polynomial degree discretizations presents an attractive alternative to methods based on bijective mappings (e.g., [2,13,3]).Constructing these sparse spectral methods for surface PDEs on spherical triangles is future work, and has applications in weather prediction [15], though it is not yet clear how to directly construct the necessary orthogonal polynomials.
The next stage is to develop an orthogonal basis for the tangent space of the spherical cap (or band), and obtain sparse differential operators for gradient, divergence etc.On the complete sphere, the vector spherical harmonics that form the orthogonal basis are simply the gradients and perpendicular gradients of the scalar spherical harmonics [1] which has been used effectively for solving PDEs on the sphere [17,6] -however, we do not have that luxury for the spherical cap or band, and hence the choice of basis will not be as straightforward.

Section 5 :
We demonstrate the proposed technique for solving differential equations on the spherical cap such as the Poisson equation, variable coefficient Helmholtz equation, and Biharmonic equation.2019-144 "Constructive approximation theory on and inside algebraic curves and surfaces".

where Y 0 := √ 2 2
and T k , U k−1 are the standard Chebyshev polynomials on the interval [−1, 1].The {Y k,i } are orthonormal with respect to the inner product p, q Y := 1 π 2π 0

)
is a normalising constant.Denote the two-parameter family of orthonormal polynomials on [α, β] by {R (a,b) n , y, z) through Q (a) n (x, y, z) in optimal O(n 2 ) complexity.Definition 5.The recurrence coefficient matrices associated with the OPs {Q (a) n,k,i } are given by the matrices A n , B n , C n , D n for n = 0, 1, 2, . . .defined above.
It is possible to recursively obtain the recurrence coefficients for the {R (a,b) n ,i (−x j , −y j , z j ) where the quadrature nodes and weights are those from Theorem 2, andM = M 1 M 2 with 2M 1 − 1 ≥ N, M 2 − 1 ≥ N (i.e.we can choose M 2 := N + 1 and M 1 := N +1 2