Fundamental Solutions and Gegenbauer Expansions of Helmholtz Operators in Riemannian Spaces of Constant Curvature

We perform global and local analysis of oscillatory and damped spherically symmetric fundamental solutions for Helmholtz operators $\big({-}\Delta\pm\beta^2\big)$ in $d$-dimensional, $R$-radius hyperbolic ${\mathbf H}_R^d$ and hyperspherical ${\mathbf S}_R^d$ geometry, which represent Riemannian manifolds with positive constant and negative constant sectional curvature respectively. In particular, we compute closed-form expressions for fundamental solutions of $\big({-}\Delta \pm \beta^2\big)$ on ${\mathbf H}_R^d$, $\big({-}\Delta+\beta^2\big)$ on ${\mathbf S}_R^d$, and present two candidate fundamental solutions for $\big({-}\Delta-\beta^2\big)$ on ${\mathbf S}_R^d$. Flat-space limits, with their corresponding asymptotic representations, are used to restrict proportionality constants for these fundamental solutions. In order to accomplish this, we summarize and derive new large degree asymptotics for associated Legendre and Ferrers functions of the first and second kind. Furthermore, we prove that our fundamental solutions on the hyperboloid are unique due to their decay at infinity. To derive Gegenbauer polynomial expansions of our fundamental solutions for Helmholtz operators on hyperspheres and hyperboloids, we derive a collection of infinite series addition theorems for Ferrers and associated Legendre functions which are generalizations and extensions of the addition theorem for Gegenbauer polynomials. Using these addition theorems, in geodesic polar coordinates for dimensions greater than or equal to three, we compute Gegenbauer polynomial expansions for these fundamental solutions, and azimuthal Fourier expansions in two-dimensions.


Introduction
In this paper we derive associated Legendre and Ferrers function expressions for fundamental solutions of Helmholtz operators −∆ ± β 2 , β 2 > 0, in Riemannian spaces of constant cur-This paper is a contribution to the Special Issue on Orthogonal Polynomials, Special Functions and Applications (OPSFA14). The full collection is available at https://www.emis.de/journals/SIGMA/OPSFA2017.html vature, namely in the d-dimensional R-radius hyperboloid H d R and hyperspherical S d R models with negative and positive sectional curvatures respectively, where R > 0. We also compute eigenfunction expansions for these fundamental solutions of Helmholtz operators in geodesic polar coordinates in these Riemannian manifolds. In particular, we derive Gegenbauer polynomial expansions for spherically symmetric fundamental solutions of Helmholtz operators on Riemannian manifolds of negative-constant and positive-constant sectional curvatures. Useful background material relevant for this paper can be found in [19,28,32,35].
This paper is organized as follows. In Section 2, we introduce and give some useful properties for the special functions and orthogonal polynomials that we will use in this paper. In particular, we summarize (and derive new) large degree asymptotics for associated Legendre and Ferrers functions of the first and second kind. We also derive a collection of infinite series addition theorems for Ferrers and associated Legendre functions (generalizations and extensions of the addition theorem for Gegenbauer polynomials) which will be used to compute Gegenbauer polynomial expansions for fundamental solutions of Helmholtz operator in spaces of constant curvature. In Section 3, for the hyperboloid model of d-dimensional hyperbolic geometry and for hyperspherical geometry, we describe some of their global properties, such as their respective geodesic distance functions, geodesic polar coordinates, Helmholtz operators, and their corresponding radial harmonics. In Section 4, we show how to compute radial harmonics in a geodesic polar coordinate system and derive fundamental solutions for −∆ ± β 2 on H d R , −∆ + β 2 on S d R , and study two candidate fundamental solutions for −∆ − β 2 on S d R . In Section 5, for d ≥ 3, we compute Gegenbauer polynomial expansions in geodesic polar coordinates for fundamental solutions of Helmholtz operators on the hyperboloid and hypersphere. We also compute azimuthal Fourier expansions for these fundamental solutions in two-dimensions. In Appendix A, the proof of an addition theorem is presented.
2 Special functions, asymptotics, and notation The set R represents the real numbers and the set C represents the complex numbers. For a 1 , a 2 , a 3 , . . . ∈ C, if i, j ∈ Z and j < i then j n=i a n := 0 and j n=i a n := 1. Note that we often adopt a common notation used for fundamental solution expansions, namely if one takes a, a ∈ R, then a ≶ := min max {a, a }.

The gamma function and factorials
The (Euler) gamma function Γ : C \ −N 0 → C (see [27,Chapter 5]), which is ubiquitous in special function theory satisfies the recurrence formula Γ(z + 1) = zΓ(z), and is an important combinatoric function which generalizes the factorial function Γ(n + 1) = n!, n ∈ N 0 . The gamma function is naturally defined through Euler's integral [27, equation (5.2.1)]. The following asymptotic approximation involving the ratio of gamma functions will also be needed. where τ a−b takes its principal value.

The associated Legendre functions of the first and second kind
We also frequently use associated Legendre functions of the first and second kind P µ ν , Q µ ν : C \ (−∞, 1] → C respectively. Gauss hypergeometric representations of these functions are [27, equation (14.3.6)] where |1 − z| < 2 and [27, equation (14.3.7)] where |z| > 1. In regard to the above definition of the associated Legendre function of the first kind, note that for µ ∈ N, Γ(1 − µ) is undefined, however the function 1 Γ(c) 2 F 1 a,b c ; z is an entire function for all a, b, c ∈ C, |z| < 1. See [27, equation (15.2.2)] and the discussion in [27,Section 14.3(ii)] for µ ∈ Z.
We shall give analogous results for associated Legendre conical functions. In order to do so, we shall require so-called envelope functions for our approximants, the Bessel and Hankel functions. The reason for this is they have zeros in our domain of validity, unlike the modified Bessel functions used in (2.12), (2.13). For the Bessel function J µ (x) (µ ≥ 0) of real argument x, an envelope function env J µ (x) is given in [27, equations (2.8.32)- (2.8.34)]. This function is continuous and nonvanishing, and has the properties that does not approach zero as x → ∞. Thus env J µ (x) has the same order of magnitude as J µ (x) uniformly for x ∈ (0, ∞), but does not have any zeros (except at x = 0 when µ > 0).
We also require envelope functions for the Hankel functions H (1) µ (z), H (2) µ (z). For positive argument these are not required, since these functions have no real zeros. However they do have complex zeros (see [27,Section 10.21(ix)]), and for this reason we shall construct them for complex z such that | arg z| ≤ 1 2 π. To this end we note the limiting behaviors which hold for | arg z| ≤ 1 2 π, namely as z → 0, (2.16) and as z → ∞, (2.17) Now for H (1) µ (z) the natural choice analogous to (2.14) would appear to be env 2 , but on referring to (2.16) we see that this function would have the A similar problem occurs for such an envelope function for H (2) µ (z). With this in mind, in the following lemma we make the required modification in our definition of the Hankel envelope functions. Note that we omit the proof, which is similar to the proof of Lemma 2.3.  has the following properties in the half-plane | arg z| ≤ 1 2 π: (i) env H (1,2) µ (z) has no zeros; We now state large degree asymptotic approximations for the associated Legendre conical functions.
Proof . The proof for the expansion of the associated Legendre function of the second kind is identical to that given in [6, Proposition 5.1], except one does not specialize to ν = µ. For the associated Legendre function of the first kind, instead start with [10, equation (8.3)].
We next define an odd Ferrers function f µ ν : (−1, 1) → C, by the difference We note that (2.48) vanishes if P µ ν (x) is even. This occurs if ν + µ is a non-negative even integer, or if ν −µ is a negative odd integer (see [26,Chapter 5,p. 187]). Similarly, the Ferrers function of the second kind difference (2.49) vanishes if ν +µ is a positive odd integer, or if ν −µ is a negative odd integer. For our purposes (2.48) will suffice, since in our applications this function will not vanish identically. (2.50) Proof . Let I denote the integral in (2.50), with positive sign in the argument of the Ferrers function. Insert the hypergeometric representation of the Ferrers function of the first kind (2.35), and make the substitution ξ = (1 − x)/2. This converts I to a Mellin transform, namely after re-mapping ξ → x, Using [12, equation (7.512.5)], we can evaluate the integral in terms of a 3 F 2 (1), namely uniformly for θ ∈ 1 2 π, π . Proof . For (2.58), combine (2.52), (2.56). For (2.59) replace θ by π − θ in (2.58) and observe Similar approximations for the Ferrers conical functions read as follows.
Then provided the Ferrers functions of the second kind are defined, where in all but (2.76) θ = θ . For µ = 0, see Corollary 2.23.

Proof . See Appendix A.
Remark 2.20. If one is interested in the above addition theorems which correspond to the Ferrers functions with ν = µ, then one can directly use (2.74), (2.75). If one examines the behavior of (2.76), when ν = µ, it is seen that only the n = 0 term in the sum survives. However, one can see through [27, equation (14.4.18)], that the left-hand side exactly matches the n = 0 term on the right-hand side. Note that the left-hand side of (2.77) is defined for ν = µ, but (except for n = 0) the terms in the series on the right-hand side are undefined in this case. However, the limit of the series exists as ν → µ, and is given by Corollary 2.21 below. Similar limiting results can be obtained if ν − µ ∈ N but we do not pursue this.
The addition theorems for the Ferrers function of the second kind with ν = −µ can be obtained in two different ways. These Ferrers functions appear in the study of a fundamental solution of the Laplace-Beltrami operator on the d-dimensional R-radius hypersphere. See (4.11) below and [3,5,7].
where the first equality is valid except when µ is a half odd integer, and the second equality is valid except when µ is an integer.
Proof . First start with (2.75), let ν = µ, and apply the connection relation ( In order to obtain the addition theorems for the Ferrers function of the second kind when the order is equal to the negative degree, and is either an integer or a half odd integer, one can then use the above corollary. The following formulae, specializations of Corollary 2.21, are exactly those results which are derived in [7, Theorems 3.1 and 4.1] for a fundamental solution of the Laplace-Beltrami operator on the d-dimensional R-radius hypersphere. Proof . First substitute µ = k ∈ N in (2.80) which produces (2.83). Next, substituting µ = m + 1 2 , m ∈ N 0 , in (2.81), to produce (2.84). Then take the limit as µ → 0 in (2.80) using (2.4) which produces (2.82).

Global analysis on Riemannian manifolds of constant curvature
In this paper, when we refer to a fundamental solution, it is meant to be a fundamental solution of a partial differential operator on a Riemannian manifold M . A fundamental solution for a linear partial differential operator E x on a d-dimensional Riemannian manifold M , is in general x ∈ M } → C which satisfies the following linear partial differential equation where x, x ∈ M , g is the Riemannian structure on M , and δ g is the Dirac delta distribution on M . Note that for the operators that are treated in this paper, a fundamental solution will always be function, but for other operators such as for wave operators (see, e.g., [18, equation (4.1.5)]), fundamental solutions are in general distributions. Furthermore, we will only perform global and local analysis for Riemannian manifolds M with constant curvature (hyperspherical S d R , hyperbolic H d R , and Euclidean geometry R d ), and we will only study linear elliptic partial differential operators of the form −∆ ± β 2 , where −∆ is the positive Laplace-Beltrami operator on M and β 2 ≥ 0, R > 0. In this paper, these operators are referred to as Helmholtz operators, and in the case where β = 0, they are Laplacian operators.
In this section we develop the necessary material in order to study fundamental solutions for these operators on these manifolds.

Hyperspherical geometry and the hyperboloid model of hyperbolic geometry
Hyperbolic space in d-dimensions is a fundamental example of a space exhibiting hyperbolic geometry. It was developed independently by Lobachevsky and Bolyai around 1830 (see [34]), and most likely by Gauss and Schweikart (although they never published this result), even earlier (see [21,Chapter 6]). It is a geometry analogous to Euclidean geometry, but such that Euclid's parallel postulate is no longer assumed to hold. There are several models of d-dimensional hyperbolic space including the Klein, Poincaré, hyperboloid, upper-half space and hemisphere models (see [32]). The hyperboloid model for d-dimensional hyperbolic geometry is closely related to the Klein and Poincaré models: each can be obtained projectively from the others. The upper-half space and hemisphere models can be obtained from one another by inversions with the Poincaré model (see [32,Section 2.2]). The model we will be focusing on in this paper is the hyperboloid model.
The hyperboloid model, also known as the Minkowski or Lorentz model, is a model of ddimensional hyperbolic geometry in which points are represented by the upper sheet (submanifold) S + of a two-sheeted hyperboloid embedded in the Minkowski space R d,1 . Minkowski space is a (d + 1)-dimensional pseudo-Riemannian manifold which is a real finite-dimensional vector space, with coordinates given by x = (x 0 , x 1 , . . . , x d ). It is equipped with a nondegenerate, symmetric bilinear form, the Minkowski bilinear form [·, ·] : The above bilinear form is symmetric, but not positive-definite, so it is not an inner product. It is defined analogously with the Euclidean inner product (·, ·) : 1 , using the language of [2] (see also [35, p. 504]), defines a pseudo-sphere of radius R. Points on the pseudo-sphere with zero radius coincide with the cone. Points on the pseudo-sphere with radius greater than zero lie within this cone, and points on the pseudo-sphere with purely imaginary radius lie outside the cone. The upper sheets of the positive radii pseudo-spheres are maximally symmetric, simply connected, negative-constant sectional curvature (given by −1/R 2 , see for instance [19, p. 148]), d-dimensional Riemannian submanifolds, embedded and with induced metric from the ambient Minkowski space R d,1 . For R ∈ (0, ∞), we refer to the upper sheet of this variety [ Similarly, we refer to the variety (x, x) = R 2 for R > 0 and x ∈ R d+1 , as the R-radius hypersphere S d R which is a maximally symmetric, simply connected, positive-constant sectional curvature (given by 1/R 2 ) d-dimensional Riemannian submanifold, embedded and with induced metric from the ambient Euclidean space. The Euclidean space R d equipped with the Pythagorean norm, is a space with vanishing curvature. We denote the unit radius hyperboloid by H d := H d 1 and the unit radius hypersphere by S d := S d 1 . In our discussion of a fundamental solution for the Helmholtz operator in the hyperboloid model of hyperbolic geometry, we focus on the positive radius pseudo-sphere which can be parameterized through subgroup-type coordinates, i.e., those which correspond to a maximal subgroup chain O(d, 1) ⊃ · · · (see for instance [28]). There exist separable coordinate systems which parameterize points on positive radius pseudo-spheres which can not be constructed using maximal subgroup chains, e.g., such as those which are analogous to parabolic coordinates, etc. We will no longer discuss these.
In order to do analysis on a fundamental solution of the Helmholtz equation on the hyperboloid and hypersphere, we need to describe how one computes distances in these spaces. One may naturally compare distances on the positive radius pseudo-sphere through analogy with the R-radius hypersphere. Distances on the hypersphere are simply given by arc lengths, angles between two arbitrary vectors, from the origin, in the ambient Euclidean space. We consider the d-dimensional hypersphere embedded in R d+1 . Points on the hypersphere can be parameterized using hyperspherical coordinate systems. Any parameterization of the hypersphere S d R , must have (x, x) = x 2 0 + · · · + x 2 d = R 2 , with R > 0. The geodesic distance between two points on the hypersphere d s : where the inverse hyperbolic cosine with argument x ∈ (1, ∞) is given by [ where x = x/R and x = x /R.

The Helmholtz equation in Riemannian spaces of constant curvature
Parametrizations of a submanifold embedded in either a Euclidean or Minkowski space are given in terms of coordinate systems whose coordinates are curvilinear. These are coordinates based on some transformation that converts the standard Cartesian coordinates in the ambient space to a coordinate system with the same number of coordinates as the dimension of the submanifold in which the coordinate lines are curved. On a d-dimensional Riemannian manifold M (a manifold together with a Riemannian metric g), the Laplace-Beltrami operator (Laplacian) ∆ : C p (M ) → C p−2 (M ), p ≥ 2, in curvilinear coordinates ξ = ξ 1 , . . . , ξ d is given by where |g| = | det(g ij )|, the infinitesimal distance is given by For a Riemannian submanifold, the relation between the metric tensor in the ambient space and g ij of (3.7), (3.8) is On H d R the ambient space is Minkowski, and therefore G = diag(1, −1, . . . , −1). On S d R , G = diag (1, 1, . . . , 1).
In any of the geodesic polar coordinate systems, the global geodesic distance between any two points on the hyperboloid and hypersphere are given by (cf. (3.4), (3.6)) where γ is the unique separation angle defined on the unit radius submanifold S d−1 which is embedded in both d-dimensional Riemannian manifolds H d R and S d R for a fixed geodesic radius. For instance, one may write the separation angle in standard geodesic polar coordinates (3.1) or (3.2) as follows sin θ j sin θ j . (3.11) Corresponding separation angle formulae for any geodesic polar coordinate system can be computed using (3.3), (3.4), and the associated formulae for the appropriate inner-products. Note that by making use of the isometry group SO(d, 1) to map x to the origin, then the geodesic distance ρ as measured from the origin to a point x ∈ H d R with curvilinear coordinate r is ρ = Rr. Hence if R = 1 which corresponds to H d (the unit radius hyperboloid), then on this Riemannian manifold, the geodesic distance is ρ = r, and there is no distinction between the global geodesic distance and the r-parameter in a geodesic polar coordinate system.
The infinitesimal distance in a geodesic polar coordinate system on the submanifold H d R is given by and on S d R by where an appropriate expression for γ in a curvilinear coordinate system is given by (3.11). If one combines (3.1) or (3.2), (3.7), (3.11) and (3.12) or (3.13), then in a particular geodesic polar coordinate system, the Helmholtz equation on H d R is given by and on S d R by where ∆ S d−1 is the corresponding Laplace-Beltrami operator on the unit radius hypersphere S d−1 .

Homogeneous solutions of the Helmholtz equation in geodesic polar coordinates
Geodesic polar coordinate systems partition H d R into a family of (d − 1)-dimensional hyperbolic hyperspheres, each with a geodesic radius Rr with r ∈ (0, ∞) on which all possible hyperspherical coordinate systems for S d−1 may be used (see for instance [35]). One then must also consider the limiting case for r = 0 to fill out all of H d R . In subgroup-type coordinate systems, one can compute the normalized hyperspherical harmonics in that space by solving the Laplace equation using separation of variables. This results in a general procedure which is given explicitly in [15,16]. These angular harmonics are given as general expressions involving trigonometric functions, Gegenbauer polynomials and Jacobi polynomials.
The harmonics in geodesic polar coordinate systems are given in terms of a radial solution multiplied by the angular harmonics. The angular harmonics are eigenfunctions of the Laplace-Beltrami operator on S d−1 which satisfy the following eigenvalue problem are normalized hyperspherical harmonics, l ∈ N 0 is the angular momentum quantum number, and K stands for the set of (d − 2)-quantum numbers identifying degenerate harmonics for each l. The degeneracy (2l + d − 2)(d − 3 + l)!/(l!(d − 2)!) (see [35, equation (9.2.11)]), tells you how many linearly independent solutions exist for a particular l value and dimension d. The hyperspherical harmonics can optionally be normalized so that where dω is the Riemannian (volume) measure (see for instance [13, Section 3.4]) on S d−1 which is invariant under the isometry group SO(d) (cf. (4.7)), and for x + iy = z ∈ C, z = x − iy, represents complex conjugation. The generalized Kronecker delta δ K K ∈ {0, 1} is defined such that it equals 1 if all of the (d − 2)-quantum numbers identifying degenerate harmonics for each l coincide, and equals zero otherwise.

A fundamental solution of −∆ ± β 2 in Riemannian spaces of constant curvature
Due to the fact that the spaces H d R and S d R are homogeneous with respect to their isometry groups, the pseudo-orthogonal group SO(d, 1) of H d R and the orthogonal group O(d) of S d R , and therefore isotropic manifolds, we expect that there exists a fundamental solution of the Helmholtz equation on each space with spherically symmetric dependence. We specifically expect these solutions to be given on H d R in terms of the associated Legendre function of the second kind with argument given by cosh r. This associated Legendre function naturally fits our requirements because it is singular at r = 0 and vanishes at infinity, whereas the associated Legendre function of the first kind, with the same argument, is regular at r = 0 and singular at infinity. One also might expect a fundamental solution of the Helmholtz equation on hyperspheres to be a Ferrers function of the second kind with argument given by cos θ, since this is the result we have found previously for Laplace's equation (see (4.11) below). However on hyperspheres for the Helmholtz equation, there is a twist, as we will see in the next subsection.

Properties of fundamental solutions of inhomogeneous
where h and s are Riemannian metrics on H d R and S d R respectively, and δ h (x, x ) and δ s (x, x ) are the corresponding Dirac delta distributions on those Riemannian manifolds. Applying the divergence theorem to solutions of the inhomogeneous Helmholtz equations Hence for the Helmholtz equation, the integral of the density distribution over the entire manifold is not required to vanish. Instead, it must be equal to the total integral of the solution multiplied by the wavenumber ±β 2 . Note that in the case where the total integral of the density distribution vanishes, this implies that the total integral of the solution must also vanish. This is the case for a fundamental solution of Laplace's equation on S d R , in which case a fundamental solution must be an odd function of the geodesic distance as measured from the equator θ = 1 2 π of the hypersphere S d R . We will also study odd solutions for Helmholtz operators on these manifolds and compare the β → 0 + limit for Laplace operators for the opposite antipodal fundamental solutions whose total integral over the manifold must vanish. Another important requirement for a fundamental solution of the −∆ − β 2 Helmholtz operator in Euclidean space is the Sommerfeld radiation condition. Sommerfeld defined this radiation condition for a scalar field satisfying the Helmholtz equation as [30, p. 189] (see also [29, p. 396]) ". . . the sources must be sources, not sinks of energy. The energy which is radiated from the sources must scatter to infinity; no energy may be radiated from infinity into the prescribed singularities of the field." On R d , this radiation condition is expressed mathematically as satisfying this condition will be given in Theorem 4.2 below.

A fundamental solution of the Helmholtz equation in spaces of constant curvature
The Dirac delta distribution with metric g is defined for an open set where dvol g is the Riemannian (volume) measure, invariant under the isometry group SO(d, and on S d R by where dω is the Riemannian volume measure on S d−1 . Notice that as r → 0 + and θ → 0 + , dvol h and dvol s go to the Euclidean volume measure, invariant under the Euclidean motion group E(d), in standard geodesic polar coordinates. Therefore in standard geodesic polar coordinates, we have the following: sin θ 2 sin 2 θ 3 · · · sin d−2 θ d−1 .
In general since we can add any function that satisfies the homogeneous Helmholtz equation to a fundamental solution of the same equation and still have a fundamental solution, we will use this freedom to make our fundamental solution as simple as possible. It is reasonable to expect that there exists a particular spherically symmetric fundamental solution on each and constant angular dependence (invariant under rotations centered about the origin), due to the influence of the point-like nature of the Dirac delta distribution. For a spherically symmetric solution of the Helmholtz equation, the corresponding ∆ S d−1 term vanishes since only the l = 0 term survives. In other words, we expect there to exist a fundamental solution of the Helmholtz equation such that, aside from a multiplicative constant which depends on R, β, and d, In the remainder of this section, we derive fundamental solutions for −∆ ± β 2 on H d R (Section 4.3), −∆ + β 2 on S d R (Section 4.4.1), and study two candidate fundamental solutions for −∆ − β 2 on S d R (Section 4.4.2). We will require a fundamental solution of Helmholtz operators in Euclidean space R d . These are well-known and can be found in, for instance, [4, p. 139]. (4.9) then G d,± β respectively are fundamental solutions for −∆ ± β 2 in Euclidean space R d , where ∆ is the Laplace operator on R d .
Note most authors only present the above theorem for the case d ≥ 2 but it is easily verified to also be valid for the case d = 1 as well. Fundamental solutions for −∆ on H d R , S d R (opposite antipodal) are respectively (see [6, (4.11) where the geodesic distances are given by (3.9), (3.10). These fundamental solutions will be useful to constrain the possible homogeneous solutions which are to be used to obtain our fundamental solutions in the β → 0 + limit (when they exist, see Section 4.1).
Due to the fact that the spaces H d R and S d R are homogeneous with respect to their isometry groups SO(d, 1) and O(d) respectively, and therefore isotropic manifolds, without loss of generality, we are free to map the point x ∈ H d R or x ∈ S d R to the origin. In this case, the global geodesic distance function ρ h coincides with the radial parameter r in geodesic polar coordinates, and ρ s with θ; and therefore we may interchange r with ρ h and θ with ρ s accordingly (cf. (3.9) with r = 0) in our representation of a fundamental solution for the Helmholtz equation on each manifold. Notice that we can add any homogeneous solution to a fundamental solution of the Helmholtz equations (3.14), (3.15), respectively, and still have a fundamental solution of Helmholtz operators. This is because a fundamental solution of the Helmholtz equation must satisfy

Fundamental solutions of Helmholtz operators on H
is the geodesic distance between x and x on the pseudo-sphere of unit radius H d , with x = x/R, x = x /R. Then H d,+ R,β is a fundamental solution for the Helmholtz operator −∆ + β 2 on H d R . Remark 4.4. It has been brought to our attention by one of the referees that a fundamental solution for the Helmholtz equation H d,+ 1,β (x, x ) on the unit-radius hyperboloid H d is given in [24,Theorem 3.3]. Note that the half-space model of hyperbolic geometry is used there, rather than the hyperboloid model which is used in this paper.
Proof . Our derivation for a fundamental solution of the Helmholtz equation H d,+ R,β (x, x ) on the R-radius hyperboloid H d R is as follows. By starting with the spherically symmetric solution 1 which is singular at r = 0, to the homogeneous the Helmholtz equation, we have The constant c is obtained by matching to a Euclidean fundamental solution of the Helmholtz equation in the flat-space limit R → ∞ (see [7,Section 2.4]) with geodesic distance Rρ h and thus where we have used sinh ρ h ∼ ρ h as ρ h → 0. By equating (4.12) with G d,+ β and solving for c, we where n is an odd-integer, so that the associated Legendre function of the second kind becomes a toroidal harmonic. All toroidal harmonics can be obtained, for instance, by using [ is the geodesic distance between x and x on the pseudo-sphere of unit radius H d , x = x/R, x = x /R. Then H d,− R,β is a fundamental solution for the Helmholtz operator −∆−β 2 on H d R . (2π) which proves the full behavior of H d,− R,β for 4β 2 R 2 > (d − 1) 2 . For the case 4β 2 R 2 ≤ (d − 1) 2 , we match up to a fundamental solution of −∆ on H d R (4.10) as β → 0 which requires a fundamental solution with functional dependence as follows Matching the constants c and g at β = (d − 1)/(2R), completes the proof.
where d h (x, x ) is the geodesic distance between two points x, x ∈ H d R .
Proof . Existence is clear. For uniqueness, suppose H ± andH ± are two such functions. Let Since H d R is locally Euclidean, one has by local elliptic regularity that h ± can be extended to a C ∞ -functionĥ ± : H d R → R. It follows from (4.15) for H ± andH ± that lim d h,s (x,x )→∞ĥ ± (x) = 0. (4.16) The strong elliptic maximum/minimum principle on a Riemannian manifold for a bounded domain Ω states that if u is a homogeneous solution to an elliptic operator, then the supremum/infimum of u in Ω coincides with the supremum/infimum of u on the boundary ∂Ω. By using a compact exhaustion sequence Ω k in a non-compact connected Riemannian manifold and passing to a subsequence x k ∈ ∂Ω k such that x k → ∞, the strong elliptic maximum/minimum principle can be extended to non-compact connected Riemannian manifolds with boundary conditions at infinity (see for instance [13,Section 8.3.2]). Taking Ω k ⊂ H d R , the strong elliptic maximum/minimum principle for non-compact connected Riemannian manifolds implies using (4.16) thatĥ ± = 0. Therefore h ± = 0 and H ± (x, x ) =H ± (x, x ) for all x ∈ H d R \{x }.
By Proposition 4.7, for d ≥ 2, the functions H d,± R,β are the unique fundamental solutions for these Helmholtz operators which satisfy the vanishing decay (4.15).  Section 3.3), one should study the behavior of these homogeneous solutions as x → −1 + . In fact, P −µ ν (−x) is the unique solution to the associated Legendre differential equation with the required behavior as x → −1 + (see Section 2.4.1). In particular, (sin ρ s ) −µ P −µ ν (−x) has the required spherically symmetric regular behavior on the opposite pole.

Fundamental solutions of Helmholtz operators
Furthermore, as was seen in Section 4.1, a fundamental solution of the Helmholtz equation on S d R must also satisfy a total integral normalization property (4.4), due to applying the divergence theorem. An alternative integration constraint due to the divergence theorem on a hypersphere can be considered as follows. Consider the compact manifold with boundary by a ball of radius R , 0 < 1, embedded within S d R with center at the origin θ = 0. Let X = ∇S d,± R,β (x, x ) and through (4.2), one has Now use the divergence theorem (4.3), and one obtains what we refer to as the -ball integral constraint, namely which can be verified directly. This is because for some differentiable function u(θ) on S d R , θ ∈ (0, π), the normal derivative in standard geodesic polar coordinates (3.2) is Apart from these local behaviors, fundamental solutions of Helmholtz operators on S d R may or may not have the desired limiting behaviors as β → 0 or in the flat-space limit R → ∞. For instance, there does not exist a fundamental solution of Laplace's equation on S d R with non-vanishing total integral. Hence for fundamental solutions of a single Dirac delta distribution S d,± R,β (4.2), the β → 0 limit should not exist. However, for an opposite antipodal fundamental solution A d,± R,β (4.5), the β → 0 limit should exist, namely (4.11) One would also expect that since all manifolds are locally Euclidean, that a fundamental solution of Helmholtz operators should look locally like a Euclidean fundamental solutions of Helmholtz operators (4.8), (4.9), namely Since the Sommerfeld radiation condition requirement does not apply to the Helmholtz operator −∆ + β 2 , for this operator, problems associated with this requirement should not arise. In fact, we will see the above requirements for this operator is easily satisfied.
A fundamental solution with a single Dirac delta distribution at the origin for −∆ + β 2 on S d R is given in the following theorem.
ρ s := cos −1 (( x, x )) is the geodesic distance between x and x on the hypersphere of unit radius S d , x = x/R, x = x /R. Then S d,+ R,β is a fundamental solution for the Helmholtz operator −∆ + β 2 on S d R . Proof . Apart from a constant in ν, µ, we guess a form of a fundamental solution for this Helmholtz operator for a single Dirac delta distribution at the origin on S d R . Hence our ansatz will be , where x = cos ρ s , and c(ν, µ) is constant in x. This constant, for instance, can be determined through the integral normalization requirement (4.4). We must verify Starting with the left-hand side of (4.20), and using (2.51), noting the standard volume integral after simplifying, one requires that All that remains is to demonstrate that (4.18) produces (4.8) in the flat-space limit, R → ∞. In the flat-space limit ρ s ∼ ϑ R and ν ∼ − 1 2 +iβR, with ϑ → x−x . Using the uniform asymptotics given by (2.64), the demonstration is validated. This completes the proof of Theorem 4.8.
Remark 4.9. Note that (2.47) gives the strength of the cusp on the opposite pole (at θ = π) for fundamental solutions on hyperspheres, such as that given in Theorem 4.8.
A fundamental solution with two opposite antipodal Dirac delta distributions, one at the origin and another on the opposite pole of S d R , for −∆ + β 2 is given in the following theorem. (4.22) where µ and ν are given by (4.19), ρ s = cos −1 (( x, x )) is the geodesic distance between x and x on the hypersphere of unit radius S d , x = x/R, x = x /R, and f µ ν is the odd Ferrers function defined by (2.48). Then A d,+ R,β is an opposite antipodal fundamental solution for the Helmholtz operator −∆ + β 2 on S d R .
Proof . Since the opposite antipodal fundamental solution of the Helmholtz operator −∆+β 2 must be an odd function, it should be given through (2.48). Starting with (4.22) we then perform the flat-space limit R → ∞ by using the asymptotics given by (2.70), along with sin(π(ν − µ)) ∼ − 1 2 e πτ e ±iπµ ; as in the proof of Theorem 4.8, we must consider that ρ s ∼ ϑ/R, ϑ ∼ x−x , x, x ∈ R d . Note that the asymptotic contribution due to −P −µ ν (cos ρ s ) is negligible as compared to P −µ ν (− cos ρ s ). One can then easily see that (4.22) correctly approaches (4.8). Furthermore, since (4.22) is an odd function about θ = π 2 on the hypersphere S d R , the integral normalization requirement (4.4) is clearly satisfied. This completes the proof.
Remark 4.13. If one considers the β ∼ 0 approximation for (4.22), we see that as β → 0, the opposite antipodal fundamental solution of Laplace's equation on the hypersphere (4.11). In order to demonstrate this limit, we have used the binomial, small angle approximations for trigonometric functions, as well as (2.42).
Remark 4.14. The -ball integral constraint (4.17) for A d,+ R,β is also satisfied. This can be verified in an identical fashion to Remark 4.11 while also noting that P −µ ν (x) and its derivative vanish as x → 1 − , µ > 0.
where  Szmytkowski (2007) [31, equation (3.23)] in terms of Gegenbauer functions (cf. (2.37)) for unit radius hyperspheres. Note that in order to convert the wavenumber given in [31], so that it is the same as given here, one must take β 2 R 2 = λ(λ + d − 1), and therefore take with the plus sign chosen to match the formulae given here.
Remark 4.19. The -ball integral constraint (4.17) for S d,− R,β is also satisfied. This can be verified by noting that as ρ s → 0 + , , and using (4.21).
Remark 4.20. If one considers the small β approximation for S d,− R,β (x, x ) we see that as β → 0, In order to obtain this we have used [27, equation (14.5.18)], the binomial and small angle approximations for trigonometric functions. Hence lim β→0 S d,− R,β (x, x ) = ∞, and as one would expect (see Section 4.1) the β → 0 (Laplace) limit of S d,− R,β (x, x ) does not exist.
Remark 4.21. The flat-space limit R → ∞ of the function S d,− R,β should not exist. In fact, it oscillates wildly in this limit, as it approaches the following function Remark 4.27. The opposite antipodal candidate function A d,− R,β (x, x ) approaches the Laplace fundamental solution as β → 0 + , namely (4.11) and in the flat-space limit it approaches   17-4.20. Apart from having the correct flat-space limit, S d,− R,β does not satisfy the integral normalization requirement (4.4) which is not allowed, as shown in Section 4.1. The only critical property that S d,− R,β satisfies is the correct Euclidean flat-space limit. The flat-space behavior of S d,− R,β , displayed in Remark 4.21, may well be due to imposing incorrect boundary conditions since S d R is a compact manifold without boundary, whereas Euclidean space is actually a noncompact manifold with boundary (at infinity). Furthermore, S d,− R,β is partially 'antipodal' to the source point (diverges at the source point as well as at the antipole), which is not a desirable property. It was also pointed out by one of the referees that due to the Γ(µ − ν) factor, S d,− R,β will diverge when β = 0, as is necessary (see Remark 4.19), and at an infinite number of points in the underdamped (oscillatory) regime which correspond to the wavenumber being one of the eigenvalues of the Laplace-Beltrami operator on S d R . Similarly, A d,− R,β has all the right properties that an opposite antipodal fundamental solution of the Helmholtz operator −∆−β 2 on S d R should satisfy, with the exception of the correct flat-space limit. Given a resolution of the flat-space limit for S d,− R,β , the same problem will be solved for A d,− R,β . The ultimate resolution of this Conjecture 4.29 must also be associated with a correct and consistent formulation of the Sommerfeld radiation condition (4.6) on S d R .

Gegenbauer expansions in geodesic polar coordinates
For convenience, define the following degrees In the boundary case, (d − 1) 2 − 4β 2 R 2 = 0, or equivalently, the wavenumber only has dependence on the dimension of the space d and the radius of curvature of the manifold R. In this case, the Ferrers and associated Legendre functions reduce to elementary or elementary transcendental functions. For instance, if µ = d 2 − 1 ∈ Z (d even), then the fundamental solutions are given in terms of complete elliptic integrals (the Legendre functions are toroidal harmonics, and the equivalent for the Ferrers functions). Alternatively, if µ is a half odd integer (d odd), then the Ferrers functions and associated Legendre functions reduce to trigonometric and hyperbolic functions respectively. Given starting points for various degrees and orders, a lattice of degrees and orders which satisfy these properties can be obtained using recurrence relations for these functions, namely [27, equations (14.10 In fact, recent results due to Maier [23, Theorem 6.1] have shown that the space of Ferrers and associated Legendre functions is much more rich than was previously expected. In fact if the orders are integers and the degrees differ by ±1/r (r = 2, 3, 4, 6) from an integer, then the resulting Legendre functions are also given in terms of complete elliptic integrals of the first and second kind! Define the constants

Azimuthal Fourier series in two dimensions
It is interesting to consider how we can obtain azimuthal Fourier expansions for our fundamental solutions. For the particular case where a fundamental solution is given in terms of an associated Legendre function of the second kind, we can use the addition theorem (2.33) for the twodimensional case. One has the following azimuthal Fourier expansions in d = 2.
Similar expansions are possible for the other fundamental solutions. Two-dimensions is the only case where we are able to compute azimuthal Fourier expansions of our fundamental solutions whose coefficients are given without a sum. For dimensions greater than two, the approach used previously for Laplace's equation [6,Section 4] does not easily generalize for the Helmholtz operator.
Lemma A.4. Let l ∈ N 0 . Then Proof . The first equality of (A.14) comes from (A.13) and noting that (−l) n = 0 for n ≥ l + 1. Next, use the addition theorem for Gegenbauer polynomials [27, equation (18.18.8)] with degree l, which is a finite sum over n = 0, 1, 2, . . . , l. The Gegenbauer polynomial on the lefthand side and the first two Gegenbauer polynomials appearing in the sum are expressed in terms of Ferrers functions of the first kind using (2.37). Performing these replacements and simplifying produces the second equality of (A.14).
In what follows we shall make use the following uniqueness theorem of complex analysis, due to Carlson (cf. [33,Section 5.81]).
Remark A.6. The order restriction as stated is only required for z = 0, and can be relaxed to F(z) being of (arbitrary) exponential type for z > 0; however the above weaker version of Carlson's theorem suffices for our purposes.
In order to apply the above theorem, we need estimates on A(z; λ, θ, θ ) for large complex z.