Addition Formula and Related Integral Equations for Heine–Stieltjes Polynomials

Abstract

It is shown that symmetric products of Heine–Stieltjes quasi-polynomials satisfy an addition formula. The formula follows from the relationship between Heine–Stieltjes quasi-polynomials and spaces of generalized spherical harmonics, and from the known explicit form of the reproducing kernel of these spaces. In special cases, the addition formula is written out explicitly and verified. As an application, integral equations for Heine–Stieltjes quasi-polynomials are found.

1. Introduction

Let S be a set and $\mathcal{H}$ a finite dimensional vector space of real-valued functions defined on S. Suppose there is an inner product $\langle f,g\rangle$ defined for $f,g\in \mathcal{H}$ which turns $\mathcal{H}$ into a Hilbert space. In this situation, there exists a unique reproducing kernel $K(x,y)$ for $\mathcal{H}$, see Aronszajn ([1] §1.1–1.3). This kernel is a real-valued function defined on $S\times S$ with the property that $K(x,y)$ as a function of y belongs to $\mathcal{H}$ for every $x\in S$, and

$$f(x)=\langle f,K(x,·)\rangle \mathrm{for}\mathrm{every}f\in \mathcal{H},x\in S.$$

(1)

In fact, the existence of K is easy to see. We choose an orthonormal basis $e_1,e_2,\cdots ,e_d$ of $\mathcal{H}$, $d=dim\mathcal{H}$. Then,

$$K(x,y)=\sum _{i=1}^d{e}_i(x)e_i(y)$$

(2)

is a reproducing kernel. In some applications, the basis functions $e_i$ are represented by special functions. Then, Equation (2) can be called an addition formula for these special functions if the kernel admits an explicit representation.

As an example, we consider $S={\mathbf{S}}^k$, where ${\mathbf{S}}^k$ is the unit sphere in ${\mathbb{R}}^{k+1}$ and k is a positive integer. For a non-negative integer m, let $\mathcal{H}={\mathcal{H}}_m$ be the space of spherical (surface) harmonics of degree m. A spherical harmonic of degree m is a harmonic polynomial in $k+1$ variables $x_0,x_1,\cdots ,x_k$, homogeneous of degree m. The inner product of ${\mathcal{H}}_m$ is given by

$$\langle f,g\rangle ={\int }_{{\mathbf{S}}^k}f(\mathbf{x})g(\mathbf{x})dS(\mathbf{x}),$$

where we normalize the surface measure of ${\mathbf{S}}^k$ to one:

$$1={\int }_{{\mathbf{S}}^k}dS.$$

In this case, there is an explicit formula for the reproducing kernel:

$$K_m(\mathbf{x},\mathbf{y})=\frac{m+\frac{k-1}{2}}{\frac{k-1}{2}}{C}_m^{(\frac{k-1}{2})}(x_0{y}_0+\cdots +x_k{y}_k),$$

(3)

where $C_m^{(\lambda )}$ denotes the Gegenbauer polynomial of degree m as shown by Hochstadt ([2] §6.3 (25)). If $k=2$, then $K_m=(2m+1)P_m(x_0{y}_0+x_1{y}_1+x_2{y}_2)$, where $P_m$ denotes the Legendre polynomial. Using spherical coordinates, we construct the standard basis of ${\mathcal{H}}_m$ employing associated Legendre functions, and we obtain the classical addition theorem for Legendre polynomials ([2] §5.5).

In this paper, we treat another special case of Equation (2). In this instance, the special functions appearing on the right-hand side of Equation (2) are Heine–Stieltjes quasi-polynomials—see Section 2 for their definition. Heine–Stieltjes quasi-polynomials are solutions $E(t)$ of a Fuchsian equation in the form $E(t)=g(t)p(t)$, where $g(t)$ is an explicitly known function and $p(t)$ is a polynomial in one variable t. In the special case that $g(t)=1$ they are also called Heine–Stieltjes polynomials. These polynomials were introduced by Stieltjes [3] based on previous work by Heine ([4] pp. 445–479). There exists a considerable amount of research on Heine–Stieltjes polynomials in mathematics and physics, for example, see [5,6,7,8]. Heine–Stieltjes polynomials also appear in the Digital Library of Mathematical Functions [9] in Section 31.15 (they are called Stieltjes polynomials in [9]). The kernel on the left-hand side of Equation (2) involves an integral over a Gegenbauer polynomial, see Section 4. It is the reproducing kernel for a space ${\mathcal{H}}_m$ of generalized spherical harmonics, see Section 3 for the definition of ${\mathcal{H}}_m$. This reproducing kernel was found by Xu ([10] Theorem 3.3). This result is related to a product formula for Jacobi polynomials due to Dijksma and Koornwinder [11].

In this way, we obtain an addition theorem for Heine–Stieltjes quasi-polynomials. In general, Heine–Stieltjes quasi-polynomials cannot be represented explicitly but they can be computed numerically as shown in the numerical example in Section 5. However, the kernel K has the explicit representation Equation (25). Therefore, we observe the remarkable fact that although Heine–Stieltjes quasi-polynomials do not allow an explicit representation, a certain combination of them does admit an explicit representation. In the theory of Heine–Stieltjes quasi-polynomials, such formulas of an explicit nature are very rare.

In the special case $k=2$, Heine–Stieltjes polynomials reduce to Heun polynomials ([12] A 3.6) and we find an addition formula for Heun quasi-polynomials which is still new. Lamé polynomials are special cases of Heun polynomials. The corresponding addition formula for Lamé polynomials can be found in Hobson ([13] p. 475).

In Section 6, we apply equation Equation (1) to obtain nonlinear integral equations for Heine–Stieltjes polynomials. In the special case $k=2$, we obtain integral equations for Heun polynomials comparable to but different from equations given by Arscott [14] and Sleeman [15]. The reader is also referred to papers by Kalnins and Miller [16,17] which contain related results.

Finally, we make a few remarks concerning notation. We denote the non-negative integers by $\mathbb{N}$ and the real numbers by $\mathbb{R}$. Boldface letters denote vectors or multi-indices. Throughout, we use $\mathbf{n}=(n_1,n_2,\cdots ,n_k)$ and $\mathbf{p}=(p_0,p_1,\cdots ,p_k)$, where k is a given positive integer. $\mathbf{n}$ is an “oscillation multi-index” where $n_j\in \mathbb{N}$ for all $j=1,2,\cdots ,k$. It counts the number of zeros of a Heine–Stieltjes quasi-polynomial in k disjoint intervals. $\mathbf{p}$ is a “parity multi-index”, where $p_j\in \{0,1\}$ for every $j=0,1,\cdots ,k$. It determines the parity of a function $f:{\mathbb{R}}^{k+1}\to \mathbb{R}$. The function $f(x_0,x_1,\cdots ,x_k)$ is said to have parity $\mathbf{p}$ if it is an even function of $x_j$ when $p_j=0$ and an odd function of $x_j$ when $p_j=1$ for every $j=0,1,\cdots ,k$. We use the standard notation $|\mathbf{n}|=n_1+n_2+\cdots +n_k$, $|\mathbf{p}|=p_0+p_1+\cdots +p_k$.

2. Heine–Stieltjes Quasi-Polynomials

Throughout the paper, two sets of parameters are given as

$$a_0<a_1<\cdots <a_k$$

(4)

and

$${\alpha }_0,{\alpha }_1,\cdots ,{\alpha }_k\in (-\frac{1}{2},\infty ).$$

(5)

Consider multi-indices $\mathbf{n}=(n_1,n_2,\cdots ,n_k)\in {\mathbb{N}}^k$ of non-negative integers and $\mathbf{p}=(p_0,p_1,\cdots ,p_k)\in {\{0,1\}}^{k+1}$. We define the Heine–Stieltjes quasi-polynomial $E_{\mathbf{n},\mathbf{p}}$ by

$$E_{\mathbf{n},\mathbf{p}}(t):=d_{\mathbf{n},\mathbf{p}}\prod _{j=0}^k{|t-a_j|}^{p_j/2}\prod _{\ell =1}^{|\mathbf{n}|}(t-{\theta }_{\ell }),t\in \mathbb{R},$$

(6)

where $d_{\mathbf{n},\mathbf{p}}>0$ is a normalization factor to be determined by Equation (22), and

$${\theta }_1<{\theta }_2<\cdots <{\theta }_{|n|},|\mathbf{n}|:=n_1+\cdots +n_k,$$

with the first $n_1$ of $\theta$’s lying in $(a_0,a_1)$, then the next $n_2$ of $\theta$’s lying in $(a_1,a_2)$, and so on until the last $n_k$ of $\theta$’s lying in $(a_{n_{k-1}},a_{n_k})$. The $\theta$’s are uniquely determined by the condition that $E_{\mathbf{n},\mathbf{p}}$ is a solution of the Fuchsian Equation (2.1) [18]

$$\prod _{j=0}^k(t-a_j)\left[E^{″}+\sum _{j=0}^k\frac{{\alpha }_j+\frac{1}{2}}{t-a_j}{E}^{\prime }\right]+\left[-\frac{1}{2}\sum _{j=0}^k\frac{p_j{\alpha }_j{A}_j}{t-a_j}+\sum _{i=0}^{k-1}{\lambda }_i{t}^i\right]E=0$$

(7)

for some suitable values of ${\lambda }_0,\cdots ,{\lambda }_{k-1}$ (which are also uniquely determined.) The constants $A_j$ in Equation (7) are defined by

$$A_j:=\prod _{\genfrac{}{}{0pt}{}{i=0}{i\ne j}}^k(a_j-a_i)\mathrm{for}j=0,1,\cdots ,k.$$

(8)

Equivalently, the $\theta$’s are determined by the system of equations

$$\sum _{\genfrac{}{}{0pt}{}{q=1}{q\ne \ell }}^{|\mathbf{n}|}\frac{2}{{\theta }_{\ell }-{\theta }_q}+\sum _{j=0}^k\frac{{\alpha }_j+\frac{1}{2}+p_j}{{\theta }_{\ell }-a_j}=0,\ell =1,2,\cdots ,|\mathbf{n}|.$$

(9)

For the existence and uniqueness statements, see Szegö ([19] Section 6.8). Apart from the constant factor $d_{\mathbf{n},\mathbf{p}}$, the definition of $E_{\mathbf{n},\mathbf{p}}$ agrees with the one given in ([18] §2).

3. Generalized Spherical Harmonics

By definition, a generalized spherical harmonic of degree m ([18] §4) is a polynomial $f(x_0,x_1,\cdots ,x_k)$, homogeneous of degree m, satisfying the equation

$${\Delta }_{\mathbf{\alpha }}f:=\sum _{j=0}^k{\mathcal{D}}_j^{2}f=0,\mathbf{\alpha }=({\alpha }_0,{\alpha }_1,\cdots ,{\alpha }_k),$$

(10)

introduced by Dunkl [20] and also by Dunkl and Xu [21]. In Equation (10) we use the generalized partial derivatives

$${\mathcal{D}}_{j}f(\mathbf{x}):=\frac{\partial }{\partial x_j}f(\mathbf{x})+{\alpha }_j\frac{f(\mathbf{x})-f({\sigma }_j\mathbf{x})}{x_j},\mathbf{x}=(x_0,x_1,\cdots ,x_k),$$

(11)

where ${\sigma }_j$ is the reflection at the jth coordinate plane:

$${\sigma }_j(x_0,x_1,\cdots ,x_k)=(x_0,x_1,\cdots ,x_{j-1},-x_j,x_{j+1},\cdots ,x_k).$$

(12)

Equation (10) contains the given parameters ${\alpha }_0,{\alpha }_1,\cdots ,{\alpha }_k$ from Equation (5). If ${\alpha }_j=0$ for all j then the equation reduces to the Laplace equation, and we are in the classical case.

Let $\mathcal{H}={\mathcal{H}}_m$ denote the linear space of generalized spherical harmonics of degree m considered as functions on ${\mathbf{S}}^k$. Its dimension is

$$dim{\mathcal{H}}_m=\left(\genfrac{}{}{0pt}{}{m+k-1}{k-1}\right)+\left(\genfrac{}{}{0pt}{}{m+k-2}{k-1}\right).$$

We introduce a weighted inner product on ${\mathcal{H}}_m$ by defining

$$\langle f,g\rangle :={\int }_{{\mathbf{S}}^k}w(\mathbf{x})f(\mathbf{x})g(\mathbf{x})dS(\mathbf{x}),$$

(13)

where the weight function is given by

$$w(\mathbf{x}):=M|x_0{|}^{2{\alpha }_0}|x_1{|}^{2{\alpha }_1}\cdots {|x_k|}^{2{\alpha }_k}$$

(14)

with the constant M chosen such that

$${\int }_{{\mathbf{S}}^k}w(\mathbf{x})dS(\mathbf{x})=1.$$

Since we assume that ${\alpha }_j>-\frac{1}{2}$, for each $j=0,1,\cdots ,k$, $w(\mathbf{x})$ is integrable on ${\mathbf{S}}^k$.

Let $\mathbf{x}=(x_0,x_1,\cdots ,x_k)\in {\mathbf{S}}_{+}^k$, where

$${\mathbf{S}}_{+}^k:=\{\mathbf{x}\in {\mathbf{S}}^k:x_j>0\mathrm{for}\mathrm{all}j=0,1,\cdots ,k\}.$$

(15)

Its sphero-conal coordinates $(s_1,\cdots ,s_k)$ ([22] §1.3), ([18] §4) lie in

$$Q:=(a_0,a_1)\times (a_1,a_2)\times \cdots \times (a_{k-1},a_k)$$

(16)

and they are determined by the equations

$$\sum _{j=0}^k\frac{x_j^2}{s_i-a_j}=0\mathrm{for}i=1,\cdots ,k.$$

(17)

In Equations (16) and (17), we use the parameters in Equation (4).

This defines a bijective map from ${\mathbf{S}}_{+}^k$ to Q with inverse

$$x_j^2=\frac{{\prod }_{i=1}^k(s_i-a_j)}{{\prod }_{\genfrac{}{}{0pt}{}{i=0}{i\ne j}}^k(a_i-a_j)}.$$

(18)

Let $E_{\mathbf{n},\mathbf{p}}$ be a Heine–Stieltjes quasi-polynomial. We introduce the generalized sphero-conal harmonic ([18] §4)

$$G_{\mathbf{n},\mathbf{p}}(x_0,x_1,\cdots ,x_k):=E_{\mathbf{n},\mathbf{p}}(s_1)E_{\mathbf{n},\mathbf{p}}(s_2)\cdots E_{\mathbf{n},\mathbf{p}}(s_k),$$

(19)

where $s_1,\cdots ,s_k$ denote sphero-conal coordinates for $\mathbf{x}\in {\mathbf{S}}_{+}^k$. It is shown in [18] that

$$G_{\mathbf{n},\mathbf{p}}(\mathbf{x})={(d_{\mathbf{n},\mathbf{p}})}^k{c}_{\mathbf{n},\mathbf{p}}{x}_0^{p_0}{x}_1^{p_1}\cdots x_k^{p_k}\prod _{\ell =1}^{|\mathbf{n}|}\sum _{j=0}^k\frac{x_j^2}{{\theta }_{\ell }-a_j},$$

(20)

when $E_{\mathbf{n},\mathbf{p}}$ is written in the form Equation (6),

$$c_{\mathbf{n},\mathbf{p}}:={(-1)}^{|\mathbf{n}|}\prod _{j=0}^k{|A_j|}^{p_j/2}\prod _{\ell =1}^{|\mathbf{n}|}\prod _{i=0}^k(a_i-{\theta }_{\ell }),$$

(21)

and $A_j$ is defined in Equation (8). Moreover, $G_{\mathbf{n},\mathbf{p}}\in {\mathcal{H}}_m$ is a generalized spherical harmonic of degree $m=2|\mathbf{n}|+|\mathbf{p}|$, $|\mathbf{p}|:=p_0+p_1+\cdots +p_k$, and it has parity $\mathbf{p}$, that is

$$G_{\mathbf{n},\mathbf{p}}({\sigma }_j\mathbf{x})={(-1)}^{p_j}{G}_{\mathbf{n},\mathbf{p}}(\mathbf{x}).$$

By choice of $d_{\mathbf{n},\mathbf{p}}$ in Equation (6), we have

$$\langle G_{\mathbf{n},\mathbf{p}},G_{\mathbf{n},\mathbf{p}}\rangle =1,$$

(22)

using the inner product in ${\mathcal{H}}_m$, defined by Equation (13). The normalization can be carried out based on the formula ([23] Appendix B)

$${\int }_{{\mathbf{S}}^k}w(\mathbf{x})x_0^{2q_0}{x}_1^{2q_1}\cdots x_k^{2q_k}dS(\mathbf{x})=\frac{{({\rho }_0)}_{q_0}\cdots {({\rho }_k)}_{q_k}}{{(\rho )}_{q_0+\cdots +q_k}},$$

(23)

where ${(a)}_q=a(a+1)\cdots (a+q-1)$ (${(a)}_0:=1$) denotes the Pochhammer symbol, and here and in the following we abbreviate

$${\rho }_j:={\alpha }_j+\frac{1}{2},\rho :=\sum _{j=0}^k{\rho }_j.$$

(24)

The following result is known from ([18] Theorem 3).

Theorem 1.

Let $m\in \mathbb{N}$. The system of all generalized sphero-conal harmonics $G_{\mathbf{n},\mathbf{p}}$ of degree $m=2|\mathbf{n}|+|\mathbf{p}|$ forms an orthonormal basis for ${\mathcal{H}}_m$ with respect to the inner product Equation (13).

4. Addition Formula

For given m, k, and $\mathbf{\alpha }=({\alpha }_0,{\alpha }_1,\cdots ,{\alpha }_k)$, we consider the space ${\mathcal{H}}_m$ of generalized spherical harmonics equipped with the inner product Equation (13). Let $K_m$ be its reproducing kernel. If all ${\alpha }_j$ are positive, Xu ([10] Theorem 3.3) derived the following beautiful formula for $K_m$:

$$\begin{array}{c}{K}_m(\mathbf{x},\mathbf{y})=\frac{m+\rho -1}{\rho -1}\underset{{[-1,1]}^{k+1}}{\int }{C}_m^{(\rho -1)}(x_0{y}_0{u}_0+\cdots +x_k{y}_k{u}_k)\\ \times \prod _{i=0}^{k}c({\alpha }_i)(1+u_i){(1-u_i^2)}^{{\alpha }_i-1}d{u}_0\cdots d{u}_k.\end{array}$$

(25)

Here, $C_n^{(\lambda )}$ denotes the Gegenbauer polynomial, $\rho$ is from Equation (24), $c(\mu )$ is defined by

$$\frac{1}{c(\mu )}={\int }_{-1}^1{(1-u^2)}^{\mu -1}du=\frac{\Gamma (\frac{1}{2})\Gamma (\mu )}{\Gamma (\mu +\frac{1}{2})},\mu >0,$$

and

$$\mathbf{x}=(x_0,x_1,\cdots ,x_k)\in {\mathbf{S}}^k,\mathbf{y}=(y_0,y_1,\cdots ,y_k)\in {\mathbf{S}}^k.$$

The kernel $K_m$ can be calculated without carrying out a numerical integration. When we expand $C_m^{(\rho -1)}(x_0{y}_0{u}_0+\cdots +x_k{y}_k{u}_k)$, we obtain a linear combination of monomials ${(x_0{y}_0{u}_0)}^{q_0}{(x_1{y}_1{u}_1)}^{q_1}\cdots {(x_k{y}_k{u}_k)}^{q_k}$, where $q_0,q_1,\cdots ,q_k$ are non-negative integers. To carry out the integration in Equation (25), we use that, for all $q\in \mathbb{N}$ and $\mu >0$,

$$c(\mu ){\int }_{-1}^1{u}^q(1+u){(1-u^2)}^{\mu -1}du=\frac{{(\frac{1}{2})}_r}{{(\mu +\frac{1}{2})}_r},$$

(26)

where

$$r=\lfloor \frac{1}{2}(q+1)\rfloor$$

denotes the largest integer less than or equal to $\frac{1}{2}(q+1)$. The formula in Equation (25) can be used only for ${\alpha }_j>0$. However, using Equation (26), we see that $K_m$ is a well-defined polynomial in the variables $x_0,\cdots ,x_k,y_0,\cdots ,y_k$ whenever ${\alpha }_j>-\frac{1}{2}$ for all $j=0,1,\cdots ,k$. Therefore, the function $K_m$ defined by Equation (25) can be extended analytically to the parameter domain ${\rho }_j>0$, $j=0,1,\cdots ,k$, and this extended function is the reproducing kernel for ${\mathcal{H}}_m$.

For illustration, let us compute $K_2$. We start with the Gegenbauer polynomial

$$C_2^{(\rho -1)}(z)=(\rho -1)(2\rho z^2-1).$$

We substitute $z={\sum }_{j=0}^k{x}_j{y}_j{u}_j$ and use Equation (26). We obtain

$$K_2(\mathbf{x},\mathbf{y})=\rho (\rho +1)\left[\sum _{j=0}^k\frac{x_j^2{y}_j^2}{{\rho }_j}+\sum _{i=0}^k\sum _{j=i+1}^k\frac{x_i{x}_j{y}_i{y}_j}{{\rho }_i{\rho }_j}\right]-\rho -1.$$

(27)

We see that this kernel is well-defined when ${\rho }_j>0$ for all $j=0,1,\cdots ,k$. If $\mathbf{\alpha }=\mathbf{0}$, $K_m$ reduces to the classical kernel Equation (3) through a limiting process using that

$$\underset{\alpha \to 0^{+}}{lim}\frac{{\int }_{-1}^{1}g(t){(1-t^2)}^{\alpha -1}dt}{{\int }_{-1}^1{(1-t^2)}^{\alpha -1}dt}=\frac{1}{2}(g(-1)+g(1))$$

for every continuous function $g(t)$ defined on $[-1,1]$.

We now obtain the following addition formula for generalized sphero-conal harmonics.

Theorem 2.

For every $m\in \mathbb{N}$ and all $\mathbf{x},\mathbf{y}\in {\mathbf{S}}^k$ we have

$$K_m(\mathbf{x},\mathbf{y})=\sum _{2|\mathbf{n}|+|\mathbf{p}|=m}{G}_{\mathbf{n}\mathbf{p}}(\mathbf{x})G_{\mathbf{n},\mathbf{p}}(\mathbf{y}),$$

(28)

where the summation extends over all pairs $\mathbf{n},\mathbf{p}$ with $2|\mathbf{n}|+|\mathbf{p}|=m$.

Proof. 

The space ${\mathcal{H}}_m$ equipped with the inner product Equation (13) has the reproducing kernel $K_m(\mathbf{x},\mathbf{y})$ given by Equation (25). By Theorem 1, the system of generalized sphero-conal harmonics, $G_{\mathbf{n},\mathbf{p}},2|\mathbf{n}|+|\mathbf{p}|=m$, is an orthonormal basis of ${\mathcal{H}}_m$. Then, Equation (28) follows from Equation (2).  □

We may decompose the kernel $K_m$ as

$$K_m(\mathbf{x},\mathbf{y})=\sum _{\mathbf{p}}{K}_{m,\mathbf{p}}(\mathbf{x},\mathbf{y}),$$

(29)

where the sum extends over all $\mathbf{p}\in {\{0,1\}}^{k+1}$ and $K_{m,\mathbf{p}}(\mathbf{x},\mathbf{y})$ is the sum of all monomials $c{x}_0^{q_0}{x}_1^{q_1}\cdots x_k^{q_k}{y}_0^{q_0}{y}_1^{q_1}\cdots y_k^{q_k}$ appearing in $K_m(\mathbf{x},\mathbf{y})$ that have parity $\mathbf{p}=(p_0,p_1,\cdots ,p_k)$, that is, $p_j+q_j$ is even for every $j=0,1,\cdots ,k$. Therefore, using the reflections from Equation (12) we have

$$K_{m,\mathbf{p}}({\sigma }_j\mathbf{x},\mathbf{y})=K_{m,\mathbf{p}}(\mathbf{x},{\sigma }_j\mathbf{y})={(-1)}^{p_j}{K}_{m,\mathbf{p}}(\mathbf{x},\mathbf{y})\mathrm{for}\mathrm{all}j=0,1,\cdots ,k.$$

We note that $K_{m,\mathbf{p}}$ is the reproducing kernel of the subspace ${\mathcal{H}}_{m,\mathbf{p}}$ of ${\mathcal{H}}_m$ consisting of functions with parity $\mathbf{p}$. For example, Equation (27) yields

$$K_{2,\mathbf{0}}(\mathbf{x},\mathbf{y})=\rho (\rho +1)\sum _{j=0}^k\frac{x_j^2{y}_j^2}{{\rho }_j}-\rho -1.$$

(30)

Since $G_{\mathbf{n},\mathbf{p}}$ has parity $\mathbf{p}$, Theorem 2 implies

$$K_{m,\mathbf{p}}(\mathbf{x},\mathbf{y})=\sum _{2|\mathbf{n}|=m-|\mathbf{p}|}{G}_{\mathbf{n}\mathbf{p}}(\mathbf{x})G_{\mathbf{n},\mathbf{p}}(\mathbf{y})\mathrm{for}\mathbf{x},\mathbf{y}\in {\mathbf{S}}^k,$$

(31)

where m and $\mathbf{p}$ are given and the sum extends over all $\mathbf{n}$ satisfying $2|\mathbf{n}|=m-|\mathbf{p}|$.

Let $\mathbf{s}=(s_1,\cdots ,s_k)\in Q$, $=(t_1,\cdots ,t_k)\in Q$ be sphero-conal coordinates for $\mathbf{x},\mathbf{y}\in {\mathbf{S}}_{+}^k$, respectively. Let $L_m(\mathbf{s},\mathbf{t})$ be the reproducing kernel $K_m$ written in sphero-conal coordinates. The transformation formulas are given by Equation (18).

Then Theorem 2 and Equation (19) imply the following addition formula for Heine–Stieltjes quasi-polynomials.

Theorem 3.

For every $m\in \mathbb{N}$ and all $\mathbf{s},\mathbf{t}\in Q$ we have

$$L_m(\mathbf{s},\mathbf{t})=\sum _{2|\mathbf{n}|+|\mathbf{p}|=m}\prod _{i=1}^k{E}_{\mathbf{n},\mathbf{p}}(s_i)E_{\mathbf{n},\mathbf{p}}(t_i).$$

(32)

Similarly, let $L_{m,\mathbf{p}}(\mathbf{s},\mathbf{t})$ be the kernel $K_{m,\mathbf{p}}$ written in sphero-conal coordinates. Then, using Equations (19) and (31), we obtain for any m and $\mathbf{p}$

$$L_{m,\mathbf{p}}(\mathbf{s},\mathbf{t})=\sum _{2|\mathbf{n}|=m-|\mathbf{p}|}\prod _{i=1}^k{E}_{\mathbf{n},\mathbf{p}}(s_i)E_{\mathbf{n},\mathbf{p}}(t_i).$$

(33)

In the special case $\mathbf{p}=\mathbf{0}$, Equation (33) is an addition formula for Heine–Stieltjes polynomials and the kernel $L_{m,\mathbf{p}}$ is a polynomial in $s_1,\cdots ,s_k$, $t_1,\cdots ,t_k$.

For illustration, let us write out the addition formula in detail if $m=2$ when

$$dim{\mathcal{H}}_2=\frac{k(k+1)}{2}+k.$$

We begin by finding all pairs $\mathbf{n},\mathbf{p}$ with $2|\mathbf{n}|+|\mathbf{p}|=2$. There are two cases: (1) $|\mathbf{n}|=0$, $|\mathbf{p}|=2$; and (2) $|\mathbf{n}|=1$, $|\mathbf{p}|=0$. In the first case, we have $\mathbf{n}=\mathbf{0}$ and there are $\frac{k(k+1)}{2}$ different $\mathbf{p}$’s with $|\mathbf{p}|=2$. One of these $\mathbf{p}$ is $\mathbf{p}=(p_0,p_1,\cdots ,p_k)$ with $p_i=p_j=1$, $i<j$ and all other components of $\mathbf{p}$ equal to zero. Then, Equation (20) gives

$$G_{\mathbf{0},\mathbf{p}}(\mathbf{x})={(d_{0,\mathbf{p}})}^k|A_i{|}^{1/2}{|A_j|}^{1/2}{x}_i{x}_j.$$

(34)

From Equations (22) and (23), we obtain

$$d_{\mathbf{0},\mathbf{p}}={\left(\frac{\rho (\rho +1)}{{\rho }_i{\rho }_j|A_i||A_j|}\right)}^{1/(2k)}$$

so

$$G_{\mathbf{0},\mathbf{p}}(\mathbf{x})={\left(\frac{\rho (\rho +1)}{{\rho }_i{\rho }_j}\right)}^{1/2}{x}_i{x}_j,$$

(35)

and, by Equation (6),

$$E_{\mathbf{0},\mathbf{p}}(t)={\left(\frac{\rho (\rho +1)}{{\rho }_i{\rho }_j|A_i||A_j|}\right)}^{1/(2k)}|t-a_i{|}^{1/2}{|t-a_j|}^{1/2}.$$

(36)

From Equation (27), we find

$$K_{2,\mathbf{p}}(\mathbf{x},\mathbf{y})=\frac{\rho (\rho +1)}{{\rho }_j{\rho }_j}{x}_i{x}_j{y}_i{y}_j.$$

Now, Equation (35) verifies Equation (31) because the sum on the right-hand side of Equation (31) contains only one term.

In the second case, we have $\mathbf{p}=\mathbf{0}$ and there are k different $\mathbf{n}$’s with $|\mathbf{n}|=1$. One of these $\mathbf{n}$ is $\mathbf{n}=(n_1,\cdots ,n_k)$ with $n_h=1$ and all other components of $\mathbf{n}$ equal to zero. Let us abbreviate $E_{\mathbf{n},\mathbf{p}}$ as $E_h$ and $G_{\mathbf{n},\mathbf{p}}$ as $G_h$. Then,

$$E_h(t)=d_h(t-{\theta }_h),$$

where $d_h=d_{\mathbf{n},\mathbf{0}}$, and, by Equation (9), ${\theta }_h$ is the unique solution of the equation

$$\sum _{j=0}^k\frac{{\alpha }_j+\frac{1}{2}}{{\theta }_h-a_j}=0$$

(37)

lying in the interval $(a_{h-1},a_h)$. By Equation (20), we obtain

$$G_h(\mathbf{x})=e_h\sum _{j=0}^k\frac{x_j^2}{{\theta }_h-a_j},$$

(38)

where

$$e_h=-d_h^k\prod _{i=0}^k(a_i-{\theta }_h).$$

(39)

Using Equation (23), we calculate

$$\begin{array}{c}{\int }_{{\mathbf{S}}^k}w(\mathbf{x})G_h{(\mathbf{x})}^{2}dS(\mathbf{x})=e_h^2\sum _{i,j=0}^k\frac{1}{({\theta }_h-a_i)({\theta }_h-a_j)}{\int }_{{\mathbf{S}}^k}w(\mathbf{x})x_i^2{x}_j^{2}dS(\mathbf{x})\hfill \\ =\frac{e_h^2}{\rho (\rho +1)}\left(\sum _{j=0}^k\frac{{\rho }_j({\rho }_j+1)}{{({\theta }_h-a_j)}^2}+\sum _{\genfrac{}{}{0pt}{}{i,j=0}{i\ne j}}^k\frac{{\rho }_i{\rho }_j}{({\theta }_h-a_i)({\theta }_h-a_j)}\right).\hfill \end{array}$$

Applying Equation (37), we simplify to

$${\int }_{{\mathbf{S}}^k}w(x)G_h{(\mathbf{x})}^{2}dS(\mathbf{x})=\frac{e_h^2}{\rho (\rho +1)}\sum _{j=0}^k\frac{{\rho }_j}{{({\theta }_h-a_j)}^2}.$$

(40)

We introduce the polynomials

$$f(t):=\prod _{i=0}^k(t-a_i),g(t):=\prod _{i=1}^k(t-{\theta }_i).$$

Then, we have

$$g(t)=\left(\frac{1}{\rho }\sum _{j=0}^k\frac{{\rho }_j}{t-a_j}\right)f(t).$$

(41)

In order to prove Equation (41), we note that both sides are polynomials in t of degree k with leading coefficient 1, and, by Equation (37), both sides have the same zeros ${\theta }_1,\cdots ,{\theta }_k$. Therefore, we have the partial fraction expansion

$$\frac{g(t)}{f(t)}=\frac{1}{\rho }\sum _{j=0}^k\frac{{\rho }_j}{t-a_j}.$$

(42)

We first differentiate Equation (42), with respect to t, and then substitute $t={\theta }_h$ to write Equation (40) as

$${\int }_{{\mathbf{S}}^k}w(\mathbf{x})G_h{(\mathbf{x})}^{2}dS(\mathbf{x})=-\frac{e_h^2}{\rho +1}\frac{g^{\prime }({\theta }_h)}{f({\theta }_h)}\mathrm{for}h=1,2,\cdots ,k.$$

(43)

By Equation (22), the integral on the left-hand side is equal to 1, so

$$e_h^2=-(\rho +1)\frac{f({\theta }_h)}{g^{\prime }({\theta }_h)}\mathrm{for}h=1,2,\cdots ,k.$$

(44)

It follows from Equation (39) and $d_h>0$ that ${(-1)}^{h+1}{e}_h>0$ for every $h=1,2,\cdots ,k$. Therefore, Equations (38) and (44) give

$$G_h(\mathbf{x})={(-1)}^{h+1}{\left((\rho +1)\frac{|f({\theta }_h)|}{|g^{\prime }({\theta }_h)|}\right)}^{1/2}\sum _{j=0}^k\frac{x_j^2}{{\theta }_h-a_j}\mathrm{for}h=1,2,\cdots ,k.$$

(45)

From Equation (39), we find

$$d_h={\left(\frac{|e_h|}{|f({\theta }_h)|}\right)}^{1/k}\mathrm{for}h=1,2,\cdots ,k,$$

so, by Equation (6),

$$E_h(t)={\left(\frac{\rho +1}{|f({\theta }_h)||g^{\prime }({\theta }_h)|}\right)}^{1/(2k)}(t-{\theta }_h)\mathrm{for}h=1,2,\cdots ,k.$$

(46)

If we replace the constant term $-(\rho +1)$ on the right-hand side of Equation (30) by $-(\rho +1){\sum }_{j=0}^k{x}_j^2{\sum }_{j=0}^k{y}_j^2$ and compare coefficients, the addition formula Equation (31) with $m=2$ and $\mathbf{p}=\mathbf{0}$ holds, provided that

$$\frac{f^{\prime }(a_j)}{g(a_j)}-1=-\sum _{h=1}^k\frac{f({\theta }_h)}{g^{\prime }({\theta }_h)}\frac{1}{{({\theta }_h-a_j)}^2}\mathrm{for}j=0,1,\cdots ,k$$

(47)

and

$$-1=-\sum _{h=1}^k\frac{f({\theta }_h)}{g^{\prime }({\theta }_h)}\frac{1}{({\theta }_h-a_i)({\theta }_h-a_j)}\mathrm{for}0\le i\ne j\le k.$$

(48)

In fact, we may derive the formulas in Equations (47) and (48) from suitable partial fraction expansions. Therefore, combining both cases, we verified the addition formula Equation (28) for $m=2$, and the reader might imagine the formidable task that would be required to verify the addition formula in this way for general m.

If $k=2$, the addition formula Equation (32) can be written in a slightly different form using sphero-conal coordinates involving Jacobi elliptic functions ([9] Ch. 22]), see ([24] p. 24) and ([9] 29.18.2).

These coordinates have the advantage that they are valid on the entire sphere ${\mathbf{S}}^2$. We choose a modulus $\kappa \in (0,1)$ and set

$$\begin{array}{ccc}\hfill x_0& =& \kappa \mathrm{s}n{\beta }_1\mathrm{s}n{\beta }_2,\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill x_1& =& i\frac{\kappa }{{\kappa }^{\prime }}\mathrm{c}n{\beta }_1\mathrm{c}n{\beta }_2,\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill x_2& =& \frac{1}{{\kappa }^{\prime }}\mathrm{d}n{\beta }_1\mathrm{d}n{\beta }_2,\hfill \end{array}$$

(51)

where $\mathrm{s}nz=\mathrm{s}n(z,\kappa )$, $\mathrm{c}nz=\mathrm{c}n(z,\kappa )$, and $\mathrm{d}nz=\mathrm{d}n(z,\kappa )$ denote Jacobi elliptic functions corresponding to the modulus $\kappa$ and ${\kappa }^{\prime }=\sqrt{1-{\kappa }^2}$ is the complementary modulus. If $K=K(\kappa )$ and $K^{\prime }=K^{\prime }(\kappa )$ denote elliptic integrals ([9] Ch. 19), then we obtain a coordinate system in ${\mathbf{S}}^2$ by letting ${\beta }_1$ vary in $[0,4K]$ and ${\beta }_2^{\prime }$ in $[0,2K^{\prime }]$, where ${\beta }_2=K+i{\beta }_2^{\prime }$. These coordinates are connected to sphero-conal coordinates in algebraic form by $s_1={\mathrm{s}n}^2{\beta }_1$, $s_2={\mathrm{s}n}^2{\beta }_2$ when ${\beta }_1\in (0,K)$, ${\beta }_2^{\prime }\in (0,K^{\prime })$, with $a_0=0$, $a_1=1$, $a_2={\kappa }^{-2}>1$.

If ${\beta }_1\in (0,K)$, then

$$|{\mathrm{s}n}^2{\beta }_1-a_j{|}^{1/2}=\left\{\begin{array}{cc}\mathrm{s}n{\beta }_1\hfill & \mathrm{if}j=0,\hfill \\ \mathrm{c}n{\beta }_1\hfill & \mathrm{if}j=1,\hfill \\ \frac{1}{\kappa }\mathrm{d}n{\beta }_1\hfill & \mathrm{if}j=2,\hfill \end{array}\right.$$

so $E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_1)$ becomes a polynomial in $\mathrm{s}n{\beta }_1$, $\mathrm{c}n{\beta }_1$, $\mathrm{d}n{\beta }_1$ of degree $2|\mathbf{n}|+|\mathbf{p}|$. We define $E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_1)$ to be this polynomial for all values of ${\beta }_1$ allowing for a slight abuse of notation. Similarly, if ${\beta }_2^{\prime }\in (0,K^{\prime })$ then

$$|{\mathrm{s}n}^2{\beta }_2-a_j{|}^{1/2}=\left\{\begin{array}{cc}\mathrm{s}n{\beta }_2\hfill & \mathrm{if}j=0,\hfill \\ i\mathrm{c}n{\beta }_2\hfill & \mathrm{if}j=1,\hfill \\ \frac{1}{\kappa }\mathrm{d}n{\beta }_2\hfill & \mathrm{if}j=2,\hfill \end{array}\right.$$

so $E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_2)$ becomes a polynomial in $\mathrm{s}n{\beta }_2$, $\mathrm{c}n{\beta }_2$, $\mathrm{d}n{\beta }_2$ of degree $2|\mathbf{n}|+|\mathbf{p}|$.

If $\mathbf{x}=(x_0,x_1,x_2)\in {\mathbf{S}}^2$, $\mathbf{y}=(y_0,y_1,y_2)\in {\mathbf{S}}^2$ are represented by coordinates $({\beta }_1,{\beta }_2)$, $({\gamma }_1,{\gamma }_2)$, respectively, we have

$$G_{\mathbf{n},\mathbf{p}}(\mathbf{x},\mathbf{y})=E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_1)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_2).$$

(52)

If we set

$${\tilde{L}}_m({\beta }_1,{\beta }_2,{\gamma }_1,{\gamma }_2)=K_m(\mathbf{x}({\beta }_1,{\beta }_2),\mathbf{y}({\gamma }_1,{\gamma }_2)),$$

(53)

then we obtain the addition formula Equation (32) in the form

$${\tilde{L}}_m({\beta }_1,{\beta }_2,{\gamma }_1,{\gamma }_2)=\sum _{2|\mathbf{n}|+|\mathbf{p}|=m}{E}_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_1)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_2)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\gamma }_1)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\gamma }_2).$$

(54)

In particular, if we choose ${\alpha }_0={\alpha }_1={\alpha }_2=0$ then Equations (3) and (53) imply

$${\tilde{L}}_m({\beta }_1,{\beta }_2,{\gamma }_1,{\gamma }_2)=(2m+1)P_m(H),$$

(55)

where $P_m$ is the Legendre polynomial of degree m and

$$H={\kappa }^2\mathrm{s}n{\beta }_1\mathrm{s}n{\beta }_2\mathrm{s}n{\gamma }_1\mathrm{s}n{\gamma }_2-\frac{{\kappa }^2}{{\kappa }^{\prime 2}}\mathrm{c}n{\beta }_1\mathrm{c}n{\beta }_2\mathrm{c}n{\gamma }_1\mathrm{c}n{\gamma }_2+\frac{1}{{\kappa }^{\prime 2}}\mathrm{d}n{\beta }_1\mathrm{d}n{\beta }_2\mathrm{d}n{\gamma }_1\mathrm{d}n{\gamma }_2$$

and Equation (54) reduces to the addition formula for Lamé quasi-polynomials given by Hobson ([13] p. 475):

$$(2m+1)P_m(H)=\sum _{2|\mathbf{n}|+|\mathbf{p}|=m}{E}_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_1)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_2)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\gamma }_1)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\gamma }_2).$$

(56)

5. A Numerical Example

Heine–Stieltjes quasi-polynomials $E_{\mathbf{n},\mathbf{p}}$ do not admit an explicit representation. However, they can be computed numerically, then we can also verify the addition formula numerically in some examples. We consider the example $k=2$, $m=4$, and

$$a_0=0,a_1=1,a_2=3,{\alpha }_0=\frac{1}{2},{\alpha }_1=\frac{1}{3},{\alpha }_2=\frac{4}{3}.$$

The space ${\mathcal{H}}_4$ has dimension 9. To calculate the 9 corresponding Heine–Stieltjes quasi-polynomials $E_{\mathbf{n},\mathbf{p}}$, we start by finding all $\mathbf{n}=(n_1,n_2)$ and $\mathbf{p}=(p_0,p_1,p_2)$ that satisfy $2|\mathbf{n}|+|\mathbf{p}|=4$. There are three pairs $\mathbf{n},\mathbf{p}$ with $|\mathbf{n}|=2$, $|\mathbf{p}|=0$ and six pairs with $|\mathbf{n}|=1,|\mathbf{p}|=2$. These pairs are listed explicitly in

Table 1. Zeros and norming constants for Heine–Stieltjes polynomials.

n	p	${\mathit{\theta }}_1$	${\mathit{\theta }}_2$	${\mathit{c}}_{\mathbf{n},\mathbf{p}}$	${\mathit{d}}_{\mathbf{n},\mathbf{p}}$
(2,0)	(0,0,0)	0.1962642033	0.7837291906	0.1661415463	2.3125743278
(1,1)	(0,0,0)	0.3317225234	2.1827247789	−1.2479976882	1.7110564475
(0,2)	(0,0,0)	1.3576619901	2.3831914313	1.6215061837	1.3837750463
(1,0)	(1,1,0)	0.4745185575		1.5425185368	2.1309937152
(0,1)	(1,1,0)	2.2313637955		−5.1731280029	1.5747131736
(1,0)	(1,0,1)	0.6061553109		2.4246057510	1.5175072801
(0,1)	(1,0,1)	1.7467858656		−6.9358117249	1.1668545986
(1,0)	(0,1,1)	0.2790107946		1.8961260741	1.6279553075
(0,1)	(0,1,1)	1.8974597937		−6.5038818353	1.1962339394

. The next step is to solve Equation (9) numerically. If $|\mathbf{n}|=2$ there will be two equations for ${\theta }_1,{\theta }_2$ and three pairs of solutions—the first $a_0<{\theta }_1<{\theta }_2<a_1$ for $\mathbf{n}=(2,0)$, the second $a_0<{\theta }_1<a_1<{\theta }_2<a_2$ for $\mathbf{n}=(1,1)$, and the third $a_1<{\theta }_1<{\theta }_2<a_2$ for $\mathbf{n}=(0,2)$. The equations can be solved easily with mathematical software (we use Maple). The three different solutions pairs can be found by choosing suitable initial values for the unknowns ${\theta }_1,{\theta }_2$. If $|\mathbf{n}|=1$, there is only one equation, Equation (9), that we can also solve numerically (or exactly). Our results are listed in Table 1. The table also contains the values of the constants $c_{\mathbf{n},\mathbf{p}}$, $d_{\mathbf{n},\mathbf{p}}$. The computation of $c_{\mathbf{n},\mathbf{p}}$ from Equation (21) is trivial. The computation of $d_{\mathbf{n},\mathbf{p}}$ is a little more involved but also easy. We are using Equations (20), (22) and (23).

Using Table 1 we can numerically evaluate the right-hand side of the addition formula in Equation (28) (or Equation (32)). The left-hand side of Equation (28) is given by

$$\begin{array}{cc}\hfill K_4(\mathbf{x},\mathbf{y})& =\frac{1}{2}\rho (\rho +3)-\rho (\rho +1)(\rho +3)\left[\sum _{i=0}^2\frac{x_i^2{y}_i^2}{{\rho }_i}+\sum _{i=0}^1\sum _{j=i+1}^2\frac{x_i{y}_i{x}_j{y}_j}{{\rho }_i{\rho }_j}\right]\hfill \\ & +\rho (\rho +1)(\rho +2)(\rho +3)[\frac{1}{2}\sum _{i=0}^2\frac{x_i^4{y}_i^4}{{\rho }_i({\rho }_i+1)}+\sum _{\genfrac{}{}{0pt}{}{i,j=0}{i\ne j}}^2\frac{x_i^3{y}_i^3{x}_j{y}_j}{{\rho }_i({\rho }_i+1){\rho }_j}\hfill \\ & +\sum _{i=0}^1\sum _{j=i+1}^2\frac{x_i^2{y}_i^2{x}_j^2{y}_j^2}{{\rho }_i{\rho }_j}+\frac{x_0{y}_0{x}_1{y}_1{x}_2{y}_2}{{\rho }_0{\rho }_1{\rho }_2}\sum _{i=0}^2{x}_i{y}_i].\hfill \end{array}$$

It should be noted that the kernel $K_4$ is given explicitly. Choosing various points $\mathbf{x},\mathbf{y}\in {\mathbf{S}}^2$, one can verify numerically that the addition formula Equation (28) holds. Of course, there will be a small difference between the left-hand and right-hand sides of Equation (28) because the right-hand side cannot be calculated exactly.

6. Integral Equations

By Equations (1) and (13), we have, for $f\in {\mathcal{H}}_m$,

$$f(\mathbf{x})={\int }_{{\mathbf{S}}^k}w(\mathbf{y})K_m(\mathbf{x},\mathbf{y})f(\mathbf{y})dS(\mathbf{y})\mathrm{for}\mathbf{x}\in {\mathbf{S}}^k,$$

where $w(\mathbf{y})$ is defined by Equation (14). Since $G_{\mathbf{n},\mathbf{p}}\in {\mathcal{H}}_m$, when $m=2|\mathbf{n}|+|\mathbf{p}|$, we obtain the following integral equation for generalized spherical harmonics.

Theorem 4.

If $m=2|\mathbf{n}|+|\mathbf{p}|$, then

$$G_{\mathbf{n},\mathbf{p}}(\mathbf{x})={\int }_{{\mathbf{S}}^k}w(\mathbf{y})K_m(\mathbf{x},\mathbf{y})G_{\mathbf{n},\mathbf{p}}(\mathbf{y})dS(\mathbf{y})for\mathbf{x}\in {\mathbf{S}}^k,$$

(57)

which holds regardless of whether $G_{\mathbf{n},\mathbf{p}}$ is normalized or not.

Noting that

$$K_{m,\mathbf{p}}(\mathbf{x},{\sigma }_j\mathbf{y})G_{\mathbf{n},\mathbf{p}}({\sigma }_j\mathbf{y})=K_{m,\mathbf{p}}(\mathbf{x},\mathbf{y})G_{\mathbf{n},\mathbf{p}}(\mathbf{y})\mathrm{for}\mathrm{all}j=0,1,\cdots ,k,$$

with the reflection ${\sigma }_j$ defined by Equation (12) and $G_{\mathbf{n},\mathbf{p}}\in {\mathcal{H}}_{m,\mathbf{p}}$ when $m=2|\mathbf{n}|+|\mathbf{p}|$, it follows that

$$G_{\mathbf{n},\mathbf{p}}(\mathbf{x})=2^{k+1}{\int }_{{\mathbf{S}}_{+}^k}w(\mathbf{y})K_{m,\mathbf{p}}(\mathbf{x},\mathbf{y})G_{\mathbf{n},\mathbf{p}}(\mathbf{y})dS(\mathbf{y})$$

(58)

with $S_{+}^k$ defined in Equation (15). We transform the integral on the right-hand side of Equation (58) to sphero-conal coordinates $t_1,t_2,\cdots ,t_k$ by using that

$$2^{k+1}w(\mathbf{y}(t))dS(\mathbf{y}(t))=v(\mathbf{t})d\mathbf{t},$$

(59)

where

$$v(\mathbf{t}):=\frac{\Gamma (\rho )}{{\prod }_{j=0}^k\Gamma ({\rho }_j)}\prod _{j=0}^k|A_j{|}^{-{\alpha }_j}\prod _{j=0}^k\prod _{i=1}^k{|t_i-a_j|}^{{\alpha }_j-1/2}\prod _{1\le i<j\le k}(t_j-t_i).$$

This transformation formula is derived from the metric tensor for sphero-conal coordinates given in ([22] Equation (18)). Transforming the integral in Equation (58) to sphero-conal coordinates and using Equation (19) we obtain the following integral equation for Heine–Stieltjes quasi-polynomials.

Theorem 5.

If $m=2|\mathbf{n}|+|\mathbf{p}|$, then

$$\prod _{j=1}^k{E}_{\mathbf{n},\mathbf{p}}(s_j)={\int }_{Q}v(\mathbf{t})L_{m,\mathbf{p}}(\mathbf{s},\mathbf{t})\prod _{j=1}^k{E}_{\mathbf{n},\mathbf{p}}(t_j)d\mathbf{t}for\mathbf{s}\in Q$$

(60)

with Q from Equation (16).

As a special case, consider $k=2$, $\mathbf{p}=(0,0,0)$; and $a_0=0$, $a_1=1$, and $a_2=a>1$. Then, the polynomial $E_{\mathbf{n},\mathbf{0}}$ with ${\theta }_l$ from Equation (9) is called a Heun polynomial ([12] A 3.6). For example, for $m=2$ we employ the kernel

$$\begin{array}{cc}\hfill L_{2,\mathbf{0}}(s_1,s_2,t_1,t_2)& =\frac{\rho (\rho +1)s_1{s}_2{t}_1{t}_2}{{\rho }_0{a}^2}+\hfill \\ & +\frac{\rho (\rho +1)(s_1-1)(s_2-1)(t_1-1)(t_2-1)}{{\rho }_1{(1-a)}^2}+\hfill \\ & +\frac{\rho (\rho +1)(s_1-a)(s_2-a)(t_1-a)(t_2-a)}{{\rho }_2{a}^2{(1-a)}^2}-(\rho +1).\hfill \end{array}$$

(61)

Sleeman ([15] (4.35)) gave an integral equation for Heun polynomials in the form

$$E_{\mathbf{n},\mathbf{0}}=\lambda {\int }_{Q}K(u,t_1,t_2)E_{\mathbf{n},\mathbf{0}}(t_1)E_{\mathbf{n},\mathbf{0}}(t_2)d{t}_{1}d{t}_2.$$

(62)

This integral equation is similar to Equation (60) with $k=2$. The difference is that Equation (62) contains an eigenvalue $\lambda$ on the right-hand side and only one Heun polynomial on the left-hand side. The kernel functions are different because we use sphero-conal coordinates in ${\mathbf{S}}^2$ while Sleeman ([15] (4.1)) uses ellipsoidal coordinates in ${\mathbb{R}}^3$.

If $k=2$, we may also employ the kernel ${\tilde{L}}_m$ introduced in Section 4. This allows us to use the same kernel for all Heun quasi-polynomials $E_{\mathbf{n},\mathbf{p}}$ with $m=2|\mathbf{n}|+|\mathbf{p}|$. Using the sphero-conal coordinates in the form Equations (49)–(51), we obtain

$$w(\mathbf{y})dS(\mathbf{y})=\tilde{v}({\gamma }_1,{\gamma }_2)d{\gamma }_{1}d{\gamma }_2^{\prime },$$

where

$$\begin{array}{cc}\hfill \tilde{v}({\gamma }_1,{\gamma }_2):& =\frac{\Gamma (\rho )}{2\Gamma ({\rho }_0)\Gamma ({\rho }_1)\Gamma ({\rho }_2)}\times \hfill \\ & {\left|\kappa \mathrm{s}n{\gamma }_1\mathrm{s}n{\gamma }_2\right|}^{2{\alpha }_0}{\left|\frac{\kappa }{{\kappa }^{\prime }}\mathrm{c}n{\gamma }_1\mathrm{c}n{\gamma }_2\right|}^{2{\alpha }_1}{\left|\frac{1}{{\kappa }^{\prime }}\mathrm{d}n{\gamma }_1\mathrm{d}n{\gamma }_2\right|}^{2{\alpha }_2}{\kappa }^2({\mathrm{s}n}^2{\gamma }_2-{\mathrm{s}n}^2{\gamma }_1).\hfill \end{array}$$

Then, Theorem 4 and Equation (52) yield the following integral equation for Heun quasi-polynomials in Jacobi form.

Theorem 6.

Let $k=2$ and $m=2|\mathbf{n}|+|\mathbf{p}|$. Then,

$$\begin{array}{cc}\hfill & E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_1)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\beta }_2)\hfill \\ =& {\int }_0^{2K^{\prime }}{\int }_0^{4K}\tilde{v}({\gamma }_1,{\gamma }_2){\tilde{L}}_m({\beta }_1,{\beta }_2,{\gamma }_1,{\gamma }_2)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\gamma }_1)E_{\mathbf{n},\mathbf{p}}({\mathrm{s}n}^2{\gamma }_2)d{\gamma }_{1}d{\gamma }_2^{\prime }.\hfill \end{array}$$

(63)

In the special case ${\alpha }_0={\alpha }_1={\alpha }_2=0$, this is an integral equation for Lamé quasi-polynomials ([1] Ch. IX) involving

$$\tilde{v}({\gamma }_1,{\gamma }_2)=\frac{{\kappa }^2}{4\pi }({\mathrm{s}n}^2{\gamma }_2-{\mathrm{s}n}^2{\gamma }_1)$$

(64)

and the kernel ${\tilde{L}}_m$ given in Equation (55).

Arscott ([14] (5.6)) obtained a similar integral equation for Lamé polynomials, but of the type Equation (62). Comparison of the kernels shows that Arscott’s kernel is more complicated than ${\tilde{L}}_m$. Moreover, our integral equation does not involve unknown constants. It is interesting to note that the kernel ${\tilde{L}}_m$ from Equation (55) also appears in linear integral equation for Lamé polynomials. In this setting, the kernel enters the theory as the Riemann function of a partial differential equation, see [25].

7. Conclusions

We showed in Theorem 3 that Heine–Stieltjes quasi-polynomials satisfy an addition formula. The proof is simple because it is based on two known results, namely, Theorem 1 and Equation (3.1) for a reproducing kernel. However, it is a nontrivial task to define Heine–Stieltjes quasi-polynomials and the corresponding generalized sphero-conal harmonics. The addition formula is one of very few results of an explicit nature in the theory of Heine Stieltjes quasi-polynomials. We verified the correctness of the formula in special cases exactly ($m=2$) and numerically ($m=4$). As a mathematical application, we found integral relations for Heine–Stieltjes polynomials that appear to be new.