The Associated Legendre Polynomials are implicitly defined as the solution to a second-order differential equation, but most practical uses require an efficient means of explicitly evaluating the functions for any degree, order, and argument. In this article I introduce the implementation used in Legendre.jl which is based on evaluating a series of recurrence relations.

## Introduction¶

In the first part of this series,
we saw that the Associated Legendre Polynomials are implicitly defined as
the solution to a second-order differential equation arising from solving
Laplace’s equation in spherical coordinates.
This doesn’t help much if you want to evaluate the functions, though, so
we require a mathematical definition which is easier to utilize than solving
a differential equation on every use.
The method used by Legendre.jl is to evaluate a few recurrence relations.
I was first introduced to this particular method in Limpanuparb &
Milthorpe^{1}, but several other well-known
sources^{2}^{,}^{3} mention or use the same underlying
principle.

## Recurrence Relations¶

Our method of calculation makes use of the following three recurrence
relations:
`\begin{align} P_{\ell+1}^{\ell+1}(x) &= -(2\ell+1)\sqrt{1-x^2} P_\ell^\ell(x) \label{eqn:std_rr_1term_lm} \\ P_{\ell+1}^\ell(x) &= (2\ell+1) x P_\ell^\ell(x) \label{eqn:std_rr_1term_l} \\ (\ell - m + 1)P_{\ell+1}^m(x) &= (2\ell+1)xP_\ell^m(x) - (\ell+m)P_{\ell-1}^m(x) \label{eqn:std_rr_2term} \end{align}`

When combined with an appropriate initial condition, the Associated Legendre
Polynomials can be computed for any degree, order, and argument.
In standard form, the required initial condition is `\(P_0^0(x)=1\)`

for all `\(x\)`

.
We also assume that `\(m \ge 0\)`

for this and all future postings for a couple
of reasons.
First, the negative-order Associated Legedre Polynomials are related to the
positive orders via the identity
`\begin{align} P_\ell^{-m}(x) = (-1)^m \frac{(\ell-m)!}{(\ell+m)!} P_\ell^m(x) \end{align}`

and therefore the negative orders can be easily computed from only the postive
ones.
Second, the CMB is a real quantity, so the harmonic coefficients for the
negative orders must be related to the coefficients of the postive orders via
a symmetry relation.
Therefore in practical use, it is unnecessary (and may in fact be wasteful) to
explicitly calculate the negative order terms.

The recurrence relations above can now be separated into two groups —
the first (Eq. `\(\ref{eqn:std_rr_1term_lm}\)`

) changes both degree
and order, while the latter two
(Eqs. `\(\ref{eqn:std_rr_1term_l}\)`

–`\(\ref{eqn:std_rr_2term}\)`

)
only change the degree.
Also note that the first two are single-term recurrence relations that require
only one prior value to evaluate, but the last equation is a two-term
recurrence relation that requires two prior values.

With a bit of introspection, an algorithm which can reach any pair
`\((\ell,m)\)`

using only these three relations can be constructed.
The method is shown schematically below in Fig. 2.1.

- Starting from
`\(P_0^0\)`

, boost along the diagonal to the desired order`\(m\)`

with the first recurrence relation. - Use the second recurrence relation to boost from
`\(P_m^m\)`

to`\(P_\ell^m\)`

where`\(\ell = m+1\)`

. - Finally, loop using the third recurrence relation until the desired
degree
`\(\ell\)`

is reached.

We now have everything necessary to actually compute the Associated Legendre Polynomials, so in the next section we’ll implement a basic calculation implementation.

## Implementing the recurrence relations¶

Let us now translate the algorithm into working Julia code.
We’ll call the function `legendre`

, taking in as arguments the degree `ℓ`

,
the order `m`

, and an argument `x`

.
The initial condition is simply the constant value `1`

, so if `ℓ == 0`

and `m == 0`

, the function returns immediately.

```
function legendre(ℓ, m, x)
P = one(x)
ℓ == 0 && m == 0 && return P
```

The next step in the outlined algorithm above is to move along the diagonal
where `\(\ell = m\)`

at each step until the desired order `\(m\)`

is reached.
Again, the function can return early if the requested degree and order
was on the diagonal.

```
# Use recurrence relation 1 to increase to the target m
y = sqrt(one(x) - x^2)
for m′ in 0:m-1
P = -(2m′ + 1) * y * P
end
ℓ == m && return P
```

Now the remaining steps involve increasing in degree until the requested
`ℓ`

is reached.
The two-term recurrence relation requires two lower degree terms, so we
need to keep track of a little bit more state at each step.
At the beginning, the on-diagonal term corresponds to `\(P_{\ell-1}^m\)`

in the 2-term recurrence (`\(\ref{eqn:std_rr_2term}\)`

), and we use
the 1-term recurrence (`\(\ref{eqn:std_rr_1term_l}\)`

) to generate the
second “prior” term `\(P_{\ell}^m\)`

(returning early if the calculation is done).

```
# Setup for 2-term recurrence, using recurrence relation 2 to set
# Pmin1
Pmin2 = P
Pmin1 = (2m + 1) * x * Pmin2
ℓ == m+1 && return Pmin1
```

Finally, the 2-term recurrence proceeds to the desired `\(\ell\)`

, with the
final iteration having the requested value of `\(P_\ell^m(x)\)`

in `P`

:

```
# Use recurrence relation 3
for ℓ′ in m+1:ℓ-1
P = ((2ℓ′ + 1) * x * Pmin1 - (ℓ′ + m) * Pmin2) / (ℓ′ - m + 1)
Pmin2, Pmin1 = Pmin1, P
end
return P
end
```

That’s all there is to a naive implementation for calculating the polynomials. We’ll see in future parts of this series that there are some issues which should be improved upon, but for now this serves as a useful and very easily understood base.

## Plots of the polynomials¶

We’ll conclude this posting by using our new function to make a couple of
summary plots.
The first one we’d like to make reproduces a figure often shown in summaries
of the polynomials — the functions on the interval `\(x \in [-1, 1]\)`

for the
first few degrees (here `\(0 \le \ell \le 4\)`

):

```
using PyPlot
x = range(-1, 1, length=1000)
for ℓ in 0:4
plot(x, legendre.(ℓ, 0, x), label="\$P_$ℓ^0\$")
end
axhline.([-1, 0, 1], color="gray", linewidth=0.5)
axvline.([-1, 0, 1], color="gray", linewidth=0.5)
legend(loc = "lower right")
title("Legendre Polynomials (\$m = 0\$)")
```

Note that for `\(m = 0\)`

, these are all truly polynomials of degree `\(\ell\)`

, and
we can qualitatively see expected characteristics of polynomials such as each
function having `\(\ell\)`

roots.

Despite their name, the Associated Legendre Polynomials are not actually
polynomials for odd `\(m\)`

.
For instance, we can plot the first 5 “polynomials” for order `\(m = 1\)`

, and
we qualitatively see asymptoting behavior as `\(x \rightarrow \pm 1\)`

.

```
for m in 1:5
plot(x, legendre.(m, 1, x), label="\$P_{$m}^1\$")
end
axhline.([-1, 0, 1], color="gray", linewidth=0.5)
axvline.([-1, 0, 1], color="gray", linewidth=0.5)
legend(loc = "lower right")
title("Associated Legendre Polynomials (\$m = 1\$)")
```

## Growth and numeric overflow¶

Another observation to note is that the `\(m = 1\)`

functions have peak amplitudes
outside of the range `\([-1, 1]\)`

.
In general, the extrema grow in maximum (absolute) amplitude for increasing
`\(\ell\)`

/ `\(m\)`

, and the numerical computation quickly overflows for only
moderate degrees and orders that are critical to CMB analysis.
Specifically, the envelope of `\(P_\ell^m(x)\)`

which bounds the local extrema for
all values of `\(x\)`

is
`\begin{align} \left| P_\ell^m(x) \right| \le \frac{Γ(\ell+m+1)}{Γ(\ell+\frac{3}{2})} \left( \frac{2}{\pi\sqrt{1 - x^2}} \right)^{1/2} \end{align}`

(see Abramowitz and Stegun^{4}), which asymptotically approaches
for large `\(\ell\)`

,
`\begin{align} \DeclareMathOperator*{\env}{env} \env_{\ell\rightarrow\infty}\left(P_\ell^m\right) &\propto \ell^{m-1/2} \end{align}`

With a bit of trial and error, we can find that `\(\ell,m \sim 150\)`

is near the
point where the computations overflow.
For instance, plotting `\(P_{151}^{150}(x)\)`

and `\(P_{152}^{150}(x)\)`

for
`\(x \in [-1, 1]\)`

:

```
plot(x, legendre.(151, 150, x) ./ 1e308, label="\$P_{151}^{150}\$")
plot(x, legendre.(152, 150, x) ./ 1e308, label="\$P_{152}^{150}\$")
ylim((-1, 1) .* floatmax(Float64) ./ 1e308)
ylabel("\$\\times 10^{308}\$")
legend()
```

The break in the orange line from `\(-0.25\)`

to `\(0.25\)`

is due to the
calculations overflowing the range of 64-bit floating point numbers.

The spherical harmonics require calculating the polynomials to degrees and
orders much greater than `\(150\)`

, so this poses a substantial problem to
resolve.
In the next part of this series, we’ll confront this issue and show how
the normalization term in the spherical harmonics can be combined with the
recurrence relations to produce a new pre-normalized function which no
longer overflows.

T. Limpanuparb and J. Milthorpe. “Associated Legendre Polynomials and Spherical Harmonics Computation for Chemistry Applications” In:

*Proceedings of the 40th Congress on Science and Technology of Thailand*(Dec 2014) arXiv: 1410.1748 ↩︎“Legendre and Related Functions” (Chapter 14) In:

*NIST Digital Library of Mathematical Functions*↩︎M. Reinecke. “Libpsht - algorithms for efficient spherical harmonic transforms" In:

*Astronomy & Astrophysics*(Feb 2011) arXiv: 1010.2084 ↩︎Eq. 8.10.7 (p336) of Abramowitz and Stegun, In:

*Handbook of Mathematical Functions*10th printing (1972) ↩︎