HEALPix conventions

A bandlimited function $f$ on the sphere can be expanded in spherical harmonics, $Y_{\ell m}$, as

$$f({ \gamma}) = \sum_{\ell =0}^{\ell_{\mathrm{max}}}\sum_{m}a_{\el... ... \gamma})&\myequal&\sum_{\ell =0}^{\lmax}\sum_{m}a_{\ell m}Y_{\ell m}(\gamma),$$ (4)

where ${{\gamma}}$ denotes a unit vector pointing at polar angle $\theta\in[0,\pi]$ and azimuth $\phi\in[0,2\pi)$. Here we have assumed that there is insignificant signal power in modes with $\ell>\ell_{\mathrm{max}}$ and introduce the notation that all sums over $m$ run from $-\ell_{\mathrm{max}}$ to $\ell_{\mathrm{max}}$ but all quantities with index ${\ell m}$ vanish for $m>\ell$. Our conventions for $Y_{\ell m}$ are defined in subsection A.4 below. The spherical harmomics coefficients are then

$$a_{\ell m} \equiv \int d\gamma Y^\ast_{\ell m}(\gamma) f(\gamma),$$ (5)

where, the integral is done over the whole sphere, and the superscript star denotes complex conjugation.

Pixelating $f({\gamma})$ corresponds to sampling it at $N_{\mathrm{pix}}$ locations $\gamma_{p}$, $p\in[0,N_{\mathrm{pix}}-1]$. The sample function values $f_p$ can then be used to estimate $a_{\ell m}$. A straightforward estimator is

$$\hat{a}_{\ell m} = \frac{4\pi}{N_{\mathrm{pix}}}\sum_{p=0}^{N_{\mathrm{pix}}-1} Y^\ast_{\ell m}(\gamma_p) f(\gamma_p),$$ (6)

where an equal weight was assumed for each pixel. This zeroth order estimator, as well as higher order estimators, are implemented in various Fortran90, C++, IDL, and python facilities included in the package, such as anafast or anafast_cxx.

Angular power spectrum conventions

These $\hat{a}_{\ell m}$ can be used to compute estimates of the angular power spectrum $\hat{C}_\ell$ as

$$\hat{C}_\ell =\frac{1}{2\ell +1}\sum_{m} \vert\hat{a}_{\ell m}\vert^2.$$ (7)

Equations (6) and (7) above do not consider the impact of a pixel masking or weighting $f(\gamma_p) \longrightarrow f(\gamma_p) w(\gamma_p)$ on the power spectrum estimation of $f$, which is described in Wandelt, Hivon & Górski (2001) and addressed in Hivon et al. (2002), Chon et al. (2004), Tristram et al. (2005), Rocha et al. (2009) and Planck 2015-XI (2015) among others.

The HEALPix package contains the Fortran90 facility synfast, which takes as input a power spectrum $C_\ell$ and generates a realisation of $f(\gamma_p)$ on the HEALPix grid. The convention for power spectrum input into synfast is straightforward: each $C_\ell$ is just the expected variance of the $a_{\ell m}$ at that $\ell$.

Example: The spherical harmonic coefficient $a_{00}$ is the integral of the $f(\gamma)/\sqrt{4 \pi}$ over the sphere. To obtain realisations of functions which have $a_{00}$ distributed as a Gaussian with zero mean and variance 1, set $C_0$ to 1. The value of the synthesised function at each pixel will be Gaussian distributed with mean zero and variance $1/(4\pi)$. As required, the integral of $f(\gamma)$ over the full $4\pi$ solid angle of the sphere has zero mean and variance $4\pi$.

Note that this definition implies the standard result that the total power at the angular wavenumber $\ell$ is $(2\ell+1)C_\ell$, because there are $2\ell+1$ modes for each $\ell$.

This defines unambiguously how the $C_\ell$ have to be defined given the units of the physical quantity $f$. In cosmic microwave background research, popular choices for simulated maps are

• $\Delta T/T$, a dimensionless quantity measuring relative fluctuations about the average CMB temperature.
• The absolute quantity $\Delta T$ in $\mu K$ or $K$.

HEALPix and Boltzmann codes

CMBFAST

A widely used solver of the Boltzmann equations for the computation of theoretical predictions of the spectrum of CMB anisotropy used to be CMBFAST (https://lambda.gsfc.nasa.gov/toolbox/tb_cmbfast_ov.cfm).

CMBFAST made its outputs in ASCII files, which instead of $C_{X,\ell}$ contain quantities defined as

$$D_{X,\ell} = \frac{\ell(\ell+1)}{(2\pi)T_{CMB}^2}C_{X,\ell},$$ (8)

where $T_{CMB}=2.726K$ is the temperature of the CMB today and $X$ stands for T, E, B or C (see § A.3).

The version 4.0 of CMBFAST also created a FITS file containing the power spectra $C_{X,\ell}$, designed for interface with HEALPix. The spectra for polarization were renormalized to match the normalization used in HEALPix 1.1, which was different from the one used by CMBFAST and by HEALPix 1.2 (see § A.3.2 for details).

A later version of CMBFAST (4.2, released in Feb. 2003) generated FITS files containing $C_{X,\ell}$, with the same convention for polarization as the one used internally. It therefore matches the convention adopted by HEALPix in its version 1.2.

For backward compatibility, we provide an IDL code (convert_oldhpx2cmbfast) to change the normalization of existing FITS files created with CMBFAST 4.0. When created with the correct normalization (with CMBFAST 4.2) or set to the correct normalization (using convert_oldhpx2cmbfast), the FITS file will include a specific keyword (POLNORM = CMBFAST) in their header to identify them. The map simulation code synfast will issue a warning if the input power spectrum file does not contain the keyword POLNORM, but no attempt will be made to renormalize the power spectrum. If the keyword is present, it will be inherited by the simulated map.

CAMB and CLASS

Newer and actively maintained Boltzmann codes currently include camb and class:

• camb (https://camb.info) is written in Fortran 90 with a python wrapper, and can optionally output into FITS files the $C_X(\ell)$ power spectra in [K]$^2$ in a format directly usable by HEALPix;
  
• class (http://class-code.net) is written in C and C++, and only outputs $\frac{\ell(\ell+1)}{2\pi}C_X(\ell)$ in plain text files (optionally in [$\mu$K]$^2$ and in a order of columns for polarized spectra matching the one of camb).

Both codes are parallelized for faster computations and provide fine control of the output accuracy.

Polarisation convention

Figure 4: Orthographic projection of a fake full sky for temperature (color coded) and polarization (represented by the rods). All the input Spherical Harmonics coefficients are set to 0, except for $a_{21}^{TEMP}=\ -\ a_{2-1}^{TEMP}=1$ and $a_{21}^{GRAD}=\ -\ a_{2-1}^{GRAD}=1$

Internal convention

Starting with version 1.20 (released in Feb 2003),HEALPix uses the same conventions as CMBFAST for the sign and normalization of the polarization power spectra, as exposed below (adapted from Zaldarriaga (1998)). How this relates to what was used in previous releases is exposed in A.3.2.

The CMB radiation field is described by a $2\, \times \, 2$ intensity tensor $I_{ij}$ (Chandrasekhar, 1960). The Stokes parameters $Q$ and $U$ are defined as $Q=(I_{11}-I_{22})/4$ and $U=I_{12}/2$, while the temperature anisotropy is given by $T=(I_{11}+I_{22})/4$. The fourth Stokes parameter $V$ that describes circular polarization is not necessary in standard cosmological models because it cannot be generated through the process of Thomson scattering. While the temperature is a scalar quantity $Q$ and $U$ are not. They depend on the direction of observation $\textbf{n}$ and on the two axis $(\textbf{e}_{1}, \textbf{e}_{2})$ perpendicular to $\textbf{n}$ used to define them. If for a given $\textbf{n}$ the axes $(\textbf{e}_{1}, \textbf{e}_{2})$ are rotated by an angle $\psi$ such that ${\textbf{e}_{1}}^{\prime}=\cos \psi \ {\textbf{e}_{1}}+\sin\psi \ {\textbf{e}_{2}}$ and ${\textbf{e}_{2}}^{\prime}=-\sin \psi \ {\textbf{e}_{1}}+\cos\psi \ {\textbf{e}_{2}}$ the Stokes parameters change as

$$Q^{\prime} = \cos 2\psi \ Q + \sin 2\psi \ U$$

$$U^{\prime} = -\sin 2\psi \ Q + \cos 2\psi \ U$$ (9)

To analyze the CMB temperature on the sky, it is natural to expand it in spherical harmonics. These are not appropriate for polarization, because the two combinations $Q\pm iU$ are quantities of spin $\pm 2$ (Goldberg, 1967). They should be expanded in spin-weighted harmonics $\, _{\pm2}Y_l^m$ (Seljak & Zaldarriaga, 1997; Zaldarriaga & Seljak, 1997),

$$T(\textbf{n}) = \sum_{lm} a_{T,lm} Y_{lm}(\textbf{n})$$

$$(Q+iU)(\textbf{n}) = \sum_{lm} a_{2,lm}\;_2Y_{lm}(\textbf{n})$$

$$(Q-iU)(\textbf{n}) = \sum_{lm} a_{-2,lm}\;_{-2}Y_{lm}(\textbf{n}).$$ (10)

To perform this expansion, $Q$ and $U$ in equation (10) are measured relative to $(\textbf{e}_{1}, \textbf{e}_{2})=(\textbf{e}_\theta , \textbf{e}_\phi )$, the unit vectors of the spherical coordinate system. Where $\textbf{e}_\theta$ is tangent to the local meridian and directed from North to South, and $\textbf{e}_\phi$ is tangent to the local parallel, and directed from West to East. The coefficients $_{\pm 2}a_{lm}$ are observable on the sky and their power spectra can be predicted for different cosmological models. Instead of $_{\pm 2}a_{lm}$ it is convenient to use their linear combinations

$$a_{E,lm} = -(a_{2,lm}+a_{-2,lm})/2$$

$$a_{B,lm} = -(a_{2,lm}-a_{-2,lm})/2i,$$ (11)

which transform differently under parity. Four power spectra are needed to characterize fluctuations in a gaussian theory, the autocorrelation between $T$, $E$ and $B$ and the cross correlation of $E$ and $T$. Because of parity considerations the cross-correlations between $B$ and the other quantities vanish and one is left with

$$\langle a_{X,lm}^{*} a_{X,lm^\prime}\rangle = \delta_{m,m^\prime}C_{Xl},$$

$$\langle a_{T,lm}^{*}a_{E,lm}\rangle = \delta_{m,m^\prime}C_{Cl},$$ (12)

where $X$ stands for $T$, $E$ or $B$, $\langle\cdots \rangle$ means ensemble average and $\delta_{i,j}$ is the Kronecker delta.

We can rewrite equation (10) as

$$T(\textbf{n}) = \sum_{lm} a_{T,lm} Y_{lm}(\textbf{n})$$

$$Q(\textbf{n}) = -\sum_{lm} a_{E,lm} X_{1,lm} +i a_{B,lm}X_{2,lm}$$

$$U(\textbf{n}) = -\sum_{lm} a_{B,lm} X_{1,lm}-i a_{E,lm} X_{2,lm}$$ (13)

where we have introduced $X_{1,lm}(\textbf{n})=(\;_2Y_{lm}+\;_{-2}Y_{lm})/2$ and $X_{2,lm}(\textbf{n})=(\;_2Y_{lm}-\;_{-2}Y_{lm})/ 2$. They satisfy $Y^{*}_{lm} = (-1)^m Y_{l-m}$, $X^{*}_{1,lm}=(-1)^m X_{1,l-m}$ and $X^*_{2,lm}=(-1)^{m+1}X_{2,l-m}$ which together with $a_{T,lm}=(-1)^m a_{T,l-m}^*$, $a_{E,lm}=(-1)^m a_{E,l-m}^*$ and $a_{B,lm}=(-1)^m a_{B,l-m}^*$ make $T$, $Q$ and $U$ real.

In fact $X_{1,lm}(\textbf{n})$ and $X_{2,lm}(\textbf{n})$ have the form, ${X_{1,lm}(\textbf{n})=\sqrt{(2l+1) / 4\pi} F_{1,lm}(\theta)\ e^{im\phi}}$ and ${X_{2,lm}(\textbf{n})=\sqrt{(2l+1) / 4\pi} F_{2,lm}(\theta)\ e^{im\phi}}$, ${F_{(1,2),lm}(\theta)}$ can be calculated in terms of Legendre polynomials (Kamionkowski et al., 1997)

$$F_{1,lm}(\theta) = N_{lm} \left[ -\left({l-m^2 \over \sin^2\theta} +{1 \over 2}l(l... ...s \theta) +(l+m) {\cos \theta \over \sin^2 \theta} P_{l-1}^m(\cos\theta)\right]$$

$$F_{2,lm}(\theta) = N_{lm}{m \over \sin^2 \theta} [ -(l-1)\cos \theta P_l^m(\cos \theta)+(l+m) P_{l-1}^m(\cos\theta)],$$ (14)

where

$$N_{lm}(\theta) = 2 \sqrt{(l-2)!(l-m)! \over (l+2)!(l+m)!}.$$ (15)

Note that $F_{2,lm}(\theta)=0$ if $m=0$, as it must to make the Stokes parameters real.

The correlation functions between 2 points on the sky (noted 1 and 2) separated by an angle $\beta$ can be calculated using equations (12) and (13). However, as pointed out in Kamionkowski et al. (1997), the natural coordinate system to express the correlations is one in which $\textbf{e}_{1}$ vectors at each point are tangent to the great circle connecting these 2 points, with the $\textbf{e}_{2}$ vectors being perpendicular to the $\textbf{e}_{1}$ vectors. With this choice of reference frames, and using the addition theorem for the spin harmonics (Hu & White, 1997),

$$\sum_m \;_{s_1} Y_{lm}^*(\textbf{n}_1) \;_{s_2} Y_{lm}(\textbf{n}_2) = \sqrt{2l+1 \over 4 \pi} \;_{s_2} Y_{l-s_1}(\beta,\psi_1)e^{-is_2\psi_2}$$ (16)

we have (Kamionkowski et al., 1997)

$$\langle T_1T_2 \rangle = \sum_l {2l+1 \over 4 \pi} C_{Tl} P_l(\cos \beta)$$

$$\langle Q_{r}(1)Q_{r}(2) \rangle = \sum_l {2l+1 \over 4 \pi} [C_{El} F_{1,l2}(\beta)-C_{Bl} F_{2,l2}(\beta)]$$

$$\langle U_{r}(1)U_{r}(2) \rangle = \sum_l {2l+1 \over 4 \pi} [C_{Bl} F_{1,l2}(\beta)-C_{El} F_{2,l2}(\beta) ]$$

$$\langle T(1)Q_{r}(2) \rangle = - \sum_l {2l+1 \over 4 \pi} C_{Cl} F_{1,l0}(\beta)$$

$$\langle T(1)U_{r}(2) \rangle = 0.$$ (17)

The subscript $r$ here indicate that the Stokes parameters are measured in this particular coordinate system. We can use the transformation laws in equation (9) to write $(Q,U)$ in terms of $(Q_r,U_r)$.

Using the fact that, when $\beta \rightarrow 0$, $P_\ell(\cos\beta) \rightarrow 1$ and $P_\ell^2(\cos \beta) \rightarrow \sin^2 \beta \frac{(\ell+2)!}{8 (\ell-2)!}$, the definitions above imply that the variances of the temperature and polarization are related to the power spectra by

$$\langle TT \rangle = \sum_\ell {2\ell+1 \over 4 \pi} C_{T\ell}$$

$$\langle QQ \rangle + \langle UU\rangle = \sum_l {2\ell+1 \over 4 \pi} \left(C_{E\ell} +C_{B\ell}\right)$$

$$\langle TQ\rangle = \langle TU\rangle = 0.$$ (18)

It is also worth noting that with these conventions, the cross power $C_{C\ell}$ for scalar perturbations must be positive at low $\ell$, in order to produce at large scales a radial pattern of polarization around cold temperature spots (and a tangential pattern around hot spots) as it is expected from scalar perturbations (Crittenden et al., 1995).

Note that Eq. (13) implies that, if the Stokes parameters are rotated everywhere via

$$\left(\begin{array}{c} Q'\\ U' \end{array}\right) = \left(\begin{... ...end{array} \right) \left(\begin{array}{c} Q\\ U \end{array} \right),%\nonumber$$ (19)

then the polarized $a_{\ell m}$ coefficients are submittted to the same rotation

$$\left(\begin{array}{c} a_{E,\ell m}'\\ a_{B,\ell m}' \end{array}\... ...ft(\begin{array}{c} a_{E,\ell m}\\ a_{B,\ell m} \end{array} \right).%\nonumber$$ (20)

Finally, with these conventions, a polarization with ($Q>0,U=0$) will be along the North–South axis, and ($Q=0,U>0$) will be along a North-West to South-East axis (see Fig. 5)

Relation to previous releases

Even though it was stated otherwise in the documention, HEALPix used a different convention for the polarization in its previous releases. The tensor harmonics approach (Kamionkowski et al. (1997), hereafter KKS) was used, instead of the current spin weighted spherical harmonics. These two approaches differ by the normalisation and sign of the basis functions used, which in turns change the normalisation of the power spectra. Table 1 summarizes the relations between the CMB power spectra in the different releases. See § A.2.1 about the interface between HEALPix and CMBFAST.

Table 1: Relation between CMB power spectra conventions used in HEALPix, CMBFAST and KKS. The power spectra on the same row are equal.

Component	HEALPix $\ge$ 1.21	CMBFAST	KKS	HEALPix $\le$ 1.12
Temperature	$C_{\ell}^{\mathrm{TEMP}}$	$C_{\mathrm{T},\ell}$	$C_{\ell}^{\mathrm{T}}$	$C_{\ell}^{\mathrm{TEMP}}$
Electric or Gradient	$C_{\ell}^{\mathrm{GRAD}}$	$C_{\mathrm{E},\ell}$	$2C_{\ell}^{\mathrm{G}}$	$2C_{\ell}^{\mathrm{GRAD}}$
Magnetic or Curl	$C_{\ell}^{\mathrm{CURL}}$	$C_{\mathrm{B},\ell}$	$2C_{\ell}^{\mathrm{C}}$	$2C_{\ell}^{\mathrm{CURL}}$
Temp.-Electric cross correlation	$C_{\ell}^{\mathrm{T-GRAD}}\rule[.3cm]{0cm}{.2cm}$	$C_{\mathrm{C},\ell}$	$-\sqrt{2}$ $C_{\ell}^{\mathrm{TG}}$	$\sqrt{2}C_{\ell}^{\mathrm{T-GRAD}}$

¹ Version 1.2 (Feb 2003) or more recent of HEALPix package

² Version 1.1 or older of HEALPix package

Introducing the matrices

$$M_{\ell m} = \left( \begin{array}{cc} X_{1,\ell m} & i X_{2,\ell m} \\ -i X_{2,\ell m} & X_{1,\ell m} \end{array}\right)$$ (21)

where the basis functions $X_1$ and $X_2$ have been defined in Eqs. (14) and above, the decomposition in spherical harmonics coefficients (13) of a given map of the Stokes parameter $Q$ and $U$ can be written in the case of HEALPix 1.2 as

$$\phantom{1.2}{ \left( \begin{array}{c} Q \rule[.3cm]{0cm}{.2c... ...le[.3cm]{0cm}{.2cm}\rule[-.3cm]{0cm}{.2cm}\end{array}\right) }.$$ (22)

For KKS, with the same definition of $M$, the decomposition reads

$${ \left( \begin{array}{c} Q \rule[.3cm]{0cm}{.2cm}\rule[-.3cm... ...le[.3cm]{0cm}{.2cm}\rule[-.3cm]{0cm}{.2cm}\end{array}\right) },$$ (23)

whereas in HEALPix 1.1 it was

$$\phantom{1.1}{ \left( \begin{array}{c} Q \rule[.3cm]{0cm}{.2c... ...le[.3cm]{0cm}{.2cm}\rule[-.3cm]{0cm}{.2cm}\end{array}\right) }.$$ (24)

The difference between KKS and 1.1 was due to an error of sign on one the basis functions.

Relation with IAU convention

In a cartesian referential with axes $x$ and $y$, the Stokes parameters for linear polarisation are defined such that $+Q$ is aligned with $+x$, $-Q$ with $+y$ and $+U$ with the bisectrix of $+x$ and $+y$. Although this definition is universally accepted, some confusion may still arise from the relation of this local cartesian system to the global spherical one, as described below (Hamaker & Leahy, 2003), and as illustrated in Fig. 5.

Figure 5: Coordinate conventions for HEALPix (lhs panels) and IAU (rhs panels). The upper panels illustrate how the spherical coordinates are measured, and the lower panel how the $Q$ and $U$ Stokes parameters are identified in the tangential plan.

The polarization conventions defined by the International Astronomical Union (IAU, 1974) are summarized in Hamaker & Bregman (1996). They define at each point on the celestial sphere a cartesian referential with the $x$ and $y$ axes pointing respectively toward the North and East, and the $z$ axis along the line of sight pointing toward the observer (ie, inwards) for a right-handed system.

On the other hand, following the mathematical and CMB litterature tradition, HEALPix defines a cartesian referential with the $x$ and $y$ axes pointing respectively toward the South and East, and the $z$ axis along the line of sight pointing away from the observer (ie, outwards) for a right-handed system. The Planck CMB mission follows the same convention (Ansari et al., 2003).

The consequence of this definition discrepency is a change of sign of $U$, which, if not accounted for, jeopardizes the calculation of the Electric and Magnetic CMB polarisation power spectra.

How HEALPix deals with these discrepancies: POLCCONV keyword

The FITS keyword POLCCONV has been introduced in HEALPix 2.0 to describe the polarisation coordinate convention applied to the data contained in the file. Its value is either 'COSMO' for files following the HEALPix/CMB/Planck convention (default for sky map synthetized with HEALPix routine synfast) or 'IAU' for those following the IAU convention, as defined above. Absence of this keyword is interpreted as meaning 'COSMO' (as it is the case for WMAP maps).

Starting with HEALPix 3.40, when dealing with a polarized (full-sky or cut-sky) signal map,
– the F90 subroutine input_map in its default mode,
– the F90 facilities calling it and dealing with the $I$, $Q$ and $U$ Stokes parameters as a whole, ie anafast and smoothing,
– as well as their IDL wrappers ianafast and ismoothing,
– the IDL visualisation routines azeqview, cartview, gnomview, mollview and orthview called with Polarization=2 or 3,
– and all C++ facilities (and the input routine read_Healpix_map_from_fits)
will all
– issue an error message and crash if POLCCONV is explicitely set to a value different from 'COSMO' and 'IAU',
– issue a warning (except in C++), and swap the sign of the $U$ polarisation stored into memory if the FITS file being read contains POLCCONV='IAU',
– issue a warning (except in C++) if the keyword POLCCONV is totally absent, and then carry on with the original data,
– or work silently with the original data if POLCCONV='COSMO'.
On the other hand, and as in previous releases, routines treating or showing each of $I$, $Q$ and $U$ fields separately, such as the F90 facilities median_filter, ud_grade, or map2gif as well as their IDL counterparts median_filter, ud_grade, or mollview et al run with Polarization=0 or 1 will ignore the value of POLCCONV (copying it unchanged into their output files, when applicable) and preserve the sign of $U$.

Finally, the IDL subroutine change_polcconv.pro and the Python facility change_polcconv.py are provided to add the POLCCONV keyword or change/update its value and swap the sign of the $U$ Stokes parameter, when applicable, in an existing FITS file.

Spherical harmonic conventions

The Spherical Harmonics are defined as

$$Y_{\ell m}(\theta,\phi) = \lambda_{\ell m}(\cos\theta) e^{{i} m\phi}$$ (25)

where

$$\lambda_{\ell m}(x) = \sqrt{ \frac{2\ell+1}{4\pi} \frac{(\ell-m)!}{(\ell+m)!} } P_{\ell m}(x), \quad\textrm{for~} m\ge 0$$ (26)

$$\lambda_{\ell m} = (-1)^m \lambda_{\ell \vert m\vert}, \quad\textrm{for~} m < 0,$$

$$\lambda_{\ell m} = 0, \quad\textrm{for}\, \vert m\vert > \ell.$$

Introducing $x\equiv\cos\theta$, the associated Legendre Polynomials $P_{\ell m}$ solve the differential equation

$$(1-x^2)\frac{d^2}{dx^2}P_{\ell m} - 2x \frac{d}{dx}P_{\ell m} + \left(\ell(\ell+1) - \frac{m^2}{1-x^2}\right) P_{\ell m} = 0.$$ (27)

They are related to the ordinary Legendre Polynomials $P_\ell$ by

$$P_{\ell m} = (-1)^m (1-x^2)^{m/2} \frac{d^m}{dx^m} P_{\ell}(x),$$ (28)

which are given by the Rodrigues formula

$$P_{\ell}(x) = \frac{1}{2^\ell \ell!}\frac{d^\ell}{dx^\ell} (x^2-1)^\ell.$$ (29)

Note that our $Y_{\ell m}$ are identical to those of Edmonds (1957), even though our definition of the $P_{\ell m}$ differ from his by a factor $(-1)^m$ (a.k.a. Condon-Shortley phase).