# Ultraspherical/Gegenbauer polynomials to
unify 2D/3D Ambisonic directivity designs

Franz Zotter

###### Abstract

This report on axisymmetric ultraspherical/Gegenbauer polynomials [Geg77] and
their use in Ambisonic directivity design in 2D and 3D presents an alternative
mathematical formalism to what can be read in, e.g., my and Matthias Frank’s
book on Ambisonics or Jérôme Daniel’s thesis, Gary Elko’s differential array
book chapters, or Boaz Rafaely’s spherical microphone array book.

My original—but discarded—intention was to already include ultraspherical
polynomials in the Ambisonics book to pursue a simple thought: retrieving
suitable axisymmetric continuous-direction functions that serve the discrete-
direction metrics for experimental psychoacoustic findings in spatial audio
uniformly, for any number $\mathrm{D}$ of space dimensions, of course
typically $\mathrm{D}=2,3$. But the cost for this natural way to understand
spherical/circular polynomials and their directional sampling was too high.
This would have required carrying out both the derivations of two entire
formalisms from scratch, the ultraspherical/Gegenbauer polynomials and the
circular/spherical harmonics, of which only the latter ones define the
Ambisonic format.

Nevertheless, ultraspherical/Gegenbauer polynomials are highly valuable when
designing axisymmetric beams and understanding spherical $t$ designs so that
this report will shed some light on what circular, spherical, and
ultraspherical axisymmetric polynomials are. While mathematically interesting
by themselves already, they can be useful in spherical beamforming as
described in the literature on spherical and differential microphone arrays,
e.g. [Raf19, Elk00, Elk04]. What are polynomials bases for axisymmetric
functions, say axisymmetric harmonics (or ultraspherical or Gegenbauer)
polynomials? How can they be utilized to define Ambisonic order weightings or
spherical beamformers?

In this report, these ultraspherical/Gegenbauer polynomials will be used to
uniformly derive for arbitrary dimensions $\mathrm{D}$ the various directivity
designs or Ambisonic order weightings known from literature: max-DI/basic
[Dan01], max-$r_{\mathrm{E}}$ [Dan01], supercardioid [Raf19, Elk00, Elk04],
cardioid/inphase [Dan01]. Is there a way to relate higher-order cardioids and
supercardioids? How could one define directivity patterns with an on-axis
flatness constraint?

###### Contents

  1. 1 Axisymmetric directivity functions
    1. 1.1 Axisymmetric harmonics: Gegenbauer polynomials
    2. 1.2 Continuous axisymmetric directivity/panning patterns in D dimensions
    3. 1.3 Optimal discretization: spherical $t$ design
    4. 1.4 Continuous circular patterns in Chebyshev polynomials
    5. 1.5 Addition theorem of the circular harmonics
    6. 1.6 Optimal circular discretization: circular $t$ design
    7. 1.7 Continuous spherical patterns in Legendre polynomials
    8. 1.8 Addition theorem of the spherical harmonics
    9. 1.9 Optimal spherical discretization: spherical $t$-design
  2. A Laplacian in axisymmetric coordinates
  3. B Jacobi determinant of axisymmetric coordinates
  4. C Solving the ultraspherical differential equation: Gegenbauer polynomials
  5. D Orthogonality of axisymmetric harmonics
  6. E Orthogonal polynomials
    1. E.1 The transform integral: a benefit of orthogonality
    2. E.2 Symetries and interval $-1\leq x\leq 1$.
    3. E.3 Three-term recurrence of orthogonal polynomials
    4. E.4 Explicit three-term recurrence
    5. E.5 Squared norm of orthogonal polynomials
    6. E.6 Christoffel-Darboux formula
    7. E.7 Polynomials $C_{n}^{(\alpha)}(x)$ and $\mathcal{P}_{n}(x)$ at $x=1$
    8. E.8 Derivative of the polynomials $\mathcal{P}_{n}(x)$ at $x=1$
    9. E.9 Value of the polynomials $C^{(\alpha)}_{n}(x)$ and $\mathcal{P}_{n}(x)$ at $x=0$
    10. E.10 Rodrigues’ formula for Gegenbauer polynomials.
    11. E.11 Norm of Gegenbauer polynomials
    12. E.12 Integral recurrence of the Gegenbauer polynomials
    13. E.13 Additional recurrence relations
  7. F Designs and metrics $P$, $E$, $Q$, $FBR$, $\bm{r}_{\mathrm{V}}$, and $\bm{r}_{\mathrm{E}}$ for axisymmetric directivity/panning patterns
    1. F.1 Supercardioid: Gram matrix for half interval
    2. F.2 Basic pattern with maximum DI/maximum $\bm{r}_{\mathrm{V}}$ vector
    3. F.3 Pattern with maximum $\bm{r}_{\mathrm{E}}$ vector
    4. F.4 Particular max-$r_{\mathrm{E}}$ solutions for $\mathrm{D}=3$
    5. F.5 (Ultra-)Spherical cap function
    6. F.6 Inphase/higher-order cardioid pattern
    7. F.7 Max-flat pattern: generalized higher-order cardioid
  8. G MATLAB Code examples
    1. G.1 gegenbauer_me(x,N,D)
    2. G.2 chebyshev(x,N)
    3. G.3 legendre_u(x,N)
    4. G.4 maxre(N,D)
    5. G.5 one_over_gegenbauer_norm_square(N,D)
    6. G.6 supercard(N)
    7. G.7 inphase(N)
    8. G.8 inphase_maxflat(N,L,D)
    9. G.9 spherecap(N,D)
    10. G.10 eval_PEQrVrEFBR(a,D)
    11. G.11 test_PEQrVrEFBR(a,D)
  9. H 2D Ambisonic panning as Fourier series
    1. H.1 Panning function as Fourier series
    2. H.2 Finite-order Fourier series
    3. H.3 Finite-order series with maximum $\bm{r}_{\mathrm{E}}$ vector length
    4. H.4 Optimally arranged surround loudspeaker ring

##  1 Axisymmetric directivity functions

In the seventies, Cooper [CS72] described the idea of using angular amplitude
functions to control the auditory event, which could be thought of as an
extension to mid/side decomposition, or as recording with figure-of-eight and
omnidirectional directivity patterns. Stemming from this approach and from
this time, the notion and technology of Ambisonics was developed by Felgett
[Fel74], Gerzon [Ger75], and Craven [CG77]. Since then, the term Ambisonics
became common for technology using angular harmonic functions to express
surround sound. Read more about all of this in the open Ambisonics book at
Springer https://doi.org/10.1007/978-3-030-17207-7.

Appendix H shows how 2D Ambisonic panning can be expressed in terms of a
finite-order Fourier series using the exponentials $e^{\mathrm{i}m\phi}$ with
real coefficients, i.e., a finite cosine $\cos m\phi$ series in the azimuthal
distance $\phi=\varphi-\varphi_{s}$ between array angles
$\varphi=\upvarphi_{l}$ and the given angle $\varphi_{s}$. For this symmetric
finite-resolution function, a regular circular polygon of sufficiently many
loudspeaker directions would automatically yield correct aiming of the vector
measures $\bm{r}_{\mathrm{E}}$ and $\bm{r}_{\mathrm{V}}$ towards the panning
direction $\bm{\theta}_{\mathrm{s}}$ and panning-invariant loudness and width
magnitudes $E$, $P$, $\|\bm{r}_{\mathrm{E}}\|$, $\|\bm{r}_{\mathrm{V}}\|$, and
a certain directivity factor $Q$.

Instead, we desire defining axisymmetric functions around the panning
direction $\bm{\theta}_{\mathrm{s}}$ in a space of arbitrarily many dimensions
$\mathrm{D}$, and their decomposition into suitable limited degree harmonic
functions. The projection of the variable direction vector $\bm{\theta}$ on
given direction $\bm{\theta}_{\mathrm{s}}$ yields the cosine of the enclosed
angle to the axis $\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}=\cos\phi$,
and hereby expresses the only dimension of the axisymmetric function,
$x=\cos\phi$. As the $m^{\mathrm{th}}$ power thereof, $\cos^{m}\phi$, can be
expressed in terms of a cosine series up to the $m^{\mathrm{th}}$ degree, a
power series in $x$ is a reasonable representation of the panning functon on
the sphere.

To do this properly, we will first learn about the term harmonics, solve for
axisymmetric harmonics, establish the panning function on the circle using
these harmonics, and then finally come to the spherical panning function.

###  1.1 Axisymmetric harmonics: Gegenbauer polynomials

The Laplacian is defined in the $\mathrm{D}$-dimensional Cartesian space as

$\displaystyle\bigtriangleup$
$\displaystyle=\sum_{j=1}^{\mathrm{D}}\frac{\partial^{2}}{\partial
x_{j}^{2}}.$ (1)

Solving the eigenvalue problem $\bigtriangleup=-m^{2}\,f=0$ for one dimension
$\frac{\partial^{2}}{\partial x^{2}}f=-m^{2}f$ yields with the exponential
guess $f=e^{ax}$, inserted and differentiated to $a^{2}\,f=-m^{2}\,f$,
_harmonic functions_ $f=e^{\pm\mathrm{i}mx}$ with $m\in\mathbb{Z}$ if the
function should be periodic $f(x)=f(x+2\pi\,k)$, $k\in\mathbb{Z}$. If only a
limited number of eigenvalues $m\leq K$ is used, this yields a limited-degree
expansion into _harmonics_ , a band-limited Fourier series. It is useful to
limit the resolution of a represented function.

For 2 dimensions Appendix H defines panning function with limited $m$, defined
using a series of real part of this complex exponential
$e^{\mathrm{i}m\phi}=\cos(m\phi)$ with the polar angle
$\phi=\arctan\frac{y}{x}$ enclosed with the $x$ axis. For 3 dimensions, these
trigonometric or complex-exponential angular harmonics are not useful anymore,
and lose some beneficial properties, for instance, $\cos(m\phi)$ has no
uniformly limited resolution and is not analytic around the poles of the
sphere.

Generally, homogeneous solutions to the Laplacian $\bigtriangleup g=0$ are
called harmonics. In case of a spherical coordinate system with a unit-
direction vector $\bm{\theta}$ and a radius $r$ and a separation approach
$g=f(\bm{\theta})\,h(r)$, we obtain $\bigtriangleup
g=f(\bm{\theta})\bigtriangleup h(r)+h(r)\bigtriangleup f(\bm{\theta})=0$, for
which there needs to be a separation constant $\lambda=\frac{\bigtriangleup
h}{h}$ defining Laplacian eigensolutions $\bigtriangleup f=-\lambda\,f$ for
the directional part. For suitable Eigenvalues $-\lambda$, these _directional
harmonics_ span an orthogonal set of basis functions as solutions.

Chosing the last Cartesian coordinate axis as panning direction
$\bm{\theta}_{\mathrm{s}}=[\dots,\,0,\,1]^{\mathrm{T}}$, and writing the
direction vector as
$\bm{\theta}=\frac{[x_{1},\,\dots,\,x_{\mathrm{D}}]^{\mathrm{T}}}{\sqrt{\sum_{j}x_{j}^{2}}}=\frac{[x_{1},\,\dots,\,x_{\mathrm{D}}]^{\mathrm{T}}}{r}$,
i.e.
$\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}=\frac{x_{\mathrm{D}}}{r}=\cos\phi$,
a purely axisymmetric panning function defined on $x:=\cos\phi$ becomes

$\displaystyle f$ $\displaystyle=f(\cos\phi)=f(x).$

By this choice, the _directional harmonics_ become _axisymmetric harmonics_
with regard to the enclosed angle $\phi$ to some directional axis. To find
their formal solution, we re-write the Laplacian to the axial dimension
$-1\leq x\leq 1$. Such harmonics are useful to represent limited-resolution
patterns of rotational symmetry with regard to their axis. As derived in the
Appendix A, the Laplacian with regard to the axial coordinate $x=\cos\phi$
becomes:

$\displaystyle
r^{2}\bigtriangleup_{x}=(1-x^{2})\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}-(\mathrm{D}-1)\,x\,\frac{\mathrm{d}}{\mathrm{d}x}.$
(2)

#### The Ggegenbauer/ultraspherical differential equation: Laplacian
eigenvalue problem for axisymmetric functions.

With the above, we establish the ultraspherical differential equation
$(r^{2}\bigtriangleup_{x}+\lambda)f=0$ with the eigenvalue $-\lambda$ to the
Laplacian

$\displaystyle(1-x^{2})f^{\prime\prime}-(\mathrm{D}-1)\,x\,f^{\prime}+\lambda\,f=0,$
(3)

in the hope that if a limited set of discrete eigenvalues $\lambda_{n}\leq K$
is used, then we obtain harmonics orthogonal on the interval $-1\leq x\leq 1$
that contain only limited powers of $x^{n}$ with $n\leq\mathrm{N}$, given a
suitable of $K$ and $\lambda_{n}$. We know trigonometrically that such
functions can only contain oscillations up to $\cos(\mathrm{N}\phi)$, thus
have limited angular resolution.

We solve the equation with the power series approach $f=\sum_{k}c_{k}\,x^{k}$,
see Appendix C, and get an interrelation of its coefficients
$c_{k+2}=\frac{k^{2}+k(\mathrm{D}-2)-\lambda}{(k+1)(k+2)}\;c_{k}$. Choosing
$\lambda=\lambda_{n}=n^{2}+(\mathrm{D}-2)n$, the power series terminates at
$k=n$, so that $c_{k>n}=0$, hence it becomes an $n^{\mathrm{th}}$-degree
polynomial $f=\mathcal{P}_{n}$, and the relation between the coefficients of
$\mathcal{P}_{n}$ remains

$\displaystyle
c_{k+2}^{(n)}=\frac{(k-n)(k+n+\mathrm{D}-2)}{(k+1)(k+2)}\;c_{k}^{(n)}.$ (4)

Also solutions with different $\lambda$ could be acceptable, but their
infinite series may diverge, e.g., at $x=1$ or $x=-1$, so we disregard. We re-
write the above equation accordingly:

$\displaystyle(1-x^{2})\mathcal{P}_{n}^{\prime\prime}-(\mathrm{D}-1)\,x\,\mathcal{P}_{n}^{\prime}+n(n+\mathrm{D}-2)\,\mathcal{P}_{n}=0.$
(5)

The above equation is called the _ultraspherical_ or _Gegenbauer_ differential
equation [Geg77], and solutions can be constructed from
$\mathcal{P}_{0}(x)=1$, and $\mathcal{P}_{1}(x)=x$ with the recurrence derived
in subsection E.3, eq. (87) with the standardization $\mathcal{P}_{n}(1)=1$ of
the _ultraspherical_ polynomials, therefore differing from Gegenbauer’s
typical definition
$C_{n}^{(\alpha)}=\frac{(2\alpha)^{\overline{n}}}{n!}\,\mathcal{P}_{n}(x)$
[Geg77], cf. subsection E.7 eq. (97),

$\displaystyle\mathcal{P}_{n+1}(x)=\frac{2n+\mathrm{D}-2}{n+\mathrm{D}-2}\;x\,\mathcal{P}_{n}(x)-\frac{n}{n+\mathrm{D}-2}\mathcal{P}_{n-1}(x).$
(6)

Figure 1 shows these ultraspherical polynomials for $\mathrm{D}=2,3,4$,
$\alpha=\frac{D-2}{2}$. For $\mathrm{D}=2,3$, they become the Chebyshev and
the Legendre polynomials.

Figure 1: Ultraspherical polynomials
$\mathcal{P}_{n}(x)=\frac{n!}{(2\alpha)^{\overline{n}}}C_{n}^{(\alpha)}(x)\big{|}_{1}=1$
for $\mathrm{D}=2,3,4$, $\alpha=\frac{D-2}{2}$, i.e. Chebyshev polynomials
$T_{n}(x)$ for $\mathrm{D}=2$, and Legendre polynomials $P_{n}(x)$ for
$\mathrm{D}=3$. They are orthogonal wrt. the weight
$w(x)=\sqrt{1-x^{2}}^{\mathrm{D}-3}$.

#### Orthogonality.

The ultraspherical polynomials $\mathcal{P}_{n}$ are orthogonal according to
the integral
$\iiint\,\mathcal{P}_{n}\,\mathcal{P}_{m}\,\delta(r-1)\prod_{j}\mathrm{d}x_{j}=N_{n}^{2}\delta_{nm}$.
We re-write the integration elements to get the orthogonality relation in $x$,
see Appendix B, yielding

$\displaystyle\int_{-1}^{1}\mathcal{P}_{n}(x)\mathcal{P}_{m}(x)\,\sqrt{1-x^{2}}^{\mathrm{D}-3}\,\mathrm{d}x=N_{n}^{2}\,\delta_{nm},$
(7)

with the squared norm

$\displaystyle N_{n}^{2}$
$\displaystyle=\int_{-1}^{1}\mathcal{P}_{n}^{2}(x)\;\sqrt{1-x^{2}}^{\mathrm{D}-3}\,\mathrm{d}x=\frac{n!}{(2n+\mathrm{D}-2)(n+\mathrm{D}-3)!}\frac{S_{\mathrm{D}-1}}{S_{\mathrm{D}-2}}$
(8)

of the polynomial, see subsection E.11 eq. (109). Orthogonality can be shown
by integrating eq. (5) over $\mathcal{P}_{m}$ using the corresponding
integration element, and by integration by parts, as shown in Appendix D.

Orthogonality has the benefit of enabling a transform integral to find the
coefficients $\gamma_{n}$ of a function decomposed into the polynomials
$\mathcal{P}_{n}$, see Appendix E,

$\displaystyle g(x)$
$\displaystyle=\sum_{n=0}^{\infty}\gamma_{n}\,\mathcal{P}_{n}(x),$
$\displaystyle\gamma_{n}$
$\displaystyle=\frac{1}{N_{n}^{2}}\int_{-1}^{1}g(x)\,\mathcal{P}_{n}(x)\,\sqrt{1-x^{2}}^{\mathrm{D}-3}\,\mathrm{d}x.$
(9)

###  1.2 Continuous axisymmetric directivity/panning patterns in D dimensions

The axisymmetric directivity/panning function is based on the decomposition of
a Dirac delta pointing towards the panning direction
$\bm{\theta}_{\mathrm{s}}$,

$\displaystyle\delta(1-\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})$
$\displaystyle=\begin{cases}\infty,&\text{for
$\bm{\theta}=\bm{\theta}_{\mathrm{s}}$}\\\ 0,&\text{otherwise}\end{cases},$
$\displaystyle\text{with
}\int_{\mathbb{S}^{\mathrm{D}-1}}\delta(1-\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})\,\mathrm{d}\bm{\theta}$
$\displaystyle=1,$ (10)

into infinitely many axisymmetric harmonics

$\displaystyle\delta(1-\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})$
$\displaystyle=\frac{1}{S_{\mathrm{D}-2}}\sum_{n=0}^{\infty}\gamma_{n}\,\mathcal{P}_{n}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}),$
(11)

where the factor $\frac{1}{S_{\mathrm{D}-2}}$ accounts for normalization with
regard to the integrals over the independent, axisymmetric coordinates, whose
value corresponds to the surface of the next lesser dimensional unit sphere
$\mathbb{S}^{\mathrm{D}-2}$.

The desired axisymmetric panning function is based on the representation of
the Dirac delta, but with limited degree and involving a sidelobe-suppressing
weight $a_{n}$,

$\displaystyle g(\bm{\theta})$
$\displaystyle=\frac{1}{S_{\mathrm{D}-2}}\sum_{n=0}^{\mathrm{N}}a_{n}\,\gamma_{n}\,\mathcal{P}_{n}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}).$
(12)

#### Simple $\gamma_{n}$ for axisymmetric Dirac delta at $x=x_{0}$.

Given a Dirac delta located at $x=x_{0}$ as
$\frac{\delta(x_{0}-x)}{\sqrt{1-x^{2}}^{\mathrm{D}-3}}$ which is correctly
normalized to the axial coordinate to fulfill
$\int_{-1}^{1}\delta(1-x)\;\,\mathrm{d}x=1$, we use the transform integral eq.
(9) to find its coefficients $\gamma_{m}$:

$\displaystyle\frac{\delta(x_{0}-x)}{\sqrt{1-x^{2}}^{\mathrm{D}-3}}$
$\displaystyle=\sum_{n=0}^{\infty}\gamma_{n}\,\mathcal{P}_{n}(x),$
$\displaystyle\gamma_{n}=\frac{1}{N_{n}^{2}}\;\int_{-1}^{1}\delta(1-x)\,\mathcal{P}_{n}(x)\;\mathrm{d}x=\frac{\mathcal{P}_{n}(x_{0})}{N_{n}^{2}}.$
(13)

The resulting expression
$\frac{\delta(x_{0}-x)}{w(x)}=\sum_{n=0}^{\infty}\frac{\mathcal{P}_{n}(x_{0})\mathcal{P}_{n}(x)}{N_{n}^{2}}$
may be simplified by the Christoffel-Darboux identity eq. (94) for summation
over a limited degree $n\leq\mathrm{N}$.

When located at the pole $x_{0}=1$, the ultraspherical/Gegenbauer polynomials
normalized to $\mathcal{P}_{n}(1)=1$ yield a simple expansion coefficient, cf.
eq. (8), with $a_{n}=\mathcal{P}_{n}(x_{0})=1$,

$\displaystyle\gamma_{n}$ $\displaystyle=\frac{1}{N_{n}^{2}},$ $\displaystyle
g(x)$
$\displaystyle=\frac{1}{S_{\mathrm{D}-2}\,}\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\,\mathcal{P}_{n}(x),$
(14)

and $a_{n}=1$ results from eq. (13), while $a_{n}\neq 1$ can be chosen to
suppress sidelobes.

#### Metrics $P$, $E$, $Q$, $\bm{r}_{\mathrm{V}}$, $\bm{r}_{\mathrm{V}}$,
$\mathrm{FBR}$ for finite-degree axisymmetric patterns.

We still need to find $a_{n}$, and these criteria were considered helpful. As
derived in Appendix F, the metrics of the continuous, limited-degree
axisymmetric panning function can be found by solving integrals over the
continuous loudspeaker direction $\bm{\theta}$. The measures only depend on
the weight $a_{n}$; with the integration measure abbreviated as
$w(x)=\sqrt{1-x^{2}}^{\mathrm{D}-3}$,

$\displaystyle P$
$\displaystyle=S_{\mathrm{D}-2}\,\int_{-1}^{1}g(x)\,w(x)\mathrm{d}x=a_{0},$
(15) $\displaystyle E$
$\displaystyle=S_{\mathrm{D}-2}\int_{-1}^{1}g^{2}(x)\,w(x)\mathrm{d}x=\sum_{n=0}^{\mathrm{N}}\frac{|a_{n}|^{2}}{S_{\mathrm{D}-2}N_{n}^{2}},$
(16) $\displaystyle Q$
$\displaystyle=\frac{S_{D-1}\,g^{2}(1)}{E}=\frac{S_{D-1}\left(\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{S_{D-2}N_{n}^{2}}\right)^{2}}{E},$
(17) $\displaystyle\bm{r}_{\mathrm{V}}$
$\displaystyle=\frac{\int_{-1}^{1}g(x)\,x\;w(x)\mathrm{d}x}{P}\,\bm{\theta}_{\mathrm{s}}=\frac{a_{1}}{a_{0}}\,\bm{\theta}_{\mathrm{s}},$
(18) $\displaystyle\bm{r}_{\mathrm{E}}$
$\displaystyle=\frac{\int_{-1}^{1}g^{2}(x)\,x\;w(x)\mathrm{d}x}{E}\,\bm{\theta}_{\mathrm{s}}=\frac{\sum_{n=0}^{\mathrm{N}}\frac{2\beta_{n+1}}{N_{n}^{2}}\,a_{n}a_{n+1}}{\sum_{n=0}^{\mathrm{N}}\frac{1}{N_{n}^{2}}|a_{n}|^{2}}\,\bm{\theta}_{\mathrm{s}}.$
(19) $\displaystyle\mathrm{FBR}$
$\displaystyle=\textstyle\frac{\int_{0}^{1}g^{2}(x)\,w(x)\,\mathrm{d}x}{\int_{-1}^{0}g^{2}(x)\,w(x)\,\mathrm{d}x}=\frac{\sum_{n=0}^{\mathrm{N}}\sum_{m=0}^{\mathrm{N}}a_{n}g_{nm}a_{m}}{\sum_{n=0}^{\mathrm{N}}\sum_{m=0}^{\mathrm{N}}a_{n}g_{nm}a_{m}(-1)^{m}}.$
(20)

Obviously, $P$ and $\bm{r}_{\mathrm{V}}$ only depend on the weights $a_{0}$
and $a_{1}$, while $E$, $Q$, and $\bm{r}_{\mathrm{E}}$ involve the weights
$a_{n}$ of all the degrees; integrals of the FBR require the entries of the
Gram matrix $g_{nm}$ as derived in subsection F.1.

#### Basic/maximum-DI weights.

With the weights $a_{n}=1$, there is no sidelobe suppression but the narrowest
main lobe and highest directivity that can be represented using the limited
degrees $n\leq\mathrm{N}$, cf. subsection F.2. For this choice, the
directivity factor is highest
$Q=\sum_{n=0}^{\mathrm{N}}\frac{1}{S_{D-2}N_{n}^{2}}$, which is
$\sum_{n=0}^{\mathrm{N}}(2-\delta_{m})$ or $\sum_{n=0}^{\mathrm{N}}(2n+1)$ and
so $Q=\\{2\mathrm{N}+1,\,(\mathrm{N}+1)^{2}\\}$ for $D=\\{2,3\\}$. Because
there is an equivalent least-squares problem
$\int_{-1}^{1}[\delta(1-x)-g(x)]^{2}w(x)\mathrm{d}x\rightarrow\min$ to the
defining eq. (14) with neutral weights $a_{n}=1$, the squared pattern
$g^{2}(x)$ gets minimized everywhere except at $x=1$, so that a maximum
$\mathrm{DI}=10\,\lg Q$ is reached. Fig. 2(a) shows resulting patterns for
$\mathrm{D}=2,3$ and $\mathrm{N}=1,2,5$.

Its simple, closed-form expression found in [Raf19] is based on the
Christoffel-Darboux identity eq. (96) subsection E.6 and for any number of
dimensions with $x=\cos\vartheta$, $y=1$

$\displaystyle g(\cos\vartheta)$
$\displaystyle=\frac{\mathrm{D}-1}{2\mathrm{N}+\mathrm{D}-1}\frac{\mathcal{P}_{\mathrm{N}+1}(\cos\vartheta)-\mathcal{P}_{\mathrm{N}}(\cos\vartheta)}{\cos\vartheta-1}.$
(21)

(a) basic (b) max-$r_{\mathrm{E}}$

Figure 2: Basic and max-$r_{\mathrm{E}}$, $\mathrm{D}=2,3$ (solid,dashed),
$\mathrm{N}=1,2,5$ (green,red,blue).

#### Max$-\bm{r}_{\mathrm{E}}$ weights.

To shape the finite-degree axisymmetric panning function, we choose
max-$\bm{r}_{\mathrm{E}}$ weights [DRP99] $a_{n}$ that maximize the magnitude
$r_{\mathrm{E}}$ of the $\bm{r}_{\mathrm{E}}$ vector. In great mathematical
detail, subsection F.3 shows the solution to be an axisymmetric Dirac delta,
so $a_{n}=\mathcal{P}_{n}$, evaluated at the greatest root
$r_{\mathrm{E,max}}$ of $\mathcal{P}_{\mathrm{N}+1}$ in the next higher
degree,

$\displaystyle a_{n}$ $\displaystyle=\mathcal{P}_{n}(r_{\mathrm{E,max}}),$
$\displaystyle r_{\mathrm{E,max}}$
$\displaystyle=\mathrm{arg\,max}_{r_{\mathrm{E}}}\\{\mathcal{P}_{\mathrm{N}+1}(r_{\mathrm{E}})=0\\}.$
(22)

In general, this requires a root finder or approximation of the greatest root.
For $D=\\{2,3\\}$, the value is
$r_{\mathrm{E,max}}=\\{\cos\frac{\pi}{2(\mathrm{N}+1)},\,\cos\frac{\pi}{180^{\circ}}\frac{137.9^{\circ}}{\mathrm{N}+1.51}\\}$,
where the value for $D=3$ is an approximation [ZF12]. Fig. 2(b) shows
max-$r_{\mathrm{E}}$ patterns for $\mathrm{D}=2,3$ and $\mathrm{N}=1,2,5$.

The Christoffel-Darboux identity eq. (94) subsection E.6 simplifies by
$\mathcal{P}_{\mathrm{N}+1}(r_{\mathrm{E}})=0$ and yields the simple, closed-
form expression for the max-$r_{\mathrm{E}}$ directivity normalized to
$g(1)=1$, with removable singularity at $x=r_{\mathrm{E,max}}$

$\displaystyle g(\cos\vartheta)$
$\displaystyle=\frac{(1-r_{\mathrm{E,max}})\mathcal{P}_{\mathrm{N}+1}(\cos\vartheta)}{\cos\vartheta-
r_{\mathrm{E,max}}}.$ (23)

(a) supercardioid (b) inphase

Figure 3: Supercardioid and in-phase, $\mathrm{D}=2,3$ (solid,dashed),
$\mathrm{N}=1,2,5$ (blue,red).

#### Supercardioid/Maximum front-to-back ratio weights.

For the supercardioid functions, the front-to-back energy ratio (FBR) is
maximized. Fig. 3(a) shows the resulting functions. The resulting problem is
an eigenvector problem

$\displaystyle\bm{a}^{\top}\bm{G}\bm{a}$ $\displaystyle\rightarrow\max$ (24)
$\displaystyle\text{s.t. }\|\bm{a}\|$ $\displaystyle=1.$

The Gram matrix $\bm{G}$ contains
$g_{nn^{\prime}}=\int_{0}^{1}\frac{\mathcal{P}_{n}(x)\mathcal{P}_{n^{\prime}}(x)}{N_{n}^{2}N_{n^{\prime}}^{2}}w(x)\mathrm{d}x$
and its entries are derived in subsection F.1. The vector $\bm{a}$ is a
normalized vector with the largest-possible eigenvalue. Apparently, the
resulting FBR yields an exponential improvement of the optimal front-to-back
energy ratios for the circle and sphere, yielding after regression for
$N=1,\dots,5$,

$\displaystyle
10\lg\mathrm{FBR_{max}}=13.75\,\mathrm{N}-\begin{cases}3.6\thinspace{\mathrm{dB}},&\mathrm{D=2},\\\
3.0\thinspace{\mathrm{dB}},&\mathrm{D=3}.\end{cases}$ (25)

As obtaining $a_{n}$ by eigendecomposition can be impractical, the similarity
to the weights $a_{n,\mathrm{inphase}}$ below is helpful: $a_{n}$ can be
approximated by $a_{n,\mathrm{inphase}}^{\beta}$, modified by an approximate
exponent $0.57\leq\beta\leq 0.7$ determined numerically for $\mathrm{D}=2,3$
and $\mathrm{N}=1,\dots,10$,

$\displaystyle a_{n}$ $\displaystyle\approx a_{n,\mathrm{inphase}}^{\beta},$
$\displaystyle\beta$
$\displaystyle=\frac{0.73\,\mathrm{N}+0.67\mathrm{D}-1.11}{\mathrm{N}+1.11\mathrm{D}-1.5}.$
(26)

Its absolute approximation error is limited
$|a_{n,\mathrm{inphase}}^{\beta}-a_{n}|<-44\thinspace\mathrm{dB}$.

#### In-phase/Higher-order cardioid.

According to [DRP99], it can be helpful to define weights that entirely
suppress the side lobes by causing an $\mathrm{N}$-fold zero at $x=-1$, which
is often referred to as higher-order cardioid, and which is a max-flat design
at $x=-1$,

$\displaystyle g(x)$ $\displaystyle\propto(1+x)^{\mathrm{N}},$ (27)

and subsection F.6 complements the original derivation of [Dan01] to formally
define the in-phase weighting for arbitrary dimensions

$\displaystyle
a_{n}=\frac{\mathrm{N}!}{(\mathrm{N}-n)!}\frac{(\mathrm{N}+\mathrm{D}-2)!}{(\mathrm{N}+n+\mathrm{D}-2)!}.$
(28)

Fig. 3(b) shows in-phase patterns for $\mathrm{D}=2,3$ and $\mathrm{N}=1,2,5$.

(a) $\mathrm{D}=2$ (b) $\mathrm{D}=3$

Figure 4: Examples for maximally flat designs with $\mathrm{L}+1$ flatness
degrees on axis, and a zero with $\mathrm{M}+1$ flatness degrees for
$\mathrm{N}=\mathrm{L}+\mathrm{M}+1$ and for $\mathrm{D}=2$ and
$\mathrm{D}=3$.

#### Other max-flat/Butterworth designs

allow to employ some of the polynomial degrees to flatten the directivity on
axis, as inspired by the FIR filter design firflatLp(N,K,M) from the book
[Sch10, p.92]. The recursive algorithm is defined in subsection F.6 in the
equations (152) and (153). The total oder $\mathrm{N}=\mathrm{L}+\mathrm{M}+1$
is composed of the degree of flatness $\mathrm{M}$ for the zero at $x=-1$ and
$\mathrm{L}$ for unity at $x=1$. Examples for $\mathrm{N}=5$, varying
$\mathrm{L}$ and $\mathrm{D}=2$ or $\mathrm{D}=3$ are shown in Fig. 4; angular
flatness at $\phi=0,\pi$ is twice as high.
_Inconsistency warning: Butterworth as used by Devaraju[Dev15] refers to the
weights $a_{n}=1/{\sqrt{1+(n/n_{c})^{2k}}}$ discretizing the continuous
spectral roll-off of a Butterworth low pass curve. Despite the inconsistency
with Devaraju, I chose to keep the name Butterworth synonymous with the
existence of a flatness constraint of a continuous function, which is
applicable to the directional pattern across its continuous angle._

Figure 5: Cap functions for $\mathrm{D}=2,3$ (solid,dashed), and
$\mathrm{N}=7,30$ (blue,red).

#### Cap functions.

The weights $a_{n}$ for a spherical cap function from $x_{0}$ to $1$ are
derived in subsection F.5 to be

$\displaystyle
a_{n}=\frac{w(x_{0})}{2n+2\alpha}\,\bigg{[}\mathcal{P}_{n-1}(x_{0})-\mathcal{P}_{n+1}(x_{0})\bigg{]}.$
(29)

with the exception at $n=0$, where $a_{0}=\upvarphi_{0}$ for $\mathrm{D}=2$
and $a_{0}=1-x_{0}$ for $\mathrm{D}=3$. Fig. 5 shows the cap functions for an
$80^{\circ}$ cap represented with the degrees $\mathrm{N}=\\{7,\,30\\}$ for
$\mathrm{D}=2,3$.

These weights are interesting to model loudspeakers on rigid spheres [ZZFK17].
In [ZRF22] a trapezoidal function was used to approximate the VBAP kernel, and
it consisted of the convolution of two differently sized cap functions,
obtained by multiplying $a_{n}^{(1.375\alpha)}\,a_{n}^{(0.75\alpha)}$ with an
average spacing angle $\alpha$, defining either
$x_{0}=\cos\frac{1.375\alpha}{2}$ or $x_{0}=\cos\frac{0.75\alpha}{2}$ for the
two caps.

#### FIGURES WITH METRICS

The Figs. 6,7,8,9 evaluate the patterns presented above in terms of the $Q$,
$r_{\mathrm{V}}$, $r_{\mathrm{E}}$, $\mathrm{FBR}$ metrics, with
implementation example and verification in Appendices G.10 and G.11.

Clearly, in Fig. 6, the basic/max-DI pattern outperforms the others, followed
by max-$r_{\mathrm{E}}$, supercardioid, and inphase as the pattern with the
worst directivity factor $Q$.
In Fig. 7, the conversion of $\|\bm{r}_{\mathrm{V}}\|$ to a spread angle by
$\arccos\|\bm{r}_{\mathrm{V}}\|$ clearly indicates that the basic/max-DI is
the pattern with zero spread or max-$r_{\mathrm{V}}$ that is unity, followed
by max-$r_{\mathrm{E}}$, supercardioid, and at last inphase exhibiting the
largest spread; similarly as in Fig. 6 about the directivity factor Q.
This changes in Fig. 8, where the conversion of $\|\bm{r}_{\mathrm{E}}\|$ to a
spread angle by $\arccos\|\bm{r}_{\mathrm{E}}\|$ clearly leads the smallest
spread for max-rE, followed by supercardioid, and inphase en par with basic
are worst and exhibit the largest spread.
The front-to-back ratio in Fig. 9 is clearly best for supercardioid, followed
by inphase at high degrees and max-rE at low degrees, and basic exhibits the
worst FBR.

This verifies that by definition, basic/max-DI optimally isolates a direction
within a diffuse sound field, max-$r_{\mathrm{E}}$ optimally defines the
perceived direction of panned signals in surround sound, inphase is clearly
the best in suppressing everything around the opposite-axis direction, and
supercardioid is best in isolating the on-axis halfspace from the opposite-
axis one.

Dolph-Chebyshev designs were presented by Rafaely [Raf19] and are omitted here
for brevity, also some designs mentioned on https://spaudiopy.readthedocs.io
the notable repository of Chris Hold.

Figure 6: $Q$ values for different pattern designs. Figure 7: $\arccos
r_{\mathrm{V}}$ values for different pattern designs. Figure 8: $\arccos
r_{\mathrm{E}}$ values for different pattern designs. Figure 9:
$10\lg\mathrm{FBR}$ values for different pattern designs.

###  1.3 Optimal discretization: spherical $t$ design

Assume that we sample a polynomial in
$\left(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}\right)^{n}$ limited in
degree by $n\leq t$, $t\in\mathbb{N}$. The discretization uses the directions
$\\{\bm{\uptheta}_{l}\\}$ at $\mathrm{L}$ points and that this is called $t$
design whenever discrete summation and integral are equivalent, i.e., discrete
summation accurately determines the direct component (DC) of the polynomial
despite the variable orientation of $\bm{\theta}_{s}$, according to the
definition in [DGS77],

$\displaystyle\int_{\mathbb{S}_{\mathrm{D}-1}}\mathcal{P}_{n}\left(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}\right)\,\mathrm{d}\bm{\theta}$
$\displaystyle={\textstyle\frac{S_{\mathrm{D}-1}}{\mathrm{L}}}\sum_{l=1}^{\mathrm{L}}\mathcal{P}_{n}\left(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\uptheta}_{l}\right).$
(30)

This simplistic criterion of just estimating the DC correctly is relevant for
all models:

  * $\bullet$

In the calculation of $P$, an integral over the
$\mathrm{N}^{\mathrm{th}}$-degree polynomial pattern $g(\bm{\theta})$ is
evaluated, which is why at least a $t\geq\mathrm{N}$ design should be used to
get correct results for $P$ by summing up the discretized pattern
$g(\bm{\uptheta}_{l})$ times $\mathrm{d}\Omega=\frac{S_{D-1}}{\mathrm{L}}$.

  * $\bullet$

In the calculation of $\bm{r}_{\mathrm{V}}$, the
$\mathrm{N}^{\mathrm{th}}$-degree pattern $g(\bm{\theta})$ is multiplied by
$\bm{\theta}^{\mathrm{T}}_{\mathrm{s}}\bm{\theta}$, causing an integral over
an $(\mathrm{N}+1)^{\mathrm{th}}$-degree polynomial, therefore discretization
requires at least a $t\geq\mathrm{N}+1$ design for the numerator of
$\bm{r}_{\mathrm{V}}$ to sum up
$\bm{\theta}^{\mathrm{T}}_{\mathrm{s}}\bm{\theta}\,g(\bm{\uptheta}_{l})$
correctly.

  * $\bullet$

For calculating $E$, the squared continuous pattern $g(\bm{\theta})^{2}$ is
integrated, a $2\mathrm{N}^{\mathrm{th}}$-degree polynomial, thus
discretization should at least consider a $t\geq 2\mathrm{N}$ design so that
summing up the discrete pattern $g(\bm{\uptheta}_{l})^{2}$ times
$\mathrm{d}\Omega=\frac{S_{D-1}}{\mathrm{L}}$ correctly estimates ${E}$.

  * $\bullet$

In the calculation of $\bm{r}_{\mathrm{E}}$, the squared continuous function
$g(\bm{\theta})^{2}$ multiplied by
$\bm{\theta}^{\mathrm{T}}_{\mathrm{s}}\bm{\theta}$ is integrated, a
$(2\mathrm{N}+1)^{\mathrm{th}}$-degree polynomial, so discretization should
consider a $t\geq 2\mathrm{N}+1$ design for correct results for the numerator
of $\bm{r}_{\mathrm{E}}$ that sums up
$\bm{\theta}^{\mathrm{T}}_{\mathrm{s}}\bm{\theta}\,g(\bm{\uptheta}_{l})^{2}$.

In conclusion, using at least a $t\geq 2\mathrm{N}+1$ design ensures that $E$,
$P$, $\bm{r}_{\mathrm{V}}$, and $\bm{r}_{\mathrm{E}}$ obtained by from
discrete patterns $g(\bm{\uptheta}_{l})$ complies with results for the
continuous function $g(\bm{\theta})$.

_Note that the FBR measure is an exception because its lower integration limit
may exclude a varying number of discretization nodes, depending on orientation
and is therefore more accurate when solved with the analytic values $g_{nm}$
of the Gram matrix, see subsection F.1._

#### Solutions in 2 and 3 dimensions.

For $\mathrm{D}=2$ we get the Chebyshev differential equation and the
Chebyshev polynomials $T_{n}$ as solutions, and for $\mathrm{D}=3$ the
Legendre differential equation with Legendre polynomials $P_{n}$ as solutions

$\displaystyle\mathrm{D}$ $\displaystyle=2:$
$\displaystyle(1-x^{2})T_{n}^{\prime\prime}-x\,T_{n}^{\prime}+n^{2}\,T_{n}$
$\displaystyle=0,$ (31) $\displaystyle\mathrm{D}$ $\displaystyle=3:$
$\displaystyle(1-x^{2})P_{n}^{\prime\prime}-2\,x\,P_{n}^{\prime}+n(n+1)\,P_{n}$
$\displaystyle=0.$ (32)

These polynomials and their three-term recurrences are listed in Tables 1,2,
yielding polynomials normalized at $x=1$, i.e., $T_{n}(1)=1$ and $P_{n}(1)=1$
for all $n$, and they fulfill the orthogonality relations

$\displaystyle\mathrm{D}$ $\displaystyle=2:$
$\displaystyle\int_{-1}^{1}T_{n}(x)T_{m}(x)\,\frac{\mathrm{d}x}{\sqrt{1-x^{2}}}$
$\displaystyle=\frac{2}{2-\delta_{m,0}}\,\delta_{nm},$
$\displaystyle\text{with }N_{n}^{2}$ $\displaystyle=\frac{2-\delta_{m,0}}{2},$
(33) $\displaystyle\mathrm{D}$ $\displaystyle=3:$
$\displaystyle\int_{-1}^{1}P_{n}(x)P_{m}(x)\,{\mathrm{d}x}$
$\displaystyle=\frac{2}{2n+1}\,\delta_{nm},$ $\displaystyle\text{with
}N_{n}^{2}$ $\displaystyle=\frac{2n+1}{2}.$ (34)

and recurrences

$\displaystyle\mathrm{D}$ $\displaystyle=2:$ $\displaystyle T_{m+1}(x)$
$\displaystyle=(2-\delta_{m,0})x\,T_{m}(x)-T_{m-1}(x),$ (35)
$\displaystyle\mathrm{D}$ $\displaystyle=3:$ $\displaystyle P_{n+1}(x)$
$\displaystyle=\frac{2n+1}{n+1}xP_{n}(x)-\frac{n}{n+1}P_{n-1}(x).$ (36)

The Dirac delta $\delta(1-\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})$
pointing to $\bm{\theta}_{\mathrm{s}}$ is represented by the coefficients

$\displaystyle\mathrm{D}=2:\quad\gamma_{n}$
$\displaystyle=\frac{2-\delta_{n,0}}{2\pi},$
$\displaystyle\mathrm{D}=3:\quad\gamma_{n}$ $\displaystyle=\frac{2n+1}{4\pi}.$
(37)

###  1.4 Continuous circular patterns in Chebyshev polynomials

On the circle, the equivalence of a Dirac delta evaluated at the projection of
the variable direction $\bm{\theta}$ on the panning direction
$\bm{\theta}_{\mathrm{s}}$ is given in terms of an infinite sum of Chebyshev
polynomials, see Tab. 1, cf. [ZF12],

$\displaystyle\sum_{m=0}^{\infty}\frac{2-\delta_{m,0}}{2\pi}\,T_{m}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})=\delta(1-\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}),$

with
$\int_{\mathbb{S}^{1}}\delta(1-\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})\,\mathrm{d}\bm{\theta}=1$,
and
$\delta(1-\bm{\theta}^{\mathrm{T}}_{\mathrm{s}}\bm{\theta})=\begin{cases}\infty,&\text{for
$\bm{\theta}=\bm{\theta}_{s}$},\\\ 0,&\text{otherwise}.\end{cases}$
The limited-degree panning function should involve the sidelobe-suppressing
weights $a_{m}$

$\displaystyle
g(\bm{\theta})=\sum_{m=0}^{\mathrm{N}}\frac{2-\delta_{m,0}}{2\pi}\,a_{m}\,T_{m}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}),$

for which we may insert $\mathrm{max}-\bm{r}_{\mathrm{E}}$ weights [DRP99],
$a_{m}=T_{m}\Bigl{(}\cos\frac{90^{\circ}}{\mathrm{N}+1}\Bigr{)}$, as derived
in subsection F.3.

Table 1: Coefficients of the Chebyshev polynomials that are orthogonal for integration on the interval $[-1,1]$ with the integration element $\frac{\mathrm{d}x}{\sqrt{1-x^{2}}}$. The Chebyshev polynomials obey the recurrence $T_{m+1}(x)+T_{m-1}(x)-(2-\delta_{m,0})x\,T_{m}(x)=0$ for $m>1$.  $T_{m}(x)=\sum_{k=0}^{m}c_{k}\,x^{k}$ | Coefficients of the Chebyshev polynomials |  |
---|---|---|---
$m=$ | $c_{9}$ | $c_{8}$ | $c_{7}$ | $c_{6}$ | $c_{5}$ | $c_{4}$ | $c_{3}$ | $c_{2}$ | $c_{1}$ | $c_{0}$
$0$ |  |  |  |  |  |  |  |  |  | $1$
$1$ |  |  |  |  |  |  |  |  | $1$ |
$2$ |  |  |  |  |  |  |  | $2$ |  | $-1$
$3$ |  |  |  |  |  |  | $4$ |  | $-3$ |
$4$ |  |  |  |  |  | $8$ |  | $-8$ |  | $1$
$5$ |  |  |  |  | $16$ |  | $-20$ |  | $5$ |
$6$ |  |  |  | $32$ |  | $-48$ |  | $18$ |  | $-1$
$7$ |  |  | $64$ |  | $-112$ |  | $56$ |  | $-7$ |
$8$ |  | $128$ |  | $-256$ |  | $160$ |  | $-32$ |  | $1$
$9$ | $256$ |  | $-576$ |  | $432$ |  | $-120$ |  | $9$ |

#### Measures.

The function is of limited angular resolution $\mathrm{N}$ and allows to be
sampled correctly. The $E$, $P$, $Q$, $\bm{r}_{\mathrm{V}}$, and
$\bm{r}_{\mathrm{E}}$ measures of the continuous function can be derived, see
Appendix F, and for verification with the Fourier transform Appendix H,

$\displaystyle P$ $\displaystyle=a_{0},$ $\displaystyle\bm{r}_{\mathrm{V}}$
$\displaystyle=\frac{a_{1}}{a_{0}}\;\bm{\theta}_{\mathrm{s}},$ (38)
$\displaystyle E$
$\displaystyle=\sum_{n=0}^{\mathrm{N}}\frac{|a_{n}|^{2}}{2\pi},$
$\displaystyle Q$
$\displaystyle=\frac{\left|\sum_{n=0}^{\mathrm{N}}(2-\delta_{n})a_{n}\right|^{2}}{\sum_{n=0}^{\mathrm{N}}(2-\delta_{n})|a_{n}|^{2}},$
$\displaystyle\bm{r}_{\mathrm{E}}$
$\displaystyle=\frac{\sum_{n=0}^{\mathrm{N}}2\Re\\{a_{n}a_{n+1}^{*}\\}}{\sum_{n=0}^{\mathrm{N}}(2-\delta_{n})|a_{n}|^{2}}\;\bm{\theta}_{\mathrm{s}}.$
(39)

###  1.5 Addition theorem of the circular harmonics

The 2D Ambisonic format represents source directions at the azimuth angle
$\varphi_{\mathrm{s}}$ along the continuous azimuth $\varphi$ using the
circular harmonics. These are the typical terms of a real-valued Fourier
transform of $2\pi$-periodic functions, here in the azimuth angle $\varphi$,

$\displaystyle\Phi_{m}(\bm{\theta})=\sqrt{\frac{2-\delta_{m}}{2\pi}}\begin{cases}\cos|m|\varphi,&m\geq
0\\\ \sin|m|\varphi,&m<0,\end{cases}$ (40)

The addition theorem
$\cos(\alpha-\beta)=\cos\alpha\cos\beta+\sin\alpha\sin\beta$ of the circular
harmonics

$\displaystyle\Phi_{m}(\bm{\theta})\Phi_{m}(\bm{\theta}_{\mathrm{s}})+(1-\delta_{m})\Phi_{-m}(\varphi)\Phi_{-m}(\varphi_{\mathrm{s}})={\frac{2-\delta_{m}}{2\pi}}\,T_{m}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})$
(41)

relates them to the polynomials $T_{m}$, so that we find the equivalence for
$g(\bm{\theta})$

$\displaystyle
g(\bm{\theta})=\sum_{m=-\mathrm{N}}^{\mathrm{N}}a_{|m|}\Phi_{m}(\bm{\theta}_{\mathrm{s}})\Phi_{m}(\bm{\theta})$
$\displaystyle=\sum_{n=0}^{\mathrm{N}}{\frac{2-\delta_{m}}{2\pi}}\,a_{m}\,T_{m}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}).$
(42)

This relation shows for single source directions $\bm{\theta}_{\mathrm{s}}$
that axisymmetric functions represented by the 2D Ambisonic format are
equivalent to the polynomial representation, here.

###  1.6 Optimal circular discretization: circular $t$ design

Discretization uses the directions $\\{\bm{\uptheta}_{l}\\}$ at $\mathrm{L}$
points that this is called $t$ design equivalently with subsection 1.3
whenever discrete summation of the directional function and integral are
equivalent. For $\mathrm{D}=2$, this involves the Chebyshev polynomials

$\displaystyle\int_{-\pi}^{\pi}T_{m}\left(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}\right)\,\mathrm{d}\bm{\theta}$
$\displaystyle={\textstyle\frac{2\pi}{\mathrm{L}}}\sum_{l=1}^{\mathrm{L}}T_{m}\left(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\uptheta}_{l}\right).$
(43)

On the circle, every equiangular sampling
$\upvarphi_{l}=\frac{2\pi}{\mathrm{L}}(l-1)$ is a $t$-design if the number of
samples fulfills the sampling criterion $t\leq\mathrm{L}+1$. There should be
equally many sampling points than zero crossings of $\cos(m\varphi)$ or
$\sin(m\varphi)$, or two samples per oscillation (equivalent to
Whittaker/Kotelnikov/Shannon sampling theorem).

Our knowledge about $t$ designs, equiangular sampling on the circle, and the
polynomial degrees in the various metrics allow to draw the following
conclusion: We need

$\displaystyle\mathrm{L}\geq 2\mathrm{N}+2$ (44)

equiangularly spaced nodes on the horizon to obtain an optimally placed
Ambisonic playback system of the degree $\mathrm{N}$. The
$\bm{r}_{\mathrm{E}}$ vector is the strictest of the criteria.

###  1.7 Continuous spherical patterns in Legendre polynomials

On the sphere, the equivalence of a Dirac delta evaluated at the projection of
the variable direction $\bm{\theta}$ on the panning direction
$\bm{\theta}_{\mathrm{s}}$ is given in terms of an infinite sum of Legendre
polynomials, see Tab. 2, cf. [ZF12],

$\displaystyle\sum_{n=0}^{\infty}\frac{2n+1}{4\pi}\,P_{n}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})=\delta(1-\bm{\theta}^{\mathrm{T}}_{\mathrm{s}}\bm{\theta}),$

for which
$\int_{\mathbb{S}^{2}}\delta(1-\bm{\theta}^{\mathrm{T}}_{\mathrm{s}}\bm{\theta}),\,\mathrm{d}\bm{\theta}=1$
and
$\delta(1-\bm{\theta}^{\mathrm{T}}_{\mathrm{s}}\bm{\theta})=\begin{cases}\infty,&\text{for
$\bm{\theta}=\bm{\theta}_{s}$},\\\ 0,&\text{otherwise}.\end{cases}$
The limited-degree summation of the panning function involves the coefficients
$a_{n}$

$\displaystyle g(\bm{\theta})$
$\displaystyle=\sum_{n=0}^{\mathrm{N}}\frac{2n+1}{4\pi}\,a_{n}\,P_{n}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}),$
(45)

and it may use $\mathrm{max}-\bm{r}_{\mathrm{E}}$ weights according to [ZF12,
DRP99], as derived in subsection F.3,

$\displaystyle a_{n}$
$\displaystyle=P_{n}\Bigl{(}\cos\frac{137.9^{\circ}}{\mathrm{N}+1.51}\Bigr{)}.$
(46)

Table 2: Coefficients of the Legendre polynomials that are orthogonal for the integration on the interval $[-1,1]$ with the integration element $\mathrm{d}x$. The recurrence to obtain these coefficients is $(n+1)P_{n+1}(x)-(2n+1)xP_{n}(x)+nP_{n-1}(x)=0$.  $P_{n}(x)=\sum_{k=0}^{n}c_{k}\,x^{i}$ | Coefficients of the Legendre polynomials |  |
---|---|---|---
$n=$ | $c_{9}$ | $c_{8}$ | $c_{7}$ | $c_{6}$ | $c_{5}$ | $c_{4}$ | $c_{3}$ | $c_{2}$ | $c_{1}$ | $c_{0}$
$0$ |  |  |  |  |  |  |  |  |  | $1$
$1$ |  |  |  |  |  |  |  |  | $1$ |
$2$ |  |  |  |  |  |  |  | $\frac{3}{2}$ |  | $-\frac{1}{2}$
$3$ |  |  |  |  |  |  | $\frac{5}{2}$ |  | $-\frac{3}{2}$ |
$4$ |  |  |  |  |  | $\frac{35}{8}$ |  | $-\frac{30}{8}$ |  | $\frac{3}{8}$
$5$ |  |  |  |  | $\frac{63}{8}$ |  | $-\frac{70}{8}$ |  | $\frac{15}{8}$ |
$6$ |  |  |  | $\frac{231}{16}$ |  | $-\frac{315}{16}$ |  | $\frac{105}{16}$ |  | $-\frac{5}{16}$
$7$ |  |  | $\frac{429}{16}$ |  | $-\frac{693}{16}$ |  | $\frac{315}{16}$ |  | $-\frac{35}{16}$ |
$8$ |  | $\frac{6435}{128}$ |  | $-\frac{12012}{128}$ |  | $\frac{6930}{128}$ |  | $-\frac{1260}{128}$ |  | $\frac{35}{128}$
$9$ | $\frac{12155}{128}$ |  | $-\frac{25740}{128}$ |  | $\frac{18018}{128}$ |  | $-\frac{4620}{128}$ |  | $\frac{315}{128}$ |

#### Measures.

The function is of limited angular resolution $\mathrm{N}$ and allows to be
sampled correctly. The $E$, $Q$, $P$, $\bm{r}_{\mathrm{V}}$, and
$\bm{r}_{\mathrm{E}}$ measures of the continuous function can be derived, see
Appendix F,

$\displaystyle P$ $\displaystyle=a_{0},$ $\displaystyle\bm{r}_{\mathrm{V}}$
$\displaystyle=\frac{a_{1}}{a_{0}}\;\bm{\theta}_{\mathrm{s}},$ (47)
$\displaystyle E$
$\displaystyle=\sum_{n=0}^{\mathrm{N}}\frac{2n+1}{4\pi}|a_{n}|^{2},$
$\displaystyle Q$
$\displaystyle=\frac{\left|\sum_{n=0}^{\mathrm{N}}(2n+1)a_{n}\right|^{2}}{\sum_{n=0}^{\mathrm{N}}(2n+1)|a_{n}|^{2}},$
$\displaystyle\bm{r}_{\mathrm{E}}$
$\displaystyle=\frac{\sum_{n=0}^{\mathrm{N}}(n+1)2\Re\\{a_{n}a_{n+1}^{*}\\}}{\sum_{n=0}^{\mathrm{N}}(2n+1)|a_{n}|^{2}}\;\bm{\theta}_{\mathrm{s}}.$
(48)

###  1.8 Addition theorem of the spherical harmonics

The 3D Ambisonic format is based on representing source directions using the
spherical harmonics. These involve the azimuth angle $\varphi$ and the zenith
angle $\vartheta$ of the direction vector $\bm{\theta}$, and the associated
Legendre functions $P_{n}^{m}(x)$, which are equivalent to the Legendre
polynomials for $m=0$,

$\displaystyle
Y_{n}^{m}(\bm{\theta})=(-1)^{|m|}\sqrt{\frac{(2-\delta_{m})}{4\pi}\frac{(2n+1)\,(n-|m|)!}{(n+|m|)!}}\,P_{n}^{|m|}(\cos\vartheta)\,\begin{cases}\cos|m|\varphi,&m\geq
0\\\ \sin|m|\varphi,&m<0.\end{cases}$ (49)

Their addition theorem relates them to the Legendre polynomials $P_{n}(x)$

$\displaystyle\sum_{m=-n}^{n}Y_{n}^{m}(\bm{\theta})Y_{n}^{m}(\bm{\theta}_{\mathrm{s}})={\frac{2n+1}{4\pi}}\,P_{n}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}),$
(50)

so that we find the equivalence for $g(\bm{\theta})$

$\displaystyle
g(\bm{\theta})=\sum_{n=0}^{\mathrm{N}}\sum_{m=-n}^{n}a_{n}Y_{n}^{m}(\bm{\theta})Y_{n}^{m}(\bm{\theta}_{\mathrm{s}})$
$\displaystyle=\sum_{n=0}^{\mathrm{N}}{\frac{2n+1}{4\pi}}\,a_{n}\,P_{n}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta}).$
(51)

Single source directions $\bm{\theta}_{\mathrm{s}}$ represented within 3D
Ambisonics are therefore equivalent to the Legendre-polynomial representation.

###  1.9 Optimal spherical discretization: spherical $t$-design

Spherical $t$-designs are $\mathrm{L}$ sampling nodes
$\\{\bm{\uptheta}_{l}\\}$ on the sphere that allow to replace the integration
over any spherical polynomial of limited degree $n\leq t$ to be evaluated
correctly using a summation over the sampled polynomials. Note that the
direction $\bm{\theta}_{\mathrm{s}}$ stays arbitrary, and a $t$-design
establishes the exact equivalence

$\displaystyle\int_{\mathbb{S}^{2}}P_{n}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\theta})\,\mathrm{d}\bm{\theta}={\textstyle\frac{4\pi}{\mathrm{L}}}\sum_{l=1}^{\mathrm{L}}P_{n}(\bm{\theta}_{\mathrm{s}}^{\mathrm{T}}\bm{\uptheta}_{l}).$
(52)

The measures $P$, $E$, $\bm{r}_{\mathrm{V}}$, and $\bm{r}_{\mathrm{E}}$ of the
continuous functions, as shown in the Appendix F, require integration over
$\mathrm{N}^{\mathrm{th}}$, $(\mathrm{N}+1)^{\mathrm{th}}$,
$2\mathrm{N}^{\mathrm{th}}$, and $(2\mathrm{N}+1)^{\mathrm{th}}$ degree
polynomials, respectively. Exact equivalence between continuous and
discretized functions therefore require the discrete sampling nodes
$\\{\bm{\uptheta}_{l}\\}$ to be at least a

$\displaystyle t\geq 2\mathrm{N}+1\text{ design.}$ (53)

Apart from the exception of 5 regular polyhedra (tetrahedron with
$\mathrm{L}=4$ and $t=2$, octahedron with $\mathrm{L}=6$ and $t=3$, cube with
$\mathrm{L}=8$ and $t=3$, icosahedron with $\mathrm{L}=12$ and $t=5$,
dodecahedron with $\mathrm{L}=20$ and $t=5$), spherical designs are otherwise
generally found by optimization, as, e.g., by Hardin and Sloane [HS96] and
Gräf and Potts [GP11], as available on the following websites
_http://neilsloane.com/sphdesigns/dim3/_ and
_http://homepage.univie.ac.at/manuel.graef/quadrature.php_
(Chebyshev-type Quadratures on $\mathbb{S}^{2}$).

##  Appendix A Laplacian in axisymmetric coordinates

A function in the last, normalized Cartesian coordinate
$x:=\frac{x_{\mathrm{D}}}{r}=\cos\vartheta$ is axisymmetric with regard to
$\mathbf{i}_{\mathrm{D}}=[0,\dots,0,1]^{\top}$, and the generalized chain rule
yields for the Laplacian

$\displaystyle\bigtriangleup=\sum_{j}\frac{\partial^{2}}{\partial
x_{j}^{2}}=\sum_{j}\frac{\partial}{\partial x_{j}}\left(\frac{\partial
x}{\partial
x_{j}}\frac{\mathrm{d}}{\mathrm{d}x}\right)=\sum_{j}\frac{\partial^{2}x}{\partial
x_{j}^{2}}\frac{\mathrm{d}}{\mathrm{d}x}+\sum_{j}\left(\frac{\partial
x}{\partial x_{j}}\right)^{2}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}.$ (54)

The derivative of $\frac{1}{r}$ keeps re-appearing. It is with regard to any
of the coordinates

$\displaystyle\frac{\partial}{\partial
x_{j}}\frac{1}{r}=-\frac{1}{r^{2}}\frac{\partial r}{\partial
x_{j}}=-\frac{1}{r^{2}}\frac{\partial[\sum_{j}x_{j}^{2}]^{\frac{1}{2}}}{\partial
x_{j}}=-\frac{1}{r^{2}}\frac{[\sum_{j}x_{j}^{2}]^{-\frac{1}{2}}}{2}\frac{\partial\sum_{j}x_{j}^{2}}{\partial
x_{j}}=-\frac{1}{r^{3}}\frac{1}{2}2x_{j}=-\frac{x_{j}}{r^{3}}.$ (55)

The derivatives of $x:=\frac{x_{\mathrm{D}}}{r}$ are for $j\neq 1$:

$\displaystyle\frac{\partial x}{\partial x_{j}}$
$\displaystyle=\frac{\partial}{\partial
x_{j}}\frac{x_{\mathrm{D}}}{r}=-\frac{x_{\mathrm{D}}\,x_{j}}{r^{3}},$ (56)
$\displaystyle\frac{\partial^{2}x}{\partial^{2}x_{j}}$
$\displaystyle=-\frac{\partial}{\partial
x_{j}}\frac{x_{\mathrm{D}}\,x_{j}}{r^{3}}=-\frac{x_{\mathrm{D}}}{r^{3}}+3\frac{x_{\mathrm{D}}x_{j}^{2}}{r^{5}}=\frac{3x_{\mathrm{D}}x_{j}^{2}-r^{2}x_{\mathrm{D}}}{r^{5}}.$
(57)

The first derivative of $x$ to $x_{\mathrm{D}}$ involves one more term
depending on $x_{\mathrm{D}}$ that we could skip in the product rule before:

$\displaystyle\frac{\partial x}{\partial x_{\mathrm{D}}}$
$\displaystyle=\frac{\partial}{\partial
x_{\mathrm{D}}}\frac{x_{\mathrm{D}}}{r}=-\frac{x_{\mathrm{D}}^{2}}{r^{3}}+\frac{1}{r}=\frac{r^{2}-x_{\mathrm{D}}^{2}}{r^{3}}=\frac{\sum_{j\backslash\mathrm{D}}x_{j}^{2}}{r^{3}}=\frac{r_{\mathrm{D}-1}^{2}}{r^{3}},$
(58) $\displaystyle\frac{\partial^{2}x}{\partial x_{\mathrm{D}}^{2}}$
$\displaystyle=\frac{\partial}{\partial
x_{\mathrm{D}}}\frac{r_{\mathrm{D}-1}^{2}}{r^{3}}=\frac{-3\,x_{\mathrm{D}}\,r_{\mathrm{D}-1}}{r^{5}}.$
(59)

We obtain for the factors involved in front of the second and first-order
derivatives of the axisymmetric Laplacian:

$\displaystyle\sum_{j}\left(\frac{\partial x}{\partial x_{j}}\right)^{2}$
$\displaystyle=\frac{r_{\mathrm{D}-1}^{4}}{r^{6}}+\frac{x_{\mathrm{D}}^{2}\sum_{j\backslash\mathrm{D}}x_{j}^{2}}{r^{6}}=\frac{(r^{2}-x_{\mathrm{D}}^{2})r_{\mathrm{D}-1}^{2}}{r^{6}}+\frac{x_{\mathrm{D}}^{2}r_{\mathrm{D}-1}^{2}}{r^{6}}=\frac{r^{2}r_{\mathrm{D}-1}^{2}}{r^{6}}=\frac{r_{\mathrm{D}-1}^{2}}{r^{4}},$
(60) $\displaystyle\sum_{j}\frac{\partial^{2}\phi}{\partial x_{j}^{2}}$
$\displaystyle=\frac{-3x_{\mathrm{D}}\,r_{\mathrm{D}-1}}{r^{5}}+\frac{3x_{\mathrm{D}}\sum_{j\backslash\mathrm{D}}x_{j}^{2}-\sum_{j\backslash\mathrm{D}}r^{2}x_{\mathrm{D}}}{r^{5}}=-(\mathrm{D}-1)\frac{x_{\mathrm{D}}}{r^{3}}.$
(61)

Considering $\frac{x_{\mathrm{D}}}{r}=x$ and
$\frac{r_{\mathrm{D}-1}}{r}=\sqrt{1-x^{2}}$, i.e.,
$\frac{x_{\mathrm{D}}}{r^{3}}=\frac{x}{r^{2}}$ and
$\frac{r_{\mathrm{D}-1}^{2}}{r^{4}}=\frac{1-x^{2}}{r^{2}}$, the axisymmetric
Laplacian $\bigtriangleup_{x}$ becomes

$\displaystyle
r^{2}\bigtriangleup_{x}=(1-x^{2})\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}-(\mathrm{D}-1)\,x\,\frac{\mathrm{d}}{\mathrm{d}x}.$
(62)

We can re-write the axisymmetric Laplacian (which can be proven using the
chain rule), involving the definition $\alpha=\frac{\mathrm{D}}{2}-1$, and
$w(x)=(1-x^{2})^{\alpha-\frac{1}{2}}$, see Appendix B

$\displaystyle r^{2}\bigtriangleup_{x}\,f$
$\displaystyle=(1-x^{2})f^{\prime\prime}-(2\alpha+1)\,x\,f^{\prime}$ (63)
$\displaystyle=\frac{1}{(1-x^{2})^{\alpha-\frac{1}{2}}}\left[(1-x^{2})^{\alpha+\frac{1}{2}}f^{\prime}\right]^{\prime}=\frac{1}{w(x)}\left[(1-x^{2})\,w(x)\;f^{\prime}\right]^{\prime}.$
(64)

##  Appendix B Jacobi determinant of axisymmetric coordinates

For the axisymmetric coordinate $x$, we re-write the integration constants
$\prod_{j=1}^{\mathrm{D}}\mathrm{d}x_{j}=|\bm{J}|\,\mathrm{d}x\,\prod_{j=1}^{\mathrm{D}-1}\mathrm{d}\tilde{x}_{j}$
with the determinant of the _Jacobian_ matrix that describes how the surface
element is stretched. The Jacobi matrix is filled with derivatives of the old
coordinates with regard to the new ones. One of the new coordinates is $x$
that defines $x_{\mathrm{D}}=r\,x$. The other $\mathrm{D}-1$ coordinates are
just copies of $x_{j}$ but normalized to preserve the information about the
radius as the new radius $r=\sqrt{\sum_{i}\tilde{x}_{i}^{2}}$ of the
$(\mathrm{D}-1)$-dimensional subspace. This defines

$\displaystyle x_{\mathrm{D}}=x\,r,$ $\displaystyle x_{j}$
$\displaystyle=\tilde{x}_{j}\,\sqrt{1-x^{2}},$ (65)

and herewith the derivatives for $\bm{J}$ are

$\displaystyle\bm{J}$ $\displaystyle=\begin{bmatrix}[\frac{\partial
x_{j}}{\partial\tilde{x}_{i}}]_{i,j<\mathrm{D}},&[\frac{\partial
x_{\mathrm{D}}}{\partial\tilde{x}_{i}}]_{i<\mathrm{D}}\\\ [\frac{\partial
x_{j}}{\partial x}]_{j<\mathrm{D}}&\frac{\partial x_{\mathrm{D}}}{\partial
x}\end{bmatrix}=\begin{bmatrix}\mathbf{I}\sqrt{1-x^{2}},&\bigl{[}\frac{x\,\tilde{x}_{i}}{r}\bigr{]}_{i<\mathrm{D}}\\\
\Bigl{[}-\frac{\tilde{x}_{j}\,x}{\sqrt{1-x^{2}}}\Bigr{]}_{j<\mathrm{D}}&r\end{bmatrix},$
(66)

which makes the determinant particularly simple to calculate: (i) the
determinant of a product of matrices is the product of the determinants, (2)
for a block-triangular matrix, the product of the diagonal blocks would be the
determinant. A block-matrix $\bm{E}$ can be made block-triangular by LU
decomposition

$\displaystyle\bm{J}$ $\displaystyle=\begin{bmatrix}\bm{A}&\bm{B}\\\
\bm{C}&\bm{D}\end{bmatrix}=\begin{bmatrix}\bm{A}&\bm{0}\\\
\bm{C}&\bm{I}\end{bmatrix}\begin{bmatrix}\bm{I}&\bm{A}^{-1}\bm{B}\\\
\bm{0}&\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}\end{bmatrix},$ $\displaystyle\text{so
that }\det\bm{J}$
$\displaystyle=\det\bm{A}\,\det(\bm{D}-\bm{C}\bm{A}^{-1}\bm{B}).$ (67)

We get the weight of the integration constant $|\bm{J}|\,\mathrm{d}x$:

$\displaystyle|\bm{J}|$
$\displaystyle=\sqrt{1-x^{2}}^{\mathrm{D}-1}\;\left(r+\frac{1}{\sqrt{1-x^{2}}}\,\sum_{j}\frac{x\,\tilde{x}_{j}}{r}\frac{x\,\tilde{x}_{j}}{\sqrt{1-x^{2}}}\right)$
(68)
$\displaystyle=\frac{\sqrt{1-x^{2}}^{\mathrm{D}-1}}{r}\left(1+\frac{x^{2}r^{2}}{r^{2}(1-x^{2})}\right)=\frac{\sqrt{1-x^{2}}^{\mathrm{D}-3}}{r}=\frac{w(x)}{r}.$
(69)

##  Appendix C Solving the ultraspherical differential equation:
Gegenbauer polynomials

With $\alpha=\frac{\mathrm{D}}{2}-1$, the Laplacian eigenfunction problem
$r^{2}\bigtriangleup y(x)=-\lambda\,y(x)$ is the ultraspherical or Gegenbauer
differential equation for the purely axisymmetric function $y$

$\displaystyle(1-x^{2})y^{\prime\prime}-(2\alpha+1)\,x\,y^{\prime}+\lambda\,y=0.$

It is solved by a power series (Frobenius Ansatz)

$\displaystyle y=\sum_{k=0}^{\infty}c_{k}\;x^{k}.$

Inserted into the differential equation, we can simplify, shift the indices
$k$ to fit the powers $x^{k}$, and find by comparison of coefficients.

$\displaystyle(1-x^{2})\,\sum_{k=2}^{\infty}c_{k}\,k\,(k-1)\,x^{k-2}-(2\alpha+1)\,x\,\sum_{k=1}^{\infty}c_{k}\,k\,x^{k-1}+\lambda\sum_{k=0}^{\infty}c_{k}\,x^{k}$
$\displaystyle=0$
$\displaystyle\sum_{k=2}^{\infty}c_{k}\,k\,(k-1)\,x^{k-2}-\sum_{k=2}^{\infty}c_{k}\,k\,(k-1)\,x^{k}-(2\alpha+1)\,\sum_{k=1}^{\infty}c_{k}\,k\,x^{k}+\lambda\sum_{k=0}^{\infty}c_{k}\,x^{k}$
$\displaystyle=0$
$\displaystyle\sum_{k=0}^{\infty}c_{k+2}\,(k+1)\,(k+1)\,x^{k}-\sum_{k=2}^{\infty}c_{k}\,k\,(k-1)\,x^{k}-(2\alpha+1)\,\sum_{k=1}^{\infty}c_{k}\,k\,x^{k}+\lambda\sum_{k=0}^{\infty}c_{k}\,x^{k}$
$\displaystyle=0$
$\displaystyle\sum_{k=0}^{\infty}\Bigl{\\{}\bigl{[}\lambda-k\,(k-1)\,u[k-2]-(2\alpha+1)\,k\,u[k-1]\bigr{]}\,c_{k}+(k+1)\,(k+2)\,c_{k+2}\Bigr{\\}}x^{k}$
$\displaystyle=0$ $\displaystyle
c_{k+2}=\frac{(2\alpha+1)\,k+k\,(k-1)-\lambda}{(k+1)(k+2)}\;c_{k}.\quad\text{for
$k\geq 2$}$ $\displaystyle\Longrightarrow
c_{k+2}=\frac{k(k+2\alpha)-\lambda}{(k+1)(k+2)}\;c_{k}.\quad\text{for $k\geq
2$}$

Choosing $\lambda=\lambda_{n}=n(n+2\alpha)$, the series terminates at $k=n$.
The recurrence of the coefficients for a certain $n$ can be re-written to

$\displaystyle
c_{k+2}^{(n)}=\frac{k^{2}+2\alpha\,k-n^{2}-2\alpha\,n}{(k+1)(k+2)}\;c_{k}=-\frac{(n-k)(n+k+2\alpha)}{(k+1)(k+2)}\;c_{k}^{(n)}.$

For even $n=2m$, odd-degree coefficients vanish $c_{2l+1}^{(2m)}=0$ and for
odd $n=2m+1$, coefficients of even powers vanish $c_{2l}^{(2m+1)}$. The non-
zero coefficient series are composed factorials and double factorials
(containing every second factor) when applying the recurrence repeatedly, and
they are either based on a starting value $c_{0}^{(n)}$ or $c_{1}^{(n)}$.

We may replace the constants $k+2=n-2l$ and arrive at

$\displaystyle c_{n-2l}^{(n)}$
$\displaystyle=-4(l+1)(n-l-1+\alpha)\frac{c_{n-2l-2}^{(n)}}{(n-2l)(n-2l-1)}$
(70)

Starting from $l=0$, re-iterate the recurrence until arriving at
$c_{n-2\lfloor\frac{n}{2}\rfloor}^{(n)}$ after $2l\leq n$ iterations, or with
round-off operator $l\leq\lfloor\frac{n}{2}\rfloor$. For
$l=\lfloor\frac{n}{2}\rfloor$, the term $n-l$ becomes
$n-\lfloor\frac{n}{2}\rfloor=\lceil\frac{n}{2}\rceil$. We use rising
factorials or Gamma functions
$(\lceil\frac{n}{2}\rceil+\alpha-1)^{\overline{\lfloor\frac{n}{2}\rfloor-l}}=\frac{\Gamma(n-l+\alpha)}{\Gamma(\lceil\frac{n}{2}\rceil+\alpha)}=(\lceil\frac{n}{2}\rceil+\alpha-1)\dots(n-l+\alpha-2)(n-l+\alpha-1)$
to generalize the ratio of factorials to real arguments

$\displaystyle c_{n-2l}^{(n)}$
$\displaystyle=\frac{(-1)^{l}\,2^{n-2l}\lfloor\frac{n}{2}\rfloor!(\lceil\frac{n}{2}\rceil+\alpha-1)^{\overline{\lfloor\frac{n}{2}\rfloor-l}}}{l!(n-2l)!}c_{n-2\lfloor\frac{n}{2}\rfloor}^{(n)}$
$\displaystyle=\frac{(-1)^{l}\,2^{n-2l}\lfloor\frac{n}{2}\rfloor!}{l!(n-2l)!}\frac{\Gamma(n-l+\alpha)}{\Gamma(\lceil\frac{n}{2}\rceil+\alpha)}c_{n-2\lfloor\frac{n}{2}\rfloor}^{(n)}.$

The typical choice for the coefficient
$c_{n-2\lfloor\frac{n}{2}\rfloor}^{(n)}=\frac{\Gamma(\lceil\frac{n}{2}\rceil+\alpha)}{\lfloor\frac{n}{2}\rfloor!\Gamma(\alpha)}$
removes clutter depending on integer division of $n$ and keeps the values of
the Gamma functions $\frac{\Gamma(n-l+\alpha)}{\Gamma(\alpha)}$ finite, or
more slim with the factorial $(\alpha)^{\overline{n-l}}$ rising from
$(\alpha)$ to $(\alpha+n-l-1)$

$\displaystyle c_{n-2l}^{(n)}$
$\displaystyle=\frac{(-1)^{l}\,2^{n-2l}}{l!(n-2l)!}\frac{\Gamma(n-l+\alpha)}{\Gamma(\alpha)}=\frac{(-1)^{l}\,2^{n-2l}(\alpha)^{\overline{n-l}}}{l!(n-2l)!},$
(71)

defining the Gegenbauer polynomials [Geg77]

$\displaystyle C_{n}^{(\alpha)}(x)$
$\displaystyle=\sum_{l=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{l}\,2^{n-2l}(\alpha)^{\overline{n-l}}}{l!(n-2l)!}\,x^{n-2l},$
(72)

and for $\alpha=0$ the rising factorial simplifies to
$(\alpha)^{\overline{n-l}}=(n-l-1)!$.

##  Appendix D Orthogonality of axisymmetric harmonics

The axisymetric Laplacian differential equation
$(1-x^{2})\mathcal{P}^{\prime\prime}_{n}-(\mathrm{D}-1)\,x\,\mathcal{P}^{\prime}_{n}+\lambda_{n}\,\mathcal{P}_{n}=0$
with $\lambda_{n}=n(n+2\alpha)$ can be re-written compactly as

$\displaystyle\frac{1}{w(x)}\left[(1-x^{2})\,w(x)\,\mathcal{P}^{\prime}_{n}\right]^{\prime}+\lambda_{n}\mathcal{P}_{n}=0.$
(73)

Multiplied by another solution $\mathcal{P}_{m}$ and integrated over
$w(x)\,\mathrm{d}x$ it becomes

$\displaystyle\int_{-1}^{1}\left[(1-x^{2})\,w(x)\,\mathcal{P}^{\prime}_{n}\right]^{\prime}\mathcal{P}_{m}\,\mathrm{d}x+\lambda_{n}\,\int_{-1}^{1}\mathcal{P}_{n}\,\mathcal{P}_{m}\,w(x)\,\mathrm{d}x=0.$

The first expression can be dealt with using partial integration $\int
f^{\prime}g=fg-\int fg^{\prime}:$

$\displaystyle\int_{-1}^{1}\left[(1-x^{2})\,w(x)\,\mathcal{P}^{\prime}_{n}\right]^{\prime}\mathcal{P}_{m}\,\mathrm{d}x$
$\displaystyle=(1-x^{2})\,w(x)\,\mathcal{P}^{\prime}_{n}\,\mathcal{P}_{m}\Big{|}_{-1}^{1}-\int_{-1}^{1}(1-x^{2})\,w(x)\,\mathcal{P}^{\prime}_{n}\,\mathcal{P}^{\prime}_{m}\,\mathrm{d}x.$

To show that the solutions $\mathcal{P}_{n}$ and $\mathcal{P}_{m}$ with $n\neq
m$ are orthogonal, we subtract a similar equation, the one with the roles of
$\mathcal{P}_{n}$ and $\mathcal{P}_{m}$ interchanged, and get:

$\displaystyle(1-x^{2})\,w(x)[\mathcal{P}_{n}^{\prime}\mathcal{P}_{m}-\mathcal{P}_{m}^{\prime}\mathcal{P}_{n}]\bigg{|}_{-1}^{1}+(\lambda_{n}-\lambda_{m})\,\int_{-1}^{1}\mathcal{P}_{n}\,\mathcal{P}_{m}\,w(x)\,\mathrm{d}x$
$\displaystyle=0,$ (74)

because the two integrals
$\int(1-x^{2})\,w(x)\,\mathcal{P}_{n}^{\prime}\mathcal{P}_{m}^{\prime}\mathrm{d}x$
are identical. With the integration boundaries $-1\leq x\leq 1$, the factor
$(1-x^{2})\,w(x)|_{\pm 1}=0$ in front of the mixed derivatives vanishes, which
happens for any $-\frac{1}{2}\geq\alpha$, i.e., $\mathrm{D}\geq 1$. We get:

$\displaystyle(\lambda_{n}-\lambda_{m})\,\int_{-1}^{1}\mathcal{P}_{n}\,\mathcal{P}_{m}\,w(x)\,\mathrm{d}x=0.$
(75)

If the indices match $n=m$, the first term vanishes and the integral needs not
vanish. For different indices $n\neq m$ our eigenvalues are unique as
$\lambda_{n}=n(n+\mathrm{D}-2)$, hence $\lambda_{n}\neq\lambda_{m}$, so the
first term won’t vanish, consequently the integral must vanish:

Any pair of solutions $\mathcal{P}_{n}$ and $\mathcal{P}_{m}$ with different
index $n\neq m$ is orthogonal under the integration weight
$w(x)=(1-x^{2})^{\alpha-\frac{1}{2}}$.

##  Appendix E Orthogonal polynomials

A set of orthogonal polynomials $\mathcal{P}_{n}(x)$ sorted by their degree
$n$, up to which they contain powers of the variable
$\mathcal{P}_{n}=\sum_{k=0}^{n}c_{n,k}\,x^{k}$ can naturally be used to re-
expand any other limited-degree polynomial
$\mathscr{P}_{n}=\sum_{k=0}^{n}c_{k}\,x^{k}$:

$\displaystyle\mathscr{P}_{n}(x)=\sum_{m=0}^{n}\gamma_{m}\,\mathcal{P}_{m}(x).$
(76)

This is not surprising. What makes orthogonal polynomials useful is their
orthogonality on the interval $[a,b]$:

$\displaystyle\int_{a}^{b}\mathcal{P}_{n}(x)\mathcal{P}_{m}(x)\;w(x)\,\mathrm{d}x=\begin{cases}0\qquad\text{for
$n\neq m$},\\\ N_{n}^{2}\qquad\text{for $n=m$}.\end{cases}$ (77)

Here, $N_{n}^{2}=\int_{a}^{b}\mathcal{P}_{n}^{2}(x)\;w(x)\,\mathrm{d}x$ is the
squared magnitude of each function integrated over the entire interval.

Note that by the ability of limited-degree polynomials $\mathcal{P}_{m<n}$ to
express any $m^{\mathrm{th}}$-degree polynomial with $m<n$, the orthogonality
automatically extendeds to _any_ lower-degree polynomial
$\mathscr{P}_{m<n}=\sum_{l=0}^{m}\gamma_{l}\,\mathcal{P}_{l}$

$\displaystyle\int_{a}^{b}\mathcal{P}_{n}(x)\mathscr{P}_{m<n}(x)\;w(x)\,\mathrm{d}x=0.$
(78)

###  E.1 The transform integral: a benefit of orthogonality

The orthogonality allows to decompose any given function $g(x)$ into
coefficients $\gamma_{m}$ of the orthogonal polynomials rather easily by an
integral transformation. We write down the decomposition that we desire and
integrate over a test function $\mathcal{P}_{n}(x)w(x)\mathrm{d}x\int$.
Orthogonality simplifies the expression that we get

$\displaystyle g(x)$
$\displaystyle=\sum_{m=0}^{\infty}\gamma_{m}\,\mathcal{P}_{m}(x)\qquad\qquad\Big{|}\;\mathcal{P}_{n}(x)w(x)\mathrm{d}x\int_{a}^{b}$
$\displaystyle\int_{a}^{b}g(x)\,\mathcal{P}_{n}(x)\;w(x)\,\mathrm{d}x$
$\displaystyle=\sum_{m=0}^{\infty}\gamma_{m}\,\underbrace{\int_{a}^{b}\mathcal{P}_{m}(x)\mathcal{P}_{n}(x)\;w(x)\,\mathrm{d}x}_{N_{n}^{2}\delta_{nm}}=N_{n}^{2}\;\gamma_{n},$

We hereby found the expansion coefficients by the integral without matrix
inversion. The transform pair is:

$\displaystyle g(x)$
$\displaystyle=\sum_{m=0}^{\infty}\gamma_{m}\,\mathcal{P}_{m}(x)\qquad\qquad\gamma_{n}=\frac{1}{N_{n}^{2}}\;\int_{a}^{b}g(x)\,\mathcal{P}_{n}(x)\;w(x)\,\mathrm{d}x.$
(79)

As infinitely many orthogonal polynomials consist of any power $x^{k}$ with
$k=0,1,2,\dots,\infty$, the representation is able to represent any function
that can be described by an infinite power series.

###  E.2 Symetries and interval $-1\leq x\leq 1$.

If $-a=b=1$ and $w(x)=w(-x)$, the integration range and integration weight is
symmetrical. We can re-write the orthogonality condition

$\displaystyle\int_{-1}^{1}\mathcal{P}_{n}\,\mathcal{P}_{m}\;w\,\mathrm{d}x$
$\displaystyle=\phantom{-}\int_{-1}^{0}\mathcal{P}_{n}(x)\,\mathcal{P}_{m}(x)\;w(x)\,\mathrm{d}x+\int_{0}^{1}g(x)\,\mathcal{P}_{n}(x)\;w(x)\,\mathrm{d}x$
$\displaystyle=-\int_{0}^{-1}\mathcal{P}_{n}(x)\,\mathcal{P}_{m}(x)\;w(x)\,\mathrm{d}x+\int_{0}^{1}\mathcal{P}_{n}(x)\,\mathcal{P}_{m}(x)\;w(x)\,\mathrm{d}x$
$\displaystyle=-\int_{0}^{1}\mathcal{P}_{n}(-x)\,\mathcal{P}_{m}(-x)\;w(-x)\,\mathrm{d}(-x)+\int_{0}^{1}\mathcal{P}_{n}(x)\,\mathcal{P}_{m}(x)\;w(x)\,\mathrm{d}x$
$\displaystyle=\phantom{-}\int_{0}^{1}\mathcal{P}_{n}(-x)\,\mathcal{P}_{m}(-x)\;w(x)\,\mathrm{d}x+\int_{0}^{1}\mathcal{P}_{n}(x)\,\mathcal{P}_{m}(x)\;w(x)\,\mathrm{d}x$
$\displaystyle=\phantom{-}\int_{0}^{1}\underbrace{\Bigl{[}\mathcal{P}_{n}(-x)\,\mathcal{P}_{m}(-x)+\mathcal{P}_{n}(x)\,\mathcal{P}_{m}(x)\Bigr{]}}_{=0\text{
for $n\neq m$}}\;w(x)\,\mathrm{d}x=N_{n}^{2}\,\delta_{nm}.$ (80)

The expression in brackets
$[\mathcal{P}_{n}(-x)\,\mathcal{P}_{m}(-x)+\mathcal{P}_{n}(x)\,\mathcal{P}_{m}(x)]$
vanishes whenever exclusively one of the two functions is odd, either
$\mathcal{P}_{n}(-x)=-\mathcal{P}_{n}(x)$ or
$\mathcal{P}_{m}(-x)=-\mathcal{P}_{m}(x)$, while the other one is even. The
highest power in a polynomial corresponds to its degree, odd powers of
$x^{2k+1}$ to odd function, even powers of $x^{2k}$ to even functions, so
consequently orthogonal polynomials for the interval $-1\leq x\leq 1$ and
weight $w(-x)=w(x)$ alternate between even and odd in $m$ to obtain
orthogonality: Odd-degree orthogonal polynomials contain only odd powers,
even-degree orthogonal polynomials only even powers.

###  E.3 Three-term recurrence of orthogonal polynomials

Normally, we define $\mathcal{P}_{0}(x)=1$ and any higher-degree polynomial by
a simple recurrence relation, hence based on the lower degree polynomials.
Finding the recurrence relation requires some consideration.

Orthogonality yields that we don’t need an orthogonal polynomial of the
$n^{\mathrm{th}}$ degree $\mathcal{P}_{n}(x)$ to represent an arbitrary lower-
degree polynomial $\mathscr{P}_{m<n}(x)$. This becomes clear when expanding
the polynomial
$\mathscr{P}_{m<n}(x)=\sum_{k=0}^{m<n}\beta_{k}\,\mathcal{P}_{k}(x)$ and
considering the transform integral with $\mathcal{P}_{n}$ vanishing
$\int\mathcal{P}_{k<n}\mathcal{P}_{n}=0$ for every term because of
orthogonality. We may write

$\displaystyle\int_{a}^{b}\mathscr{P}_{m<n}(x)\,\mathcal{P}_{n}(x)\;w(x)\,\mathrm{d}x$
$\displaystyle=0.$ (81)

After multiplication with $x$, an orthogonal polynomial of degree $n-1$
becomes an $n^{\mathrm{th}}$-degree polynomial

$\displaystyle x\,\mathcal{P}_{n-1}(x)=\mathscr{P}_{n}(x).$ (82)

We can find the expansion coefficients, $\beta_{m}$, thereof into orthogonal
polynomials by integrating $x\,\mathcal{P}_{n-1}(x)$ over
$\mathcal{P}_{m}(x)$. By applying the knowledge from above, we see which
coefficients must vanish:

$\displaystyle\beta_{m<n-2}$
$\displaystyle=\frac{1}{N_{m<n-2}^{2}}\int_{a}^{b}\overbrace{\mathcal{P}_{m<n-2}(x)\,x}^{\mathscr{P}_{m<n-1}(x)}\,\mathcal{P}_{n-1}(x)\;w(x)\,\mathrm{d}x=0,$
$\displaystyle\beta_{m>n}$
$\displaystyle=\frac{1}{N_{m>n}^{2}}\int_{a}^{b}\mathcal{P}_{m>n}(x)\,\underbrace{x\,\mathcal{P}_{n-1}(x)}_{\mathscr{P}_{n}(x)}\;w(x)\,\mathrm{d}x=0.$

Obviously, $x\,\mathcal{P}_{n-1}(x)$ only consists of $\mathcal{P}_{m}(x)$
with $m=n-2,n-1,n$: _Any set of orthogonal polynomials obeys a three-term
recurrence:_

$\displaystyle x\,\mathcal{P}_{n-1}(x)$
$\displaystyle=\beta_{n}^{(n)}\,\mathcal{P}_{n}(x)+\beta_{n-1}^{(n)}\,\mathcal{P}_{n-1}(x)+\beta_{n-2}^{(n)}\,\mathcal{P}_{n-2}(x)$
(83)

that can be used as
$\mathcal{P}_{n}(x)=\frac{1}{\beta_{n}^{(n)}}\Bigl{[}\bigl{(}x-\beta_{n-1}^{(n)}\bigr{)}\,\mathcal{P}_{n-1}(x)-\beta_{n-2}^{(n)}\,\mathcal{P}_{n-2}(x)\Bigr{]}$
to calculate higher-degree polynomials based on the lower-degree ones. The
recurrence reduces for $-a=b=1$ and $w(-x)=w(x)$ as this requires
$\mathcal{P}_{n}$ to only contain even or odd powers of $x$,

$\displaystyle
x\,\mathcal{P}_{n-1}(x)=\beta_{n}^{(n)}\,\mathcal{P}_{n}(x)+\beta_{n-2}^{(n)}\,\mathcal{P}_{n-2}(x).$
(84)

###  E.4 Explicit three-term recurrence

For the Gegenbauer polynomials
$C_{n}^{(\alpha)}(x)=\sum_{l=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{l}\,2^{n-2l}(\alpha)^{\overline{n-l}}}{l!(n-2l)!}\,x^{n-2l}$
of eq. (72), we can obtain an explicit three-term recurrence by comparing the
highest powers; divided by common factors
$\frac{2^{n}(\alpha)^{\overline{n}}}{0!n!}$ and
$\frac{2^{n-1}(\alpha)^{\overline{n-1}}}{0!(n-2)!}$ in the columns on the
right

$\displaystyle C_{n+1}^{(\alpha)}$
$\displaystyle=\frac{2^{n+1}(\alpha)^{\overline{n+1}}}{0!(n+1)!}x^{n+1}$
$\displaystyle-\frac{2^{n-1}(\alpha)^{\overline{n}}}{1!(n-1)!}x^{n-1}$
$\displaystyle\dots\bigg{|}$ $\displaystyle\frac{2(\alpha+n)}{n+1}$
$\displaystyle x^{n+1}\bigg{|}$ $\displaystyle-(\alpha+n-1)$ $\displaystyle
x^{n-1}$ $\displaystyle xC_{n}^{(\alpha)}$
$\displaystyle=\frac{2^{n}(\alpha)^{\overline{n}}}{0!n!}x^{n+1}$
$\displaystyle-\frac{2^{n-2}(\alpha)^{\overline{n-1}}}{1!(n-2)!}x^{n-1}$
$\displaystyle\dots\bigg{|}$ $\displaystyle x^{n+1}\bigg{|}$
$\displaystyle-\frac{n-1}{2}$ $\displaystyle x^{n-1}$ $\displaystyle
C_{n-1}^{(\alpha)}$ $\displaystyle=$
$\displaystyle\frac{2^{n-1}(\alpha)^{\overline{n-1}}}{0!(n-2)!}x^{n-1}$
$\displaystyle\dots\bigg{|}$ $\displaystyle\phantom{x^{n+1}}\bigg{|}$
$\displaystyle x^{n-1}$

Combining $(n+1)C_{n+1}^{(\alpha)}-2(\alpha+n)xC_{n}^{(\alpha)}$ obviously
cancels the power $x^{n+1}$ and leaves
$-(n+1)(\alpha+n-1)+(n-1)(\alpha+n)=-2\alpha-n+1$ to cancel by
$C_{n-1}^{(\alpha)}$. We get

$\displaystyle
2(n+\alpha)\,x\,C_{n}^{(\alpha)}=(n+1)C_{n+1}^{(\alpha)}+(n+2\alpha-1)C_{n-1}^{(\alpha)}.$
(85)

#### For standardization to $\mathcal{P}_{n}(1)=1$

If we standardize the Gegenbauer polynomials to ultraspherical ones
$\mathcal{P}_{n}(1)=1$ for all $n$, eq. (84) imposes the condition that
$\beta_{n}^{(n)}+\beta_{n-2}^{(n)}=1$, so that
$\beta_{n-2}^{(n)}=1-\beta_{n}^{(n)}$. We skip the superscript and write
$\beta_{n}$ for the recurrence coefficient,

$\displaystyle
x\,\mathcal{P}_{n-1}(x)=\beta_{n}\,\mathcal{P}_{n}(x)+(1-\beta_{n})\,\mathcal{P}_{n-2}(x).$
(86)

We obtain the recurrence coefficient $\beta_{n}$ from comparing coefficients
of the highest powers in $\mathcal{P}_{n}$ and $\mathcal{P}_{n-1}$. This is
done based on the coefficient recurrence for $c_{n}^{(n)}$,

$\displaystyle\beta_{n}\,c_{n}^{(n)}$
$\displaystyle=c_{n-1}^{(n-1)}\Longrightarrow\beta_{n}=\frac{c_{n-1}^{(n-1)}}{c_{n}^{(n)}}.$

From above we know about the existence and nature of three-term recurrences.
For the polynomials with $\mathcal{P}_{n}(1)=1$ in eqs.(72) (98), we get a
ratio of the $l=0$ coefficients for $n$ and $n+1$ as
$\beta_{n+1}=\frac{c_{n}^{(n)}}{c_{n+1}^{(n+1)}}=\frac{2^{n}\,\Gamma(n+\alpha)}{\Gamma(n+2\alpha)}\frac{\Gamma(n+2\alpha+1)}{2^{n+1}\,\Gamma(n+\alpha+1)}=\frac{n+2\alpha}{2(n+\alpha)}$
and $1-\beta_{n+1}=\frac{n}{2(n+\alpha)}$ for a recurrence as in eq. (86)

$\displaystyle 2(n+\alpha)\,x\,\mathcal{P}_{n}$
$\displaystyle=(n+2\alpha)\mathcal{P}_{n+1}+n\mathcal{P}_{n-1}.$ (87)

###  E.5 Squared norm of orthogonal polynomials

According to Law and Sledd [LS75], the recurrence coefficient $\beta_{n}$ can
also be used to evaluate the normalization
$N_{n}^{2}=\int_{-1}^{1}[\mathcal{P}_{n}(x)]^{2}\,w(x)\,\mathrm{d}x$ of the
polynomial recursively. This is done by re-formulating the three-term
recurrence to

$\displaystyle\mathcal{P}_{n}$
$\displaystyle=(A_{n}x+B_{n})\mathcal{P}_{n-1}-C_{n}\mathcal{P}_{n-2}.$ (88)

Law and Sledd’s derivation is based on explicitly evaluating only the first
term of the recurrence until $\mathcal{P}_{0}=1$, while the rest must stay
some polynomial of the degree $n-1$

$\mathcal{P}_{n}=A_{n}x\mathcal{P}_{n-1}+\mathscr{A}_{n-1}=A_{n}A_{n-1}x^{2}\mathcal{P}_{n-2}+\mathscr{B}_{n-1}=A_{n}A_{n-1}\dots
A_{1}x^{n}+\mathscr{C}_{n-1},$

or more neatly

$\displaystyle x^{n}=(A_{n}A_{n-1}\dots
A_{1})^{-1}\mathcal{P}_{n}+\mathscr{P}_{n-1}.$ (89)

Integration over $\mathcal{P}_{n}\,w\,\mathrm{d}x\,\int_{a}^{b}$ yields

$\displaystyle\int_{a}^{b}x^{n}\mathcal{P}_{n}\,w\,\mathrm{d}x$
$\displaystyle=(A_{n}A_{n-1}\dots
A_{1})^{-1}\underbrace{\int_{a}^{b}(\mathcal{P}_{n})^{2}\,w\,\mathrm{d}x}_{=N_{n}^{2}}+\underbrace{\int_{a}^{b}\mathscr{P}_{n-1}\,\mathcal{P}_{n}\,w\,\mathrm{d}x}_{=0}$
$\displaystyle\int_{a}^{b}x^{n}\mathcal{P}_{n}\,w\,\mathrm{d}x$
$\displaystyle=(A_{n}A_{n-1}\dots A_{1})^{-1}\,N_{n}^{2}.$ (90)

$x^{n}\mathcal{P}_{n}$ is re-expressed by re-arranging eq. (88) to
$x\mathcal{P}_{n}=A_{n+1}^{-1}[\mathcal{P}_{n+1}-B_{n+1}\mathcal{P}_{n}+C_{n+1}\mathcal{P}_{n-1}]$
and multiplication with $x^{n-1}=(A_{n-1}\dots
A_{1})^{-1}[\mathcal{P}_{n-1}+\mathscr{P}_{n-2}]$ from eq. (89),

$\displaystyle\int_{a}^{b}x^{n}\mathcal{P}_{n}\,w\,\mathrm{d}x$
$\displaystyle=\frac{\int_{a}^{b}[\mathcal{P}_{n+1}-B_{n+1}\mathcal{P}_{n}+C_{n+1}\mathcal{P}_{n-1}][\mathcal{P}_{n-1}+\mathscr{P}_{n-2}]\,w\,\mathrm{d}x}{A_{n+1}(A_{n-1}\dots
A_{1})}$ $\displaystyle=\frac{C_{n+1}}{A_{n+1}(A_{n-1}\dots
A_{1})}\underbrace{\int_{a}^{b}(\mathcal{P}_{n-1})^{2}\,w\,\mathrm{d}x}_{=N_{n-1}^{2}}.$
(91)

Together with eq. (90), we get for

$\displaystyle N_{n}^{2}=\frac{C_{n+1}A_{n}}{A_{n+1}}N_{n-1}^{2},$ (92)

which is in our $\mathcal{P}_{n}(1)=1$ standardized case with
$\beta_{n}\mathcal{P}_{n}(x)=x\mathcal{P}_{n-1}-(1-\beta_{n})\mathcal{P}_{n-2}$
hence $A_{n}=\frac{1}{\beta_{n}}$, and $C_{n}=\frac{1-\beta_{n}}{\beta_{n}}$,
as well as $\mathcal{P}_{0}=1$ and
$\int_{-1}^{1}(\mathcal{P}_{0})^{2}\,w(x)\,\mathrm{d}x=\frac{S_{\mathrm{D}-1}}{S_{\mathrm{D}-2}}$,
and so

$\displaystyle N_{n}^{2}$
$\displaystyle=\frac{1-\beta_{n+1}}{\beta_{n}}\,N_{n-1}^{2}=\frac{S_{\mathrm{D}-1}}{S_{\mathrm{D}-2}}\prod_{k=1}^{n}\frac{1-\beta_{k+1}}{\beta_{k}};$
(93)

see Appdx. E.11 for more details.

###  E.6 Christoffel-Darboux formula

From the general three-term recurrence

$\displaystyle A_{n+1}x\mathcal{P}_{n}(x)$
$\displaystyle=\mathcal{P}_{n+1}(x)-B_{n+1}\mathcal{P}_{n}(x)+C_{n+1}\mathcal{P}_{n-1}(x)$

one finds the Cristoffel-Darboux formula by expanding with
$\mathcal{P}_{n}(y)$ and $(x-y)$

$\displaystyle A_{n+1}(x-y)\mathcal{P}_{n}(x)\mathcal{P}_{n}(y)$
$\displaystyle=[\mathcal{P}_{n+1}(x)\mathcal{P}_{n}(y)-\mathcal{P}_{n}(x)\mathcal{P}_{n+1}(y)]$
$\displaystyle+C_{n+1}[\mathcal{P}_{n-1}(x)\mathcal{P}_{n}(y)-\mathcal{P}_{n}(x)\mathcal{P}_{n-1}(y)].$

Multiplying with $\frac{1}{A_{n+1}N_{n}^{2}}$ yields with the normalization
recurrence $\frac{A_{n+1}}{C_{n+1}}N_{n}^{2}=A_{n}\,N_{n-1}^{2}$ from eq. (92)

$\displaystyle\frac{\mathcal{P}_{n}(x)\mathcal{P}_{n}(y)}{N_{n}^{2}}$
$\displaystyle=\underbrace{\frac{\mathcal{P}_{n+1}(x)\mathcal{P}_{n}(y)-\mathcal{P}_{n}(x)\mathcal{P}_{n+1}(y)}{(x-y)\,A_{n+1}N_{n}^{2}}}_{\mathcal{F}_{n}(x,y)}+\underbrace{\frac{\mathcal{P}_{n-1}(x)\mathcal{P}_{n}(y)-\mathcal{P}_{n}(x)\mathcal{P}_{n-1}(y)}{(x-y)\,A_{n}N_{n-1}^{2}}}_{-\mathcal{F}_{n-1}(x,y)}.$

Adding all expressions of smaller $n$ until $n=0$, where
$\mathcal{F}_{-1}(x,y)=0$, yields the Christoffel-Darboux formula decomposing
an axisymmetric Dirac delta $\delta(x-y)$

$\displaystyle\sum_{n=0}^{\mathrm{N}}\frac{\mathcal{P}_{n}(x)\mathcal{P}_{n}(y)}{N_{n}^{2}}$
$\displaystyle=\frac{\mathcal{P}_{\mathrm{N}+1}(x)\mathcal{P}_{\mathrm{N}}(y)-\mathcal{P}_{\mathrm{N}}(x)\mathcal{P}_{\mathrm{N}+1}(y)}{(x-y)\,A_{\mathrm{N}+1}N_{\mathrm{N}}^{2}}.$
(94)

If $y=x$, the fraction exhibits a pole and the confluent Christoffel-Darboux
identity is obtained by $y=x+\Delta x$ and taking the limit $\lim_{\Delta
x\rightarrow 0}$

$\displaystyle\sum_{n=0}^{\mathrm{N}}\frac{[\mathcal{P}_{n}(x)]^{2}}{N_{n}^{2}}$
$\displaystyle=\lim_{\Delta x\rightarrow
0}\frac{\mathcal{P}_{\mathrm{N}+1}(x)\mathcal{P}_{\mathrm{N}}(x+\Delta
x)-\mathcal{P}_{\mathrm{N}}(x)\mathcal{P}_{\mathrm{N}+1}(x+\Delta x)}{-\Delta
x\,A_{\mathrm{N}+1}N_{\mathrm{N}}^{2}}$
$\displaystyle=\frac{\mathcal{P}_{\mathrm{N}}(x)\mathcal{P}_{\mathrm{N}+1}^{\prime}(x)-\mathcal{P}_{\mathrm{N}+1}(x)\mathcal{P}_{\mathrm{N}}^{\prime}(x)}{A_{\mathrm{N}+1}N_{\mathrm{N}}^{2}},$
(95)

which is obviously finite. This equation is useful to simplify the finite-
degree representation of a Dirac delta at $y=1$, as, e.g., found in Rafaely’s
book [Raf19]. For normalization at $x=1$, all factors can be disregarded and
using the derivative at $x=1$ in eq. (99), we find
$\mathcal{P}_{\mathrm{N+1}}^{\prime}(1)-\mathcal{P}_{\mathrm{N}}^{\prime}(1)=\frac{(\mathrm{N}+1)(\mathrm{N}+\mathrm{D}-1)-\mathrm{N}(\mathrm{N}+\mathrm{D}-2)}{\mathrm{D}-1}=\frac{2\mathrm{N}+\mathrm{D}-1}{\mathrm{D}-1}$
so that $g(1)=1$ and elsewhere eq. (94) holds with the newly obtained factor

$\displaystyle g(x)$
$\displaystyle=\frac{\mathrm{D}-1}{2\mathrm{N}+\mathrm{D}-1}\frac{\mathcal{P}_{\mathrm{N}+1}(x)-\mathcal{P}_{\mathrm{N}}(x)}{x-1}.$
(96)

###  E.7 Polynomials $C_{n}^{(\alpha)}(x)$ and $\mathcal{P}_{n}(x)$ at $x=1$

To find the values $C_{n}^{(\alpha)}(1)$, we regard the three-term recurrence
eq. (85) for $x=1$

$\displaystyle(n+1)C_{n+1}^{(\alpha)}-2(\alpha+n)\,C_{n}^{(\alpha)}+(2\alpha+n-1)C_{n-1}^{(\alpha)}=0.$

As there is multiplication of $C_{n+1}^{(\alpha)}$ with $(n+1)$, which is a
hint for $\frac{1}{n!}$, and multiplication of $C_{n-2}^{(\alpha)}$ with
$(2\alpha+n-1)$, hinting towards $(2\alpha)^{\overline{n}}$, we test a
solution $C_{n}^{(\alpha)}=\frac{(2\alpha)^{\overline{n}}}{n!}$

$\displaystyle\frac{(2\alpha)^{\overline{n}}}{(n-1)!}-2(\alpha+n-1)\,\frac{(2\alpha)^{\overline{n-1}}}{(n-1)!}+\frac{(2\alpha)^{\overline{n-1}}}{(n-2)!}=0$
$\displaystyle(2\alpha+n-1)(2\alpha)^{\overline{n-1}}-(2\alpha+2n-2)\,(2\alpha)^{\overline{n-1}}+(n-1)(2\alpha)^{\overline{n-1}}=0.\qquad\blacksquare$

Therefore the value of the Gegenbauer polynomials at $x=1$ become

$\displaystyle C_{n}^{(\alpha)}(1)$
$\displaystyle=\frac{(2\alpha)^{\overline{n}}}{n!},$ $\displaystyle\text{and
}C_{n}^{(\alpha)}(x)$
$\displaystyle=\frac{(2\alpha)^{\overline{n}}}{n!}\mathcal{P}_{n}^{(\alpha)}(x).$
(97)

Division of eq. (72) by $\frac{(2\alpha)^{\overline{n}}}{n!}$ lets us also get
an explicit power-series expansion of $\mathcal{P}_{n}$

$\displaystyle\mathcal{P}_{n}(x)$
$\displaystyle=\sum_{l=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{l}2^{n-2l}\,n!}{l!(n-2l)!}\frac{(\alpha)^{\overline{n-l}}}{(2\alpha)^{\overline{n}}}\,x^{n-2l},$
(98)

which yields numerical difficulties at $\alpha\rightarrow 0$. We re-formulate
the rising factorials as
$(\alpha)^{\overline{n-l}}=\Gamma(\alpha+n-l)|_{0}/\Gamma(\alpha)=(n-l-1)!/\Gamma(\alpha)$
and
$(2\alpha)^{\overline{n}}=\Gamma(2\alpha+n)|_{0}/\Gamma(2\alpha)=(n-1)!/\Gamma(2\alpha)$,
use duplication [OBC00, 5.5.E5]
$\Gamma(2z)/\Gamma(z)=2^{2z-1/2}\Gamma(z+1/2)/\sqrt{2\pi}$, and get with
$\Gamma(1/2)=\sqrt{\pi}$ the limit $\lim_{\alpha\rightarrow
0}\frac{\Gamma(2\alpha)}{\Gamma(\alpha)}=1/2$, and so $\lim_{\alpha\rightarrow
0}\frac{(\alpha)^{\overline{n-l}}}{(2\alpha)^{\overline{n}}}=\frac{(n-l-1)!}{(n-1)!}$,

$\displaystyle\mathcal{P}_{n}(x)$ $\displaystyle=\lim_{\alpha\rightarrow
0}\frac{\Gamma(2\alpha)}{\Gamma(\alpha)}\sum_{l=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{l}2^{n-2l}\,n!}{l!(n-2l)!}\frac{(n-l-1)!}{(n-1)!}\,x^{n-2l},$
$\displaystyle=\sum_{l=0}^{\lfloor\frac{n}{2}\rfloor}\frac{(-1)^{l}\,2^{n-2l-1}n!\;(n-l-1)!}{l!(n-2l)!}\,x^{n-2l}.$

###  E.8 Derivative of the polynomials $\mathcal{P}_{n}(x)$ at $x=1$

The ultraspherical differential equation
$(1-x^{2})\,\mathcal{P}_{n}^{\prime\prime}-(2\alpha+1)\,x\,\mathcal{P}_{n}^{\prime}+n(n+2\alpha)\,\mathcal{P}_{n}=0$
becomes at $x=1$ with $(1-x^{2})=0$ and $\mathcal{P}_{n}(1)=1$

$\displaystyle\mathcal{P}_{n}^{\prime}(1)$
$\displaystyle=\frac{n(n+2\alpha)}{2\alpha+1}=\frac{n(n+\mathrm{D}-2)}{\mathrm{D}-1}.$
(99)

###  E.9 Value of the polynomials $C^{(\alpha)}_{n}(x)$ and
$\mathcal{P}_{n}(x)$ at $x=0$

By the explicit series expansions, we can evaluate the Gegenbauer polynomials
at $x=0$. Only such with even degree have a non-zero value there:

$\displaystyle C_{2m}^{(\alpha)}(0)$
$\displaystyle=\sum_{l=0}^{m}\frac{(-1)^{l}(\alpha)^{\overline{2m-l}}}{l!(2m-2l)!}\underbrace{(2x)^{n-2l}\bigg{|}_{0}}_{\neq
0:l=m}=\frac{(-1)^{m}\,(\alpha)^{\overline{m}}}{m!}.$ (100)

With the ultraspherical polynomials $\mathcal{P}_{n}$, this becomes by
division with $\frac{(2\alpha)^{\overline{n}}}{n!}$

$\displaystyle\mathcal{P}_{2m}^{(\alpha)}(0)$
$\displaystyle=\frac{(-1)^{m}\,(2m)!(\alpha)^{\overline{m}}}{m!(2\alpha)^{\overline{2m}}}.$
(101)

###  E.10 Rodrigues’ formula for Gegenbauer polynomials.

The French mathematician Rodrigues provided a formula for the Legendre
polynomials that are orthogonal using the integral measure $w(x)=1$. A general
Rodrigues formula for a more arbitrary measure
$w(x)=(1-x^{2})^{\alpha-\frac{1}{2}}$ can be given to specify the
corresponding orthogonal polynomials $\mathcal{P}_{n}$.

The ultraspherical (Gegenbauer) differential equation for axisymmetric
harmonic solutions in $\alpha=\frac{\mathrm{D}}{2}-1$ dimensions is

$\displaystyle(1-x^{2})\mathcal{P}_{n}^{\prime\prime}(x)-(2\alpha+1)x\mathcal{P}_{n}^{\prime}(x)+n(n+2\alpha)\mathcal{P}_{n}(x)$
$\displaystyle=0,$ (102)

and it is solved by the generalized Rodrigues formula

$\displaystyle\mathcal{P}_{n}(x)=e_{n}^{-1}(1-x^{2})^{-\alpha+{\textstyle\frac{1}{2}}}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}(1-x^{2})^{\alpha-{\textstyle\frac{1}{2}}+n}.$

#### Proof.

To show this, we start with the inner-most function

$f_{n}=(1-x^{2})^{\alpha-{\textstyle\frac{1}{2}}+n}$

that the solution is constructed of and derive it

$f_{n}^{\prime}=(n+\alpha-{\textstyle\frac{1}{2}})(1-x^{2})^{\alpha-\frac{1}{2}+n-1}(-2x)=-2(n+\alpha-{\textstyle\frac{1}{2}})x(1-x^{2})^{-1}f_{n},$

yielding a differential equation after multiplication with $(1-x^{2})$:

$(1-x^{2})f_{n}^{\prime}+(2n+2\alpha-1)xf_{n}=0.$

Deriving again yields a second-order differential equation

$\displaystyle(1-x^{2})f_{n}^{\prime\prime}-2xf_{n}^{\prime}+(2n+2\alpha-1)xf_{n}^{\prime}+(2n+2\alpha-1)f_{n}$
$\displaystyle=0$
$\displaystyle(1-x^{2})f_{n}^{\prime\prime}+(2n+2\alpha-3)xf_{n}^{\prime}+(2n+2\alpha-1)f_{n}$
$\displaystyle=0.$

Now this second-order differential equation in
$f_{n}=(1-x^{2})^{\alpha-\frac{1}{2}+n}$ is derived $n$ times requiring for
the first two expressions the Leibniz rule for repeated differentiation of
products
$\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}(f\,g)=\sum_{k=0}^{n}\frac{n!}{k!(n-k)!}f^{(k)}\,g^{(n-k)}$,

$\displaystyle(1-x^{2})f_{n}^{(n+2)}-2nxf_{n}^{(n+1)}-n(n-1)f_{n}^{(n)}$
$\displaystyle+(2n+2\alpha-3)xf_{n}^{(n+1)}+n(2n+2\alpha-3)f_{n}^{(n)}$
$\displaystyle+(2n+2\alpha-1)f_{n}^{(n)}$ $\displaystyle=0$
$\displaystyle(1-x^{2})f_{n}^{(n+2)}+(2\alpha-3)xf_{n}^{(n+1)}+\left[n(n+2\alpha)+(2\alpha-1)\right]f_{n}^{(n)}$
$\displaystyle=0.$

The above differential equation is second-order for
$f_{n}^{(n)}=\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}(1-x^{2})^{n}$, which is
related to the Rodrigues solution by
$f_{n}^{(n)}=e_{n}(1-x^{2})^{\alpha-\frac{1}{2}}\,\mathcal{P}_{n}$, which upon
insertion and omission of $e_{n}$ becomes

$\displaystyle(1-x^{2})[(1-x^{2})^{\alpha-\frac{1}{2}}\,\mathcal{P}_{n}]^{\prime\prime}+(2\alpha-3)x[(1-x^{2})^{\alpha-\frac{1}{2}}\,\mathcal{P}_{n}]^{\prime}$
$\displaystyle+\left[n(n+2\alpha)+(2\alpha-1)\right](1-x^{2})^{\alpha-\frac{1}{2}}\,\mathcal{P}_{n}$
$\displaystyle=0$
$\displaystyle(1-x^{2})[(1-x^{2})^{\alpha-\frac{1}{2}}\,\mathcal{P}_{n}^{\prime\prime}-2(2\alpha-1)x(1-x^{2})^{\alpha-\frac{3}{2}}\mathcal{P}_{n}^{\prime}\quad$
$\displaystyle-(2\alpha-1)(1-x^{2})^{\alpha-\frac{3}{2}}\mathcal{P}_{n}+(2\alpha-1)(2\alpha-3)x^{2}(1-x^{2})^{\alpha-\frac{5}{2}}\mathcal{P}_{n}]$
$\displaystyle+(2\alpha-3)x[(1-x^{2})^{\alpha-\frac{1}{2}}\,\mathcal{P}_{n}^{\prime}-(2\alpha-1)x(1-x^{2})^{\alpha-\frac{3}{2}}\,\mathcal{P}_{n}]$
$\displaystyle+\left[n(n+2\alpha)+(2\alpha-1)\right](1-x^{2})^{\alpha-\frac{1}{2}}\,\mathcal{P}_{n}(x)$
$\displaystyle=0\quad\Big{|}(1-x^{2})^{-\alpha+\frac{1}{2}}$
$\displaystyle(1-x^{2})[\mathcal{P}_{n}^{\prime\prime}-2(2\alpha-1)x(1-x^{2})^{-1}\mathcal{P}_{n}^{\prime}\quad$
$\displaystyle-\cancel{(2\alpha-1)(1-x^{2})^{-1}\mathcal{P}_{n}}+\cancel{(2\alpha-1)(2\alpha-3)x^{2}(1-x^{2})^{-2}\mathcal{P}_{n}}]$
$\displaystyle+(2\alpha-3)x[\mathcal{P}_{n}^{\prime}-\cancel{(2\alpha-1)x(1-x^{2})^{-1}\,\mathcal{P}_{n}}]$
$\displaystyle+\left[n(n+2\alpha)+\cancel{(2\alpha-1)}\right]\,\mathcal{P}_{n}(x)$
$\displaystyle=0$
$\displaystyle(1-x^{2})\mathcal{P}_{n}^{\prime\prime}(x)-(2\alpha+1)x\mathcal{P}_{n}^{\prime}(x)+n(n+2\alpha)\mathcal{P}_{n}(x)$
$\displaystyle=0.$

As this is the ultraspherical (Gegenbauer) differential equation, and we
obtained it by the $(n+2)^{\mathrm{th}}$ derivative of
$f_{n}=(1-x^{2})^{\alpha-\frac{1}{2}+n}$ and its substitution
$f_{n}^{(n)}=e_{n}(1-x^{2})^{\alpha-\frac{1}{2}}\mathcal{P}_{n}$, we hereby
proved Rodrigues’ solution
$\mathcal{P}_{n}=e_{n}^{-1}(1-x^{2})^{-\alpha+\frac{1}{2}}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}(1-x^{2})^{\alpha-\frac{1}{2}+n}$.

To find the normalization factor $e_{n}$ for standardization of
$\mathcal{P}_{n}(1)$, here for $\mathcal{P}_{n}(1)=1$, we apply Leibniz’s rule
for repeated differentiation

$\displaystyle\mathcal{P}_{n}(x)$
$\displaystyle=e_{n}^{-1}(1-x^{2})^{-\alpha+\frac{1}{2}}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}[(1-x)^{\alpha-\frac{1}{2}+n}(1+x)^{\alpha-\frac{1}{2}+n}]$
$\displaystyle=e_{n}^{-1}\sum_{k=0}^{n}\frac{n!}{k!(n-k)!}\frac{(1+x)^{-\alpha+\frac{1}{2}}\mathrm{d}^{k}(1+x)^{\alpha-\frac{1}{2}+n}}{\mathrm{d}x^{k}}\frac{(1-x)^{-\alpha+\frac{1}{2}}\mathrm{d}^{n-k}(1-x)^{\alpha-\frac{1}{2}+n}}{\mathrm{d}x^{n-k}}$
$\displaystyle=e_{n}^{-1}\sum_{k=0}^{n}\frac{n!}{k!(n-k)!}\frac{\Gamma(\alpha+n+\frac{1}{2})(1+x)^{n-k}}{\Gamma(\alpha+n-k+\frac{1}{2})}\underbrace{\frac{\Gamma(\alpha+n+\frac{1}{2})(-1)^{n-k}(1-x)^{k}}{\Gamma(\alpha+k+\frac{1}{2})}}_{\neq
0:x=1:k=0}$ $\displaystyle\mathcal{P}_{n}(x)\Big{|}_{1}$
$\displaystyle=e_{n}^{-1}\,\frac{\Gamma(\alpha+n+\frac{1}{2})(2)^{n}(-1)^{n}}{\Gamma(\alpha+\frac{1}{2})}=1$
$\displaystyle e_{n}$
$\displaystyle=\frac{(2)^{n}\Gamma(\alpha+n+\frac{1}{2})}{(-1)^{n}\Gamma(\alpha+\frac{1}{2})}=\frac{2^{n}(\alpha+\frac{1}{2})^{\overline{n}}}{(-1)^{n}}.$
(103)

Hereby, Rodrigues’ formula for Gegenbauer polynomials standardized to
ultraspherical ones $\mathcal{P}_{n}(1)=1$ becomes:

$\displaystyle\mathcal{P}_{n}(x)=\frac{(-1)^{n}}{2^{n}(\alpha+\frac{1}{2})^{\overline{n}}}(1-x^{2})^{-\alpha+\frac{1}{2}}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}(1-x^{2})^{\alpha-\frac{1}{2}+n},$
(104)

and alternatively we could write with the integration weight
$w(x)=(1-x^{2})^{\alpha-\frac{1}{2}}$ inserted
$\mathcal{P}_{n}(x)=\frac{(-1)^{n}}{(2\alpha+1)^{\overline{n}}\,w(x)}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}\left[(1-x^{2})^{n}\,w(x)\right]$.

The typical Gegenbauer standardization is
$C_{n}^{(\alpha)}(1)=\frac{\Gamma(n+2\alpha)}{n!\Gamma(2\alpha)}=\frac{(2\alpha)^{\overline{n}}}{n!}$
as in eq. (97), and its Rodrigues’ formula is obtained by multiplying the
above with $\frac{(2\alpha)^{\overline{n}}}{n!}$, cf. eq. (97),

$\displaystyle
C_{n}^{(\alpha)}(x)=\frac{(-1)^{n}(2\alpha)^{\overline{n}}}{n!2^{n}(\alpha+\frac{1}{2})^{\overline{n}}}(1-x^{2})^{-\alpha+\frac{1}{2}}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}(1-x^{2})^{\alpha-\frac{1}{2}+n}.$
(105)

###  E.11 Norm of Gegenbauer polynomials

For the transform integral, the squared norm
$N_{n}^{2}=\int_{-1}^{1}|\mathcal{P}_{n}(x)|^{2}\,w(x)\,\mathrm{d}x$ plays an
important role for normalization. As the zeroth-degree polynomial is
$C_{0}^{(\alpha)}(x)=\mathcal{P}_{0}(x)=1$, its norm
$N_{0}^{2}=\int_{-1}^{1}w(x)\,\mathrm{d}x$ should be easiest to find. The
others follow from eq. (92).

#### Finding $N_{0}^{2}$.

Using all coordinates of the $\mathrm{D}$-dimensional space, the integral over
$x$ the last one of multiple integrals. In the sequence of the
$(\mathrm{D}-1)$-fold integration over $1$ to determine the surface
$S_{\mathrm{D}-1}$ of the $\mathbb{S}^{\mathrm{D}-1}$ unit sphere in a
$\mathrm{D}$-dimensional space, the integral over $x$ delivers the
distinguishing factor between the surfaces $S_{\mathrm{D}-1}$ and
$S_{\mathrm{D}-2}$, the one of the next lower-dimensional unit sphere, cf.
[HC37, vol.2, p.387],

$\displaystyle S_{\mathrm{D}-1}$
$\displaystyle=\underbrace{\iiint\,\mathrm{d}S_{\mathrm{D}-2}}_{=S_{\mathrm{D}-2}}\cdot\,\underbrace{\int_{-1}^{1}\,\sqrt{1-x^{2}}^{\mathrm{D}-3}\,\mathrm{d}x}_{=N_{0}^{2}}.$
(106)

With
$S_{\mathrm{D}-1}=\frac{2\,\sqrt{\pi}^{\mathrm{D}}}{\Gamma(\frac{\mathrm{D}}{2})}$
[HC37, vol.2, p.387], we conclude that

$\displaystyle N_{0}^{2}$
$\displaystyle=\frac{S_{\mathrm{D}-1}}{S_{\mathrm{D}-2}}=\frac{\sqrt{\pi}\,\Gamma(\frac{\mathrm{D}-1}{2})}{\Gamma(\frac{\mathrm{D}}{2})}.$
(107)

For 2 dimensions $\mathrm{D}=2$ we get $N_{0}^{2}=\frac{2\pi}{2}=\pi$, and for
3 dimensions $\mathrm{D}=3$ we get $N_{0}^{2}=\frac{4\pi}{2\pi}=2$, for both
kinds $\mathcal{P}_{n}$ and $C_{n}^{(\alpha)}$.

#### Squared norm for any degree $N_{n}^{2}$ for $C_{n}^{(\alpha)}$.

eq. (92) obtains the norm $N_{n}^{2}=\frac{C_{n+1}A_{n}}{A_{n+1}}N_{n-1}^{2}$
of higher-degree polynomials recursively from $N_{0}^{2}$. We use the
recurrence
$C_{n}^{(\alpha)}=\underbrace{\textstyle\frac{2(n+\alpha-1)}{n}}_{A_{n}}x\,C_{n-1}^{(\alpha)}-\underbrace{\textstyle\frac{n+2\alpha-2}{n}}_{C_{n}}\,C_{n-2}^{(\alpha)}$
to find the squared norm, and $\alpha=\frac{\mathrm{D}}{2}-1$,

$\displaystyle N_{n}^{2}$
$\displaystyle=\frac{C_{n+1}A_{n}}{A_{n+1}}N_{n-1}^{2}=\frac{\frac{n+2\alpha-1}{n+1}\frac{2(n+\alpha-1)}{n}}{\frac{2(n+\alpha)}{n+1}}N_{n-1}^{2}=\frac{(n+2\alpha-1)(2n+2\alpha-2)}{n(2n+2\alpha)}N_{n-1}^{2}$
$\displaystyle=\frac{(n+2\alpha-1)\cancel{(2n+2\alpha-2)}}{n(2n+2\alpha)}\frac{(n+2\alpha-2)(2n+2\alpha-4)}{(n-1)\cancel{(2n+2\alpha-2)}}N_{n-2}^{2}$
$\displaystyle=\frac{(n+2\alpha-1)!}{n!(2n+2\alpha)}N_{0}^{2}=\frac{(n+\mathrm{D}-3)!}{n!(2n+\mathrm{D}-2)}\frac{S_{\mathrm{D}-1}}{S_{\mathrm{D}-2}},$
(108)

what yields for $C_{n}^{(0)}$ the norm $N_{n}^{2}=\pi(2-\delta_{n,0})$ when
$\mathrm{D}=2$ and $N_{n}^{2}=\frac{2}{2n+1}$ for $\mathrm{D}=3$.

#### Squared norm for any degree $N_{n}^{2}$ for $\mathcal{P}_{n}$.

For the $\mathcal{P}_{n}(1)=1$ ultraspherical polynomials with
$\mathcal{P}_{n}=\underbrace{\textstyle\frac{2(n+\alpha-1)}{n+2\alpha-1}}_{A_{n}}x\,\mathcal{P}_{n-1}-\underbrace{\textstyle\frac{n-1}{n+2\alpha-1}}_{C_{n}}\,\mathcal{P}_{n-2}$
this is for $n>0$

$\displaystyle N_{n}^{2}$
$\displaystyle=\frac{C_{n+1}A_{n}}{A_{n+1}}N_{n-1}^{2}=\frac{\frac{n}{n+2\alpha}\frac{2(n+\alpha-1)}{n+2\alpha-1}}{\frac{2(n+\alpha)}{n+2\alpha}}N_{n-1}^{2}=\frac{n(n+\alpha-1)}{(n+\alpha)(n+2\alpha-1)}N_{n-1}^{2}$
$\displaystyle=\frac{n!(\alpha)^{\overline{n}}N_{0}^{2}}{(\alpha+1)^{\overline{n}}(2\alpha)^{\overline{n}}}=\frac{n!}{(2n+\mathrm{D}-2)\;(n+\mathrm{D}-3)!}\,\frac{S_{\mathrm{D}-1}}{S_{\mathrm{D}-2}},$
(109)

yielding $N_{n}^{2}=\frac{\pi}{2-\delta_{n,0}}$ for $\mathrm{D}=2$ and
$N_{n}^{2}=\frac{2}{2n+1}$ for $\mathrm{D}=3$;
$N_{0}^{2}=\frac{S_{D-1}}{S_{D-2}}$ can be omitted.

###  E.12 Integral recurrence of the Gegenbauer polynomials

With the power-series representation of the Gegenbauer polynomials
$C_{n}^{(\alpha)}$ in eq. (72), we find derivative recurrences after deriving
$C_{n+1}^{(\alpha)}$ and $C_{n-1}^{(\alpha)}$. We separating the factors
distinguishing the series to the coefficients of $C_{n}^{(\alpha)}$, after an
index shift in $l\rightarrow l-1$ for the derivative of $C_{n-1}^{(\alpha)}$:

$\displaystyle C_{n+1}^{(\alpha)^{\prime}}$
$\displaystyle=\sum_{l=0}^{\lfloor\frac{n+1}{2}\rfloor}\frac{(-1)^{l}2^{n-2l+1}(\alpha)^{\overline{n-l+1}}\cancel{(n-2l+1)}x^{n-2l}}{l!(n-2l\cancel{+1})!}$
$\displaystyle=\sum_{l=0}^{\lfloor\frac{n+1}{2}\rfloor}2(\alpha+n-l)\,\frac{(-1)^{l}2^{n-2l}(\alpha)^{\overline{n-l}}}{l!(n-2l)!}\,x^{n-2l}$
$\displaystyle C_{n-1}^{(\alpha)^{\prime}}$
$\displaystyle=\sum_{l=0}^{\lfloor\frac{n-1}{2}\rfloor}\frac{(-1)^{l}2^{n-2l-1}(\alpha)^{\overline{n-l-1}}\cancel{(n-2l-1)}x^{n-2l-2}}{l!(n-2l-\cancel{1}^{2})!}$
$\displaystyle=\sum_{l=\cancel{1}^{0}}^{\lfloor\frac{n+1}{2}\rfloor}(-1)2l\frac{(-1)^{l}2^{n-2l}(\alpha)^{\overline{n-l}}x^{n-2l}}{l!(n-2l)!}.$

Combining $C_{n+1}^{(\alpha)^{\prime}}-C_{n-1}^{(\alpha)^{\prime}}$ obviously
leaves

$\displaystyle
2(\alpha+n)\,\sum_{l=0}^{\lfloor\frac{n+1}{2}\rfloor}\frac{(-1)^{l}2^{n-2l}(\alpha)^{\overline{n-l}}}{l!(n-2l)!}\,x^{n-2l}=2(\alpha+n)\,C_{n}^{(\alpha)},$

and hence

$\displaystyle 2(\alpha+n)C_{n}^{(\alpha)}$
$\displaystyle=[C_{n+1}^{(\alpha)}-C_{n-1}^{(\alpha)}]^{\prime}.$ (110)

Similarly, we find for the Gegenbauer polynomials standardized to
$\mathcal{P}_{n}(1)=1$ eq. (98)

$\displaystyle\mathcal{P}^{\prime}_{n+1}$
$\displaystyle=\sum_{l=0}^{\lfloor\frac{n+1}{2}\rfloor}\frac{2(n+1)(\alpha+n-l)}{2\alpha+n}\,\frac{(-1)^{l}2^{n-2l}\,n!}{l!(n-2l)!}\frac{(\alpha)^{\overline{n-l}}}{(2\alpha)^{\overline{n}}}\,x^{n-2l}$
(111) $\displaystyle\mathcal{P}^{\prime}_{n-1}$
$\displaystyle=\sum_{l=\cancel{1}^{0}}^{\lfloor\frac{n+1}{2}\rfloor}\frac{(-1)2l(2\alpha+n-1)}{n}\,\frac{(-1)^{l}2^{n-2l}\,n!}{l!(n-2l)!}\frac{(\alpha)^{\overline{n-l}}}{(2\alpha)^{\overline{n}}}\,x^{n-2l},$

leaving

$\displaystyle
2(\alpha+n)\mathcal{P}_{n}=\Big{[}\frac{2\alpha+n}{n+1}\mathcal{P}_{n+1}-\frac{n}{2\alpha+n-1}\mathcal{P}_{n-1}\Big{]}^{\prime},$
(112)

which leaves for $\alpha=\frac{1}{2}$ and Legendre polynomials
$(2n+1)P_{n}=[P_{n+1}-P_{n-1}]^{\prime}$, and $\alpha=0$ with Chebyshev
polynomials $2T_{n}=[\frac{1}{n+1}T_{n+1}-\frac{1}{n-1}T_{n-1}]^{\prime}$. By
this, we get the recurrence for integration

$\displaystyle
2(n+\alpha)\int\mathcal{P}_{n}\mathrm{d}x=\frac{n+2\alpha}{n+1}\mathcal{P}_{n+1}-\frac{n}{n+2\alpha-1}\mathcal{P}_{n-1}.$
(113)

###  E.13 Additional recurrence relations

From the definition of the differential equation (102), we directly get

$\displaystyle(1-x^{2})\,\mathcal{P}_{n}^{\prime\prime}$
$\displaystyle=(1+2\alpha)\,x\,\mathcal{P}_{n}^{\prime}-n(n+2\alpha)\,\mathcal{P}_{n}.$
(114)

From partial integration of $x\,\mathcal{P}_{n}^{\prime}$, we get via the
multiplication recurrence (87) and integration recurrence (113):

$\displaystyle\int x\,\mathcal{P}_{n}^{\prime}\,\mathrm{d}x$
$\displaystyle=x\,\mathcal{P}_{n}-\int\mathcal{P}_{n}\,\mathrm{d}x$
$\displaystyle=\big{[}{\textstyle\frac{n+2\alpha}{2(n+\alpha)}}-{\textstyle\frac{n+2\alpha}{2(n+\alpha)(n+1)}}\big{]}\mathcal{P}_{n+1}+\big{[}{\textstyle\frac{n}{2(n+\alpha)}}+{\textstyle\frac{n}{2(n+\alpha)(n+2\alpha-1)}}\big{]}\mathcal{P}_{n-1}$
$\displaystyle 2(n+\alpha)\int x\,\mathcal{P}_{n}^{\prime}\,\mathrm{d}x$
$\displaystyle=n(n+2\alpha)\bigg{[}\frac{1}{n+1}\mathcal{P}_{n+1}+\frac{1}{n+2\alpha-1}\mathcal{P}_{n-1}\bigg{]}.$
(115)

From partial integration of the differential equation, we get

$\displaystyle\int(1-x^{2})\,\mathcal{P}_{n}^{\prime\prime}\,\mathrm{d}x-(2\alpha+1)\,\int
x\,\mathcal{P}_{n}^{\prime}\,\mathrm{d}x+n(n+2\alpha)\,\int\mathcal{P}_{n}\,\mathrm{d}x$
$\displaystyle=0$
$\displaystyle\textstyle\big{[}(1-x^{2})\mathcal{P}_{n}^{\prime}+2\int
x\,\mathcal{P}_{n}^{\prime}\,\mathrm{d}x\big{]}-(2\alpha+1)\int
x\,\mathcal{P}_{n}^{\prime}\,\mathrm{d}x+n(n+2\alpha)\int\mathcal{P}_{n}\,\mathrm{d}x$
$\displaystyle=0$
$\displaystyle\textstyle(1-x^{2})\mathcal{P}_{n}^{\prime}-(2\alpha-1)\int
x\,\mathcal{P}_{n}^{\prime}\,\mathrm{d}x+n(n+2\alpha)\int\mathcal{P}_{n}\,\mathrm{d}x$
$\displaystyle=0$
$\displaystyle\textstyle(1-x^{2})\mathcal{P}_{n}^{\prime}-(2\alpha-1)\big{[}x\,\mathcal{P}_{n}-\int\mathcal{P}_{n}\,\mathrm{d}x\big{]}+n(n+2\alpha)\int\mathcal{P}_{n}\,\mathrm{d}x$
$\displaystyle=0$
$\displaystyle\textstyle(1-x^{2})\mathcal{P}_{n}^{\prime}-(2\alpha-1)\,x\,\mathcal{P}_{n}+[n(n+2\alpha)+(2\alpha-1)]\int\mathcal{P}_{n}\,\mathrm{d}x$
$\displaystyle=0$
$\displaystyle\textstyle(1-x^{2})\mathcal{P}_{n}^{\prime}-(2\alpha-1)\,x\,\mathcal{P}_{n}+(n+1)(n+2\alpha-1)\int\mathcal{P}_{n}\,\mathrm{d}x$
$\displaystyle=0$ $\displaystyle(1-x^{2})\mathcal{P}_{n}^{\prime}$
$\displaystyle=(2\alpha-1)\,x\,\mathcal{P}_{n}-(n+1)(n+2\alpha-1)\int\mathcal{P}_{n}\,\mathrm{d}x$
$\displaystyle=\big{[}{\textstyle\frac{(2\alpha-1)(n+2\alpha)}{2(n+\alpha)}}-{\textstyle\frac{\cancel{(n+1)}(n+2\alpha-1)(n+2\alpha)}{2(n+\alpha)\cancel{(n+1)}}}\big{]}\mathcal{P}_{n+1}+\big{[}{\textstyle\frac{n(2\alpha-1)}{2(n+\alpha)}}+{\textstyle\frac{n(n+1)\cancel{(n+2\alpha-1)}}{2(n+\alpha)\cancel{(n+2\alpha-1)}}}\big{]}\mathcal{P}_{n-1}$
$\displaystyle 2(n+\alpha)(1-x^{2})\mathcal{P}_{n}^{\prime}$
$\displaystyle=n(n+2\alpha)\big{[}-\mathcal{P}_{n+1}+\mathcal{P}_{n-1}\big{]}$
(116)

##  Appendix F Designs and metrics $P$, $E$, $Q$, $FBR$,
$\bm{r}_{\mathrm{V}}$, and $\bm{r}_{\mathrm{E}}$ for axisymmetric
directivity/panning patterns

The panning function or directivity pattern

$\displaystyle g(x)$
$\displaystyle=\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{S_{\mathrm{D}-2}N_{n}^{2}}\,\mathcal{P}_{n}(x)$
(117)

is rotationally symmetric around the $x$ axis of a $\mathrm{D}$-dimensional
space. For the scalar measures, we integrate over $x$ with the surface weight
$w(x)=\sqrt{1-x^{2}}^{\mathrm{D}-3}$, and we use the surface $S_{D-2}$ of the
unit sphere in the $(\mathrm{D}-1)$-dimensional space to denote the result of
the integrals independent of $x$:

$\displaystyle P$
$\displaystyle=S_{\mathrm{D}-2}\,\int_{-1}^{1}\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{S_{\mathrm{D}-2}\,N_{n}^{2}}\,\mathcal{P}_{n}(x)\,w(x)\mathrm{d}x=\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\overbrace{\int_{-1}^{1}\mathcal{P}_{n}(x)\,w(x)\mathrm{d}x}^{N_{n}^{2}\delta_{n0}}=a_{0}$
(118) $\displaystyle E$
$\displaystyle=S_{\mathrm{D}-2}\sum_{n,n^{\prime}}\frac{a_{n^{\prime}}^{*}}{S_{\mathrm{D}-2}N_{n^{\prime}}^{2}}\frac{a_{n}}{S_{\mathrm{D}-2}N_{n}^{2}}\overbrace{\int_{-1}^{1}\mathcal{P}_{n}(x)\mathcal{P}_{n^{\prime}}(x)\,w(x)\mathrm{d}x}^{N_{n}^{2}\delta_{nn^{\prime}}}=\sum_{n=0}^{\mathrm{N}}\frac{|a_{n}|^{2}}{S_{\mathrm{D}-2}N_{n}^{2}},$
(119) $\displaystyle Q$
$\displaystyle=\frac{S_{D-1}\left(\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{S_{D-2}\,N_{n}^{2}}\right)^{2}}{E},\qquad\quad
FBR=\frac{\sum_{n^{\prime}=0}^{\mathrm{N}}\sum_{n=0}^{\mathrm{N}}g_{nn^{\prime}}a_{n}a_{n^{\prime}}^{*}}{\sum_{n^{\prime}=0}^{\mathrm{N}}\sum_{n=0}^{\mathrm{N}}g_{nn^{\prime}}a_{n}a_{n^{\prime}}^{*}(-1)^{n+n^{\prime}}},$
(120)

with the Gram matrix $g_{nn^{\prime}}$ of the polynomials integrated over the
frontal half-space, only

$\displaystyle g_{nn^{\prime}}$
$\displaystyle=\frac{1}{N_{n}N_{n^{\prime}}}\int_{0}^{1}\mathcal{P}_{n}(x)\mathcal{P}_{n^{\prime}}(x)\,w(x)\,\mathrm{d}x,$
(121)

see subsection F.1, also [Raf19, Elk00, Elk04] for more information.

For rotationally symmetric shapes in space, the vector models always indicate
the direction of axial symmetry. The remaining integral for
$\bm{r}_{\mathrm{V}}$ and $\bm{r}_{\mathrm{E}}$ is a one-dimensional
calculation along $x$ and determines the length.

$\displaystyle r_{\mathrm{V}}$
$\displaystyle=\frac{S_{\mathrm{D}-2}\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}S_{\mathrm{D}-2}}\overbrace{\int_{-1}^{1}\mathcal{P}_{n}(x)\,x\,w(x)\mathrm{d}x}^{N_{1}^{2}\delta_{n,1}}}{P}=\frac{a_{1}}{a_{0}}.$
(122)

For the $\bm{r}_{\mathrm{V}}$, the expression below the brace was evaluated by
inserting the recurrence of the Gegenbauer polynomials
$x\mathcal{P}_{n}(x)=\dots$,

$\displaystyle\int_{-1}^{1}[(1-\beta_{n+1})P_{n-1}(x)+\beta_{n+1}P_{n+1}(x)]w(x)\mathrm{d}x=(1-\beta_{2})N_{0}^{2}\delta_{n,1}=N_{1}^{2}\delta_{n,1},$
(123)

and by noticing that one can re-write $(1-\beta_{2})N_{0}^{2}$ to
$\beta_{1}N_{1}^{2}$ because of the recurrence of the normalization in eq.
(92), $N_{1}^{2}=\frac{1-\beta_{2}}{\beta_{1}}N_{0}^{2}$, and further to
$N_{1}^{2}$ given that $\beta_{1}=1$.

For $\bm{r}_{\mathrm{E}}$, the expression with the integral yields

$\displaystyle r_{\mathrm{E}}$
$\displaystyle=\frac{S_{\mathrm{D}-2}\,\sum_{n,n^{\prime}}\frac{a_{n}}{S_{\mathrm{D}-2}N_{n}^{2}}\frac{a_{n^{\prime}}^{*}}{S_{\mathrm{D}-2}N_{n^{\prime}}^{2}}\overbrace{\int_{-1}^{1}\mathcal{P}_{n}(x)\mathcal{P}_{n^{\prime}}(x)\,x\,w(x)\,\mathrm{d}x}}{E}=\frac{\sum_{n}\frac{\beta_{n+1}}{N_{n}^{2}}2\Re\\{a_{n}a_{n+1}^{*}\\}}{\sum_{n=0}^{\mathrm{N}-1}\frac{1}{N_{n}^{2}}|a_{n}|^{2}}.$
(124)

To evaluate the expression below the brace, the recurrence was used

$\displaystyle\textstyle\int_{-1}^{1}\mathcal{P}_{n}(x)\bigl{[}(1-\beta_{n^{\prime}+1})\mathcal{P}_{n^{\prime}-1}(x)$
$\displaystyle+\beta_{n^{\prime}+1}P_{n^{\prime}+1}(x)\bigr{]}w(x)\mathrm{d}x$
$\displaystyle=N_{n}^{2}(1-\beta_{n^{\prime}+1})\delta_{n+1,n^{\prime}}+N_{n}^{2}\beta_{n^{\prime}+1}\delta_{n,n^{\prime}+1},$

and re-insertion gave

$\displaystyle r_{\mathrm{E}}$
$\displaystyle=\frac{\sum_{n=0}^{\mathrm{N}-1}\frac{a_{n}a_{n+1}^{*}}{S_{\mathrm{D}-2}N_{n+1}^{2}}(1-\beta_{n+2})+\sum_{n^{\prime}=0}^{\mathrm{N}-1}\frac{a_{n^{\prime}+1}a_{n^{\prime}}^{*}}{S_{\mathrm{D}-2}N_{n^{\prime}+1}^{2}}\beta_{n^{\prime}+1}}{E}=\frac{\sum_{n^{\prime}=0}^{\mathrm{N}-1}\frac{a_{n}a_{n+1}^{*}+a_{n}^{*}a_{n+1}}{S_{\mathrm{D}-2}N_{n^{\prime}}^{2}}\beta_{n+1}}{E},$
(125)

for which eq. (92) was used to re-write
$\frac{1-\beta_{n+2}}{N_{n+1}^{2}}=\frac{\beta_{n+1}}{N_{n}^{2}}$ in the first
numerator term, and a shift of $n$ by one in the second one.

###  F.1 Supercardioid: Gram matrix for half interval

For the design of a supercardioid pattern, the entries of the Gram matrix are
needed

$\displaystyle g_{nm}$
$\displaystyle=\int_{x_{0}}^{1}\frac{\mathcal{P}_{n}(x)}{N_{n}^{2}}\frac{\mathcal{P}_{m}(x)}{N_{m}^{2}}\,w(x)\,\mathrm{d}x.$
(126)

As
$\mathcal{P}_{n}(x)\mathcal{P}_{m}(x)=(-1)^{n+m}\mathcal{P}_{n}(-x)\mathcal{P}_{m}(-x)$
uses $(-1)^{2(l+k)}$ or $(-1)^{2(l+k)+2}$ there is no sign change between the
half intervals for matching parity. In these cases, eq. (126) multiplied by
$\frac{1}{2}$ and integrated over the full range $\int_{-1}^{1}$ is equivalent
and orthogonality preserved, therefore and the entries
$g_{2l,2k}=g_{2l+1,2k+1}=0$ for $l\neq k$. The entries of the diagonal $n=m$
is the half-interval integral of the squared function, therefore half the
squared normalization $\frac{N_{n}^{2}}{2}$ divided by $N_{n}^{4}$. Omitting
the unit-sphere ratio $\frac{S_{D-1}}{S_{D-2}}$ by setting $N_{0}^{2}=1$ in
eq. (109) for simplicity, we get

$\displaystyle g_{nn}$
$\displaystyle=\frac{1}{2\,N_{n}^{2}}=\frac{(\alpha+1)^{\overline{n}}(2\alpha)^{\overline{n}}}{2\;n!(\alpha)^{\overline{n}}}.$
(127)

eq. (74) gives a hint for $n\neq m$, using
$w(x)=(1-x^{2})^{\alpha-\frac{1}{2}}$, so
$(1-x^{2})\,w(x)=(1-x^{2})^{\alpha+\frac{1}{2}}$, which vanishes at the upper
bound $x=1$:

$\displaystyle
g_{nm}=\int_{x_{0}}^{1}\frac{\mathcal{P}_{n}(x)}{N_{n}^{2}}\frac{\mathcal{P}_{m}(x)}{N_{m}^{2}}\,w(x)\,\mathrm{d}x$
$\displaystyle=\frac{(1-x^{2})\,w(x)}{N_{n}^{2}N_{m}^{2}}\,\frac{\mathcal{P}_{n}^{\prime}\mathcal{P}_{m}-\mathcal{P}_{m}^{\prime}\mathcal{P}_{n}}{n(n+2\alpha)-m(m+2\alpha)}\bigg{|}_{x_{0}}^{1}$
(128)
$\displaystyle=\frac{-(1-x_{0}^{2})^{\alpha+\frac{1}{2}}}{N_{n}^{2}N_{m}^{2}}\,\frac{\mathcal{P}_{n}^{\prime}(x_{0})\mathcal{P}_{m}(x_{0})-\mathcal{P}_{m}^{\prime}(x_{0})\mathcal{P}_{n}(x_{0})}{n(n+2\alpha)-m(m+2\alpha)}.$

At the lower bound $x_{0}=0$, the factor $(1-x^{2})^{\alpha+\frac{1}{2}}=1$ is
unity, and we get $\mathcal{P}_{2l+1}(0)=0$ because of the odd parity and
$\mathcal{P}_{2k}^{\prime}(0)=0$ due to the even parity

$\displaystyle g_{2l+1,2k}$
$\displaystyle=\frac{-1}{N_{n}^{2}N_{m}^{2}}\,\frac{\mathcal{P}_{2l+1}^{\prime}(0)\mathcal{P}_{2k}(0)-\mathcal{P}_{2k}^{\prime}(0)\mathcal{P}_{2l+1}(0)}{(2l+1)(2l+2\alpha+1)-(2k)(2k+2\alpha)}$
$\displaystyle=\frac{1}{N_{2l+1}^{2}N_{2k}^{2}}\,\frac{-\mathcal{P}_{2l+1}^{\prime}(0)\mathcal{P}_{2k}(0)}{(2l+1)(2l+2\alpha+1)-(2k)(2k+2\alpha)}.$
(129)

With
$\mathcal{P}_{2k}(0)=\frac{(-1)^{k}(2k)!(\alpha)^{\overline{k}}}{k!\,(2\alpha)^{\overline{2k}}}$
and
$\mathcal{P}_{n+1}^{\prime}(x)=\sum_{l}\frac{(-1)^{l}2^{n-2l+1}(n+1)!(\alpha)^{\overline{n-l+1}}}{l!(n-2l)!\,(2\alpha)^{\overline{n+1}}}x^{n-2l}$
from eqs. (101) (111), with $n=2l$ and $x=0$, so
$\mathcal{P}_{2l+1}^{\prime}(0)=\frac{2\,(-1)^{l}(2l+1)!(\alpha)^{\overline{l+1}}}{l!\,(2\alpha)^{\overline{2l+1}}}$
we have with $\frac{1}{N_{n}^{2}}$,$\frac{1}{N_{m}^{2}}$

$\displaystyle g_{2l+1,2k}$
$\displaystyle=\textstyle\frac{1}{N_{0}^{2}}\frac{(\alpha+1)^{\overline{2l+1}}\cancel{(2\alpha)^{\overline{2l+1}}}}{\cancel{(2l+1)!}(\alpha)^{\overline{2l+1}}}\frac{2\,(-1)^{l}\cancel{(2l+1)!}(\alpha)^{\overline{l+1}}}{l!\,\cancel{(2\alpha)^{\overline{2l+1}}}}\frac{(\alpha+1)^{\overline{2k}}\cancel{(2\alpha)^{\overline{2k}}}}{\cancel{(2k)!}(\alpha)^{\overline{2k}}}\frac{(-1)^{k}\cancel{(2k)!}(\alpha)^{\overline{k}}}{k!\,\cancel{(2\alpha)^{\overline{2k}}}}\cdots\,$
$\displaystyle\qquad\textstyle\cdots\frac{-1}{(2l+1)(2l+2\alpha+1)-(2k)(2k+2\alpha)}$
$\displaystyle=\textstyle\frac{\Gamma(\alpha+1)}{\sqrt{\pi}\Gamma(\alpha+\frac{1}{2})}\frac{(-1)^{l+1}(\alpha+1)^{\overline{2l+1}}}{l!\,(\alpha+1+l)^{\overline{l}}}\frac{(-1)^{k}(\alpha+1)^{\overline{2k}}}{k!\,(\alpha+k)^{\overline{k}}}\,\frac{2}{(2l+1)(2l+2\alpha+1)-(2k)(2k+2\alpha)}.$
(130)

Interchanging $n\leftrightarrow m$ in eq. (126) yields the values
$g_{2l+1,2k}=g_{2k,2l+1}$ on the other side of the diagonal. _(No idea from
where the $N_{0}^{2}=\frac{S_{\mathrm{D-1}}}{S_{\mathrm{D-2}}}$ sneaked in,
but it was numerically tested to be necessary for
$\mathrm{D}=\\{2,2.1,2.5,3,3.5,4,4.5,5\\}$) to restore balance between
diagonal/offdiagonal)._

###  F.2 Basic pattern with maximum DI/maximum $\bm{r}_{\mathrm{V}}$ vector

With $a_{n}=1$, the value $r_{\mathrm{V}}=\frac{a_{1}}{a_{0}}=1$. The maximum
directivity factor $Q$ is obtained, because $a_{n}=1$ decomposes a Dirac
delta, which practically minimizes the squares of the directivity everywhere
except at $x=1$. The resulting value of the directivity factor is

$\displaystyle\mathrm{max}Q$
$\displaystyle=S_{\mathrm{D}-2}\frac{(\sum_{n=0}^{\mathrm{N}}\frac{1}{N_{n}^{2}})^{2}}{\sum_{n=0}^{\mathrm{N}}\frac{1}{N_{n}^{2}}}=\sum_{n=0}^{\mathrm{N}}\frac{S_{\mathrm{D}-2}}{N_{n}^{2}}=\begin{cases}\sum_{n=0}^{\mathrm{N}}(2-\delta_{n})=2\mathrm{N}+1,&\mathrm{D}=2\\\
\sum_{n=0}^{\mathrm{N}}(2n+1)=(\mathrm{N}+1)^{2},&\mathrm{D}=3.\end{cases}$
(131)

###  F.3 Pattern with maximum $\bm{r}_{\mathrm{E}}$ vector

The energy vector length is determined by the real-valued coefficients $a_{n}$
in terms of
$r_{\mathrm{E}}=\frac{\sum_{n=0}^{\mathrm{N}-1}\frac{2\beta_{n+1}}{N_{n}^{2}}a_{n}a_{n+1}}{\sum_{n=0}^{\mathrm{N}}\frac{1}{N_{n}^{2}}a_{n}^{2}}=\frac{A}{E}$.
We may zero its derive it with regard to $a_{n}$ to maximize its length,
yielding
$r_{\mathrm{E}}^{\prime}=\frac{A^{\prime}}{E}-\frac{A}{E^{2}}E^{\prime}=\frac{A^{\prime}-r_{\mathrm{E}}\,E^{\prime}}{E}=0$,
where $A/E=r_{\mathrm{E}}$ was re-inserted as described by Daniel [DRP99,
Dan01]. It is sufficient to zero $A^{\prime}-r_{\mathrm{E}}\,E^{\prime}=0$ and
insert the derivatives of the numerator and denominator, $A^{\prime}$,
$E^{\prime}$, respectively. For the numerator $A^{\prime}$

$\displaystyle\frac{\partial}{\partial
a_{n}}\sum_{n^{\prime}}\frac{2\beta_{n^{\prime}+1}a_{n^{\prime}}a_{n^{\prime}+1}}{N_{n^{\prime}}^{2}}=\frac{2\beta_{n+1}a_{n+1}}{N_{n}^{2}}+\frac{2\beta_{n}a_{n-1}}{N_{n-1}^{2}}=\frac{2}{N_{n}^{2}}\left[\beta_{n+1}a_{n+1}+(1-\beta_{n+1})a_{n-1}\right]$

the equivalence
$\frac{\beta_{n}}{N_{n-1}^{2}}=\frac{1-\beta_{n+1}}{N_{n}^{2}}$ from eq. (92)
helped to unify expressions. For the denominator $E^{\prime}$ we get
$\frac{\partial}{\partial
a_{n}}\sum_{n^{\prime}}\frac{a_{n^{\prime}}^{2}}{N_{n^{\prime}}^{2}}=\frac{2}{N_{n}^{2}}a_{n}$,
and we insert $A^{\prime}$ and $E^{\prime}$ into
$A^{\prime}-r_{\mathrm{E}}\,E^{\prime}=0$

$\displaystyle\beta_{n+1}a_{n+1}+(1-\beta_{n+1})a_{n-1}-r_{\mathrm{E}}\,a_{n}=0.$
(132)

We recognize the recurrence relation
$\beta_{n+1}\mathcal{P}_{n+1}+(1-\beta_{n+1})\mathcal{P}_{n-1}-x\,\mathcal{P}_{n}(x)=0$
of the orthogonal polynomials and hereby find $a_{n}=\mathcal{P}_{n}(x)$ with
$r_{\mathrm{E}}=x$ as solution scheme

$\displaystyle a_{n}=\mathcal{P}_{n}(r_{\mathrm{E}}).$ (133)

The difference equation does not consider the finite-degree summation
boundaries of $r_{\mathrm{E}}$ in eq. (125), yet. Only for an infinite sum,
$r_{\mathrm{E}}=1$ is feasible. Otherwise, feasible choices for finite-degree
$r_{\mathrm{E}}$ are found after inserting the scheme into
$r_{\mathrm{E}}=A/E$ in the form $A=r_{\mathrm{E}}E$ of eq. (125) considering
the finite range of $n\in[0;\mathrm{N}-1]$ for $A$ and $n\in[0;\mathrm{N}]$
for $E$,

$\displaystyle\bigg{[}\sum_{n=0}^{\mathrm{N}-1}\frac{\mathcal{P}_{n}(r_{\mathrm{E}})\mathcal{P}_{n+1}(r_{\mathrm{E}})}{S_{\mathrm{D}-2}N_{n+1}^{2}}(1-\beta_{n+2})+\sum_{n^{\prime}=0}^{\mathrm{N}-1}\frac{\mathcal{P}_{n^{\prime}+1}(r_{\mathrm{E}})\mathcal{P}_{n^{\prime}}(r_{\mathrm{E}})}{S_{\mathrm{D}-2}N_{n^{\prime}+1}^{2}}\beta_{n^{\prime}}\bigg{]}$
$\displaystyle=r_{\mathrm{E}}\bigg{[}\sum_{n=0}^{\mathrm{N}}\frac{\mathcal{P}_{n}(r_{\mathrm{E}})\mathcal{P}_{n}(r_{\mathrm{E}})}{S_{\mathrm{D}-2}N_{n}^{2}}\bigg{]}$

We’d nearly precisely get the same right-hand-side term after expanding one of
the $\mathcal{P}_{n}(r_{\mathrm{E}})$ terms in the sum by the factor
$r_{\mathrm{E}}$ from outside the sum, using the recurrence
$r_{\mathrm{E}}\,\mathcal{P}_{n}(r_{\mathrm{E}})=\beta_{n+1}\mathcal{P}_{n+1}(r_{\mathrm{E}})+(1-\beta_{n+1})\mathcal{P}_{n-1}(r_{\mathrm{E}})$.
And yet, there is a decisive difference: no sum cuts away the degree
$\mathrm{N}+1$ expression, which wasn’t part of the sums on the left.
Therefore, a feasible finite-degree solution for $r_{\mathrm{E}}$ requires a
vanishing term

$\displaystyle\beta_{\mathrm{N}+1}\,\mathcal{P}_{\mathrm{N}+1}(r_{\mathrm{E}})\mathcal{P}_{\mathrm{N}}(r_{\mathrm{E}})=0.$
(134)

Because $a_{\mathrm{N}}\neq 0$ should not vanish and $\beta_{\mathrm{N}+1}\neq
0$, the feasibility constraint is that

$\displaystyle\arg\max\mathcal{P}_{\mathrm{N}+1}(r_{\mathrm{E}})=0.$ (135)

Therefore, as defined in [DRP99, Dan01], the largest root of
$\mathcal{P}_{\mathrm{N}+1}$ defines the coefficients $a_{n}$ for a maximum-
length $\bm{r}_{\mathrm{E}}$ vector.

Interestingly, this corresponds to decomposing an axisymmetric Dirac delta
$\delta(x-\max r_{\mathrm{E}})$: a concentric ring (3D) or a pair of points
(2D) around the panning direction, which could be expressed in short by the
Christoffel-Darboux identity eq. (94) with $\mathcal{P}_{\mathrm{N}+1}(\max
r_{\mathrm{E}})=0$

$\displaystyle g(\cos\vartheta)$
$\displaystyle=\frac{\mathcal{P}_{\mathrm{N}+1}(\cos\vartheta)\mathcal{P}_{\mathrm{N}}(\max
r_{\mathrm{E}})}{\cos\vartheta-\max r_{\mathrm{E}}}.$ (136)

###  F.4 Particular max-$r_{\mathrm{E}}$ solutions for $\mathrm{D}=3$

A fairly good approximation for practical use has been given in [ZF12] by

$\displaystyle r_{\mathrm{E}}=\cos\frac{137.9^{\circ}}{\mathrm{N}+1.51},$
(137)

and it saves us the time for polynomial root finding.

###  F.5 (Ultra-)Spherical cap function

For the coefficients of a spherical cap function, one needs to integrate the
polynomial over $w(x)$ and the corresponding segment of the sphere $x_{0}\leq
x\leq 1$.

The coefficient needed for expanding a spherical cap function is obtained by
integrating $a_{n}=\int_{x_{0}}^{1}\mathcal{P}_{n}(x)\,w(x)\,\mathrm{d}x$,
with the weight function $w(x)=(1-x^{2})^{\alpha-\frac{1}{2}}$.

We can manipulate Rodrigues’ formula eq. (104) in the form with $n+1$

$\displaystyle\frac{2^{n+1}(\alpha+\frac{1}{2})^{\overline{n+1}}}{(-1)^{n+1}}\;w(x)\,\mathcal{P}_{n+1}(x)$
$\displaystyle=\frac{\mathrm{d}^{n+1}}{\mathrm{d}x^{n+1}}(1-x^{2})^{\alpha+n+\frac{1}{2}},$

by carrying out one of the derivatives on the right to get

$\displaystyle\frac{2^{n+1}(\alpha+\frac{1}{2})^{\overline{n+1}}}{(-1)^{n+1}}\,$
$\displaystyle
w(x)\,\mathcal{P}_{n+1}(x)=(\alpha+n+{\textstyle\frac{1}{2}})\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}(-2x)(1-x^{2})^{\alpha+n-\frac{1}{2}},$
$\displaystyle=-(2\alpha+2n+1)\bigg{[}\underbrace{x\frac{\mathrm{d}^{n}(1-x^{2})^{\alpha+n-\frac{1}{2}}}{\mathrm{d}x^{n}}}_{\frac{2^{n}(\alpha+\frac{1}{2})^{\overline{n}}}{(-1)^{n}}x\,w(x)\,\mathcal{P}_{n}(x)}+n\,\underbrace{\frac{\mathrm{d}^{n-1}(1-x^{2})^{\alpha+n-\frac{1}{2}}}{\mathrm{d}x^{n-1}}}_{\frac{2^{n}(\alpha+\frac{1}{2})^{\overline{n}}}{(-1)^{n}}\int\mathcal{P}_{n}(x)\,w(x)\,\mathrm{d}x}\bigg{]}$
$\displaystyle\frac{2\alpha+2n+1}{2\alpha+2n+1}\;w(x)\,$
$\displaystyle\mathcal{P}_{n+1}(x)=x\,w(x)\,\mathcal{P}_{n}(x)+n\,\int\mathcal{P}_{n}(x)\,w(x)\,\mathrm{d}x$
$\displaystyle\mathcal{P}_{n+1}(x)=x\,\mathcal{P}_{n}(x)+\frac{n}{w(x)}\,\int\mathcal{P}_{n}(x)\,w(x)\,\mathrm{d}x$
(138) $\displaystyle\Longrightarrow\int\mathcal{P}_{n}(x)\,w(x)\,\mathrm{d}x$
$\displaystyle=\frac{w(x)}{n}\,\bigg{[}\mathcal{P}_{n+1}(x)-x\,\mathcal{P}_{n}(x)\bigg{]}.$
(139)

The recurrence
$x\mathcal{P}_{n}(x)=\frac{n+2\alpha}{2n+2\alpha}\mathcal{P}_{n+1}(x)+\frac{n}{2n+2\alpha}\mathcal{P}_{n-1}(x)$
yields the coefficients

$\displaystyle a_{n}=\int_{x_{0}}^{1}\mathcal{P}_{n}(x)\,w(x)\,\mathrm{d}x$
$\displaystyle=\frac{w(x)}{2n+2\alpha}\,\bigg{[}\mathcal{P}_{n+1}(x)-\mathcal{P}_{n-1}(x)\bigg{]}\Bigg{|}_{x_{0}}^{1}$
(140)
$\displaystyle=\frac{w(x_{0})}{2n+2\alpha}\,\bigg{[}\mathcal{P}_{n-1}(x_{0})-\mathcal{P}_{n+1}(x_{0})\bigg{]}.$

For the exception $n=0$, the integral to solve is
$a_{0}=\int_{x_{0}}^{1}w(x)\,\mathrm{d}x$, which is

$\displaystyle
a_{0}=\int_{x_{0}}^{1}\sqrt{1-x^{2}}^{\mathrm{D}-3}\,\mathrm{d}x=\begin{cases}\arccos(x_{0})=\upvarphi_{0},&\mathrm{D}=2\\\
1-x_{0},&\mathrm{D}=3\\\
F(\frac{1}{2},\frac{2-\mathrm{D}}{2};\frac{3}{2},1)-x_{0}\,F(\frac{1}{2},\frac{2-\mathrm{D}}{2};\frac{3}{2},x_{0}^{2}),&\text{any
$\mathrm{D}$}\end{cases}$ (141)

using Gauß’ hypergeometric function $F(a,b;c,z)$.

###  F.6 Inphase/higher-order cardioid pattern

A higher-order cardioid pattern has an $\mathrm{N}$-fold zero at $x=-1$ and is
defined as

$\displaystyle g_{\mathrm{N}}$
$\displaystyle={\textstyle\frac{1}{2^{\mathrm{N}}}}\,(1+x)^{\mathrm{N}},$
(142)

and we want to find the coefficients $a_{n}$ for its expansion

$\displaystyle g_{\mathrm{N}}$
$\displaystyle=\sum_{n=0}^{\infty}\frac{a_{n}}{N_{n}^{2}}\,\mathcal{P}_{n}(x).$
(143)

Deriving
$\frac{\mathrm{d}}{\mathrm{d}x}(1+x)^{\mathrm{N}+1}=(\mathrm{N}+1)(1+x)^{\mathrm{N}}$,
and then integrating it from its zero at $x=-1$ to $x$, we se that
$g_{\mathrm{N}}$ can be obtained from the integral equation
$(1+x)\,y=(\mathrm{N}+1)\int_{-1}^{x}y\;\mathrm{d}x$. We use the integration
and multiplication recurrences:
$2(n+\alpha)\int_{-1}^{x}\mathcal{P}_{n}\mathrm{d}x=\frac{n+2\alpha}{n+1}\mathcal{P}_{n+1}-\frac{n}{n+2\alpha-1}\mathcal{P}_{n-1}$,
and
$2(n+\alpha)(1+x)\,\mathcal{P}_{n}=(n+2\alpha)\mathcal{P}_{n+1}+2(n+\alpha)\mathcal{P}_{n}+n\mathcal{P}_{n-1}$
to write $(\mathrm{N}+1)\int y\,\mathrm{d}x-(1+x)\,y=0$:

$\displaystyle\sum_{n=0}^{\mathrm{N}}{\textstyle\frac{a_{n}}{N_{n}^{2}2(n+\alpha)}}\,[{\textstyle\frac{(\mathrm{N}+1)(n+2\alpha)-(n+1)(n+2\alpha)}{n+1}}\mathcal{P}_{n+1}-2(n+\alpha)\mathcal{P}_{n}-{\textstyle\frac{(\mathrm{N}+1)n+(n+2\alpha-1)n}{n+2\alpha-1}}\mathcal{P}_{n-1}]$
$\displaystyle=0,$
$\displaystyle\sum_{n=0}^{\mathrm{N}}{\textstyle\frac{a_{n}}{N_{n}^{2}}}\,[{\textstyle\frac{(\mathrm{N}-n)(n+2\alpha)}{2(n+\alpha)(n+1)}}\mathcal{P}_{n+1}-\mathcal{P}_{n}-{\textstyle\frac{(\mathrm{N}+n+2\alpha)n}{2(n+\alpha)(n+2\alpha-1)}}\mathcal{P}_{n-1}]$
$\displaystyle=0,$

and by applying orthogonality when integrating over
$\int_{-1}^{1}\mathcal{P}_{n^{\prime}}\mathcal{P}_{n}\mathrm{d}x=N_{n}^{2}\delta_{nn^{\prime}}$
we find after using
$N_{n}^{2}N_{n-1}^{-2}=\frac{(n+\alpha-1)n}{(n+\alpha)(n+2\alpha-1)}$ and
$N_{n}^{2}N_{n+1}^{-2}=\frac{(n+\alpha+1)(n+2\alpha)}{(n+\alpha)(n+1)}$ the
difference equation

$\displaystyle{\textstyle\frac{(\mathrm{N}-n+1)(n+2\alpha-1)}{2(n+\alpha-1)n}}{\textstyle\frac{N_{n}^{2}\,a_{n-1}}{N_{n-1}^{2}}}-{\textstyle\frac{N_{n}^{2}\,a_{n}}{N_{n}^{2}}}-{\textstyle\frac{(\mathrm{N}+n+2\alpha+1)(n+1)}{2(n+\alpha+1)(n+2\alpha)}}{\textstyle\frac{N_{n}^{2}\,a_{n+1}}{N_{n+1}^{2}}}$
$\displaystyle=0,$
$\displaystyle\textstyle(\mathrm{N}-n+1)\,a_{n-1}-2(n+\alpha)\,a_{n}-(\mathrm{N}+n+2\alpha+1)\,a_{n+1}$
$\displaystyle=0,$

and we can guess that $a_{n}$ must contain $(\mathrm{N}-n)^{-1}$ to cancel the
factor $(\mathrm{N}-n+1)$ with an index shift to $a_{n-1}$ and
$(\mathrm{N}+n+2\alpha)^{-1}$ to cancel the factor $(\mathrm{N}+n+2\alpha+1)$
with an index shift to $a_{n+1}$. We guess a decaying solution of factorials
ending at these numbers as denominator, we insert and divide by the solution
to test

$\displaystyle 1/a_{n}$ $\displaystyle=(\mathrm{N}-n)!(\mathrm{N}+n+2\alpha)!$
$\displaystyle(\mathrm{N}+n+2\alpha)-2(n+\alpha)-(\mathrm{N}-n)$
$\displaystyle=0,\qquad\blacksquare$

We express the solution normalized to its first weight $a_{0}=1$ and with the
dimensionality constant written as $2\alpha=\mathrm{D}-2$ as

$\displaystyle
a_{n}=\frac{\mathrm{N}!(\mathrm{N}+\mathrm{D}-2)!}{(\mathrm{N}-n)!(\mathrm{N}+n+\mathrm{D}-2)!}.$
(144)

The $\mathrm{N}$-fold zero at $x=-1$ can be modulated to $x=1$ by sign
alteration for the odd polynomials, i.e.
$a_{n}=(-1)^{n}\frac{\mathrm{N}!(\mathrm{N}+\mathrm{D}-2)!}{(\mathrm{N}-n)!(\mathrm{N}+n+\mathrm{D}-2)!}$.

###  F.7 Max-flat pattern: generalized higher-order cardioid

This was inspired by the maxflat FIR filter firflatLp(N,K,M) in Schüßler’s
DSP2 book [Sch10, p.92]. A polynomial with $\mathrm{L}$-fold zero at $x=1$ and
$\mathrm{M}$-fold zero at $x=-1$ integrated yields a pattern that is zero with
$\mathrm{L}$-degree flatness at $x=-1$ and a positive constant with
$\mathrm{M}$-degree flatness at $x=1$, with
$\mathrm{N}=\mathrm{L}+\mathrm{M}+1$, a Butterworth design

$\displaystyle g_{\mathrm{N}}$
$\displaystyle=\frac{(\mathrm{L+M}+1)!}{2^{\mathrm{L+M}+1}\,\mathrm{L}!\,\mathrm{M}!}\int_{-1}^{x}(1-x)^{\mathrm{L}}(1+x)^{\mathrm{M}}\,\mathrm{d}x.$
(145)

The factor provides $g_{\mathrm{N}}(1)=1$ and is obtained by applying partial
integration recursively,

$\displaystyle(m+1)\int_{-1}^{1}(1-x)^{l}(1+x)^{m}\,\mathrm{d}x$
$\displaystyle=(1-x)^{l}(1+x)^{m+1}\big{|}_{-1}^{1}+l\,\int_{-1}^{1}(1-x)^{l-1}(1+x)^{m+1}\,\mathrm{d}x$
$\displaystyle(m+1)\,g_{lm}$ $\displaystyle=g_{l,m+1}^{\prime}+l\,g_{l-1,m+1}$
$\displaystyle(m+1)(m+2)\,g_{lm}$
$\displaystyle=(m+2)g_{l,m+1}^{\prime}+l\,g_{l-1,m+2}^{\prime}+l\,(l-1)\,g_{l-2,m+2}$
$\displaystyle\frac{(m+l+1)!}{m!}g_{lm}$
$\displaystyle=\underbrace{\sum_{k=0}^{l}\frac{(m+l+1)!\,l!}{(m+k+1)!(l-k)!}\,g_{l-k,m+k+1}^{\prime}}_{g_{l-k,m+k+1}^{\prime}|_{1}=0\forall
k\neq l,\quad g_{l-k,m+k+1}^{\prime}|_{-1}=0\forall k}+0\,g_{-1,m+l+1}$
$\displaystyle\frac{(m+l+1)!}{m!}\,g_{lm}\big{|}_{1}$
$\displaystyle=l!\,g_{0,m+l+1}^{\prime}\big{|}_{1}=l!\,2^{m+l+1}$
$\displaystyle g_{lm}\big{|}_{1}$
$\displaystyle=\frac{2^{l+m+1}\,l!\,m!}{(l+m+1)!}.$

We can call $g_{\mathrm{N}}=y$ and derive twice to get
$(1-x^{2})\;y^{\prime\prime}=\mathrm{M}(1-x)\;y^{\prime}-\mathrm{L}(1+x)\;y^{\prime}$,
i.e.
$(1-x^{2})\;y^{\prime\prime}+(\mathrm{M}+\mathrm{L})\,x\,y^{\prime}-(\mathrm{M}-\mathrm{L})\,y^{\prime}=0$.
We can insert
$y=\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\mathcal{P}_{n}$ and get:

$\displaystyle(1-x^{2})\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\mathcal{P}_{n}^{\prime\prime}+(\mathrm{M}+\mathrm{L})\,x\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\mathcal{P}_{n}^{\prime}+(\mathrm{M}-\mathrm{L})\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\mathcal{P}_{n}^{\prime}$
$\displaystyle=0,$ (146)

and make use of
$(1-x^{2})\mathcal{P}_{n}^{\prime\prime}=(1+2\alpha)\,x\,\mathcal{P}_{n}^{\prime}-n(n+2\alpha)\,\mathcal{P}_{n}$,

$\displaystyle\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\,\left[(\mathrm{M}+\mathrm{L}+2\alpha+1)x\mathcal{P}_{n}^{\prime}-n(n+2\alpha)\,\mathcal{P}_{n}+(\mathrm{M}-\mathrm{L})\,\mathcal{P}_{n}^{\prime}\right]=0,$
(147)

We integrate to get

$\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\,\bigg{\\{}(\mathrm{M}+\mathrm{L}+2\alpha+1)\int
x\,\mathcal{P}_{n}^{\prime}\,\mathrm{d}x-n(n+2\alpha)\,\int\mathcal{P}_{n}\,\mathrm{d}x+(\mathrm{M}-\mathrm{L})\,\mathcal{P}_{n}\bigg{\\}}=C,$
(148)

and make use of
$2(n+\alpha)\int_{-1}^{x}x\,\mathcal{P}_{n}^{\prime}\,\mathrm{d}x=n(n+2\alpha)[\frac{1}{n+1}\mathcal{P}_{n+1}+\frac{1}{n+2\alpha-1}\mathcal{P}_{n-1}]$
$2(n+\alpha)\int_{-1}^{x}\mathcal{P}_{n}\mathrm{d}x=\frac{(n+2\alpha)}{n+1}\mathcal{P}_{n+1}-\frac{n}{n+2\alpha-1}\mathcal{P}_{n-1}$

$\textstyle\sum_{n=0}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\,\bigg{\\{}\frac{n(n+2\alpha)(\mathrm{M}+\mathrm{L}-n+1)}{(n+1)2(n+\alpha)}\,\mathcal{P}_{n+1}+\frac{n(n+2\alpha)(\mathrm{M}+\mathrm{L}+n+2\alpha+1)}{(n+2\alpha-1)2(n+\alpha)\,}\,\mathcal{P}_{n-1}+(\mathrm{M}-\mathrm{L})\,\mathcal{P}_{n}\bigg{\\}}=C.$
(149)

With the orthogonality when integrating over $\mathcal{P}_{n}$

$\textstyle\frac{N_{n}^{2}a_{n-1}}{N_{n-1}^{2}}\,\frac{(n-1)(n+2\alpha-1)(\mathrm{M}+\mathrm{L}-n+2)}{2n(n+\alpha-1)}+\frac{N_{n}^{2}a_{n+1}}{N_{n+1}^{2}}\,\frac{(n+1)(n+2\alpha+1)(\mathrm{M}+\mathrm{L}+n+2\alpha+2)}{2(n+2\alpha)(n+\alpha+1)}+a_{n}\,(\mathrm{M}-\mathrm{L})=E\,\delta_{n}.$
(150)

$N_{n}^{2}N_{n-1}^{-2}=\frac{2(n+\alpha-1)n}{2(n+\alpha)(n+2\alpha-1)}$ and
$N_{n}^{2}N_{n+1}^{-2}=\frac{2(n+\alpha+1)(n+2\alpha)}{2(n+\alpha)(n+1)}$
yield the difference equation

$\textstyle
a_{n-1}\,(n-1)(\mathrm{M}+\mathrm{L}-n+2)+a_{n+1}\,(n+2\alpha+1)(\mathrm{M}+\mathrm{L}+n+2\alpha+2)\\\
+a_{n}\,(\mathrm{M}-\mathrm{L})2(n+\alpha)=E\,\delta_{n}.$ (151)

As the resulting sequence $a_{n}$ should be limited to the degree
$\mathrm{N}=\mathrm{M}+\mathrm{L}+1$, we can abbreviate the recurrence by
$\mathrm{N}$ and $\Delta\mathrm{N}$ as

$\displaystyle\textstyle(\mathrm{N}-n+1)(n-1)\,a_{n-1}+2\Delta\mathrm{N}(n+\alpha)\,a_{n}+(\mathrm{N}+n+2\alpha+1)(n+2\alpha+1)\,a_{n+1}=E\,\delta_{n}.$

We can focus on using the recurrence for $n>0$ so that we can omit
$E\,\delta_{n}$ and write

$\displaystyle
a_{n+1}=-\frac{(\mathrm{N}-n+1)(n-1)\,a_{n-1}+2\Delta\mathrm{N}(n+\alpha)\,a_{n}}{(\mathrm{N}+n+2\alpha+1)(n+2\alpha+1)},\quad
n\geq 1.$ (152)

It is useful to start from assuming $a_{1}=1$ and obtaining $a_{n>1}$ from the
above iteration with $n=1,\dots,\mathrm{N}-1$. The missing value for $a_{0}$
can be obtained as the offset at $x=-1$, where $g_{\mathrm{N}}(-1)=0$ should
have a zero, and where $\mathcal{P}_{n}(-1)=(-1)^{n}$ consequently yields

$\displaystyle a_{0}$
$\displaystyle=-N_{0}^{2}\sum_{n=1}^{\mathrm{N}}\frac{a_{n}}{N_{n}^{2}}\mathcal{P}_{n}(-1)=-\sum_{n=1}^{\mathrm{N}}\frac{(-1)^{n}(2n+\mathrm{D}-2)\;(n+\mathrm{D}-3)!}{n!}\,a_{n}.$
(153)

Scaling can be fixed to unity by dividing all the $a_{n}\leftarrow a_{n}/b$ by
$b=g_{\mathrm{N}}(1)-g_{\mathrm{N}}(-1)$,

$\displaystyle b$
$\displaystyle=\sum_{n=1}^{\mathrm{N}}\frac{[1-(-1)^{n}](2n+\mathrm{D}-2)\;(n+\mathrm{D}-3)!}{n!}\,a_{n}$
$\displaystyle=\sum_{k=0}^{\lfloor\frac{\mathrm{N}-1}{2}\rfloor}\frac{(4k+\mathrm{D})\;(2k+\mathrm{D}-2)!}{(2k+1)!}\,a_{2k+1}.$
(154)

For a non-recurrent, closed-form solution as for the non-generalized higher-
order cardioid, only a simplified setting $\mathrm{M}=\mathrm{L}$ helps so
that the $a_{n}$ term vanishes. We get

$\displaystyle\textstyle(n-2)(\mathrm{N}-n+2)\,a_{n-2}+(n+2\alpha)(\mathrm{N}+n+2\alpha)\,a_{n}=0,$
(155) $\displaystyle\textstyle
a_{n}=\frac{(-1)^{\frac{n}{2}}(n-2)!!(2\alpha)!!(\mathrm{N}+2\alpha)!!\mathrm{N}!!}{(n+2\alpha)!!(\mathrm{N}+n+2\alpha)!!(\mathrm{N}-n)!!}$
(156)
$\displaystyle\textstyle\frac{(n-2)(\mathrm{N}-n+2)}{(n+2\alpha)(\mathrm{N}+n+2\alpha)}\frac{(n-4)!!}{(n+2\alpha-2)!!(\mathrm{N}+n+2\alpha-2)!!(\mathrm{N}-n+2)!!}-\frac{(n-2)!!}{(n+2\alpha)!!(\mathrm{N}+n+2\alpha)!!(\mathrm{N}-n)!!}=0\quad\blacksquare$
(157)

##  Appendix G MATLAB Code examples

###  G.1 gegenbauer_me(x,N,D)

Ultraspherical polynomials normalized to $\mathcal{P}_{n}(1)=1$ as given by
the recurrence of eq. (6).

⬇

1function C=gegenbauer_me(x,N,D)

2C=zeros(length(x),N+1);

3C(:,1)=1;

4if N>0

5 C(:,2)=x(:);

6end

7for n=2:N

8 C(:,n+1)=(2*n-4+D)/(n-3+D)*x(:).*C(:,n)-…

9 (n-1)/(n-3+D)*C(:,n-1);

10end

###  G.2 chebyshev(x,N)

Chebyshev polynomials $T_{n}(x)$ for $\mathrm{D}=2$ as given by the recurrence
of eq. (35).

⬇

1function T=chebyshev(x,N)

2x=x(:);

3T=zeros(length(x),N+1);

4T(:,1)=1;

5if N>0

6 T(:,2)=x;

7end

8for n=2:N

9 T(:,n+1)=2*x.*T(:,n)-T(:,n-1);

10end

###  G.3 legendre_u(x,N)

Legendre polynomials $P_{n}$ as given by the recurrence of eq. (36).

⬇

1function P=legendre_u(x,N)

2x=x(:);

3P=zeros(length(x),N+1);

4P(:,1)=1;

5if N>0

6 P(:,2)=x;

7end

8for n=1:N-1

9 P(:,(n+2))=((2*n+1)*x.*P(:,n+1)-n*P(:,n))/(n+1);

10end

###  G.4 maxre(N,D)

max-rE as defined in the optimization problem of eq. (135), and solved by
Newton method, in order to get the maximum $r_{\mathrm{E}}$ used to evaluate
the max-rE weights $a_{n}=\mathcal{P}_{n}(r_{\mathrm{E}})$.

⬇

1function [a,re]=maxre(N,D)

2re=cos(90*pi/180/(N+1)); % guessing by known D=2 max-rE

3dx0=1/(N*10); % finite difference step

4dx=dx0;

5iterate=1;

6it=0;

7while iterate % solve for P_{N+1}(rE)=0 with rE->max

8 reold=re; % store previous rE estimate

9 p0=gegenbauer_me(re,N+1,D);

10 p0=p0(end); % store previous P_{N+1}(rE)

11 pp=gegenbauer_me(re+dx,N+1,D); % estimate change+

12 pp=pp(end);

13 pm=gegenbauer_me(re-dx,N+1,D); % estimate change-

14 pm=pm(end);

15 dpdx=(pp-pm)/(2*dx); % central derivative estim.

16 re=re-p0/dpdx; % Newton-Raphson improvement step

17 dx=abs(p0)*dx0;

18 if (abs(reold-re)<1e-9)||(abs(dx)<1e-11)||(abs(p0)<1e-9)

19 iterate=0; % converged/good enough

20 end

21 it=it+1;

22end

23a=gegenbauer_me(re,N,D)’;

###  G.5 one_over_gegenbauer_norm_square(N,D)

Reciprocal of normalization squared $\frac{N_{0}^{2}}{N_{n}^{2}}$ as given in
eq. (109)

⬇

1function nn=one_over_gegenbauer_norm_square(N,D);

2nn=zeros(N+1,1);

3nn(1)=1;

4a=(D-2)/2;

5for n=1:N

6 if (n==1)&&(a==0)

7 aaa=1; % dummy a->0 that cancels

8 nn(n+1)=(a+n)*(2*aaa+0)/(n*(aaa+0))*nn(n);

9 else

10 nn(n+1)=(a+n)*(2*a+n-1)/(n*(a+n-1))*nn(n);

11 end

12end

###  G.6 supercard(N)

Supercardioid as defined in the optimization problem of eq. (24) and
subsection F.1.

⬇

1function a=supercard(N,D)

2G=gram_matrix_supercard(N,D);

3[Q,V]=eig(G);

4a=Q(:,1);

5a=a/a(1);

⬇

1function G=gram_matrix_supercard(N,D)

2G=zeros(N+1,N+1);

3a=(D-2)/2;

4n=0:N;

5N0=sqrt(pi)*gamma((D-1)/2)/gamma((D)/2);

6G(1:N+2:end)=one_over_gegenbauer_norm_square(N,D)/2;

7for k=0:N/2

8 for l=0:(N-1)/2

9 terml=(-1)^(l+1)*prod(a+(1:2*l+1))/…

10 (factorial(l)*prod(a+1+(l:2*l-1)));

11 termk=(-1)^(k)*prod(a+(1:2*k))/…

12 (factorial(k)*prod(a+(k:2*k-1)));

13 termlk=2/((2*l+1)*(2*l+2*a+1)-(2*k)*(2*k+2*a));

14 g=terml*termk*termlk/N0;

15 G(2*l+2,2*k+1)=g;

16 G(2*k+1,2*l+2)=g;

17 end

18end

###  G.7 inphase(N)

Higher-order cardioid or inphase weights $a_{n}$ as defined in eq. (28) or
derived around (144).

⬇

1function a=inphase(N,D)

2n=0:N;

3a=factorial(N)./factorial(N-n);

4b=factorial(N+D-2)./factorial(N+n+D-2);

5a=a.*b;

6a=a(:);

###  G.8 inphase_maxflat(N,L,D)

Inphase weights $a_{n}$ with $\mathrm{L}$ flatness degrees on axis, as defined
in eqs. (152) (153) (154).

⬇

1function a=inphase_maxflat(N,L,D)

2a=zeros(N+1,1);

3M=N-L-1;

4dN=L-M;

5a(2)=1;

6al=(D-2)/2;

7for n=1:N-1

8 a(n+2)=-2*dN*(n+al)*a(n+1);

9 if n>1

10 a(n+2)=a(n+2) - …

11 (N-n+1)*(n-1)*a(n);

12 end

13 a(n+2)=a(n+2)/((N+n+2*al+1)*(n+2*al+1));

14end

15nn=one_over_gegenbauer_norm_square(N,D);

16g1=a’*nn;

17nn(2:2:end)=-nn(2:2:end);

18g2=a’*nn;

19a(1)=-g2;

20a=a/(g1-g2);

###  G.9 spherecap(N,D)

Spherical cap weights $a_{n}$ for $\mathrm{D}=2,3$ as defined in eq. (29) or
derived in eqs. (140), (141), which would also contain a more arbitrary
definition for other values of $\mathrm{D}$.

⬇

1function a=spherecap(N,D,x0)

2if D==2

3 a0=acos(x0);

4elseif D==3

5 a0=1-x0;

6else

7 error(”D that is not 2 or 3 is not implemented”);

8end

9a=(D-2)/2;

10w0=(1-x0.^2).^(a-1/2);

11C0=gegenbauer_me(x0,N+1,D);

12n=1:N;

13a=w0./(2*n+2*a).*(C0(1:N)-C0(3:N+2));

14a=[a0; a(:)];

15return

###  G.10 eval_PEQrVrEFBR(a,D)

This function evaluates the level metrics $P$ and $E$, the directivity factor
$Q$, the vector metrics $r_{\mathrm{V}}$ and $r_{\mathrm{V}}$ and the front-
to-back ratio $FBR$ as given in eqs. (15)(16)(17)(18)(19)(20) and derived in
eqs. (118)(119)(120)(122)(124) of Appendix F.

⬇

1function PEQrVrEFBR=eval_PEQrVrEFBR(a,D)

2a=a(:);

3N=length(a)-1;

4n=0:N;

5Nn=one_over_gegenbauer_norm_square(N,D)’;

6SD1=2*sqrt(pi)^D/gamma(D/2); % unit sphere

7% analytic estimates of P,E,Q,rV,rE

8P=a(1)*SD1;

9E=SD1*Nn*(a(:).^2);

10Q=SD1*(Nn*a(:))^2/E;

11rV=a(2)/a(1);