Bessel Functions

Laplacian Eigenvalue Problem on a Disk

$$ \nabla^2 u = \lambda u $$

On a disk, the rotational symmetry of the geometry naturally leads to polar coordinates.

$$ \nabla^2 u(r,\theta) = \lambda u(r,\theta) $$

We solve the Laplacian eigenvalue problem on a disk of radius $ R $, with vanishing boundary conditions:

$$ u(R,\theta) = 0 $$

In polar coordinates,

$$ \nabla^2 = \frac{1}{r}\frac{\partial}{\partial r} \left( r\frac{\partial}{\partial r} \right) + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2} $$ $$ \frac{1}{r}\frac{\partial}{\partial r}\left( r\frac{\partial}{\partial r} u(r,\theta) \right) + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2} u(r,\theta) = \lambda u(r,\theta) $$

For vanishing boundary conditions, the eigenvalues of the Laplacian are negative, so $ \lambda < 0 $

Assume a separated solution:

$$ u(r,\theta) = f(r) g(\theta) $$ $$ \frac{1}{f(r)} \frac{1}{r} \frac{\partial}{\partial r} \left(rf'(r)\right) + \frac{1}{r^2}\frac{g''(\theta)}{g(\theta)} = \lambda $$

The left-hand side separates into a function of $ r $ and a function of $ \theta $, so both must equal a constant.

The angular equation is

$$ g''(\theta) = -n^2 g(\theta) \quad g_n(\theta) = e^{in\theta} \quad n \in \mathbb{Z} $$

due to periodic boundary conditions (see derivation).

The radial equation is

$$ \frac{1}{r}\frac{\partial}{\partial r} \left( r\frac{\partial f}{\partial r} \right) - \frac{n^2}{r^2}f(r) = \lambda f(r) $$

Let $ \lambda = - k^2 $

$$ \frac{1}{r}\frac{\partial}{\partial r} \left( r\frac{\partial f}{\partial r} \right) - \frac{n^2}{r^2}f(r) = - k^2 f(r) $$

Unlike the angular equation, the radial equation does not reduce to simple sinusoidal solutions. This equation can be transformed into the standard Bessel equation by introducing

$$ z = kr \quad f(r) = Z(z) $$ $$ \frac{\partial}{\partial z} = \frac{1}{k} \frac{\partial}{\partial r} $$

Substituting into the radial equation gives

$$ \frac{k}{z} k \frac{\partial}{\partial z} \left( \frac{z}{k} k \frac{\partial}{\partial z} Z(z) \right) - \frac{k^2 n^2}{z^2} Z(z) = - k^2 Z(z) $$ $$ \frac{1}{z}\frac{\partial}{\partial z} \left( z\frac{\partial Z}{\partial z} \right) - \frac{n^2}{z^2}Z(z) = - Z(z) $$

Therefore,

$$ \frac{\partial^2 Z}{\partial z^2} + \frac{1}{z}\frac{\partial Z}{\partial z} + \left( 1 - \frac{n^2}{z^2} \right)Z(z) = 0 $$

This is the Bessel equation in its standard form. Bessel functions that do not blow up at $ z = 0 $ are called Bessel functions of the first kind and are denoted with $ J_n $. Bessel functions therefore play the same role in radial geometries that sine and cosine functions play in Cartesian geometries.

Since $ f(r) = Z(kr) $, the radial solutions are

$$ f(r) = J_n(kr) $$

Unlike sine and cosine functions, Bessel functions oscillate with an amplitude and spacing that vary with radius. The first few Bessel functions of the first kind are shown below.

Just as sine and cosine functions arise naturally when solving problems with translational symmetry, Bessel functions arise naturally in problems with radial symmetry. The zeros of these functions determine the allowed standing-wave modes on a disk.

The boundary condition becomes

$$ f(R) = Z(kR) = 0 $$

Therefore, $ k $ must be a solution to the equation

$$ J_n(kR) = 0 $$

For each $ n $, we obtain a discrete set of values for $ k $ by solving $ J_n(kR) = 0 $. This equation cannot be solved in closed form. We denote the solutions by $ k_{n,m} $, where $ m $ indexes the radial nodes between $ r = 0 $ and $ r = R $ (the zeros of $ J_n(kR) $ ). The boundary condition therefore quantizes the allowed radial modes.

Bessel functions therefore play the same role in radial geometries that sine and cosine functions play in Cartesian geometries.

The corresponding eigenfunctions of $ \nabla^2 $ are

$$ u_{n,m}(r,\theta) = J_n (k_{n,m} r) e^{in\theta} $$

with radial part

$$ f_{n,m}(r) = J_n (k_{n,m} r) $$

Since the Laplacian is Hermitian under the following inner product, these eigenfunctions are orthogonal

$$ \langle u_{n,m}, u_{\tilde{n},\tilde{m}} \rangle = \int_{0}^{2\pi}\!\!\int_{0}^{R} u_{n,m}^* \, u_{\tilde{n},\tilde{m}} \, r\,dr\,d\theta = \delta_{n\tilde{n}}\,\delta_{m\tilde{m}}\, \|u_{n,m}\|^2 $$

Bessel functions are orthogonal when they have the same $ n $ but different $ m $.

Taking the real part of $ e^{in\theta} $ gives cosine modes

$$ u_{n,m}(r,\theta) = J_n (k_{n,m} r) \cos(n\theta) $$

The plots below show several eigenfunctions on the unit disk. Although the problem was formulated on a disk of radius $ R $, the radial coordinate can always be rescaled according to $ \rho = r/R $, so that the domain becomes $ 0 \leq \rho \leq 1 $. The resulting plots therefore illustrate the universal geometric structure of the modes independently of the physical size of the disk. Changing $ R $ changes the physical wavelength and eigenvalues, but not the qualitative spatial pattern of the eigenfunctions after rescaling.

The integer $ n $ controls the angular structure of the mode, while $ m $ controls the number of radial oscillations between the center and the boundary. Higher modes therefore contain increasingly fine spatial structure.

These eigenfunctions are the normal modes of vibration of a circular membrane. The surface plots below show the displacement of the membrane away from equilibrium, while the nodal lines correspond to regions that remain stationary throughout the motion.

These normal modes form a complete orthogonal basis for functions on the disk. More complicated motion can therefore be constructed by superposing many such modes with different amplitudes and frequencies.

These spatial modes can then be used to construct time-dependent motion on a disk. In the notebook below, we compute the allowed wavenumbers, visualize several Bessel-function modes, examine their orthogonality, and then construct a wave-like solution by superposing modes with different frequencies.

The final part of the notebook shows how an initially localized velocity profile can be expanded in this modal basis. As more modes are included, the solution captures increasingly fine spatial structure, giving a visual demonstration of how eigenfunction expansions represent motion on a circular domain.

Open Notebook Environment

Heat Equation in Cylindrical Coordinates

Consider the heat equation in cylindrical coordinates for the temperature distribution inside a half-cylinder of radius $ R $ and height $ H $:

$$ \dot{T} = \alpha \nabla^2 T $$

where $ T = T(r,\theta,z,t). $ The domain is $ 0 \leq r \leq R, \quad 0 \leq \theta \leq \pi, \quad 0 \leq z \leq H, \quad t \geq 0. $

The boundary conditions are

$$ T(R,\theta,z,t) = 0 $$ $$ \partial_{\theta} \, T(r,0,z,t) = 0, \quad \partial_{\theta} \, T(r,\pi,z,t) = 0 $$ $$ \partial_z \, T(r,\theta,0,t) = 0, \quad \partial_z \, T(r,\theta,H,t) = 0 $$

with initial condition

$$ T(r,\theta,z,0) = T_i $$

An observant person might notice that the geometry, boundary conditions, and initial condition are all independent of $ \theta $ and $ z $. By symmetry, we therefore look for solutions depending only on $ r $ and $ t $. The insulating boundary conditions on the flat faces means the heat flux at all the boundaries is pointing radially for the entire solution. Such a solution can be extended symmetrically to an infinite cylinder, where the absence of any preferred angular or axial direction makes the $ \theta $ and $ z $ independence especially clear.

$$ T(r,t) = f(r) p(t) $$ $$ \frac{p'(t)}{p(t)} = \alpha \frac{\nabla^2 f(r)}{f(r)} = \alpha \lambda $$

The solution therefore separates into exponentially decaying time dependence and radial eigenmodes.

$$ p'(t) = \alpha \lambda p(t) \quad \Rightarrow \quad p(t) = e^{\alpha \lambda t} $$ $$ \nabla^2 f(r) = \lambda f(r) $$ $$ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial}{\partial r} f(r) \right) = \lambda f(r) $$

Let $ \lambda = - k^2 $

$$ \frac{1}{r}\frac{\partial}{\partial r} \left( r\frac{\partial f}{\partial r} \right) = - k^2 f(r) $$

Rescale $ r $ by k

$$ z = kr \quad \frac{\partial}{\partial z} = \frac{1}{k} \frac{\partial}{\partial r} $$

Then the radial equation becomes

$$ \frac{\partial^2 f}{\partial z^2} + \frac{1}{z}\frac{\partial f}{\partial z} + f = 0 $$

This is the Bessel equation with $ n = 0 $ (because there's no $ \theta $ dependence on $ T $). The solutions finite at $ r = 0 $ are of the form $ J_0(kr) $ where $ J_n $ are Bessel functions of the first kind

$$ f(r) = J_0(kr) $$ $$ f(R) = J_0(kR) = 0 $$

We know the eigenvalues of $ \nabla^2 $ are all negative, so $ \lambda < 0 $ (vanishing BC's $ f(R) = 0 $), and $ J_0(0) = 1 $ so $ k = 0 $ is not a solution of $ J_0(kR) = 0$. Hence we must solve $ J_0(kR) = 0$ for $ \lambda < 0 $, the solutions are labeled by a positive integer $ m $, $ k_m $

$$ \nabla^2 f_m(r) = \lambda_m f_m(r) \quad \quad f_m(r) = J_0(k_mr) \quad \lambda_m = - k_m^2 $$ $$ T(r,t) = \sum_{m=1}^{\infty} c_m f_m(r) e^{\alpha\lambda_m t} = \sum_{m=1}^{\infty} c_m J_0(k_mr) e^{\alpha\lambda_m t} $$

The operator $ r^{-1} \partial_r (r \partial_r ) $ is symmetric (this follows from $ \nabla^2 $ being symmetric in general) with the inner product

$$ \langle u, v \rangle = \int_{0}^{R} u(r)v(r)\,r\,dr \quad u(R) = v(R) = 0 $$

Therefore the functions $ f_m $ form an orthogonal basis, $ \langle f_m, f_{\tilde{m}} \rangle = \delta_{m\tilde{m}}\|f_m\|^2 $

Since

$$ T(r,0) = \sum_{m=1}^{\infty} c_m f_m(r) = T_i $$

the coefficients are

$$ c_m = \frac{\langle f_m, T_i \rangle}{\langle f_m, f_m \rangle} $$

Therefore,

$$ T(r,t) = T_i \sum_{m=1}^{\infty} \frac{ \int_{0}^{R} J_0(k_mr)\,r\,dr }{ \int_{0}^{R} J_0^2(k_mr)\,r\,dr } J_0(k_mr) e^{\alpha\lambda_m t} $$

where $ J_0 $ is a Bessel function and $ k_m $ is the $ m $-th positive solution of $ J_0(k_mR) = 0 $

Higher radial modes decay more rapidly, so the temperature distribution gradually approaches equilibrium.

See the full solution with angular and axial dependence.

See a similar problem.