Legendre Polynomials

Legendre polynomials are what the angular half of a spherically symmetric problem is made of. Separate the Laplacian in spherical coordinates, drop any azimuthal ( $\phi$ ) dependence, and set $x=\cos\theta$ : the polar equation collapses to $\frac{\mathrm{d}}{\mathrm{d}x}\!\big[(1-x^2)\frac{\mathrm{d}y}{\mathrm{d}x}\big]+\lambda y=0$ . Insist that the solution stay finite at the poles $x=\pm1$ — a physical demand, since the potential there is finite — and the separation constant is forced to the ladder $\lambda=\ell(\ell+1)$ , $\ell=0,1,2,\dots$ . The regular solution is the polynomial $P_\ell(x)$ ; its companion $Q_\ell(x)$ blows up logarithmically at the endpoints and is exiled from the interior of any sphere. Three faces of the same object organize everything that follows. As a Sturm–Liouville eigenfunction on $[-1,1]$ with weight unity, $P_\ell$ is automatically orthogonal to its siblings and the family is complete — any square-integrable function has a convergent Legendre series. As a Rodrigues derivative it is a compact closed form, $P_\ell=\frac{1}{2^\ell\ell!}\frac{\mathrm{d}^\ell}{\mathrm{d}x^\ell}(x^2-1)^\ell$ , ideal for proofs. And as the generating-function coefficient of $(1-2xt+t^2)^{-1/2}$ it is no abstract machinery: that function is literally the inverse-distance kernel $1/|\mathbf r-\mathbf r'|$ of Newton and Coulomb, so every Legendre index is a multipole order and every axisymmetric potential is a Legendre series. From this hub the theory branches: restore azimuthal dependence and $P_\ell$ is promoted to the associated Legendre functions $P_\ell^m$ , which dressed with $e^{im\phi}$ become the spherical harmonics $Y_\ell^m$ ; follow the radial coordinate and the same separation delivers the spherical Bessel functions. Legendre is the hub.

Key takeaways

The Legendre polynomials $P_\ell(x)$ are the polynomial solutions of $(1-x^2)y''-2xy'+\ell(\ell+1)y=0$ on $[-1,1]$ that stay finite at the poles $x=\pm1$ ; regularity is what quantises the eigenvalue to $\ell(\ell+1)$ .
They form a complete orthogonal set with weight $w(x)=1$ : $\int_{-1}^{1}P_\ell P_m\,\mathrm{d}x=\frac{2}{2\ell+1}\delta_{\ell m}$ , normalised by $P_\ell(1)=1$ .
The generating function $(1-2xt+t^2)^{-1/2}=\sum_\ell P_\ell(x)t^\ell$ is the Coulomb/Newton kernel $1/|\mathbf r-\mathbf r'|$ , so $P_\ell$ is the angular pattern of the $\ell$ -th multipole.
Two closed forms: Rodrigues $P_\ell=\frac{1}{2^\ell\ell!}\frac{\mathrm{d}^\ell}{\mathrm{d}x^\ell}(x^2-1)^\ell$ , and Bonnet's recurrence $(\ell+1)P_{\ell+1}=(2\ell+1)xP_\ell-\ell P_{\ell-1}$ for stable numerics.
Low orders: $P_0=1,\;P_1=x,\;P_2=\tfrac12(3x^2-1),\;P_3=\tfrac12(5x^3-3x),\;P_4=\tfrac18(35x^4-30x^2+3)$ ; each $P_\ell$ has exactly $\ell$ zeros in $(-1,1)$ and parity $(-1)^\ell$ .
They are the $m=0$ spherical harmonics and the building blocks of associated Legendre functions; their zeros are the optimal Gauss–Legendre quadrature nodes.
Key idea: $P_\ell(x)$ is the unique degree- $\ell$ polynomial solution of Legendre's equation regular at $x=\pm1$ and normalised to $P_\ell(1)=1$ — equivalently the coefficients in $(1-2xt+t^2)^{-1/2}=\sum_\ell P_\ell(x)t^\ell$ , which is precisely the inverse-distance kernel of potential theory.

Where this lives today

The Legendre family runs inside working instruments and codes right now. Every time a routine evaluates an integral by Gauss–Legendre quadrature it uses the zeros of $P_n$ as nodes — the default in finite-element solvers, spectral methods, and quantum-chemistry packages, because $n$ points nail polynomials up to degree $2n-1$ . In cosmology the sky is a function on a sphere: the cosmic-microwave-background temperature map is expanded in spherical harmonics and compressed into the angular power spectrum $C_\ell$ ; Planck measured it to multipole $\ell\approx2500$ , and its acoustic peaks pin down the geometry and contents of the universe. Geodesy does the same for the Earth's gravity field — the EGM2008 geoid is a spherical-harmonic series to degree and order 2190, tens of thousands of coefficients that missions like GRACE refine. In quantum scattering, the partial-wave expansion sorts an amplitude by orbital angular momentum $\ell$ , each channel weighted by $P_\ell(\cos\theta)$ , and resonances show up as bumps in a particular $\ell$ -wave. The open frontier is high-degree stability: transforms at $\ell\sim10^4$ demand recurrence algorithms that fight underflow and round-off — keeping Bonnet's 1850s recurrence a subject of active research.

The physical picture

Think of every function on $[-1,1]$ as a vector, with $\int_{-1}^1 fg\,\mathrm{d}x$ as its dot product. Feed the monomials $1,x,x^2,\dots$ through Gram–Schmidt under that inner product and rescale each so it passes through $1$ at $x=1$ : out drops exactly $P_0,P_1,P_2,\dots$ . They are the natural orthogonal axes for any quantity depending on $\cos\theta$ over a sphere, just as sines and cosines are on a circle. Higher $\ell$ means finer angular structure — $P_\ell$ crosses zero exactly $\ell$ times inside $(-1,1)$ , a Sturm oscillation property — so a Legendre series resolves a shape the way higher multipoles resolve a charge distribution. The first figure shows the whole family, pinned together at the right endpoint.

Figure

The five lowest Legendre polynomials P₀ through P₄ on the interval −1 to 1: every curve is pinned to Pₗ(1) = 1 at the right endpoint and to Pₗ(−1) = (−1)^ℓ at the left, while Pₗ crosses zero exactly ℓ times inside the interval — the Sturm oscillation that makes higher ℓ resolve finer angular structure

The generating function turns this abstract basis into physics. Drop a point charge at $\mathbf r'$ and read its potential at $\mathbf r$ ; expand $1/|\mathbf r-\mathbf r'|$ in powers of $t=r_</r_>$ and the angular shape of the $\ell$ -th term is exactly $P_\ell(\cos\gamma)$ , with $\gamma$ the opening angle between the two vectors. So $P_0$ is the monopole, $P_1$ the dipole pattern, $P_2$ the quadrupole — the Legendre index simply is the multipole order. The second figure carries a real number through it: for a source at $r'=1$ seen from $r=2$ at $\gamma=60^\circ$ , successive multipoles march the partial sum toward the exact kernel value.

Figure

The Coulomb kernel one-over-distance expanded in Legendre polynomials for r = 2, r′ = 1, γ = 60°: the monopole through octupole terms give partial sums 0.500, 0.625, 0.609, 0.582 that converge geometrically to the exact value 0.5774, because the ratio t = r′/r = 0.5 controls the tail

Two further symmetries earn their keep. Each $P_\ell$ has definite parity $(-1)^\ell$ , so even- $\ell$ terms carry the part of a function symmetric about $x=0$ and odd- $\ell$ the antisymmetric part — halving the work when a source has reflection symmetry. And the exact endpoint values $P_\ell(1)=1$ , $P_\ell(-1)=(-1)^\ell$ make the poles the cleanest places to fix boundary coefficients.

The core equations

Everything begins with the eigenvalue problem. Legendre's equation is the $\phi$ -independent angular piece of $\nabla^2\Phi=0$ after $x=\cos\theta$ ; $P_\ell$ is its regular eigenfunction and $\ell(\ell+1)$ its eigenvalue.

(1-x^2)\frac{\mathrm{d}^2 P_\ell}{\mathrm{d}x^2}-2x\frac{\mathrm{d}P_\ell}{\mathrm{d}x}+\ell(\ell+1)P_\ell=0 \qquad (1)

where: $x$ — independent variable, $x=\cos\theta\in[-1,1]$ ; $\ell$ — non-negative integer degree, physically the multipole order; $\ell(\ell+1)$ — separation constant, quantised by regularity at $x=\pm1$ . Reading it: the coefficient $(1-x^2)$ vanishes at the endpoints, making $x=\pm1$ regular singular points — why finiteness there is a real constraint. Rodrigues' formula packages the same polynomial as a single differential operation — the workhorse for explicit expressions and orthogonality proofs.

P_\ell(x)=\frac{1}{2^\ell\,\ell!}\frac{\mathrm{d}^\ell}{\mathrm{d}x^\ell}\big(x^2-1\big)^\ell \qquad (2)

where: $\mathrm{d}^\ell/\mathrm{d}x^\ell$ — $\ell$ -fold ordinary derivative; $2^\ell\ell!$ — normalisation constant fixed so that $P_\ell(1)=1$ . Reading it: differentiating a degree- $2\ell$ polynomial $\ell$ times leaves degree $\ell$ . The generating function is the identity that makes Legendre polynomials physical rather than merely orthogonal — expanding its algebraic left side in $t$ produces the whole family as coefficients.

g(x,t)=\big(1-2xt+t^2\big)^{-1/2}=\sum_{\ell=0}^{\infty}P_\ell(x)\,t^\ell,\qquad |t|<1 \qquad (3)

where: $t$ — expansion parameter, physically $r_</r_>$ ; $x=\cos\gamma$ — cosine of the angle between field and source points. Reading it: the bracket is $|\mathbf r-\mathbf r'|^2/r_>^2$ , so $g$ is the inverse-distance kernel — every statement about $P_\ell$ is one about a potential. Bonnet's three-term recurrence connects neighbouring degrees and is the numerically stable way to evaluate $P_\ell$ (upward recursion, well-conditioned for $|x|\le1$ ).

(\ell+1)P_{\ell+1}(x)=(2\ell+1)\,x\,P_\ell(x)-\ell\,P_{\ell-1}(x) \qquad (4)

where: $P_{\ell\pm1}$ — neighbouring-degree polynomials; start from $P_0=1$ , $P_1=x$ . Reading it: $\ell=1$ gives $2P_2=3x^2-1$ , i.e. $P_2=\tfrac12(3x^2-1)$ — the recurrence builds the table. Orthogonality with weight unity, plus the explicit norm, make the family a genuine basis.

\int_{-1}^{1}P_\ell(x)P_m(x)\,\mathrm{d}x=\frac{2}{2\ell+1}\,\delta_{\ell m} \qquad (5)

where: $\delta_{\ell m}$ — Kronecker delta; $2/(2\ell+1)$ — the squared $L^2$ norm of $P_\ell$ on $[-1,1]$ . Reading it: distinct degrees are orthogonal; the norm shrinks like $1/\ell$ . Completeness turns orthogonality into an expansion: any square-integrable function is a Legendre series, coefficients read off by projection.

f(x)=\sum_{\ell=0}^{\infty}c_\ell P_\ell(x),\qquad c_\ell=\frac{2\ell+1}{2}\int_{-1}^{1}f(x)P_\ell(x)\,\mathrm{d}x \qquad (6)

where: $f$ — any function with $\int_{-1}^1 f^2\,\mathrm{d}x<\infty$ ; $c_\ell$ — coefficient, from projecting $f$ onto $P_\ell$ and dividing by the norm. Reading it: the factor $(2\ell+1)/2$ is the reciprocal norm from equation (5). The multipole expansion is the generating function in physical variables — the single most used consequence of the theory.

\frac{1}{|\mathbf r-\mathbf r'|}=\sum_{\ell=0}^{\infty}\frac{r_<^{\,\ell}}{r_>^{\,\ell+1}}\,P_\ell(\cos\gamma) \qquad (7)

where: $r_<,\,r_>$ — the smaller and larger of $|\mathbf r|,|\mathbf r'|$ ; $\gamma$ — angle between $\mathbf r$ and $\mathbf r'$ . Reading it: convergence is set by $t=r_</r_><1$ , so the series converges geometrically away from $\mathbf r=\mathbf r'$ . Finally a derivative recurrence, useful for quadrature weights and term-by-term integration of a series.

(2\ell+1)P_\ell(x)=P_{\ell+1}'(x)-P_{\ell-1}'(x) \qquad (8)

where: $P_\ell'$ — first derivative with respect to $x$ . Reading it: integrating over $[-1,1]$ collapses the right side to boundary terms — the trick behind many Legendre integrals.

Derivation AFrom regularity at the poles to the quantised eigenvalue

The goal is to show why the eigenvalue is locked to the ladder $\ell(\ell+1)$ and why the eigenfunction is a polynomial rather than an infinite series. We insert a power series into equation (1), read off the two-step recurrence it forces on the coefficients, examine how the series behaves as $x\to\pm1$ , and find that finiteness at the poles — a physical demand on a potential — can only be met if the series chops off. That truncation both quantises the eigenvalue and makes the solution a polynomial of definite degree.

Step 1. Seek $y=\sum_{n\ge0}a_n x^n$ and substitute into $(1-x^2)y''-2xy'+\lambda y=0$ . Matching the coefficient of $x^n$ gives the two-step recurrence $a_{n+2}=\dfrac{n(n+1)-\lambda}{(n+1)(n+2)}\,a_n$ , which links even coefficients among themselves and odd among themselves. Why: the equation has even symmetry, so the even and odd series are independent solutions — exactly the two families we expect from a second-order ODE.

Step 2. Examine the tail. For large $n$ , $a_{n+2}/a_n\to1$ , so the series behaves like $\sum x^{2k}\sim(1-x^2)^{-1}$ , diverging as $x\to\pm1$ . A potential that blows up at the poles is unphysical, so a generic $\lambda$ is not allowed. Why: the endpoints $x=\pm1$ are the poles of the sphere, ordinary points of space where the field must stay finite.

Step 3. The only escape is termination: if $n(n+1)-\lambda$ vanishes at some $n=\ell$ , then $a_{\ell+2}=0$ and the series stops at a polynomial of degree $\ell$ . That requires $\lambda=\ell(\ell+1)$ with $\ell$ a non-negative integer, and the companion series is killed by zeroing its seed. Why: quantisation is born here — a boundary-regularity demand, not any external postulate, selects the discrete spectrum $\ell(\ell+1)$ .

Step 4. Fix the normalisation by $P_\ell(1)=1$ ; the discarded second solution, continued as an infinite series, is $Q_\ell(x)$ , diverging logarithmically at $x=\pm1$ . Why: this pins $P_\ell$ uniquely and explains why $Q_\ell$ survives only for regions that exclude the poles, such as shells. *Check: substituting $P_2=\tfrac12(3x^2-1)$ into equation (1) with $\ell=2$ gives $(1-x^2)(3)-2x(3x)+6\cdot\tfrac12(3x^2-1)=3-3x^2-6x^2+9x^2-3=0$ .*

Go deeper · Derivation B · 1 minBonnet's recurrence straight from the generating function

Derivation BBonnet's recurrence straight from the generating function

The aim is to manufacture the three-term recurrence (4) without writing down a single polynomial, using only the generating function (3). The manoeuvre is classic: differentiate $g(x,t)$ with respect to $t$ , notice the derivative obeys a tidy first-order relation in $t$ , clear the denominator, and match the coefficient of a general power $t^\ell$ on the two sides. Because $g$ encodes the entire family at once, one differentiation yields a relation valid for all $\ell$ — far cheaper than proving the recurrence polynomial by polynomial.

Step 1. Differentiate $g=(1-2xt+t^2)^{-1/2}$ in $t$ :

\frac{\partial g}{\partial t}=\frac{x-t}{(1-2xt+t^2)^{3/2}}=\frac{x-t}{1-2xt+t^2}\,g.

Hence $(1-2xt+t^2)\,\partial_t g=(x-t)\,g$ . Why: this replaces the awkward $3/2$ power with the original $g$ , giving an exact first-order relation between $g$ and its $t$ -derivative that the series must inherit.

Step 2. Insert the series $g=\sum_\ell P_\ell t^\ell$ and $\partial_t g=\sum_\ell \ell P_\ell t^{\ell-1}$ :

(1-2xt+t^2)\sum_\ell \ell P_\ell t^{\ell-1}=(x-t)\sum_\ell P_\ell t^\ell.

Why: turning the functional relation into a statement about coefficients is what converts calculus into an algebraic recurrence.

Step 3. Collect the coefficient of $t^\ell$ : the left side gives $(\ell+1)P_{\ell+1}-2x\ell P_\ell+(\ell-1)P_{\ell-1}$ , the right side $xP_\ell-P_{\ell-1}$ . Equate and simplify to

(\ell+1)P_{\ell+1}=(2\ell+1)\,x\,P_\ell-\ell\,P_{\ell-1}.

Why: matching one power of $t$ delivers Bonnet's recurrence in a line, valid for every $\ell\ge1$ .

Step 4. Seed with $g=1+xt+O(t^2)$ , giving $P_0=1$ , $P_1=x$ , and crank: $\ell=1$ yields $P_2=\tfrac12(3x^2-1)$ , then $P_3=\tfrac12(5x^3-3x)$ , then $P_4=\tfrac18(35x^4-30x^2+3)$ . Why: the recurrence plus two seeds reproduces the entire table and is the stable algorithm used in practice. *As a consistency check: upward recursion from $P_0=1,P_1=x$ reproduces $P_2,P_3,P_4$ exactly, and $P_\ell(1)=1$ at every step since $(\ell+1)\cdot1=(2\ell+1)\cdot1-\ell\cdot1$ .*

Go deeper · Derivation C · 1 minOrthogonality and the norm 2 over (2ℓ+1)

Derivation COrthogonality and the norm 2 over (2ℓ+1)

The target is the full orthonormality statement (5): distinct Legendre polynomials are orthogonal on $[-1,1]$ , and the squared norm of $P_\ell$ is $2/(2\ell+1)$ . Orthogonality follows from the self-adjoint (Sturm–Liouville) structure, with the boundary term killed by the same factor $(1-x^2)$ that vanishes at the poles. The norm is then extracted by a slick integral: square $g(x,t)$ , integrate over $x$ , keep only the diagonal by orthogonality, and recognise the elementary left-hand integral as a logarithmic series in $t$ whose coefficients read off the norm.

Step 1. Write equation (1) in self-adjoint form, $\frac{\mathrm{d}}{\mathrm{d}x}[(1-x^2)P_\ell']+\ell(\ell+1)P_\ell=0$ , i.e. $\frac{\mathrm{d}}{\mathrm{d}x}[p\,y']+\lambda w\,y=0$ with $p=1-x^2$ , weight $w=1$ . Multiply the $\ell$ -equation by $P_m$ , the $m$ -equation by $P_\ell$ , and subtract to get a total derivative:

\frac{\mathrm{d}}{\mathrm{d}x}\!\big[(1-x^2)(P_\ell P_m'-P_m P_\ell')\big]=\big[\ell(\ell+1)-m(m+1)\big]P_\ell P_m.

Why: the self-adjoint form is what makes the left side an exact derivative, so integrating it will only leave boundary data.

Step 2. Integrate over $[-1,1]$ . The boundary term carries the factor $p=1-x^2$ , which is zero at both $x=\pm1$ , so it vanishes and

\big[\ell(\ell+1)-m(m+1)\big]\int_{-1}^{1}P_\ell P_m\,\mathrm{d}x=0.

For $\ell\neq m$ the bracket is nonzero, forcing $\int_{-1}^1 P_\ell P_m\,\mathrm{d}x=0$ . Why: orthogonality is a free gift of the endpoints being the singular points $p(\pm1)=0$ — no integration by parts survives.

Step 3. For the norm, square the generating function and integrate:

\int_{-1}^{1}g(x,t)^2\,\mathrm{d}x=\int_{-1}^{1}\frac{\mathrm{d}x}{1-2xt+t^2}=\sum_{\ell,m}t^{\ell+m}\int_{-1}^{1}P_\ell P_m\,\mathrm{d}x=\sum_{\ell}t^{2\ell}\!\int_{-1}^{1}P_\ell^2\,\mathrm{d}x,

the last step using the orthogonality just proved. Why: squaring generates every cross term, but orthogonality annihilates all but the diagonal, leaving a clean sum over the norms alone.

Step 4. Do the elementary integral: $\int_{-1}^1\frac{\mathrm{d}x}{1-2xt+t^2}=\frac{1}{2t}\ln\frac{1+t}{1-t}=\sum_{\ell\ge0}\frac{2}{2\ell+1}t^{2\ell}$ . Matching the coefficient of $t^{2\ell}$ against Step 3 gives

\int_{-1}^{1}P_\ell^2(x)\,\mathrm{d}x=\frac{2}{2\ell+1}.

Why: the entire norm sequence drops out of one logarithmic expansion — no polynomial is integrated directly. *As a concrete check: the series $\tfrac12\ln\frac{1+t}{1-t}=\sum_\ell t^{2\ell+1}/(2\ell+1)$ is standard, and direct integration gives $\int_{-1}^1 P_2^2\,\mathrm{d}x=\int_{-1}^1\tfrac14(3x^2-1)^2\mathrm{d}x=2/5=2/(2\cdot2+1)$ .*

Limiting cases

Case	Result	Meaning
Endpoint values	$P_\ell(1)=1,\ \ P_\ell(-1)=(-1)^\ell$	Fixes the normalisation; read straight off $g$ at $x=\pm1$ , where $(1\mp t)^{-1}=\sum t^\ell(\pm1)^\ell$ .
Value at the origin	$P_{2k}(0)=\dfrac{(-1)^k(2k)!}{2^{2k}(k!)^2},\ \ P_{2k+1}(0)=0$	Odd polynomials vanish at $x=0$ by parity; even ones give $1,-\tfrac12,\tfrac38,\dots$ (so $P_2(0)=-\tfrac12$ , $P_4(0)=\tfrac38$ ).
Large- $\ell$ (Laplace–Heine)	$P_\ell(\cos\theta)\sim\sqrt{\dfrac{2}{\pi\ell\sin\theta}}\,\cos\!\big[(\ell+\tfrac12)\theta-\tfrac{\pi}{4}\big]$	For fixed $\theta\in(0,\pi)$ as $\ell\to\infty$ : oscillatory, amplitude $\sim\ell^{-1/2}$ , away from the poles.
Near the pole (Mehler–Heine)	$\lim_{\ell\to\infty}P_\ell\!\big(\cos\tfrac{x}{\ell}\big)=J_0(x)$	Zooming into the forward direction turns the polynomial into a Bessel function — the bridge from spherical to cylindrical geometry.
Parity	$P_\ell(-x)=(-1)^\ell P_\ell(x)$	Even index gives an even polynomial, odd index an odd one; immediate from the $t\to-t,\,x\to-x$ symmetry of $g$ .

Common trap“Legendre polynomials are normalised to unit length, so just expand and sum the coefficients.”

They are not. The universal convention is $P_\ell(1)=1$ , not unit $L^2$ norm — the squared norm is $2/(2\ell+1)$ , so the projection coefficient carries the factor $c_\ell=\frac{2\ell+1}{2}\int f P_\ell\,\mathrm{d}x$ . Forget the $(2\ell+1)/2$ and every coefficient is wrong by that factor. Three traps travel with it. First, the argument is $x=\cos\theta$ , so the interval is $[-1,1]$ — plug in $\theta$ directly and you get nonsense. Second, the second solution $Q_\ell(x)$ is singular at $x=\pm1$ , so it is dropped for the interior of a sphere but kept for shells or exterior regions that exclude the poles. Third, do not confuse the plain $P_\ell$ ( $m=0$ ) with the associated functions $P_\ell^m$ that appear once $\phi$ -dependence returns.

Worked examples

1 · Gauss–Legendre quadrature: the zeros do the integrating (confidence builder)

The two-point Gauss–Legendre rule places its nodes at the zeros of $P_2(x)=\tfrac12(3x^2-1)$ . Find the nodes and weights, then estimate $\int_{-1}^1 e^x\,\mathrm{d}x$ and compare with the exact value.

Step 1. The nodes are the roots of $P_2$ : $3x^2-1=0\Rightarrow x_{1,2}=\pm1/\sqrt3=\pm0.57735$ . By symmetry the weights are equal, and since the rule must reproduce $\int_{-1}^1 1\,\mathrm{d}x=2$ , each is $w_1=w_2=1$ . (The general formula $w_i=2/[(1-x_i^2)P_2'(x_i)^2]$ with $P_2'=3x$ gives $w_i=2/[\tfrac23\cdot3]=1$ .)

Step 2. Apply the rule: $\int_{-1}^1 e^x\,\mathrm{d}x\approx w_1 e^{x_1}+w_2 e^{x_2}=e^{-0.57735}+e^{+0.57735}=0.56138+1.78131=\mathbf{2.34270}$ .

Step 3. Compare: the exact integral is $e-e^{-1}=2.71828-0.36788=\mathbf{2.35040}$ , so the two-point rule lands within $0.0077$ , a relative error of $0.33\%$ — from just two function evaluations. That is the payoff: with nodes at $P_n$ 's zeros, the $n$ -point rule integrates every polynomial up to degree $2n-1$ exactly, so it "sees" the cubic part of $e^x$ for free. The three-point rule (nodes $0,\pm\sqrt{3/5}=\pm0.77460$ , weights $\tfrac89,\tfrac59,\tfrac59$ ) gives $2.35034$ , cutting the error to $0.003\%$ . Two zeros of $P_2$ integrate $e^x$ to 0.33%; three zeros of $P_3$ reach 0.003%.

2 · A boundary-value problem on a sphere

A sphere of radius $a$ is held at surface potential $V(\theta)=V_0\cos^2\theta$ , with no charge outside. Find the potential everywhere outside the sphere.

Step 1. Azimuthal symmetry makes the exterior solution of $\nabla^2\Phi=0$ a Legendre series in decaying powers: $\Phi(r,\theta)=\sum_\ell A_\ell\,r^{-(\ell+1)}P_\ell(\cos\theta)$ . The task reduces to matching the boundary data — pure algebra, no integrals if we recognise the combination.

Step 2. Express the boundary data in Legendre polynomials. Since $P_2(\cos\theta)=\tfrac12(3\cos^2\theta-1)$ , invert: $\cos^2\theta=\tfrac13[2P_2+1]=\tfrac13 P_0+\tfrac23 P_2$ . So $V(\theta)=V_0[\tfrac13 P_0+\tfrac23 P_2]$ — only $\ell=0$ and $\ell=2$ appear, all other $A_\ell$ vanish.

Step 3. Match at $r=a$ term by term: $A_0 a^{-1}=\tfrac13 V_0$ and $A_2 a^{-3}=\tfrac23 V_0$ , giving $A_0=\tfrac13 V_0 a$ and $A_2=\tfrac23 V_0 a^3$ . Hence

\Phi(r,\theta)=\frac{V_0}{3}\,\frac{a}{r}+\frac{2V_0}{3}\,\frac{a^3}{r^3}\,P_2(\cos\theta),\qquad r\ge a.

Step 4. Read the physics. The $\ell=0$ piece $\tfrac13 V_0\,a/r$ is a monopole — the sphere carries net charge $Q=4\pi\varepsilon_0\cdot\tfrac13 V_0 a$ , the angular average of the surface potential. The $\ell=2$ piece is a pure quadrupole, decaying as $a^3/r^3$ ; there is no dipole because $\cos^2\theta$ is even. Far away the monopole dominates; close in, the quadrupole distortion shows through. $\Phi = \tfrac13 V_0(a/r) + \tfrac23 V_0(a^3/r^3)P_2(\cos\theta)$ : a monopole plus a quadrupole, no dipole.

Physics of the portrait

Legendre polynomials are where three enterprises — the shape of the Earth, the stability of the solar system, and the theory of orthogonal functions — turned out to be the same mathematics. The man whose name they carry met them while asking a concrete question about gravity.

Origins · 1 minAdrien-Marie Legendre and how this was discovered

Adrien-Marie Legendre

1752–1833 · potential theory and celestial mechanics

Adrien-Marie Legendre introduced the polynomials that bear his name in a 1782 memoir on the gravitational attraction of spheroids — the practical problem of computing the pull of a slightly flattened, spinning body like the Earth. Expanding the Newtonian potential $1/|\mathbf r-\mathbf r'|$ in powers of the distance ratio, he found the coefficients were a new family of polynomials in $\cos\gamma$ , and proved their orthogonality and recurrence long before anyone spoke of Sturm–Liouville theory. Pierre-Simon Laplace, working the same celestial mechanics, immediately deployed them (and their $\phi$ -dependent cousins, the spherical harmonics) against planetary perturbations, so that for decades they were known simply as the "coefficients of Laplace." Legendre was one of the great systematisers of eighteenth-century analysis — elliptic integrals, the method of least squares (1805, sparking a bitter priority dispute with Gauss), number theory — yet it is these harmonics of the sphere that keep his name in daily use. The closed differential form came later from Olinde Rodrigues (1816) and the three-term recurrence from Ossian Bonnet; but the object, and the instinct that a potential is best organised by its multipoles, are Legendre's.

History

1782

— Legendre introduces the polynomials in his memoir on the attraction of spheroids

1785

— Laplace extends them to the full spherical harmonics in his lunar and planetary theory

1816

— Olinde Rodrigues publishes the differential-operator (Rodrigues) formula for Pₗ

1861

— Heine's Handbuch der Kugelfunctionen systematises the theory of Legendre and spherical functions

1885

— Stieltjes proves that the zeros of Pₗ are the equilibrium positions of charges in a logarithmic potential

Going further

Texts. Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 15 (Legendre functions) · Hassani, Mathematical Physics, Ch. on classical orthogonal polynomials and Sturm–Liouville theory · Stone & Goldbart, Mathematics for Physics, Sturm–Liouville and special-function chapters · Byron & Fuller, Mathematics of Classical and Quantum Physics, Ch. 5–6 · Jackson, Classical Electrodynamics, 3rd ed., Ch. 3 (multipole expansion) · Abramowitz & Stegun, Handbook of Mathematical Functions, Ch. 8 (tables and asymptotics). All key results can be checked numerically: Bonnet's recurrence and Rodrigues' formula reproduce $P_0$ – $P_4$ ; the norm $2/(2\ell+1)$ is confirmed by the squared-generating-function argument and by direct integration of $P_2^2$ ; the two- and three-point Gauss–Legendre estimates of $\int_{-1}^1 e^x\,\mathrm{d}x$ are $2.34270$ and $2.35034$ against the exact value $2.35040$ ; the multipole partial sums for $r=2,r'=1,\gamma=60^\circ$ march as $0.500,0.625,0.609,0.582\to0.5774$ ; and the decomposition $\cos^2\theta=\tfrac13P_0+\tfrac23P_2$ holds exactly. Laplace–Heine and Mehler–Heine asymptotics are treated in Abramowitz & Stegun Ch. 8.

Prerequisites. Orthogonal Polynomials: General Theory · Generating Functions for Special Functions · The Gamma Function. Leads to. Associated Legendre Functions · Spherical Harmonics · Recurrence Relations for Special Functions · Spherical Bessel Functions.

Ask the tutor: Want to see how the same $P_\ell(\cos\theta)$ grows azimuthal structure into the associated Legendre functions and spherical harmonics — or how its zeros power the quadrature inside your solver? Ask, and we will build it from here.

Are Legendre polynomials the same thing as spherical harmonics?

Not quite — they are the azimuthally symmetric special case. A spherical harmonic $Y_\ell^m(\theta,\phi)$ depends on both angles; setting $m=0$ removes the $\phi$ -dependence and leaves a Legendre polynomial up to a constant: $P_\ell(\cos\theta)=\sqrt{4\pi/(2\ell+1)}\,Y_\ell^0$ . The general $Y_\ell^m$ are the associated Legendre functions $P_\ell^m(\cos\theta)$ times $e^{im\phi}$ , and the plain $P_\ell$ is just $P_\ell^0$ . So whenever a problem has rotational symmetry about an axis — a uniform field, an on-axis charge — the full spherical-harmonic machinery collapses to a single Legendre series in $\cos\theta$ .

Why does the generating function have that particular square-root form?

Because it is a disguised distance. The separation between two points at radii $r_>,r_<$ with opening angle $\gamma$ is $|\mathbf r-\mathbf r'|=\sqrt{r_>^2+r_<^2-2r_>r_<\cos\gamma}$ . Factor out $r_>$ and set $t=r_</r_>$ , $x=\cos\gamma$ : the bracket becomes $r_>\sqrt{1-2xt+t^2}$ , so $1/|\mathbf r-\mathbf r'|=(1/r_>)g(x,t)$ . The generating function is not chosen for convenience — it is the Newtonian/Coulomb kernel itself, which is why its coefficients in $t$ are the multipole patterns and why $P_\ell$ shows up in every potential problem.

When do I need Q_ℓ, the second solution?

Whenever your region excludes the poles $x=\pm1$ . The two solutions are $P_\ell$ (finite on all of $[-1,1]$ ) and $Q_\ell$ (finite on the open interval but diverging logarithmically at the endpoints). For the interior of a solid sphere, or wherever the polar axis $\theta=0,\pi$ is in the domain, $Q_\ell$ is unphysical and dropped. But for a region that avoids the axis — a cone, or the space between two coaxial cones — $Q_\ell$ is regular and both solutions are needed. It is exactly like keeping the divergent $1/r$ radial solution of Laplace's equation when the origin lies outside your domain.

How do I actually compute P_ℓ(x) for large ℓ on a computer?

Use Bonnet's recurrence upward, not Rodrigues' formula. Start from $P_0=1$ , $P_1=x$ and iterate $P_{\ell+1}=[(2\ell+1)xP_\ell-\ell P_{\ell-1}]/(\ell+1)$ . For $|x|\le1$ this is numerically stable — the polynomials stay bounded by $|P_\ell(x)|\le1$ , so nothing overflows and errors do not amplify. Rodrigues' formula, by contrast, needs $\ell$ -fold differentiation and suffers catastrophic cancellation; it is a proof tool, not a computational one. For the zeros (quadrature nodes), use Newton's method seeded with $x_i\approx\cos[\pi(i-\tfrac14)/(\ell+\tfrac12)]$ , evaluating $P_\ell$ and $P_\ell'$ via the recurrence.

What is the connection to Bessel functions?

It is a limit linking spherical and cylindrical geometry. Away from the poles, a high-degree $P_\ell(\cos\theta)$ oscillates like a damped cosine (Laplace–Heine), the spherical analogue of a plane wave. Zoom into the forward direction with argument $\cos(x/\ell)$ and let $\ell\to\infty$ : the Mehler–Heine formula gives $P_\ell(\cos(x/\ell))\to J_0(x)$ — the polynomial becomes a Bessel function. Physically, a tiny patch of a huge sphere looks flat, and the harmonics on a flat patch with circular symmetry are Bessel functions. That is why the two families are neighbours in every course.

Why is P_ℓ(1) = 1 chosen instead of unit norm?

Because $x=1$ is the physically meaningful point: it is $\gamma=0$ , two directions aligned. Then $P_\ell(1)=1$ means every multipole contributes with unit weight along the axis, so the on-axis potential of a point charge is the geometric series $\sum t^\ell=1/(1-t)=1/(r_>-r_<)$ , the elementary result. Unit $L^2$ normalisation would instead scatter factors of $\sqrt{(2\ell+1)/2}$ through every formula. The price of $P_\ell(1)=1$ is the explicit norm $2/(2\ell+1)$ you must divide by for series coefficients — a small cost for clean endpoint and generating-function values.