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¹, 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:

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

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.

1. Starting from $P_0^0$, boost along the diagonal to the desired order $m$ with the first recurrence relation.
2. Use the second recurrence relation to boost from $P_m^m$ to $P_\ell^m$ where $\ell = m+1$.
3. 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

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.