Diffusion on Curved Spaces

The animation above shows a sealed, rigid container whose walls curve inward and outward; inside, a swarm of particles begins concentrated near a single point and diffuses outward. The container is a closed 3-manifold (a deformed sphere), and its curvature is not decorative. It determines how the swarm spreads. Near the poles of the container, where the wall curves sharply, the particles are reflected at steeper angles than on the flatter equatorial band. The gold-to-blue colour encodes the heat kernel value $K_t(\mathbf{x}, \mathbf{x}_0) = \exp(-|\mathbf{x}-\mathbf{x}_0|^2/4t)$ evaluated at each particle's position: gold is high concentration, deep blue is low.

That formula is the heat kernel in flat space $\mathbb{R}^3$. On a curved manifold it is only an approximation; the true kernel acquires corrections at every order in $t$ whose coefficients are curvature invariants. The first such correction, which appears at order $t$, is the scalar curvature $R(\mathbf{x})$.

This post derives that result from first principles. We begin with the classical heat equation on $\mathbb{R}^n$ and its Gaussian fundamental solution (§1). We then build the Riemannian machinery (metric tensor, covariant gradient, divergence, and Laplace–Beltrami operator) from the ground up (§2). The heat equation on a Riemannian manifold $(M,g)$ has its own fundamental solution, derived through the parametrix construction (§3). We close with the spectral representation of the heat kernel and Weyl's law (§4).

Sign convention. The Laplace–Beltrami operator is defined as $\Delta_M = \operatorname{div}_g \circ \operatorname{grad}_g$, which is negative semi-definite on functions. The heat equation is $\partial_t u = \Delta_M u$. The competing probabilistic convention replaces this by $\partial_t u = \tfrac{1}{2}\Delta_M u$, matching the generator of Brownian motion; the two conventions differ only by the substitution $t \mapsto 2t$ throughout and produce kernels related by $K_t^{\mathrm{prob}}(\mathbf{x},\mathbf{y}) = K_{t/2}(\mathbf{x},\mathbf{y})$.

1. The Heat Equation on $\mathbb{R}^n$

1.1. Derivation from conservation

Let $u \colon \mathbb{R}^n \times [0,\infty) \to \mathbb{R}$ be the concentration of a diffusing substance. Fick's first law [1] states that diffusive flux opposes the concentration gradient:

Conservation of mass over any smooth bounded region $\Omega \subset \mathbb{R}^n$ requires

where $\hat{\mathbf{n}}$ is the outward unit normal to $\partial\Omega$. Applying the divergence theorem to the right side of [eq:conservation] and substituting [eq:fick]:

Since $\Omega$ is arbitrary and all integrands are continuous, we conclude:

1.2. Solution via the Fourier transform

For $f \in L^1(\mathbb{R}^n)$ the Fourier transform and its inverse are

Differentiation in $\mathbf{x}$ becomes multiplication in $\boldsymbol{\xi}$:

so applying [eq:fourier-deriv] twice in each coordinate and summing:

Taking the Fourier transform of [eq:heat-flat] in $\mathbf{x}$ and applying [eq:fourier-laplacian]:

This is a first-order linear ODE in $t$ for each fixed $\boldsymbol{\xi}$, with solution

By the convolution theorem, this product in frequency space corresponds to convolution in physical space:

1.3. The Gaussian heat kernel

The inverse Fourier transform of $e^{-4\pi^2|\boldsymbol{\xi}|^2 t}$ factors over coordinates, reducing to a one-dimensional integral. Completing the square:

where $\eta = \xi - ix/(4\pi t)$ and the contour shift is justified by Cauchy's theorem. Multiplying over all $n$ coordinates gives

By [eq:conv-soln], the solution to [eq:heat-flat] is $u(\mathbf{x},t) = (K_t * u_0)(\mathbf{x})$, or written symmetrically:

The visualization below shows $K_t(x) = (4\pi t)^{-1/2} e^{-x^2/4t}$ in one dimension for several values of $t$: the bell flattens and widens as $\sqrt{t}$.

Proof. (1) is immediate from the definition.

Proof of (2). Substitute $\mathbf{z} = (\mathbf{y}-\mathbf{x})/\sqrt{4t}$, so $\mathrm{d}\mathbf{y} = (4t)^{n/2}\,\mathrm{d}\mathbf{z}$. Then

where the last step uses the standard Gaussian integral $\int_{\mathbb{R}^n} e^{-|\mathbf{z}|^2}\,\mathrm{d}\mathbf{z} = \pi^{n/2}$, which factors as a product of $n$ one-dimensional integrals $\int_{-\infty}^\infty e^{-z_i^2}\,\mathrm{d}z_i = \sqrt{\pi}$ [2].

(3) is immediate from $|\mathbf{x}-\mathbf{y}|^2 = |\mathbf{y}-\mathbf{x}|^2$.

Proof of (4). Write $I = \int_{\mathbb{R}^n} K_s(\mathbf{x},\mathbf{z})\,K_t(\mathbf{z},\mathbf{y})\,\mathrm{d}\mathbf{z}$. The prefactors contribute $(4\pi s)^{-n/2}(4\pi t)^{-n/2}$. The combined exponent is

Expanding and collecting powers of $\mathbf{z}$:

The coefficient of $|\mathbf{z}|^2$ is $-(s+t)/(4st)$, and the optimal point in $\mathbf{z}$ is

Completing the square, the combined exponent equals

where the residual $-|\mathbf{x}-\mathbf{y}|^2/4(s+t)$ is independent of $\mathbf{z}$ (verified by expanding both sides and cancelling). The Gaussian integral over $\mathbf{z}$ now factors:

Combining the prefactors and residual:

(5) follows from (2) and the standard theory of approximate identities [2, Ch. 5]. □

The semigroup property has a direct visual interpretation. Hold the total time $T = s + t$ fixed and vary the split $s \mapsto T - s$. The blue kernel $K_s$ and the gold kernel $K_t = K_{T-s}$ change shape continuously, yet their convolution—the muted silver curve $K_T$—is constant throughout.

2. The Riemannian Setup

2.1. Metric tensor and volume form

Write $g = \det[g_{ij}]$ for the determinant of the metric matrix. The Riemannian volume form, the natural measure on $(M,g)$, is

This is independent of the chart: under a coordinate change the Jacobian determinant from the transformation of $\mathrm{d}x^1 \wedge \cdots \wedge \mathrm{d}x^n$ exactly cancels the factor from $\sqrt{g}$. The $L^2$ inner product on functions is

2.2. Gradient and divergence

Proof of the coordinate formula. We need components $X^i$ satisfying $g_{ij} X^i V^j = (\partial_j f) V^j$ for all vectors $\mathbf{V}$. Setting $\mathbf{V} = \partial_k$ gives $g_{ik} X^i = \partial_k f$, whence $X^i = g^{ij} \partial_j f$. □

Proof of the coordinate formula. By Cartan's formula $\mathcal{L}_{\mathbf{X}} \omega = \iota_{\mathbf{X}} \mathrm{d}\omega + \mathrm{d}(\iota_{\mathbf{X}} \omega)$ applied to the closed top-degree form $\mathrm{d}V_g$ (for which $\mathrm{d}(\mathrm{d}V_g) = 0$), we get $\mathcal{L}_{\mathbf{X}}(\mathrm{d}V_g) = \mathrm{d}(\iota_{\mathbf{X}} \mathrm{d}V_g)$. In coordinates, expanding $\iota_{\mathbf{X}}(\mathrm{d}V_g)$ and taking the exterior derivative yields $\mathrm{d}(\iota_{\mathbf{X}} \mathrm{d}V_g) = \tfrac{1}{\sqrt{g}}\partial_i(\sqrt{g}\, X^i)\, \mathrm{d}V_g$. Comparing with [eq:div-def-eq] gives [eq:div-coords] [4, Ch. 16]. □

2.3. The Laplace–Beltrami operator

On $\mathbb{R}^n$ with the Euclidean metric, $g_{ij} = \delta_{ij}$, $g = 1$, and [eq:lbo-coords] reduces to the standard Laplacian $\Delta f = \sum_i \partial_i^2 f$.

Proof. Differentiating under the integral and using the symmetry of $g^{ij}$:

It remains to show that the last integral equals $-\int_M (\Delta_M f)\varphi\,\mathrm{d}V_g$. Apply the Leibniz rule for the Riemannian divergence to the vector field $\mathbf{X} = \varphi\,\operatorname{grad}_g f$:

Integrating over $M$ and applying the Riemannian divergence theorem (boundary term vanishes because $\varphi \in C_c^\infty(M)$):

Rearranging yields $\int_M g(\operatorname{grad}_g f,\operatorname{grad}_g\varphi)\,\mathrm{d}V_g = -\int_M (\Delta_M f)\varphi\,\mathrm{d}V_g$, which is exactly [eq:first-variation] [3, Ch. 2]. □

2.4. Self-adjointness: Green's identity

Proof. Apply the Leibniz rule $\operatorname{div}_g(f\mathbf{X}) = g(\operatorname{grad}_g f,\mathbf{X}) + f\operatorname{div}_g\mathbf{X}$ to $\mathbf{X} = \operatorname{grad}_g h$:

Integrating over $M$ and applying the divergence theorem (no boundary since $M$ is closed):

Rearranging:

Applying the same argument with $f$ and $h$ exchanged (the gradient inner product is symmetric):

The two left-hand sides are therefore equal, giving all three expressions in [eq:green] [5, Ch. 1]. □

Taking $f = h$ confirms $\Delta_M$ is negative semi-definite: $\langle \Delta_M f, f\rangle = -\int_M |\operatorname{grad}_g f|_g^2 \mathrm{d}V_g \leq 0$.

2.5. Examples

The round 2-sphere. In polar coordinates $(\theta, \phi)$ on $S^2$:

The eigenfunctions of $-\Delta_{S^2}$ are the spherical harmonics $Y_\ell^m$ with eigenvalues $\lambda_\ell = \ell(\ell+1)$, multiplicity $2\ell+1$ [5, App. B].

The flat torus. On $T^n = \mathbb{R}^n/\mathbb{Z}^n$ with $g_{ij} = \delta_{ij}$, the eigenfunctions are the Fourier modes $e^{2\pi i \mathbf{k} \cdot \mathbf{x}}$ for $\mathbf{k} \in \mathbb{Z}^n$, eigenvalues $-4\pi^2|\mathbf{k}|^2$, and the heat kernel is

which converges to the flat Gaussian [eq:flat-kernel] as $t \to 0^+$ (lattice images are exponentially suppressed).

3. Diffusion on a Curved Manifold

3.1. The heat equation on $(M,g)$

For compact $M$ without boundary, a smooth solution for $t > 0$ follows from the spectral theory of §4 via the heat semigroup $u = e^{t\Delta_M} u_0$. For complete non-compact $M$, existence and uniqueness were established by Dodziuk [6] and Li–Yau [7].

3.2. The heat kernel

3.3. Geodesic distance as the natural substitute

On $(\mathbb{R}^n, \delta)$ the kernel [eq:flat-kernel] depends on $x$ and $y$ only through the Euclidean distance $|\mathbf{x}-\mathbf{y}|$. On a Riemannian manifold the canonical substitute is the geodesic distance $d_g(x,y)$. The heat kernel is not simply $(4\pi t)^{-n/2} e^{-d_g(x,y)^2/4t}$, since curvature introduces corrections at every order in $t$, but this is the leading term. The parametrix ansatz [8] is

where the transport coefficients $u_k$ are smooth functions determined by substituting [eq:parametrix] into the heat equation and matching powers of $t$.

3.4. The transport equations and the short-time expansion

We derive the coefficients by inserting [eq:parametrix] into $(\partial_t - \Delta_M^x)H_t = 0$ and reading off each power of $t$.

Setup. Fix $y \in M$ and work in normal coordinates centred at $y$ (provided by the exponential map $\exp_y$). In these coordinates [3, Ch. 5]:

• $g_{ij}(y) = \delta_{ij}$ and $\partial_k g_{ij}(y) = 0$,
• $d_g(x,y)^2 = |x|^2$ to second order,
• the volume density $\Theta(x,y) := \sqrt{g}(x)$ satisfies $\Theta(x,y) = 1 - \tfrac{1}{6}R_{ij}(y)x^ix^j + O(|x|^3)$.

Set $\phi(x,y) = d_g(x,y)^2/4$, so the exponential factor is $e^{-\phi/t}$. Write the ansatz as $H_t = (4\pi t)^{-n/2} e^{-\phi/t}\Phi$ where $\Phi = \sum_k u_k t^k$.

Applying the heat operator. Using the product formula for $\Delta_M(e^\psi f)$ with $\psi = -\phi/t$:

We evaluate two key quantities. The eikonal identity for geodesic distance:

which holds because $\operatorname{grad}_g d_g = \hat{\mathbf{r}}$ is the unit radial vector and $\phi = d_g^2/4$ gives $|{\rm grad}_g \phi|^2 = d_g^2|\operatorname{grad}_g d_g|^2/4 = d_g^2/4 = \phi$. The Laplacian of $\phi$: writing $\operatorname{grad}_g\phi = \tfrac{1}{2}d_g\hat{\mathbf{r}}$ and applying the Leibniz rule,

More precisely, since $\operatorname{div}_g\hat{\mathbf{r}} = \partial_{d_g}\!\log\Theta + (n-1)/d_g$:

Assembling [eq:product-laplacian], subtracting from $\partial_t H_t$, and using [eq:laplacian-phi]:

Collecting terms at order $t^{k-1}$ and using $g(\operatorname{grad}_g\phi, \operatorname{grad}_g\,\cdot\,) = \tfrac{d_g}{2}\partial_{d_g}$:

(with $u_{-1} \equiv 0$). Multiplying through by $2d_g^{2k-1}\Theta$ shows that the left side equals $\partial_{d_g}(d_g^{2k}\Theta u_k)/d_g^{k-1}\Theta^{1/2}$, yielding the $k$-th transport ODE:

Coefficient $u_0$. Set $k = 0$ in [eq:transport-ODE]:

so $\Theta^{1/2} u_0 = c(y)$. The normalisation condition $H_t(x,y) \to \delta_y(x)$ as $t \to 0^+$ requires the leading coefficient to equal 1 at $x = y$; since $\Theta(y,y) = 1$, we get $c(y) = 1$ and

In particular, $u_0(y,y) = 1$, recovering the flat-space leading term.

Coefficient $u_1$ on the diagonal. Integrate [eq:transport-ODE] for $k=1$ from $0$ to $d_g$:

As $x \to y$ both sides vanish; applying L'Hôpital:

It remains to compute $\Delta_M(\Theta^{-1/2})$ at $x = y$. In normal coordinates, $\Theta = \sqrt{g}$ has the expansion [8, Ch. 4]:

Since $\Delta_M(x^i x^j)|_{x=0} = 2\delta^{ij}$, applying $\Delta_M$ to [eq:theta-expand] at $x = y$ and using $R_{ii} = R(y)$ (contraction of the Ricci tensor):

Substituting [eq:delta-theta] into [eq:u1-lhopital]:

The short-time diagonal expansion. Combining [eq:u0] and [eq:u1-diag] in [eq:parametrix] evaluated at $x = y$ [9]:

The higher coefficients are universal polynomials in the curvature tensor and its covariant derivatives. The next term is [10]:

3.5. Convergence of the parametrix

The parametrix [eq:parametrix] is not yet the true heat kernel; it satisfies $(\partial_t - \Delta_M^x)H_t = O(t^{N-n/2})$ but is not exactly zero on the right. To construct the true kernel, one writes $K_t = H_t + \text{correction}$ where the correction solves a Volterra integral equation involving $H_t$. For compact $M$, the correction is smooth and $O(t^{N-n/2+1})$, confirming that [eq:heat-expansion] is asymptotically exact [11, Ch. 4].

4. The Spectral Perspective

4.1. Spectral theorem for compact manifolds

The proof uses compactness of the Sobolev embedding $H^1(M) \hookrightarrow L^2(M)$ (Rellich–Kondrachov), which forces the resolvent of $-\Delta_M$ to be compact. The spectral theorem for compact self-adjoint operators then gives the discrete spectrum; elliptic regularity implies smoothness of eigenfunctions [12, Ch. 8].

4.2. Heat kernel as a spectral series

The heat semigroup $e^{t\Delta_M}$ acts on $L^2(M)$ by

which converges in $L^2$ for all $t \geq 0$ and in $C^\infty$ for $t > 0$ (exponential decay of $e^{-\lambda_k t}$ dominates polynomial growth of $\lambda_k$). The integral kernel of $e^{t\Delta_M}$ is

confirming that [eq:kernel-rep] with this kernel reproduces [eq:semigroup-def].

On the interval $[0,\pi]$ with Dirichlet boundary conditions, $\phi_n(x) = \sin(nx)$ and $\lambda_n = n^2$. The figure below animates the solution $u(x,t) = \sum_{n \in \{1,3,5,7\}} a_n \sin(nx)\, e^{-n^2 t}$: faded traces are individual modes, bright gold is the sum. Mode $n = 7$ (eigenvalue $49$) decays 49 times faster than mode $n = 1$.

4.3. Weyl's law

Proof via the heat trace. Define the heat trace

The two expressions for $Z(t)$ are equal: substituting the spectral series [eq:spectral-kernel] into the integral gives $\sum_k e^{-\lambda_k t}\int_M \phi_k^2\,\mathrm{d}V_g = \sum_k e^{-\lambda_k t}$ by orthonormality.

Step 1: heat trace asymptotics. Integrating the diagonal expansion [eq:heat-expansion] over $M$: