Bessel Functions

Bessel's equation — the cylindrical wave equation in disguise

Bessel's equation of order n is the second-order linear ODE:

x² y'' + x y' + (x² − n²) y = 0

It arises whenever the Helmholtz or wave equation ∇²u + k²u = 0 is written in cylindrical coordinates and separated. Substitute u(r, θ, z) = R(r)·Θ(θ)·Z(z) and the radial equation, after scaling r → x/k, is exactly Bessel's equation. The order n is fixed by the angular dependence Θ(θ) = e^{inθ}.

Two independent solutions are needed. The first, J_n(x), the Bessel function of the first kind, is defined by the series:

J_n(x) = ∑_{k=0}^∞ ((−1)^k / (k! · (n+k)!)) · (x/2)^(n + 2k)

This converges for all real x and is well-behaved at the origin: J_0(0) = 1, and J_n(0) = 0 for n ≥ 1. The second solution, Y_n(x), the Bessel function of the second kind (also called Weber's function), diverges like ln x at the origin for n = 0 and like x⁻ⁿ for n ≥ 1. Any solution of Bessel's equation is a linear combination y(x) = A·J_n(x) + B·Y_n(x).

Shapes of J_n — damped oscillations

For large x, both J_n and Y_n behave like √(2/πx) times a cosine or sine:

J_n(x) ~ √(2/(πx)) · cos(x − nπ/2 − π/4)
Y_n(x) ~ √(2/(πx)) · sin(x − nπ/2 − π/4)

So they look like sinusoids whose amplitude slowly decays as 1/√x. Unlike true sinusoids the zeros are not equally spaced — they cluster slightly nearer the origin and become asymptotically uniform with spacing π. The first three zeros of J_0(x) are x ≈ 2.4048, 5.5201, 8.6537; for J_1(x), 3.8317, 7.0156, 10.1735.

The drumhead — where the first zero of J_0 = 2.4048 comes from

Consider a circular drumhead of radius R, clamped at the rim. Its vertical displacement u(r, θ, t) obeys the 2D wave equation ∂²u/∂t² = c²∇²u with u = 0 at r = R. Separating variables and assuming axial symmetry (no θ dependence) yields radial equation that becomes Bessel's equation of order 0. Solutions finite at the centre must be R(r) = J_0(kr).

The clamped boundary condition R(R) = 0 forces J_0(kR) = 0, so kR must be a zero of J_0. The frequencies are:

f_k = α_k · c / (2π R), where α_k is the k-th zero of J_0
α_1 = 2.4048, α_2 = 5.5201, α_3 = 8.6537, ...

The lowest mode (fundamental) corresponds to α_1 = 2.4048 — the entire drumhead moves up and down together with a single antinode at the centre. The second axisymmetric mode (α_2) has a node ring partway out. The non-axisymmetric modes use J_n with n ≥ 1 and produce angular nodal lines as well. Compare a string: a string of length L has modes at frequencies nπc/L for n = 1, 2, 3, … — exact integer multiples, giving harmonic overtones. A drumhead's overtones are at α_k/α_1 ≈ 1, 2.295, 3.598, … — irrational ratios, which is why a drum sounds inharmonic rather than pitched like a violin.

FM sidebands — the Jacobi–Anger expansion

A pure-tone frequency-modulated signal carries a carrier of frequency ω_c modulated by a sinusoid of frequency ω_m and modulation index β:

s(t) = cos(ω_c t + β sin(ω_m t))

The Jacobi–Anger identity expands this as an infinite sum of pure cosines:

s(t) = ∑_{n=−∞}^∞ J_n(β) · cos((ω_c + n ω_m) t)

Each term is a sinusoid at frequency ω_c + n·ω_m — the n-th sideband. The amplitude of the n-th sideband is exactly J_n(β). For β = 1 the dominant components are J_0(1) ≈ 0.7652, J_±1(1) ≈ 0.4401, J_±2(1) ≈ 0.1149, and the rest are negligible. For β = 2.4048 the carrier amplitude J_0(β) = 0 — engineers calibrate FM modulators using exactly this point (the "carrier null" method). Carson's rule estimates the FM bandwidth as 2(β + 1)·f_m by counting sidebands where |J_n(β)| is above the noise floor.

Orthogonality and the Fourier–Bessel series

Bessel functions of the same order form an orthogonal family on [0, R] with weight x:

∫₀^R J_n(α_k x/R) · J_n(α_j x/R) · x dx = 0  if k ≠ j
                                          = (R²/2) · J_{n+1}(α_k)²  if k = j

where α_k is the k-th positive zero of J_n. Any reasonable function f(r) on [0, R] (with f(R) = 0) can be expanded as a Fourier–Bessel series:

f(r) = ∑_{k=1}^∞ c_k · J_n(α_k r / R)
c_k = (2 / (R² · J_{n+1}(α_k)²)) · ∫₀^R f(r) · J_n(α_k r / R) · r dr

This is the cylindrical analogue of a Fourier sine series. It is the standard tool for solving the wave equation, heat equation, or Laplace equation on a disk or cylinder. The coefficients c_k decay rapidly for smooth f and the truncated series converges in L² and pointwise for continuous f.

Worked example — heat flow in a cylinder

A long cylinder of radius R has initial temperature T₀ throughout and is suddenly cooled with surface temperature 0. The temperature T(r, t) obeys ∂T/∂t = κ ∇²T with T(R, t) = 0 and T(r, 0) = T₀. Separating yields:

T(r, t) = ∑_{k=1}^∞ c_k · J_0(α_k r / R) · exp(−κ α_k² t / R²)

The Fourier–Bessel coefficients of the constant T₀ are c_k = 2T₀/(α_k · J_1(α_k)). The longest-lived mode decays as exp(−κ · 2.4048² · t / R²) ≈ exp(−5.78 κ t / R²), so the thermal time constant is R²/(5.78 κ) — not R²/κ as a crude estimate would suggest. The factor 5.78 = α₁² comes entirely from the geometry of the first zero of J_0.

JavaScript — computing J_n(x)

// Series for J_n(x) — accurate for small to moderate x
function besselJ(n, x, terms = 50) {
  let sum = 0;
  let term = Math.pow(x / 2, n);
  // term currently equals (x/2)^n / (0! · n!) · (sign)
  let nFact = 1;
  for (let i = 2; i <= n; i++) nFact *= i;
  term = Math.pow(x / 2, n) / nFact;
  let kFact = 1;
  for (let k = 0; k < terms; k++) {
    sum += term;
    kFact *= (k + 1);
    // multiply term by −(x/2)² / ((k+1)(n+k+1))
    term *= -(x * x / 4) / ((k + 1) * (n + k + 1));
  }
  return sum;
}

console.log(besselJ(0, 0));       // 1
console.log(besselJ(0, 2.4048));  // ≈ 0 (first zero)
console.log(besselJ(1, 3.8317));  // ≈ 0 (first zero of J_1)
console.log(besselJ(0, 5));       // ≈ -0.1776

// For large x, use the asymptotic form instead:
function besselJAsymptotic(n, x) {
  return Math.sqrt(2 / (Math.PI * x)) *
         Math.cos(x - n * Math.PI / 2 - Math.PI / 4);
}

// Newton iteration to find a zero of J_0 near a guess
function findJ0Zero(guess) {
  let x = guess;
  for (let i = 0; i < 20; i++) {
    // J_0'(x) = -J_1(x)
    x = x + besselJ(0, x) / besselJ(1, x);
  }
  return x;
}
console.log(findJ0Zero(2.5));  // ≈ 2.4048
console.log(findJ0Zero(5.5));  // ≈ 5.5201

The wider Bessel family

Function	Equation	Origin behaviour	Use case
J_n(x)	x²y'' + xy' + (x² − n²)y = 0	Regular: J_0(0) = 1	Drumheads, FM sidebands, finite-domain modes
Y_n(x)	Same	Singular: ~ ln x or x⁻ⁿ	Annular regions, exterior cylindrical scattering
I_n(x)	x²y'' + xy' − (x² + n²)y = 0	Regular, exponentially growing	Modified Bessel — heat in long fins, statistics
K_n(x)	Same modified eqn	Singular, exponentially decaying	Modified — screened Coulomb, diffusion in unbounded media
H_n^{(1,2)}(x)	= J_n ± i Y_n	Singular	Hankel functions — outgoing/incoming cylindrical waves
j_n(x), y_n(x)	Spherical: extra 1/x factor	j_n regular, y_n singular	3D quantum scattering, spherical harmonics

The Hankel functions H_n^{(1)} = J_n + iY_n and H_n^{(2)} = J_n − iY_n are the cylindrical counterparts of complex exponentials e^{±ikr} — they encode outgoing and incoming radial waves and are the natural basis for scattering and radiation problems.

Where Bessel functions actually appear

• Vibrating drumheads, timpani, gongs. The radial modes are J_0(α_k r/R) with frequencies α_k·c/(2πR). The inharmonic ratios α_k/α_1 explain why drums don't sound pitched the same way strings do.
• Optical diffraction — the Airy pattern. Light through a circular aperture produces concentric rings whose intensity is (2 J_1(x)/x)². The first dark ring sits at x = 3.8317 (first zero of J_1), setting the diffraction-limited resolution of telescopes and microscopes.
• FM radio and synthesisers. Sideband amplitudes are J_n(β). Yamaha's DX7 synthesiser, the dominant FM synthesiser of the 1980s, produced complex timbres by directly exploiting J_n(β) curves as β was modulated by an envelope.
• Coaxial cables and circular waveguides. Transverse electric/magnetic modes have field profiles J_m(k_⊥ r) cos(mθ) inside the cable. Cutoff frequencies are set by zeros of J_m or J_m'.
• Heat flow in cylinders. Cooling of rods, heating of pipes, thermal time constant R²/(α₁² κ) using the first zero of J_0.
• Quantum mechanics — particle in a cylindrical box. Radial wavefunctions are J_m(k r); energy levels set by zeros of J_m. Spherical Bessel functions j_l do the same in spherical wells.
• Statistics and signal processing. The Rice distribution (envelope of narrowband noise) and the noncentral chi-square distribution use modified Bessel I_n(x).

How Bessel functions compare to sines

	Sine / cosine	Bessel J_n
Solves	y'' + y = 0 (Cartesian)	x²y'' + xy' + (x² − n²)y = 0 (cylindrical)
Amplitude	Constant	Decays as 1/√x for large x
Zero spacing	Exactly π	Approaches π for large x; first few zeros slightly closer
Domain	Real line	x ≥ 0 (radial coordinate)
Orthogonality weight	1 (Lebesgue measure)	x (cylindrical area element)
Asymptotic form	cos(x), sin(x)	√(2/πx) · cos(x − nπ/2 − π/4)
Discoverer / date	Greek antiquity	D. Bernoulli (1732), Bessel (1824)

Common mistakes

• Using Y_n where the domain contains the origin. Y_n blows up at x = 0. On a solid drumhead (which includes the centre), the Y_n coefficient must be zero. Only on an annular geometry — where the origin is excluded — can Y_n contribute.
• Reading zeros of J_n off as multiples of π. They are not. The first three zeros of J_0 are 2.4048, 5.5201, 8.6537 — differences ≈ 3.115, 3.134. They approach π asymptotically but the first gap is noticeably larger than π.
• Treating drumhead modes as harmonic. The frequency ratios α_k/α_1 are irrational. A drum is not pitched in the harmonic-overtone sense; timpani heads are deliberately tensioned to "tune" the lowest few modes by fixing α_k/α_1 to musical intervals.
• Confusing J_n with the modified I_n. J_n oscillates; I_n grows exponentially. They satisfy different equations (sign of the x² term flips). Wrong choice gives unbounded solutions where physics demands bounded ones.
• Forgetting the weight x in the orthogonality integral. The Fourier–Bessel basis is orthogonal under ∫ f(x) g(x) x dx, not ∫ f g dx. Drop the weight and the coefficients come out wrong by a factor proportional to 1/x.
• Using the series for large x. The defining series for J_n(x) loses precision dramatically for x > 20 or so due to alternating-sign cancellation. Switch to the asymptotic form or downward Miller recurrence.