---
title: Bessel's Equation, Legendre Polynomials, and Special Functions
module: Series Solutions and Special Functions
moduleNumber: 4
lessonNumber: 3
order: 403
summary: >
  Bessel's equation puts the Frobenius machinery through all three of its cases
  and produces the functions J and Y that govern anything vibrating or
  diffusing with circular symmetry. The gamma function extends the factorial so
  that Bessel functions of every order make sense; Legendre's equation, run
  through the hypergeometric form, yields the polynomials that play the same
  role in spherical geometry. Orthogonality ties both families to the
  eigenfunction expansions of Sturm–Liouville theory.
topics: [Series Solutions and Special Functions]
sources:
  - book: Boyce
    ref: "Ch. 5 §5.7 Bessel's Equation; §5.3 (Legendre problems)"
  - book: Simmons
    ref: "Ch. 8 §44 Legendre Polynomials; §45 Properties of Legendre Polynomials; §46 Bessel Functions. The Gamma Function; §47 Properties of Bessel Functions"
draft: false
---

Two equations carry the [Frobenius method](/differential-equations/series-solutions/regular-singular-frobenius) into physics. **Bessel's equation** governs waves and diffusion in circular and cylindrical geometry (the vibrating drumhead, heat in a cylinder, the hanging chain), and its indicial roots $\pm\nu$ realize every case of the three-case theorem as $\nu$ varies. **Legendre's equation** plays the same role for spheres, and its bounded solutions are polynomials. Both families carry an orthogonality relation, and that is what makes them usable: an arbitrary function expands in Bessel functions or Legendre polynomials exactly as it expands in [sines and cosines](/differential-equations/pdes-fourier-bvp/fourier-series). One prerequisite comes first, a factorial that accepts non-integer arguments.

## The gamma function

The Frobenius series for Bessel's equation of order $\nu$ has coefficients like $1/\big( n!\,(\nu + n)! \big)$, which are meaningless when $\nu$ is not an integer unless the factorial is extended.

> **Definition (Gamma function).** For $p > 0$,
> $$
> \Gamma(p) = \int_0^{\infty} t^{p-1} e^{-t}\,\d t .
> $$
> Integration by parts gives the functional equation $\Gamma(p + 1) = p\,\Gamma(p)$, and $\Gamma(1) = 1$, so
> $$
> \Gamma(n + 1) = n! \qquad (n = 0, 1, 2, \ldots) :
> $$
> the gamma function interpolates the factorial. Writing $p! = \Gamma(p+1)$ extends the factorial symbol to all $p > -1$, and the rearrangement $\Gamma(p) = \Gamma(p+1)/p$ extends $\Gamma$ to negative non-integer arguments one unit strip at a time, with poles at $0, -1, -2, \ldots$

The one non-obvious value worth memorizing comes from the Gaussian integral: substituting $t = s^2$,

$$
\Gamma\!\left(\tfrac12\right) = 2\int_0^{\infty} e^{-s^2}\,\d s = \sqrt{\pi},
\qquad\text{so}\qquad
\left(\tfrac12\right)! = \tfrac12\sqrt{\pi},
\quad
\left(n + \tfrac12\right)! = \frac{(2n+1)!}{2^{2n+1}\,n!}\sqrt{\pi}.
$$

Because $1/\Gamma$ vanishes at the poles, the convention $1/(-m)! = 0$ for positive integers $m$ is consistent, and it is the convention the Bessel series below requires.[^simmons-gamma]

$$
% caption: The gamma function interpolates the factorial: the curve is $y = \Gamma(x)$
% for $x > 0$, and the marked points are $\Gamma(n+1) = n!$ at the integers. The dip
% bottoms out near $x = 1.46$ before the factorial growth takes over.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (-0.2,0) -- (6.4,0) node[below] {$x$};
\draw[->, black] (0,-0.3) -- (0,4.3) node[left] {$y$};
\foreach \x/\l in {1.4/1, 2.8/2, 4.2/3, 5.6/4}
  \draw[black] (\x,0.05) -- (\x,-0.05) node[below] {\l};
\foreach \y/\l in {0.55/1, 1.1/2, 3.3/6}
  \draw[black] (0.05,\y) -- (-0.05,\y) node[left] {\l};
\draw[acc, very thick] plot[smooth] coordinates {(0.35,1.99) (0.42,1.65) (0.49,1.40) (0.56,1.22) (0.70,0.97) (0.84,0.82) (1.05,0.67) (1.40,0.55) (1.75,0.50) (2.10,0.49) (2.45,0.51) (2.80,0.55) (3.50,0.73) (4.20,1.10) (4.55,1.40) (4.90,1.83) (5.25,2.43) (5.60,3.30) (5.81,4.00)};
\fill[acc] (1.40,0.55) circle (2.2pt);
\fill[acc] (2.80,0.55) circle (2.2pt);
\fill[acc] (4.20,1.10) circle (2.2pt);
\fill[acc] (5.60,3.30) circle (2.2pt);
\node[anchor=north west, text=black] at (1.36,0.42) {$0!$};
\node[anchor=north west, text=black] at (2.76,0.42) {$1!$};
\node[anchor=north west, text=black] at (4.30,0.98) {$2!$};
\node[anchor=west, text=black] at (5.72,3.25) {$3!$};
\draw[acc, dashed] (0.35,1.99) .. controls (0.28,2.6) .. (0.24,4.1);
\end{tikzpicture}
$$

## Bessel's equation

> **Definition (Bessel's equation).** For a constant $\nu \ge 0$,
> $$
> x^2 y'' + x y' + (x^2 - \nu^2)\,y = 0 .
> $$
> The origin is a regular singular point with $p_0 = 1$, $q_0 = -\nu^2$, so the indicial equation is $r^2 - \nu^2 = 0$ and the exponents at the singularity are $r = \pm\nu$.

The gap between the exponents is $2\nu$, so the parameter sweeps through every case of the Frobenius theorem: equal roots at $\nu = 0$, an integer gap at $\nu = \tfrac12, 1, \tfrac32, \ldots$, and the generic case otherwise.[^boyce-bessel]

### The first solution: $J_\nu$

Running the recurrence for the larger exponent $r = \nu$ kills the odd coefficients and steps the even ones by

$$
a_{2m} = -\frac{a_{2m-2}}{2^2\,m\,(\nu + m)} ,
$$

which telescopes into factorials. With the conventional normalization $a_0 = 1/\big( 2^{\nu}\,\Gamma(\nu + 1) \big)$:

> **Definition (Bessel function of the first kind).**
> $$
> J_{\nu}(x) = \sum_{n=0}^{\infty} \frac{(-1)^n}{n!\;\Gamma(n + \nu + 1)} \left( \frac{x}{2} \right)^{2n + \nu} ,
> $$
> convergent for all $x$. In particular
> $$
> J_0(x) = 1 - \frac{x^2}{2^2} + \frac{x^4}{2^2 4^2} - \frac{x^6}{2^2 4^2 6^2} + \cdots ,
> \qquad
> J_1(x) = \frac{x}{2} - \frac{1}{1!\,2!}\left(\frac{x}{2}\right)^3 + \frac{1}{2!\,3!}\left(\frac{x}{2}\right)^5 - \cdots .
> $$

> **Worked example.** Derive the series for $J_0$. At $\nu = 0$ the recurrence $a_{2m} = -a_{2m-2}/\big( 2^2 m(\nu + m) \big)$ becomes $a_{2m} = -a_{2m-2}/(2m)^2$, and $a_1 = 0$ forces every odd coefficient to vanish. Running it from $a_0 = 1$,[^boyce-bessel]
>
> $$
> a_2 = -\frac{1}{2^2}, \qquad
> a_4 = \frac{1}{2^2 \cdot 4^2}, \qquad
> a_6 = -\frac{1}{2^2 \cdot 4^2 \cdot 6^2},
> $$
>
> and since $2 \cdot 4 \cdots (2m) = 2^m m!$, the general term is $a_{2m} = (-1)^m/\big( 2^{2m}(m!)^2 \big)$. Hence
>
> $$
> J_0(x) = \sum_{m=0}^{\infty} \frac{(-1)^m}{2^{2m}(m!)^2}\,x^{2m}
>        = 1 - \frac{x^2}{4} + \frac{x^4}{64} - \cdots,
> $$
>
> convergent for all $x$. The equal indicial roots $r_1 = r_2 = 0$ make this the only Frobenius solution; the second carries a logarithm.

$J_0$ and $J_1$ resemble damped versions of cosine and sine, and the resemblance is quantitative. Substituting $u = \sqrt{x}\,y$ transforms Bessel's equation into

$$
u'' + \left( 1 + \frac{1 - 4\nu^2}{4x^2} \right) u = 0 ,
$$

which for large $x$ is approximately $u'' + u = 0$; undoing the substitution suggests $y \approx x^{-1/2} (\text{sinusoid})$, and a more careful analysis confirms it:[^boyce-asymp]

$$
J_{\nu}(x) \;\approx\; \sqrt{\frac{2}{\pi x}}\;\cos\!\left( x - \frac{\pi}{4} - \frac{\nu\pi}{2} \right)
\qquad (x \to \infty),
$$

with error $O(x^{-3/2})$. So every $J_\nu$ oscillates forever with amplitude decaying like $x^{-1/2}$, and has infinitely many positive zeros, asymptotically $\pi$ apart.

$$
% caption: Bessel functions of the first kind. Each oscillates with slowly decaying
% amplitude; the marked zeros of $J_0$ (2.405, 5.520, 8.654, ...) set the frequencies
% of the circular drumhead, and successive orders start flatter at the origin.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\definecolor{warm}{HTML}{B0483A}
\draw[->, black] (-0.3,0) -- (10.9,0) node[below] {$x$};
\draw[->, black] (0,-1.3) -- (0,2.15) node[left] {$y$};
\foreach \x/\l in {1.7/2, 3.4/4, 5.1/6, 6.8/8, 8.5/10, 10.2/12}
  \draw[black] (\x,0.05) -- (\x,-0.05) node[below] {\l};
\draw[black] (0.05,1.8) -- (-0.05,1.8) node[left] {1};
\draw[acc, very thick] plot[smooth] coordinates {(0.00,1.80) (0.34,1.73) (0.68,1.52) (1.02,1.21) (1.36,0.82) (1.70,0.40) (2.04,0.00) (2.38,-0.33) (2.72,-0.58) (3.06,-0.71) (3.40,-0.71) (3.74,-0.62) (4.08,-0.43) (4.42,-0.20) (4.76,0.05) (5.10,0.27) (5.44,0.44) (5.78,0.53) (6.12,0.53) (6.46,0.45) (6.80,0.31) (7.14,0.12) (7.48,-0.07) (7.82,-0.25) (8.16,-0.38) (8.50,-0.44) (8.84,-0.44) (9.18,-0.37) (9.52,-0.24) (9.86,-0.08) (10.20,0.09)};
\draw[warm, thick] plot[smooth] coordinates {(0.00,0.00) (0.34,0.35) (0.68,0.66) (1.02,0.90) (1.36,1.03) (1.70,1.04) (2.04,0.94) (2.38,0.74) (2.72,0.47) (3.06,0.17) (3.40,-0.12) (3.74,-0.36) (4.08,-0.54) (4.42,-0.62) (4.76,-0.60) (5.10,-0.50) (5.44,-0.33) (5.78,-0.12) (6.12,0.10) (6.46,0.29) (6.80,0.42) (7.14,0.49) (7.48,0.48) (7.82,0.39) (8.16,0.25) (8.50,0.08) (8.84,-0.10) (9.18,-0.26) (9.52,-0.37) (9.86,-0.42) (10.20,-0.40)};
\draw[black, thick] plot[smooth] coordinates {(0.00,0.00) (0.34,0.04) (0.68,0.14) (1.02,0.29) (1.36,0.46) (1.70,0.64) (2.04,0.78) (2.38,0.86) (2.72,0.87) (3.06,0.80) (3.40,0.66) (3.74,0.45) (4.08,0.21) (4.42,-0.04) (4.76,-0.26) (5.10,-0.44) (5.44,-0.54) (5.78,-0.56) (6.12,-0.50) (6.46,-0.38) (6.80,-0.20) (7.14,-0.01) (7.48,0.18) (7.82,0.33) (8.16,0.43) (8.50,0.46) (8.84,0.42) (9.18,0.32) (9.52,0.17) (9.86,0.01) (10.20,-0.15)};
\fill[acc] (2.04,0) circle (2.4pt);
\fill[acc] (4.69,0) circle (2.4pt);
\fill[acc] (7.36,0) circle (2.4pt);
\fill[acc] (10.02,0) circle (2.4pt);
\node[acc, anchor=south west] at (0.12,1.92) {$J_0$};
\node[warm, anchor=south] at (1.53,1.12) {$J_1$};
\node[anchor=south, text=black] at (2.72,0.92) {$J_2$};
\node[acc, anchor=north, align=center] at (1.8,-0.72) {zeros of $J_0$};
\end{tikzpicture}
$$

The first few positive zeros, which recur constantly in applications:

| Function | 1st zero | 2nd zero | 3rd zero |
| --- | --- | --- | --- |
| $J_0$ | $2.405$ | $5.520$ | $8.654$ |
| $J_1$ | $3.832$ | $7.016$ | $10.174$ |
| $J_2$ | $5.136$ | $8.417$ | $11.620$ |

The zeros interlace: between consecutive positive zeros of $J_0$ there is exactly one zero of $J_1$, and conversely. This is an immediate consequence of the derivative identities

$$
\frac{\d}{\d x} J_0(x) = -J_1(x),
\qquad
\frac{\d}{\d x}\big[ x\,J_1(x) \big] = x\,J_0(x),
$$

plus Rolle's theorem.[^simmons-props]

### The second solution: three cases realized

The second solution depends on the exponent gap $2\nu$, and Bessel's equation exhibits each behavior of the general theorem.[^boyce-bessel]

- **$2\nu$ not an integer.** Both exponents give Frobenius series, and replacing $\nu$ by $-\nu$ in the series (the gamma function keeps it meaningful) yields an independent solution:
  $$
  J_{-\nu}(x) = \sum_{n=0}^{\infty} \frac{(-1)^n}{n!\,\Gamma(n - \nu + 1)}\left(\frac{x}{2}\right)^{2n-\nu},
  \qquad
  y = c_1 J_{\nu} + c_2 J_{-\nu}.
  $$
  The leading term $x^{-\nu}$ makes $J_{-\nu}$ unbounded at the origin, so independence is visible at a glance.
- **$\nu = \tfrac12$ (integer gap, no logarithm).** The exponents $\pm\tfrac12$ differ by $1$, but the recurrence at the smaller root never fails: $a_0$ _and_ $a_1$ both remain free, the constant $a$ of the theorem is $0$, and the solutions are elementary,
  $$
  J_{1/2}(x) = \sqrt{\frac{2}{\pi x}}\,\sin x,
  \qquad
  J_{-1/2}(x) = \sqrt{\frac{2}{\pi x}}\,\cos x .
  $$
  The asymptotic sinusoid-over-$\sqrt{x}$ form is exact at half-integer orders.
- **$\nu = 0$ (equal roots, logarithm forced).** The differentiated-coefficient method gives the second solution
  $$
  y_2(x) = J_0(x) \ln x + \sum_{m=1}^{\infty} \frac{(-1)^{m+1} H_m}{2^{2m} (m!)^2}\,x^{2m},
  \qquad
  H_m = 1 + \frac12 + \cdots + \frac{1}{m} .
  $$
- **$\nu = 1$ (integer gap, logarithm present).** Here $J_{-1} = -J_1$ collapses onto the first solution — for any integer $m$, in fact, $J_{-m} = (-1)^m J_m$, because the gamma factors $1/(-m + n)!$ vanish for $n < m$ — so a genuinely new solution must carry the log:
  $$
  y_2(x) = -J_1(x)\ln x + \frac{1}{x}\left( 1 - \sum_{m=1}^{\infty} \frac{(-1)^m (H_m + H_{m-1})}{2^{2m}\,m!\,(m-1)!}\,x^{2m} \right).
  $$

> **Worked example.** The order $\nu = \tfrac12$ case in full. With $r = \tfrac12$ the recurrence is $a_{2m} = -a_{2m-2}/\big( 2m(2m+1) \big)$, so from $a_0 = 1$,[^boyce-bessel]
>
> $$
> a_2 = -\frac{1}{3!}, \qquad a_4 = \frac{1}{5!}, \qquad a_{2m} = \frac{(-1)^m}{(2m+1)!}.
> $$
>
> Then
>
> $$
> y_1 = x^{1/2}\sum_{m=0}^{\infty} \frac{(-1)^m}{(2m+1)!}\,x^{2m}
>     = x^{-1/2}\sum_{m=0}^{\infty} \frac{(-1)^m}{(2m+1)!}\,x^{2m+1}
>     = x^{-1/2}\sin x,
> $$
>
> and the standard normalization $a_0 = (2/\pi)^{1/2}$ gives $J_{1/2}(x) = \sqrt{2/(\pi x)}\,\sin x$. The smaller root $r = -\tfrac12$ runs without obstruction and yields $J_{-1/2}(x) = \sqrt{2/(\pi x)}\,\cos x$, so no logarithm appears despite the integer gap.

In practice these raw $y_2$'s are replaced by a standard choice. The **Bessel function of the second kind** $Y_\nu$ is the combination

$$
Y_{\nu}(x) = \frac{J_{\nu}(x)\cos \nu\pi - J_{-\nu}(x)}{\sin \nu\pi}
\qquad (\nu \notin \mathbb{Z}),
\qquad
Y_m = \lim_{\nu \to m} Y_{\nu},
$$

chosen so that the limit at integer orders exists and so that the large-$x$ behavior pairs with $J_\nu$: $Y_\nu(x) \approx \sqrt{2/(\pi x)}\,\sin(x - \pi/4 - \nu\pi/2)$. For every $\nu \ge 0$,

$$
y = c_1 J_{\nu}(x) + c_2 Y_{\nu}(x)
$$

is the general solution.[^simmons-second] The two behave oppositely at the origin: $J_\nu$ is bounded there ($J_0(0) = 1$, $J_\nu(0) = 0$ for $\nu > 0$) while $Y_\nu$ always blows up ($Y_0$ like $(2/\pi)\ln x$, $Y_1$ like $-2/(\pi x)$). Any physical problem that requires a solution bounded on the axis of symmetry therefore discards $Y_\nu$ outright, exactly as boundedness selected solutions in the Frobenius examples.

Many equations are Bessel's equation in disguise. The change of variables $z = ax^b$, $w = y x^c$ maps Bessel's equation in $(z, w)$ onto a two-parameter family in $(x, y)$; a notable member is Airy's equation $y'' + xy = 0$, whose general solution is $y = x^{1/2}\big[ c_1 J_{1/3}\big( \tfrac23 x^{3/2} \big) + c_2 J_{-1/3}\big( \tfrac23 x^{3/2} \big) \big]$.[^simmons-props]

### The drumhead

Separating variables in the wave equation for a circular membrane of radius $1$ (the [PDE module](/differential-equations/pdes-fourier-bvp/heat-wave-laplace-equations) does this in full) reduces the radial factor to Bessel's equation of order zero in $\lambda r$. Boundedness at the center forces the solution $J_0(\lambda r)$, and clamping the rim forces $J_0(\lambda) = 0$: the admissible frequencies $\lambda_1 < \lambda_2 < \cdots$ are precisely the zeros of $J_0$. The $k$-th radial mode vibrates in $k$ annular zones separated by motionless **nodal circles** at the radii where $J_0(\lambda_k r) = 0$. Because the zeros are not integer multiples of one another ($5.520/2.405 \approx 2.295$), the drum's overtones are inharmonic, unlike a string's.[^simmons-membrane]

$$
% caption: Third radial mode of a circular drumhead. The radial cross-section is
% $J_0(\lambda_3 r)$ with $\lambda_3 = 8.654$; nodal circles sit where $J_0$ vanishes,
% at radii $2.405/8.654 = 0.28$ and $5.520/8.654 = 0.64$, and the rim is clamped.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\fill[acc!12] (0,0) circle (2.1);
\fill[white] (0,0) circle (1.34);
\fill[acc!12] (0,0) circle (0.59);
\draw[acc, very thick] (0,0) circle (2.1);
\draw[acc, thick, dashed] (0,0) circle (1.34);
\draw[acc, thick, dashed] (0,0) circle (0.59);
\node at (0,0) {up};
\node at (0,0.97) {down};
\node at (0,1.72) {up};
\node[anchor=east, text=black, align=right] at (-2.55,1.55) {nodal circles:\\adjacent zones move\\in opposite phase};
\draw[->, black] (-2.45,1.35) .. controls (-1.85,1.15) .. (-1.22,0.68);
\node[anchor=west, text=black, align=left] at (2.55,1.55) {rim clamped:\\$J_0(8.654) = 0$};
\draw[->, black] (2.5,1.3) .. controls (2.2,1.2) .. (1.78,1.32);
\begin{scope}[yshift=-3.35cm]
\draw[->, black] (-2.35,0) -- (2.6,0) node[below] {$r$};
\draw[black] (0,0.05) -- (0,-0.05) node[below] {$0$};
\draw[black] (2.1,0.05) -- (2.1,-0.05) node[below] {$1$};
\draw[acc, thick] plot[smooth] coordinates {(0.00,0.55) (0.11,0.52) (0.21,0.45) (0.32,0.34) (0.42,0.21) (0.53,0.07) (0.63,-0.05) (0.73,-0.15) (0.84,-0.21) (0.95,-0.22) (1.05,-0.20) (1.16,-0.14) (1.26,-0.06) (1.37,0.02) (1.47,0.09) (1.58,0.14) (1.68,0.16) (1.78,0.16) (1.89,0.12) (1.99,0.06) (2.10,0.00)};
\draw[acc, thick] plot[smooth] coordinates {(0.00,0.55) (-0.11,0.52) (-0.21,0.45) (-0.32,0.34) (-0.42,0.21) (-0.53,0.07) (-0.63,-0.05) (-0.73,-0.15) (-0.84,-0.21) (-0.95,-0.22) (-1.05,-0.20) (-1.16,-0.14) (-1.26,-0.06) (-1.37,0.02) (-1.47,0.09) (-1.58,0.14) (-1.68,0.16) (-1.78,0.16) (-1.89,0.12) (-1.99,0.06) (-2.10,0.00)};
\node[acc, anchor=west, align=left] at (2.6,0.55) {radial section\\$J_0(8.654\,r)$};
\end{scope}
\draw[black, dashed] (0.59,-2.05) -- (0.59,-3.41);
\draw[black, dashed] (1.34,-2.05) -- (1.34,-3.41);
\end{tikzpicture}
$$

## Legendre's equation and its polynomials

> **Definition (Legendre's equation).** For a nonnegative integer $n$,
> $$
> (1 - x^2)\,y'' - 2x\,y' + n(n+1)\,y = 0 .
> $$
> The points $x = \pm 1$ are regular singular points; every other point is ordinary.

The equation appears when Laplace's equation is separated in spherical coordinates: the substitution $x = \cos\varphi$ turns the polar-angle factor into Legendre's equation, and the physical interval $0 \le \varphi \le \pi$ becomes $-1 \le x \le 1$ with the singular points sitting at the poles of the sphere.[^boyce-legendre] As with the drumhead, the physical problem requires boundedness, this time at $x = 1$.

The [ordinary-point method](/differential-equations/series-solutions/power-series-ordinary-points) at $x = 0$ gives two series with radius $1$; useful, but the bounded solutions are found faster through the hypergeometric form. Substituting $t = \tfrac12(1 - x)$ maps $x = 1$ to $t = 0$ and produces the hypergeometric equation with $a = -n$, $b = n + 1$, $c = 1$. The exponents at $t = 0$ are both $0$ (since $1 - c = 0$), so this is the equal-roots case: one solution is analytic and the second carries $\ln t$, hence is unbounded. Up to constant multiples there is exactly one solution bounded near $x = 1$, and because $a = -n$ the hypergeometric series terminates: it is a polynomial.[^simmons-legendre]

> **Definition (Legendre polynomial).** The $n$-th **Legendre polynomial** is the bounded solution normalized by $P_n(1) = 1$:
> $$
> P_n(x) = F\!\left( -n,\; n+1,\; 1,\; \frac{1-x}{2} \right)
> = \frac{1}{2^n} \sum_{k=0}^{\lfloor n/2 \rfloor} \frac{(-1)^k\,(2n - 2k)!}{k!\,(n-k)!\,(n - 2k)!}\;x^{n - 2k} .
> $$
> Equivalently, by **Rodrigues's formula**,
> $$
> P_n(x) = \frac{1}{2^n\,n!} \frac{\d^n}{\d x^n}\big( x^2 - 1 \big)^n .
> $$

The first few, from either formula:

$$
P_0 = 1, \quad
P_1 = x, \quad
P_2 = \tfrac12(3x^2 - 1), \quad
P_3 = \tfrac12(5x^3 - 3x), \quad
P_4 = \tfrac18(35x^4 - 30x^2 + 3).
$$

Each $P_n$ has degree exactly $n$, contains only even or only odd powers (so $P_n(-x) = (-1)^n P_n(x)$), satisfies $P_n(1) = 1$ and $P_n(-1) = (-1)^n$, and has all $n$ of its zeros real, simple, and inside $(-1, 1)$.

> **Worked example.** Generate $P_2$ and $P_3$ from the power-series recurrence. Substituting $y = \sum a_m x^m$ into $(1 - x^2)y'' - 2xy' + n(n+1)y = 0$ and collecting the coefficient of $x^m$ gives
>
> $$
> a_{m+2} = \frac{m(m+1) - n(n+1)}{(m+2)(m+1)}\,a_m .
> $$
>
> At $m = n$ the numerator vanishes, so the series terminates and the surviving parity is a degree-$n$ polynomial. For $n = 2$, the even chain has $a_2 = \frac{0 - 6}{2}\,a_0 = -3a_0$ and $a_4 = 0$, giving $a_0(1 - 3x^2)$; the normalization $P_2(1) = 1$ fixes $a_0 = -\tfrac12$, so $P_2 = \tfrac12(3x^2 - 1)$. For $n = 3$, the odd chain has $a_3 = \frac{2 - 12}{6}\,a_1 = -\tfrac53 a_1$, giving $a_1\big(x - \tfrac53 x^3\big)$; $P_3(1) = 1$ fixes $a_1 = -\tfrac32$, so $P_3 = \tfrac12(5x^3 - 3x)$.

$$
% caption: The first four Legendre polynomials on $[-1, 1]$. Degree-$n$ has $n$ interior
% zeros and alternates parity; every curve passes through $(1, 1)$, the normalization
% that fixes the constant multiple.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\definecolor{warm}{HTML}{B0483A}
\draw[->, black] (-3.0,0) -- (3.3,0) node[below] {$x$};
\draw[->, black] (0,-2.0) -- (0,2.15) node[left] {$y$};
\draw[black] (2.6,0.05) -- (2.6,-0.05) node[below] {$1$};
\draw[black] (-2.6,0.05) -- (-2.6,-0.05) node[below] {$-1$};
\draw[black] (0.05,1.6) -- (-0.05,1.6);
\node[text=black, anchor=north east] at (-0.08,1.52) {$1$};
\draw[black, thick] (-2.6,1.6) -- (2.6,1.6);
\draw[black, thick] (-2.6,-1.6) -- (2.6,1.6);
\draw[acc, very thick] plot[smooth] coordinates {(-2.60,1.60) (-2.34,1.14) (-2.08,0.74) (-1.82,0.38) (-1.56,0.06) (-1.30,-0.20) (-1.04,-0.42) (-0.78,-0.58) (-0.52,-0.70) (-0.26,-0.78) (0.00,-0.80) (0.26,-0.78) (0.52,-0.70) (0.78,-0.58) (1.04,-0.42) (1.30,-0.20) (1.56,0.06) (1.82,0.38) (2.08,0.74) (2.34,1.14) (2.60,1.60)};
\draw[warm, thick] plot[smooth] coordinates {(-2.60,-1.60) (-2.34,-0.76) (-2.08,-0.13) (-1.82,0.31) (-1.56,0.58) (-1.30,0.70) (-1.04,0.70) (-0.78,0.61) (-0.52,0.45) (-0.26,0.24) (0.00,0.00) (0.26,-0.24) (0.52,-0.45) (0.78,-0.61) (1.04,-0.70) (1.30,-0.70) (1.56,-0.58) (1.82,-0.31) (2.08,0.13) (2.34,0.76) (2.60,1.60)};
\fill[black] (2.6,1.6) circle (2.4pt);
\node[anchor=west, text=black] at (2.72,1.86) {$(1,1)$};
\node[anchor=west, text=black] at (-2.55,1.85) {$P_0$};
\node[anchor=north west, text=black] at (-2.5,-1.72) {$P_1$};
\node[acc, anchor=north] at (0,-0.94) {$P_2$};
\node[warm, anchor=south] at (-1.1,0.78) {$P_3$};
\end{tikzpicture}
$$

### Orthogonality and Legendre series

> **Theorem (Orthogonality of Legendre polynomials).**
> $$
> \int_{-1}^{1} P_n(x)\,P_m(x)\,\d x =
> \begin{cases}
> 0 & n \neq m, \\[2pt]
> \dfrac{2}{2n + 1} & n = m .
> \end{cases}
> $$

The $n \neq m$ case follows from the equation itself: writing Legendre's equation in the self-adjoint form $\big( (1 - x^2)\,y' \big)' = -n(n+1)\,y$, multiplying the equations for $P_n$ and $P_m$ crosswise, integrating by parts, and subtracting kills everything but $\big[ n(n+1) - m(m+1) \big] \int P_n P_m\,\d x$. The normalization $2/(2n+1)$ comes from Rodrigues's formula and repeated integration by parts.[^simmons-orth]

> **Worked example.** Check both parts of the theorem on $P_2 = \tfrac12(3x^2 - 1)$. Orthogonality to $P_1 = x$ is immediate by parity: $\int_{-1}^1 x \cdot \tfrac12(3x^2 - 1)\,\d x = 0$, the integrand being odd. For the normalization,
>
> $$
> \int_{-1}^1 P_2^2\,\d x
> = \frac14\int_{-1}^1 (9x^4 - 6x^2 + 1)\,\d x
> = \frac14\left( 9 \cdot \frac{2}{5} - 6 \cdot \frac{2}{3} + 2 \right)
> = \frac14 \cdot \frac{8}{5}
> = \frac{2}{5},
> $$
>
> matching $2/(2n+1)$ at $n = 2$.

Orthogonality converts function expansion into coefficient extraction, exactly as with Fourier series. Since $P_0, \ldots, P_n$ span the polynomials of degree $\le n$, any reasonable $f$ on $[-1, 1]$ has a **Legendre series**

$$
f(x) = \sum_{n=0}^{\infty} a_n P_n(x),
\qquad
a_n = \frac{2n + 1}{2} \int_{-1}^{1} f(x)\,P_n(x)\,\d x,
$$

the coefficient formula obtained by multiplying by $P_m$ and integrating term by term. For polynomials the series terminates: $x^2 = \tfrac13 P_0 + \tfrac23 P_2$, $x^3 = \tfrac35 P_1 + \tfrac25 P_3$, and so on.

The same mechanism runs on the Bessel side: with respect to the weight $x$, the rescaled functions $J_0(\lambda_j x)$ for successive zeros $\lambda_j$ are orthogonal on $[0, 1]$, giving Fourier–Bessel expansions for radially symmetric data. Both families are instances of one theorem — Legendre's and Bessel's equations are [Sturm–Liouville problems](/differential-equations/pdes-fourier-bvp/sturm-liouville), and real eigenvalues with orthogonal eigenfunctions are the general property, not a coincidence of these two examples.

## The special functions, side by side

| | Bessel ($\nu = 0$) | Legendre | Hypergeometric |
| --- | --- | --- | --- |
| Equation | $x^2y'' + xy' + x^2y = 0$ | $(1-x^2)y'' - 2xy' + n(n+1)y = 0$ | $x(1-x)y'' + [c - (a{+}b{+}1)x]y' - aby = 0$ |
| Geometry | cylinder, disk | sphere | archetype (three regular singular points) |
| Singular points | $0$, $\infty$ | $\pm 1$, $\infty$ | $0$, $1$, $\infty$ |
| Exponents | $0, 0$ | $0, 0$ at $x = 1$ | $0,\,1{-}c$ at $0$; $0,\,c{-}a{-}b$ at $1$ |
| Bounded solution | $J_0(x)$ | $P_n(x)$ | $F(a,b,c,x)$ |
| Unbounded partner | $Y_0 \sim \tfrac{2}{\pi}\ln x$ | $Q_n \sim \ln(1-x)$ | $x^{1-c}F(\cdot)$ or log case |
| Orthogonality | $\int_0^1 x\,J_0(\lambda_i x) J_0(\lambda_j x)\,\d x = 0$ | $\int_{-1}^1 P_n P_m\,\d x = 0$ | — |
| Expansion | Fourier–Bessel series | Legendre series | — |

One pattern runs through both families: a second-order equation from physics, a regular singular point imposed by the geometry, boundedness there selecting a one-dimensional family of solutions, a boundary condition quantizing the eigenvalues, and orthogonality assembling arbitrary data from the eigenfunctions. The same pattern organizes [separation of variables](/differential-equations/pdes-fourier-bvp/heat-wave-laplace-equations) and [Sturm–Liouville theory](/differential-equations/pdes-fourier-bvp/sturm-liouville) for the classical PDEs.

[^simmons-gamma]: **Simmons**, _Differential Equations with Applications and Historical Notes_, §46 — the gamma function: definition, functional equation, $\Gamma(n+1) = n!$, extension to negative arguments, $\Gamma(1/2) = \sqrt{\pi}$ (Problems 3–4), and the convention $1/(-m)! = 0$.
[^boyce-bessel]: **Boyce**, _Elementary Differential Equations and Boundary Value Problems_, §5.7 — Bessel's equation of orders $0$, $\tfrac12$, $1$ as the three cases of Theorem 5.6.1: the series computations, $H_m$ coefficients, and the logarithmic second solutions; **Simmons**, §46.
[^boyce-asymp]: **Boyce**, §5.7 — the $y'' + y = 0$ comparison for large $x$ and the asymptotic forms of $J_0$ and $Y_0$; **Simmons**, §46, equation (21) via $u = \sqrt{x}\,y$.
[^simmons-props]: **Simmons**, §46 Problems 1–2, 7–8 — the derivative identities, interlacing of zeros via Rolle's theorem, and the change of variables solving Airy's equation with $J_{\pm 1/3}$.
[^simmons-second]: **Simmons**, §46 — $J_{-p}$, the collapse $J_{-m} = (-1)^m J_m$ at integer orders, the Weber definition of $Y_p$ and its limit at integer $p$, and the matched large-$x$ behavior of $J_p$ and $Y_p$; **Boyce**, §5.7, equations (11)–(13) and (33).
[^simmons-membrane]: **Simmons**, Ch. 8, Appendix B — Bessel functions and the vibrating membrane: the radial reduction to $J_0$, the frequency condition $J_0(\lambda) = 0$, and the nodal circles of the higher modes.
[^boyce-legendre]: **Boyce**, §5.3, Problems 17–23 — the Legendre equation at the ordinary point $0$, polynomial termination for integer parameter, Rodrigues's formula, and the spherical-coordinates origin ($x = \cos\varphi$) of the equation.
[^simmons-legendre]: **Simmons**, §44 — Legendre polynomials via the hypergeometric equation at $t = \tfrac12(1-x)$: the equal-exponent log analysis, boundedness selecting $F(-n, n+1, 1, t)$, the explicit coefficient formula, and Rodrigues's formula.
[^simmons-orth]: **Simmons**, §45 — orthogonality by Rodrigues's formula and integration by parts, the $2/(2n+1)$ normalization, and Legendre series with the $(2n+1)/2$ coefficient rule; **Boyce**, §5.3, Problems 22–23.