Symmetry, Integrability and Geometry: Methods and Applications Fourier, Gegenbauer and Jacobi Expansions for a Power-Law Fundamental Solution of the Polyharmonic Equation and Polyspherical Addition Theorems

We develop complex Jacobi, Gegenbauer and Chebyshev polynomial expansions for the kernels associated with power-law fundamental solutions of the polyharmonic equation on d-dimensional Euclidean space. From these series representations we derive Fourier expansions in certain rotationally-invariant coordinate systems and Gegenbauer polynomial expansions in Vilenkin's polyspherical coordinates. We compare both of these expansions to generate addition theorems for the azimuthal Fourier coefficients.


Introduction
We have developed a technique for constructing addition theorems for the azimuthal Fourier coefficients of fundamental solutions for linear homogeneous partial differential equations on d-dimensional isotropic Riemannian manifolds. For a fundamental solution, we construct azimuthal Fourier expansions and compare them with eigenfunction expansions in rotationallyinvariant coordinate systems. The construction of eigenfunction expansions for fundamental solutions in separable coordinate systems is in general non-trivial.
In connection with the Laplace operator, we have already constructed some of these addition theorems on R 3 [5,10,11], where we have treated the spherical, cylindrical, oblate spheroidal, prolate spheroidal, parabolic, bispherical and toroidal coordinate systems. One may construct addition theorems in this manner in any rotationally-invariant coordinate system which yields solutions through separation of variables for the Laplace equation. In a similar setting, addition theorems may be generated for other inhomogeneous linear partial differential equations, such as for the Helmholtz, wave and heat equations in arbitrary dimensions. Furthermore an extension of this concept is possible when working with linear partial differential operators on Riemannian manifolds, such as for the Laplace-Beltrami operator (see [9, § 5.1]). Once a Fourier expansion for a fundamental solution is obtained for a partial differential operator on a Riemannian manifold, one must construct eigenfunction expansions for a fundamental solution corresponding to that operator and identify those nested multi-summation and multi-integration eigenfunction expansions which correspond to the Fourier coefficients for that operator.

H.S. Cohl
In this paper, we apply this technique to generate addition theorems from eigenfunction expansions for a fundamental solution of the polyharmonic operator on d-dimensional Euclidean space R d in Vilenkin-Kuznetsov-Smorodinskiȋ polyspherical coordinate systems [32, § 9.5], [34], [33, § 10.5] (hereafter Vilenkin). We have computed azimuthal Fourier expansions for a fundamental solution of this operator, as well as the corresponding eigenfunction expansions in polyspherical coordinates. In each case, the comparison of these two expansions yields new addition theorems for the azimuthal Fourier coefficients.
The main results of this paper are connected with closed-form expressions of hypergeometric orthogonal polynomial expansions for a power-law fundamental solution of the polyharmonic equation on d-dimensional Euclidean space. These expansions are the fundamental building blocks for algorithms to compute the solution of inhomogeneous linear polyharmonic boundary value problems through convolution with the source distribution. These problems are ubiquitous in physics and engineering and include elasticity, electrostatics, magnetostatics, quantum direct and exchange Coulomb interactions, Newtonian gravity, and potential fluid and heat flow, just to name a few. The expansions presented in this paper are crucial for obtaining analytic polyharmonic solutions and for numerical algorithms where expansions are needed for so-called "fast algorithms". Applications of generalized Hopf coordinates, one of the polyspherical coordinate systems we study in this paper include particle physics [24], quantum field theory [2] and cosmology [22]. Furthermore, expansions in hyperspherical coordinate systems have many applications including general atomic multibody theory (see [14,23] and references therein).
From a global analytic partial differential equation perspective, the most important results of this paper are contained in Corollaries 2 and 4. These formulae represent multipole and azimuthal Fourier decompositions for arbitrary powers of the Euclidean distance between two points. They can be used to analytically and numerically solve, in a rapidly convergent fashion, the inhomogeneous linear polyharmonic equation for isolated non-axisymmetric source distributions. Rapid convergence of the Fourier expansions is provided by the fact that for each azimuthal mode there corresponds an infinite number of meridional modes, which are all summed over in our expansions. Furthermore, if pure trigonometric azimuthal dependence exists for a particular source distribution in the inhomogeneous partial differential equation, then Corollary 4 provides a solution in a finite number of terms. In the case of an axisymmetric source distribution, the inhomogeneous polyharmonic solution is obtained from a single m = 0 term in the expansion.
Corollary 2 is powerful in that it provides a mechanism for determining the multipole moments associated with unrestricted powers of the distance between two points. This is in contrast with the Gegenbauer generating function (which is generalized by Theorem 1 and Corollary 1) which provides with the addition theorem for hyperspherical harmonics (B.13) a multipole expansion for this kernel only for powers of the distance given by ν = 2 − d, for d = 3, 4, 5, . . ..
From a special function theoretic perspective, the new results presented in this paper are Theorem 1 and Corollary 1. These series expansions represent fundamental generalizations of Heine's formula [27, (14.28.2)], Gegenbauer's generating function [27, (18.12.4)], and Heine's reciprocal square root identity [8, (3.11)]. Formula (A.14) which provides a connection between the symmetric Jacobi function of the second kind and the associated Legendre function of the second kind, is also interesting. As far as the author is aware, this has not previously appeared in the literature. The addition theorems for associated Legendre functions given by Theorems 2 and 3, and their Corollaries 5, 6, 7, also appear to be new. This paper only uses eigenfunction expansions for a power-law fundamental solution of the polyharmonic equation in Vilenkin's polyspherical coordinates to obtain new addition theorems for the azimuthal Fourier coefficients. We have only treated two different types of Vilenkin's polyspherical coordinates. In higher dimensions, many more types may be considered. The azimuthal Fourier expansion presented in this paper, Corollary 4, can be used to provide new addition theorems for the azimuthal Fourier coefficients in every rotationally-invariant coordinate system which separates the polyharmonic equation on R d . These rotationally-invariant coordinate systems include those of cylindrical, parabolic, and cyclidic type.
In this paper, we take advantage of previously derived closed-form expressions for the separated eigenfunctions in Vilenkin's polyspherical coordinates (found in [19,33] and elsewhere) to derive addition theorems from a power-law fundamental solution of the polyharmonic operator in Euclidean space R d . These addition theorems separate the complicated geometrically-relevant quantity (the azimuthal Fourier coefficients) into functions of the individual variables in the problem. We study all dimensions d ≥ 3 with an emphasis on the simplicity/explicitness of the low-dimensional examples.
This paper is organized as follows. In Section 2, we describe the kernels associated with a fundamental solution of the polyharmonic equation on Euclidean space R d for d ≥ 2 and introduce rotationally-invariant coordinate systems. In Section 3, we prove several new theorems associated with Jacobi, Gegenbauer, and Chebyshev polynomial expansions for the kernels associated with power-law fundamental solutions of the polyharmonic equation on R d . In Section 4, we derive and discuss new addition theorems in Vilenkin's polyspherical coordinates for the azimuthal Fourier coefficients of a fundamental solution for the polyharmonic equation on R d . In Appendix A, we summarize the definitions and properties of the special functions and orthogonal polynomials that we use. In Appendix B, we review Vilenkin's polyspherical coordinates and the corresponding normalized hyperspherical harmonics.
Throughout this paper we rely on the following definitions. Let a 1 , a 2 , a 3 , . . . ∈ C, with C being the set of complex numbers. If i, j ∈ Z and j < i, then j n=i a n = 0 and j n=i a n = 1.
induces a norm (the Euclidean norm) · : R d → [0, ∞), on R d , given by x := (x, x). Please see Appendix A for all notations used in this paper for special functions and orthogonal polynomials.

Fundamental solution of the polyharmonic equation in rotationally-invariant and polyspherical coordinate systems
In Euclidean space R d , let the Laplacian operator ∆ : , then Φ is called polyharmonic. We use the nonnegative Laplacian −∆ ≥ 0. The inhomogeneous polyharmonic equation is given by where we take ρ to be an integrable function so that a solution to (2.2) exists. A fundamental solution for the polyharmonic equation on R d is a function G d k : where δ is the Dirac delta function (generalized function/distribution) and x ∈ R d . Note that this equation is satisfied in the sense of distributions. A fundamental solution of the polyharmonic equation is given as follows (see for instance [4], [15, p. 202], [28, p. 45 where β p,d ∈ Q is defined as β p,d : i . The gamma function Γ : C\−N 0 → C, is a natural generalization of the factorial function. Concerning the logarithmic contribution for d even, k ≥ d/2, the polynomial x − x 2k−d is polyharmonic, so any choice for the constant β p,d is valid. Our choice for this constant is given such that −∆G d k = G d k−1 is satisfied for all k > d/2, and that for k = d/2, the constant vanishes. Note that a solution of the inhomogeneous polyharmonic equation (2.2) is obtained from G d k via a convolution.
(2.4) These coordinate systems are described by d-coordinates: an angle φ ∈ [0, 2π) plus (d − 1)curvilinear coordinates (ξ 1 , . . . , ξ d−1 ). Rotationally-invariant coordinate systems parametrize points on the (d − 1)-dimensional half-hyperplane given by φ = const and R ≥ 0 using the curvilinear coordinates (ξ 1 , . . . , ξ d−1 ). A separable rotationally-invariant coordinate system transforms the polyharmonic equation into a set of d-uncoupled ordinary differential equations with separation constants m ∈ Z and k j ∈ R for 1 ≤ j ≤ d − 2. For a separable rotationally-invariant coordinate system, this uncoupling is accomplished, in general, by assuming a product solution , where the properties of the functions R and A i , for 1 ≤ i ≤ d − 1, and the constants k j for 1 ≤ j ≤ d − 2, depend on the specific separable rotationally-invariant coordinate system in question. Separable coordinate systems are divided into two different classes, those which are simply separable (R = const), and those which are R-separable (see [26]). The Euclidean distance between two points x, x ∈ R d expressed in the rotationally-invariant coordinate system described in (2.4) is where the toroidal parameter χ is The hypersurfaces χ = const are independent of coordinate system and represent hypertori of revolution. We now rewrite (2.3) in terms of the rotationally-invariant coordinate system (2.4). From (2.3) we see that, apart from multiplicative constants, the expression l d k : (R d ×R d )\{(x, x) : x ∈ R d } → R of a fundamental solution for the polyharmonic equation on R d for d even, k ≥ d/2, is given by By expressing l d k in a rotationally-invariant coordinate system (2.4) we obtain where p = k − d/2 ∈ N 0 . Similarly, when working on an even-dimensional Euclidean space Examining (2.7) and (2.8), we see that for computation of Fourier expansions about the azimuthal separation angle (φ − φ ) of l d k and h d k , all that is required is to compute the Fourier cosine series for the following three functions f χ , h χ : R → (0, ∞) and g χ : R → R defined as where p ∈ N 0 , q ∈ N and χ > 1 is a fixed parameter.
The Fourier series of f χ is given in [8] (cf. (4.4) therein) 1 , namely for p ∈ N 0 , where n ∈ {1, 2} is the Neumann factor defined by n := 2 − δ n,0 , δ n,0 ∈ {0, 1} is the Kronecker delta. The Fourier series of h χ is given in [8, (4.5)], namely for p ∈ N, In order to compute Fourier expansion of l d k (2.7) in separable rotationally-invariant coordinate systems, all that remains is to determine the Fourier series of g χ (see [6]). A discussion of Fourier cosine expansions for a logarithmic fundamental solution of the polyharmonic equation on R d (from g χ ) can be found in [7]. The corresponding Gegenbauer polynomial expansions for a logarithmic fundamental solution of the polyharmonic equation on R d can be found in [6].

Jacobi polynomial and limiting expansions for the Euler kernel
In this section we derive Jacobi, Gegenbauer and Chebyshev polynomial of the first kind series expansions of the Euler kernel (z − x) −ν . These series expansions are used to obtain azimuthal Fourier and hyperspherical harmonic expansions for a fundamental solution of the polyharmonic equation on R d .
on any ellipse with foci at ±1 and x in the interior of that ellipse. Then Note 1. It has been brought to the author's attention by Tom Koornwinder that Theorem 1 for ν = −n, n ∈ N 0 , specializes to formula (21) in [21], namely This equivalence is provided by the interesting identity which can be obtained by comparison of Gauss hypergeometric representations.
Proof . Consider the generating function for Gegenbauer polynomials (see, e.g., [27, (18.12.4)]) where we have expressed the Gegenbauer polynomial as a symmetric Jacobi polynomial using (A.7). Utilizing (A.6) in (3.2), reversing the order of the summations, and shifting the n index yields Taking advantage of standard properties such as (A.2), (A.3) produces We substitute the definition of the 3 F 2 generalized hypergeometric function (cf. (A.1)) in the sum over n, and as previously, reverse the order of the two summations and shift the summation index. It then follows using the duplication formula and (A.2)-(A.5), that one has If we apply the right-hand side of (3.3) to (A.13) noting Theorem 12.7.3 (expansion of an analytic function in terms of orthogonal polynomials) in [31] to obtain the regions of convergence, we obtain the desired result.
on any ellipse with foci at ±1 with x in the interior of that ellipse. Then Proof . Let α = β = µ − 1 2 in Theorem 1, and use (A.14) and the definition of the Gegenbauer polynomial in terms of a symmetric Jacobi polynomial (A.7). The points ν ∈ −N 0 which must be removed are singularities originating from the associated Legendre function of the second kind on the right-hand side of (3.4). This complete the proof.
Note that these singularities are removable and correspond to non-negative integer powers of the binomial z − x. See Note 1. These singularities can be removed by taking the limits as ν approaches them.

H.S. Cohl
Corollary 3. Let ν ∈ C \ −N 0 , and x, z ∈ C such that z ∈ C \ (−∞, 1] lie on any ellipse with foci at ±1 with x in the interior of that ellipse. Then Proof . Take the limit as µ → 0 on the right-hand side of (3.4) and use (A.8).

Addition theorems in Vilenkin's polyspherical coordinates
In this section we construct power-law addition theorems in Vilenkin's polyspherical coordinates.

Power-law addition theorem on R d for d ≥ 3 in standard polyspherical coordinates
In standard polyspherical coordinates (B.11) we have the following multi-summation power-law addition theorem.
Proof . If we adopt standard polyspherical coordinates (B.11) (see Fig. 4), then one can obtain an eigenfunction expansion for a power-law fundamental solution of the polyharmonic equation in standard polyspherical coordinates using the Gegenbauer expansion (3.5) with the addition theorem for hyperspherical harmonics (B.13), and the normalized standard hyperspherical harmonics (B.19), obtaining If we expand the product of polyspherical harmonics in (4.2) with (B.19) after reversing the order of the summations, we obtain By comparing the Fourier coefficients of (4.3) with (4.4), we complete the proof of this theorem.

H.S. Cohl
This is just one example of a derived multi-summation addition theorem for arbitrary dimensions. There are an unlimited number of such straightforward examples to generate. In the next section we derive another example which is valid on R d where d is given by a power of two, generalized Hopf coordinates (B.12).

4.2
Power-law addition theorem on R 2 q for q ≥ 2 in generalized Hopf coordinates In generalized Hopf coordinates (B.12) we have the following multi-summation power-law addition theorem.
Proof . If we adopt generalized Hopf coordinates (B.12) (see Fig. 5), then we can use the corresponding harmonics (B.21) in combination with the addition theorem for hyperspherical harmonics (B.13). We compare the Gegenbauer expansion for powers of the distance (3.5) with the Fourier expansion where χ > 1 is given by (4.6). Notice that χ is independent of φ 1 − φ 1 . By using the Gegenbauer expansion (3.5) and inserting the appropriate Gegenbauer polynomial using the addition theorem for hyperspherical harmonics (B.13), we obtain By expanding the product of polyspherical harmonics in (4.8) with (B.21) expressed in terms of surrogate quantum numbers and reversing the order of the sums, we obtain a multi-summation expression for the power of the Euclidean distance between two points in generalized Hopf coordinates. Through comparison of the resulting equation with the m 1 Fourier coefficients of (4.7), we derive (4.5), a multi-summation addition theorem for the associated Legendre function of the second kind with argument χ (cf. (4.6)).

Power-law addition theorems on R 3
In d = 3 there are two ways to construct polyspherical coordinates, with trees of type b a (see Fig. 2b) and ba (see Fig. 2c). We only treat the first tree since the addition theorem from the second tree is trivially obtained from the first. where χ = r 2 + r 2 − 2rr cos θ cos θ 2rr sin θ sin θ .

Power-law addition theorems on R 4
In d = 4 there are five ways to construct polyspherical coordinates, with trees of type b 2 a (see Fig. 3a), bb a (see Fig. 3b), b ba (see Fig. 3c), b 2 a (see Fig. 3d), ca 2 (see Fig. 3e). We only treat the first and the fifth trees since the addition theorems from the second, third and fourth trees are trivially obtained from the first.
If you substitute ν = −2 (a fundamental solution for the Laplacian on R 4 ) in (4.11) then the Legendre functions of the second kind reduce to elementary functions through [1, (8.6.10-11)], and one obtains the following
If you substitute ν = −2 in (4.12), then the Legendre functions of the second kind reduce to elementary functions through [1, (8.6.10-11)], and one obtains the following Note that in the addition theorems (4.12) and (4.13), that if you make the map ϑ → ϑ− π 2 , then this transformation preserves the addition theorems such that m 1 ↔ m 2 . This transformation is equivalent to swapping the position of φ 1 and φ 2 for the tree in Fig. 3e.
The associated Legendre function of the second kind Q µ ν : C \ (−∞, 1] → C, ν + µ / ∈ −N, can be defined in terms of the Gauss hypergeometric function as follows [27, (14.3.7) and § 14.21] for |z| > 1 and elsewhere in z by analytic continuation of the Gauss hypergeometric function. One may also define the associated Legendre function of the second kind using [25, entry 24, p. 161] for |1 − z| > 2. Similarly, the associated Legendre function of the first kind can be defined using the Gauss hypergeometric function [27, (14.3.6) and § 14.21(i)] where |1 − z| < 2, and elsewhere in z by analytic continuation. We can use Whipple's formula to relate the associated Legendre function of the first kind with the associated Legendre function of the second kind. It is given by [1, (8.2.7)] for Re z > 0. The Ferrers function of the first kind (associated Legendre function of the first kind on-the-cut) P µ ν : (−1, 1) → C can be defined as [27, (14.3.1)] There is a relation between certain Gegenbauer polynomials on (−1, 1) and the Ferrers function of the first kind (cf. (8.936.2) in [16]), namely where the double factorial · : {−1, 0, 1, . . .} → N is defined such that The Jacobi function of the second kind Q where α + γ, β + γ / ∈ −N. We can derive a relation between the symmetric Jacobi function of the second kind and the associated Legendre function of the second kind Q (µ−ν+1/2,µ−ν+1/2) n+ν−1 (z) = 2 µ−ν+1/2 Γ (µ + n + 1/2) e iπ(µ−ν+1/2) Γ(ν + n)(z 2 − 1) (µ−ν)/2+1/4 Q ν−µ−1/2 n+µ−1/2 (z), (A.14) where n ∈ N 0 , µ ∈ C \ − 1 2 , − 3 2 , − 5 2 , . . . , ν ∈ C \ −N 0 . The relation (A.14) can be verified by comparing (A.13) with (A.9). and c. For type a, both branches end on a leaf node. The angle corresponding to this type of branching node is φ a ∈ [0, 2π). For type b, the left branch ends on a leaf node and the right branch ends on a branching node. The angle corresponding to this type of branching node is θ b ∈ [0, π]. For type b , the left-branch ends on a branching node and the right branch ends on a leaf node. The angle corresponding to this type of branching node is θ b ∈ − π 2 , π 2 . For type c, both the left and right branches ends on branching nodes (branching nodes of type c are only possible for d ≥ 4). The angle corresponding to this type of branching node is ϑ c ∈ 0, π 2 .

B Vilenkin's polyspherical coordinates and the method of trees
Polyspherical coordinates are hyperspherical coordinates which are described by a radial coordinate r ∈ [0, ∞) plus (d − 1)-angles which together parametrize points on S d−1 r . 2 We will first discuss a general procedure for constructing polyspherical coordinate systems called the "method of trees" 3 . We will give some examples of polyspherical coordinates which will be used in the rest of the paper and then describe the harmonic separated product solutions.
In these rooted trees, there are two types of nodes, leaf nodes and branching nodes. For a coordinate system on R d , there are d-leaf nodes, each corresponding to the particular Cartesian component of an arbitrary position vector x ∈ R d . The branching nodes split into two separate branches, one up to the left and one up to the right. Each branch emanating from a branching node will end on either a leaf node or on another branching node. There are four possibilities for branching nodes (see Fig. 1).
Separation of variables in polyspherical coordinates with (d − 1)-angles, for Laplace's equation on R d produces (d − 1)-separation constants, each of which are called quantum numbers. The quantum numbers corresponding to these angles are all integers. Quantum numbers for a particular tree label the basis of separable solutions for Laplace's equation in that particular coordinate system. With each branching node of the tree, we associate a quantum number as well as an angle. The quantum number corresponding to a (2π)-periodic (azimuthal) angle is called an azimuthal quantum number. Each azimuthal angle corresponds to a branching node of type a and an azimuthal quantum number m ∈ Z. A natural consequence of the method of trees is that there must exist at least one azimuthal angle for each tree, and therefore also for each polyspherical coordinate system.
Branching nodes of type b, b and c (see Fig. 1), are associated with angles, which in turn are associated with quantum numbers which we refer to as angular momentum quantum numbers l ∈ N 0 (see for instance Chapter 10 in [14]). There is always at least one branching node, the root node, and all branching nodes correspond to a particular angle and quantum number. Let us associate each branching node with an angle and its corresponding quantum number. These trees 2 The Riemannian manifold S d r is defined as the set of all points in R d+1 such that x 2 0 + · · · + x 2 d = r 2 (r > 0), with the metric induced from that of the ambient Euclidean space. We denote the d-dimensional hypersphere of unit radius as S d := S d 1 . 3 Describing polyspherical coordinate systems in terms of rooted trees was originally developed in [34] (see also [32, § 9.5], [33, § 10.5]) and has since been used extensively by others in a variety of contexts (see for instance [18,19,20]). parametrize points in polyspherical coordinates as follows. Starting at the root node, traverse the tree upward until you reach the leaf node corresponding to x i . The parametrization for x i is given by the hyperspherical radius r multiplied by cosine or sine of each angle encountered as you traverse the tree upward until you reach the leaf node corresponding to x i . If you branch upwards to the left or upwards to the right at each branching node, multiply by the cosine or sine of the corresponding angle respectively. This procedure produces the appropriate transformation from polyspherical coordinates to Cartesian coordinates.
There are large numbers of equivalent trees and an even larger number of total trees, each with their own specific polyspherical coordinate system. The enumeration of these trees are characterized as follows. For d, b d ∈ N, let b d be the total number of trees. Then b 1 = 1 is the number of possible 1-branch trees. In our context, a 1-branch tree does not exist in isolation. The following recurrence relation gives the total number of trees for arbitrary dimension Using the recurrence relation (B.1), the first few elements of the sequence are given by , . . . , 13}) = (1,2,5,14,42,132,429,1430,4862,16796,58786,208012).
The total number of trees are given in terms of the Catalan numbers C n (see for instance Sloane integer sequence A000108 [29] or p. 200 in [30] is the binomial coefficient for k, n ∈ N 0 with 0 ≤ k ≤ n. If a d ∈ N is the total number of equivalence classes for equivalent trees for a given dimension (determined by a left-right symmetry in the topology of the trees), then a d satisfies the following recurrence relation These are given in terms of the Wedderburn-Etherington numbers (see for instance Sloane integer sequence A001190 [29]). We use a left-to-right recursive naming language for our trees based on a depth-first search (see [12, pp. 540-549]). This naming language is given by listing the types of branching nodes available in a particular tree. In a polyspherical coordinate system the Euclidean distance between two points (2.5) can also be given as where the separation angle γ ∈ [0, π] is defined through the relation using the Euclidean inner product and norm (cf. (1.1)). The method of trees constructs the cosine of the separation angle (B.4) in a direct manner. The cosine of the separation angle will be given by the sum of d-terms, each corresponding to a leaf node of the tree. There is a unique path starting from the root node to each leaf node. It is where N i is the number of branching nodes encountered from the root node to the leaf node, A i,j : R → [−1, 1] is either the trigonometric cosine or sine function depending respectively on whether the left branch or right branch is chosen respectively, and ψ i,j ∈ R is the angle corresponding to the jth branching node for each ith leaf node. The formula for the cosine of the separation angle is unique for each tree.

B.1 Examples of Vilenkin's polyspherical coordinate systems
The simplest example of a polyspherical coordinate system on R d occurs for d = 2 (polar coordinates) where there is one branching node (the root node) and two leaf nodes (see Figs. 2a and 5a). The left-branch ends on the leaf node corresponding to x 1 and the right branch ends on the leaf node corresponding to x 2 , i.e., where r ∈ [0, ∞) and φ ∈ [0, 2π). We refer to this tree as type a. The cosine of the separation angle is given by and corresponding to the angle φ is the azimuthal quantum number m ∈ Z (see Fig. 2a). In d = 3 there are two possible topological trees, each corresponding to one of two different trees. The first tree (see Fig. 2b) corresponds to the coordinate system where θ ∈ [0, π], φ ∈ [0, 2π), i.e., standard spherical coordinates. This is a tree of type ba. The cosine of the separation angle is given by cos γ = cos θ cos θ + sin θ sin θ cos(φ − φ ), and corresponding to the angles θ ∈ [0, π] and φ ∈ [0, 2π) are the quantum numbers l ∈ N 0 and m ∈ Z respectively. The second tree (see Fig. 2c) is of type b a and is equivalent to the first. In d = 4 there are five possible topological trees for polyspherical coordinates. The first tree (see Figs. 3a and 4) corresponds to the coordinate system x 3 = r sin θ 1 sin θ 2 cos φ, x 4 = r sin θ 1 sin θ 2 sin φ, (B.8) where θ 1 , θ 2 ∈ [0, π]. This tree is of type b 2 a, and the cosine of the separation angle is given by cos γ = cos θ 1 cos θ 1 + sin θ 1 sin θ 1 cos θ 2 cos θ 2 + sin θ 2 sin θ 2 cos(φ − φ ) .
Another example of a polyspherical coordinate system which is valid for large dimensions is what we will refer to as generalized Hopf coordinates (see Fig. 5). These coordinates, valid on R 2 q for q ≥ 1, generalize two-dimensional polar coordinates (B.5, type a) and four-dimensional Hopf coordinates (B.9, type ca 2 ). These coordinates are unique in that they correspond to the only trees which contain only themselves in their equivalence class (see (B.2)). These coordinate systems have separated harmonic eigenfunctions which are given in terms of complex exponentials, and for q ≥ 2, non-symmetric Jacobi polynomials (see Fig. 5).

B.2 Hyperspherical harmonics in polyspherical coordinates
The eigenfunction expansions for a power-law fundamental solution of the polyharmonic equation in polyspherical coordinates can be derived using a Gegenbauer polynomial expansion for the H.S. Cohl relevant kernel (see Corollary 2) in conjunction with the addition theorem for hyperspherical harmonics. This addition theorem is given by (for a proof see [35]; see also § 10.2.1 in [14]), where K stands for a set of (d−2)-quantum numbers identifying harmonics for a given value of n ∈ N 0 , and cos γ is the cosine of the separation angle between two arbitrary vectors x, x ∈ R d (see (B.4)). The functions Y K n : S d−1 → C are the normalized hyperspherical harmonics. Normalization of the hyperspherical harmonics is achieved through where dΩ is the Riemannian volume measure on S d−1 .
The general basis functions that one obtains by putting coordinates on the d-dimensional unit hypersphere S d−1 , can be specified as solutions to the angular part of Laplace's equation on R d . These correspond to separated solutions of Laplace's equation, using the Laplace-Beltrami operator on the hypersphere S d−1 . The following numbers are associated with each branching node m ∈ Z, l, l α , l β ∈ N 0 (see Fig. 1). The number of vertices above each branching node l α and l β are represented by S α and S β respectively.