The Heat (Diffusion) Equation

Take any conserved stuff — heat, dissolved dye, the probability of finding a jittering particle — and let it flow down its own gradient. Two ideas, stacked, give one equation. The first is bookkeeping: whatever leaves a small region must show up as flux through its walls, the continuity law $u_t+\nabla\!\cdot\mathbf{J}=0$ . The second is a constitutive guess, older than either Fourier or Fick: the flux points downhill and is proportional to the steepness, $\mathbf{J}=-\kappa\nabla u$ . Put them together and the flux drops out, leaving $u_t=\kappa\nabla^2 u$ — the heat, or diffusion, equation. Here $u$ is a temperature, a concentration, or a probability density; $\kappa>0$ is the diffusivity, measured in $\mathrm{m^2\,s^{-1}}$ , and its very units already whisper the punchline: length squared per time, so distance spreads as the square root of time. In the classification of second-order PDEs this is the canonical parabolic case — its discriminant vanishes — and that single algebraic fact fixes its whole personality. It runs one way in time, needs exactly one initial condition, smooths its data instantly, and cannot be un-run. It sits between the elliptic Laplace equation (pure steady states) and the hyperbolic wave equation (reversible, finite-speed, two initial conditions). Two solution strategies dominate below: on the open line the Fourier transform diagonalizes $\nabla^2$ and hands us a Gaussian heat kernel — the fundamental solution; on a bounded domain separation of variables turns the PDE into a Sturm–Liouville eigenproblem whose modes decay exponentially. Both tell one story — fine structure dies first, the whole relaxing toward the harmonic steady state the boundaries demand.

Key takeaways

The heat or diffusion equation $u_t=\kappa\nabla^2 u$ is the prototypical parabolic PDE: it describes the irreversible relaxation of a temperature, concentration, or probability density toward equilibrium.
It is first order in time and second order in space, so one initial condition plus boundary data make it well-posed forward in time; running it backward is ill-posed.
On the open line the fundamental solution is the Gaussian heat kernel, whose width spreads as $\sqrt{4\kappa t}$ — diffusion length grows like the square root of time, never linearly.
On finite domains, separation of variables gives eigenfunction expansions whose modes decay as $e^{-\kappa\lambda_n t}$ , so high-spatial-frequency wiggles are erased first.
Solutions turn instantly smooth (infinite signal speed), conserve total quantity under insulated walls, and obey a maximum principle — no new interior hot spots appear.
One equation unifies Fourier's law of conduction, Fick's law of diffusion, and the Fokker–Planck dynamics of a free Brownian particle.
Key idea: A conserved quantity carried by a gradient-driven flux relaxes by smoothing — $u_t=\kappa\nabla^2 u$ damps each spatial Fourier mode at rate $\kappa k^2$ , so fine structure dies fastest and the diffusion length grows as $\sqrt{4\kappa t}$ .

Where this lives today

The heat equation refuses to stay in the furnace. In deep learning it has become the engine of the most powerful generative models yet built: score-based diffusion models (Sohl-Dickstein 2015, Ho 2020) run a diffusion process forward to turn an image into pure Gaussian noise, then learn to reverse it — sidestepping the ill-posed backward problem by training a network to estimate the score at each noise level rather than integrating the equation blindly. In quantitative finance the Black–Scholes equation for an option price is, after a change of variables, exactly the heat equation, which is why option pricing inherits its Gaussian solutions. Geophysicists use its skin-depth solution to read paleoclimate from borehole profiles and to date the cooling of the lithosphere; Kelvin famously (and wrongly, for lack of radioactivity) used it to bound the age of the Earth. In engineering it governs the thermal management of battery packs and CPUs, where the $L^2/\kappa$ timescale decides whether a hotspot triggers runaway, and in medicine the bioheat equation — diffusion plus a perfusion sink — guides thermal ablation of tumors. The same parabolic smoothing reappears as Gaussian blur and scale-space in imaging and as edge-preserving anisotropic-diffusion denoising.

The physical picture

Picture temperature — or concentration — as a crowd of blind random walkers. Each takes unbiased steps, so no walker has a preferred direction and the mean position never moves; but the crowd's spread grows relentlessly. Where walkers are packed, more wander out than wander in, so density falls; where they are sparse, more arrive than leave. That imbalance is exactly Fourier's and Fick's law — flux proportional to the local gradient — and its net effect is to average every point toward the mean of its neighbours. A sharp spike cannot survive this averaging: it collapses instantly into a smooth bump that keeps flattening.

Figure

A point of heat released at one spot spreads as a Gaussian of width σ = √(2κt): the peak falls as 1/√t and the tails widen as √t while the enclosed area stays pinned at 1, so a sharp spike becomes a broad smooth bump and the finest structure disappears first

The cleanest way to see all of this at once is to expand any profile into spatial Fourier modes $e^{ikx}$ . Substituting into $u_t=\kappa\nabla^2 u$ makes each mode obey $\dot a_k=-\kappa k^2 a_k$ , so its amplitude decays as $e^{-\kappa k^2 t}$ . The damping rate scales with $k^2$ : a wiggle twice as fine dies four times as fast. That one dispersion relation contains everything — the smoothing (short wavelengths vanish first), the irreversibility (every mode only decays, so no feature can ever be born), the maximum principle (the solution stays trapped between its initial and boundary extremes), and the $\sqrt t$ spreading law. The long-time survivor is the slowest mode, the harmonic field fixed by the boundaries — the solution of Laplace's equation. A concrete, everyday fingerprint of the same physics is how a yearly temperature swing soaks into the ground.

Figure

The seasonal temperature wave in soil (κ ≈ 5×10⁻⁷ m²/s, one-year period): the annual swing decays as e^{−x/d} with skin depth d = √(2κ/ω) ≈ 2.24 m, and lags a full half-year by depth πd ≈ 7.0 m — where the ground is coldest while the surface bakes

The core equations

The equation is assembled from two ingredients that recur across all of continuum physics: a conservation law and a flux law. Everything downstream — kernels, series, decay estimates — is a way of solving it.

\frac{\partial u}{\partial t}=\kappa\,\nabla^2 u \qquad (1)

where: $u(\mathbf{x},t)$ — temperature, concentration, or probability density; $\kappa>0$ — diffusivity, units $\mathrm{m^2\,s^{-1}}$ ; $\nabla^2=\partial_x^2+\partial_y^2+\partial_z^2$ — the Laplacian; $t$ — time. Read it as: the rate of change at a point equals the diffusivity times the curvature of the profile there — where $u$ is concave down (a local peak) it cools, where concave up (a dip) it warms, which is exactly the smoothing seen above. The engineering form keeps the material constants separate, which matters when they are not uniform, and it identifies which combination actually sets the timescale.

\rho c\,\frac{\partial T}{\partial t}=\nabla\!\cdot\!\big(k\,\nabla T\big)+q \qquad (2)

where: $T$ — temperature (K); $\rho$ — mass density; $c$ — specific heat; $k$ — thermal conductivity; $q$ — volumetric heat source. For uniform $k,\rho,c$ and $q=0$ this collapses to equation (1) with the thermal diffusivity $\kappa=k/(\rho c)$ — the ratio, not $k$ alone, is what governs how fast a shape relaxes (for copper $\kappa\approx1.16\times10^{-4}\ \mathrm{m^2\,s^{-1}}$ ). On the whole line the response to a unit point source, the fundamental solution or Green's function, is a spreading Gaussian — the single most important object attached to this equation.

G(x,t)=\frac{1}{\sqrt{4\pi\kappa t}}\,\exp\!\left(-\frac{x^2}{4\kappa t}\right),\qquad t>0 \qquad (3)

where: $G$ — temperature from a unit amount of heat deposited at $x=0,\ t=0$ ; its variance is $\sigma^2=2\kappa t$ , so the width $\sigma=\sqrt{2\kappa t}$ and the full diffusion length $\ell=\sqrt{4\kappa t}$ . As $t\to0^+$ it becomes $\delta(x)$ , and $\int_{-\infty}^{\infty}G\,\mathrm dx=1$ for every $t$ — total quantity is conserved while the shape flattens. Because the equation is linear and translation-invariant, the solution for arbitrary initial data is just the kernel smeared over that data — a convolution.

u(x,t)=\int_{-\infty}^{\infty}G(x-x',t)\,u(x',0)\,\mathrm{d}x' \qquad (4)

where: $u(x',0)$ — the initial profile; $G(x-x',t)$ — the translation-invariant kernel of equation (3). Each point of the initial data blossoms into its own Gaussian and the results are added — superposition doing all the work. On a bounded interval the open-line kernel is replaced by a discrete eigenfunction expansion; for a rod $[0,L]$ with both ends held at zero the eigenfunctions are sines and every mode decays on its own clock.

u(x,t)=\sum_{n=1}^{\infty}b_n\,\sin\!\Big(\frac{n\pi x}{L}\Big)\,e^{-\kappa(n\pi/L)^2 t} \qquad (5)

where: $L$ — rod length; $b_n$ — Fourier sine coefficients of $u(x,0)$ ; $\lambda_n=(n\pi/L)^2$ — the Sturm–Liouville eigenvalues. The $n$ -th mode's lifetime is $\tau_n=1/(\kappa\lambda_n)\propto n^{-2}$ , so fine structure ( $n$ large) vanishes fastest and the $n=1$ mode is the last to survive. Finally, the irreversibility can be made quantitative: a natural "energy" of the field decreases monotonically, which is the analytic face of the second law here.

\frac{\mathrm d}{\mathrm dt}\int u^2\,\mathrm d^3x=-2\kappa\int|\nabla u|^2\,\mathrm d^3x\le 0 \qquad (6)

where: $\int u^2$ — the $L^2$ "energy"; $\int|\nabla u|^2$ — the Dirichlet energy of the gradients, which is strictly positive unless $u$ is flat. The integral can only decrease (with Dirichlet or insulated walls), pinning the arrow of time and forbidding any return to a sharper past state.

Derivation AThe heat kernel from the Fourier transform

The goal is to solve the initial-value problem on the entire real line with no boundaries in the way, by trading the differential operator $\partial_x^2$ for multiplication. The Fourier transform turns the PDE into a one-line ODE for each wavenumber, which we integrate trivially in time, specialize to a point source, and invert with a single Gaussian integral. The payoff is equation (3), the fundamental solution from which every open-line result — spreading, conservation, the $\sqrt t$ law — follows by inspection.

Step 1. Define $\hat u(k,t)=\int_{-\infty}^{\infty}u(x,t)\,e^{-ikx}\,\mathrm dx$ and transform the PDE. Since each $x$ -derivative brings down a factor $ik$ , the Laplacian $\partial_x^2\to-k^2$ , and the equation becomes $\partial_t\hat u=-\kappa k^2\,\hat u$ . Why: the transform diagonalizes the constant-coefficient spatial operator, so the modes decouple and the PDE splits into one independent ODE per wavenumber.

Step 2. Integrate that ODE in time: $\hat u(k,t)=\hat u(k,0)\,e^{-\kappa k^2 t}$ . This is the dispersion relation $\omega=-i\kappa k^2$ — purely imaginary frequency, so every mode decays without oscillating, at a rate $\kappa k^2$ that grows with the square of the wavenumber. Why: this is the smoothing property in its most naked form — the shortest wavelengths are annihilated fastest.

Step 3. Specialize to the fundamental solution by taking a point source, $u(x,0)=\delta(x)$ , whose transform is $\hat u(k,0)=1$ . Then $\hat G(k,t)=e^{-\kappa k^2 t}$ . Why: a delta contains every wavenumber with equal unit weight, so its evolution exposes the pure response of the operator, uncontaminated by the shape of any particular initial profile.

Step 4. Invert with the definition $G(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-\kappa k^2 t}\,e^{ikx}\,\mathrm dk$ , a Gaussian integral in $k$ . Complete the square in the exponent, $-\kappa t\,k^2+ikx=-\kappa t\,(k-ix/2\kappa t)^2-x^2/4\kappa t$ , and use $\int e^{-\kappa t\,s^2}\mathrm ds=\sqrt{\pi/\kappa t}$ . Why: shifting the contour is legitimate because the Gaussian decays in every direction, so the standard real integral applies unchanged.

Step 5. Collect the pieces: $G(x,t)=\dfrac{1}{\sqrt{4\pi\kappa t}}\,e^{-x^2/4\kappa t}$ . A direct check confirms it — differentiating gives $\partial_t G=\kappa\,\partial_x^2 G$ , and $\int G\,\mathrm dx=1$ for every $t>0$ . Why: the prefactor $1/\sqrt{4\pi\kappa t}$ is exactly what keeps the area at unity as the exponential widens, encoding conservation of total heat. *A direct check confirms the result: substituting into the PDE gives $\partial_t G=\kappa\,\partial_x^2 G$ , integrating over all $x$ gives $\int_{-\infty}^{\infty}G\,\mathrm dx=1$ for every $t>0$ , and computing the second moment gives variance $\sigma^2=2\kappa t$ .*

Go deeper · Derivation B · 2 minSeparation of variables for a cooling rod

Derivation BSeparation of variables for a cooling rod

The goal is to solve $u_t=\kappa u_{xx}$ on the finite interval $0<x<L$ with the ends clamped at zero, $u(0,t)=u(L,t)=0$ , and an arbitrary starting profile $u(x,0)=f(x)$ . The strategy is to look for solutions that factor into a pure function of space times a pure function of time; the boundary conditions then quantize the spatial part into a Sturm–Liouville eigenproblem, and completeness of its eigenfunctions lets us match any initial data. The result, equation (5), makes the mode-by-mode decay explicit.

Step 1. Insert the product ansatz $u(x,t)=X(x)\,T(t)$ into the PDE and divide by $\kappa XT$ to get $\dfrac{T'}{\kappa T}=\dfrac{X''}{X}=-\lambda$ . The left side depends only on $t$ , the right only on $x$ , so both must equal a constant, named $-\lambda$ in anticipation of decay. Why: separation converts one PDE into two ODEs coupled solely through the shared separation constant $\lambda$ .

Step 2. Solve the spatial equation $X''+\lambda X=0$ with $X(0)=X(L)=0$ . The general solution $X=A\sin(\sqrt\lambda\,x)+B\cos(\sqrt\lambda\,x)$ must vanish at $x=0$ , forcing $B=0$ ; vanishing at $x=L$ then forces $\sqrt\lambda\,L=n\pi$ . Why: the homogeneous boundary conditions act as a filter, admitting only a discrete set of standing shapes — exactly the physics of a clamped string, here for temperature.

Step 3. Read off the spectrum $\lambda_n=(n\pi/L)^2$ for $n=1,2,3,\dots$ , with eigenfunctions $\sin(n\pi x/L)$ . These are orthogonal and complete: $\int_0^L\sin(m\pi x/L)\sin(n\pi x/L)\,\mathrm dx=\tfrac{L}{2}\delta_{mn}$ . Why: this is a self-adjoint Sturm–Liouville problem, so its eigenfunctions form a basis and any reasonable $f(x)$ can be expanded in them.

Step 4. Solve the temporal equation $T'=-\kappa\lambda_n T$ to get $T_n(t)=e^{-\kappa\lambda_n t}$ . The lifetime $\tau_n=1/(\kappa\lambda_n)=L^2/(\kappa n^2\pi^2)$ shrinks as $1/n^2$ . Why: finer spatial modes have larger $\lambda_n$ and therefore decay faster — the series form of "wiggles die first".

Step 5. Superpose and match the initial data: $u(x,t)=\sum_n b_n\sin(n\pi x/L)\,e^{-\kappa\lambda_n t}$ , with the coefficients projected out by orthogonality, $b_n=\dfrac{2}{L}\int_0^L f(x)\sin(n\pi x/L)\,\mathrm dx$ . Why: at $t=0$ the exponentials are 1, so the coefficients are simply the Fourier sine coefficients of $f$ — the initial condition is met term by term. *One can check directly that each term satisfies the PDE and both boundary conditions; the orthogonality integral evaluates to $\int_0^L\sin^2(n\pi x/L)\,\mathrm dx=L/2$ , and convergence and completeness of the expansion follow from Sturm–Liouville theory.*

Go deeper · Derivation C · 2 minThe self-similar solution and the diffusion length

Derivation CThe self-similar solution and the diffusion length

The goal is to show that the heat equation on the line has no intrinsic length or time scale of its own, so its solutions for scale-free data must collapse onto a single universal curve when distance is measured in the natural unit $\sqrt{4\kappa t}$ . Exploiting the scaling symmetry reduces the PDE to an ordinary differential equation in one similarity variable, recovers the error-function solution for a step initial condition, and reads off the $\sqrt t$ growth of the diffusion length directly from the ansatz.

Step 1. Notice the invariance: equation (1) is unchanged under $x\to ax,\ t\to a^2 t$ . This means the only dimensionless combination of $x$ and $t$ is $\eta=x/\sqrt{4\kappa t}$ , so for scale-free initial data (a step, say) the solution can depend on $x$ and $t$ only through $\eta$ : $u(x,t)=F(\eta)$ . Why: a symmetry of the equation must be a symmetry of the scale-free solution, which is what forces the collapse to one variable.

Step 2. Apply the chain rule, $u_t=-\dfrac{\eta}{2t}F'$ and $u_{xx}=\dfrac{1}{4\kappa t}F''$ , and substitute into $u_t=\kappa u_{xx}$ . The factors of $t$ cancel, leaving the ODE $F''+2\eta F'=0$ . Why: the cancellation of $t$ is the analytic signature of self-similarity — the problem genuinely has no preferred moment.

Step 3. Integrate once using the integrating factor $e^{\eta^2}$ : $F'(\eta)=C\,e^{-\eta^2}$ . Why: the derivative of the profile is itself a Gaussian — which is no accident, since $F'\propto G$ , the same kernel as Derivation A.

Step 4. Integrate again: $F(\eta)=C\int_0^\eta e^{-s^2}\mathrm ds+D=\tfrac{\sqrt\pi}{2}C\,\mathrm{erf}(\eta)+D$ . Fixing the constants for a unit step $\big(u\to1$ as $x\to-\infty$ , $u\to0$ as $x\to+\infty\big)$ gives $u(x,t)=\tfrac12\,\mathrm{erfc}\!\big(x/\sqrt{4\kappa t}\big)$ . Why: the error function is simply the running integral of the Gaussian kernel, so the smoothed step is the kernel's cumulative distribution.

Step 5. Read off the length scale $\ell(t)=\sqrt{4\kappa t}$ . All profiles at all times lie on one master curve once $x$ is measured in units of $\ell$ ; a front therefore advances not at constant speed but as $\sqrt t$ , and equivalently the time to diffuse a distance $L$ scales as $L^2/\kappa$ . Why: this square-root law is the single most useful estimate the equation offers, and it drops straight out of the ansatz without solving anything. *As a check, differentiating $F$ reproduces the kernel ( $F'\propto G$ ), and the erfc solution can be confirmed to satisfy both the PDE and the step initial data.*

Limiting cases

Case	Result	Meaning
Long time, insulated finite domain	$u\;\xrightarrow{t\to\infty}\;\bar u=\dfrac{1}{V}\displaystyle\int_V u(\mathbf{x},0)\,\mathrm d^3x$	With Neumann (zero-flux) walls total quantity is conserved and the profile relaxes to its spatial average — the lowest, flat eigenmode.
Long time, fixed boundary values	$u\;\xrightarrow{t\to\infty}\;u_\infty,\quad \nabla^2 u_\infty=0$	The steady state is harmonic: the solution of Laplace's equation matching the imposed boundary data.
Single Fourier mode	$u_k(t)=u_k(0)\,e^{-\kappa k^2 t}$	Each spatial frequency decays independently; the dispersion relation $\omega=-i\kappa k^2$ is purely imaginary — damping, never oscillation.
Short time / vanishing diffusivity	$u\approx u_0+\kappa t\,\nabla^2 u_0$	For $\kappa t\to0$ the field barely moves from its initial data; smoothing has only just begun, driven by the local curvature.
Periodic surface drive, half-space	$u\sim e^{-x/d}\cos(\omega t-x/d),\ d=\sqrt{2\kappa/\omega}$	A boundary oscillation penetrates a skin depth $d$ , decaying and phase-lagging with depth — day/night and seasonal thermal waves.

Common trapThe heat equation can be run backward to sharpen or deblur data

It cannot, and the reason is structural, not merely practical. Every mode evolves as $e^{-\kappa k^2 t}$ , so reversing time turns this into $e^{+\kappa k^2|t|}$ , a factor that blows up as $k\to\infty$ . Any infinitesimal high-frequency noise in the data — and real data always has some — is amplified without bound, so the backward solution does not depend continuously on the input: the problem is ill-posed in the sense of Hadamard. This is why you cannot un-blur a photograph by simply integrating the heat equation in reverse, why "inverse heat conduction" (inferring a past or surface temperature from later interior measurements) demands regularization, and why the equation encodes a genuine arrow of time. A companion trap: do not expect a finite propagation speed. The Gaussian kernel is nonzero for every $x$ at any $t>0$ , so a point source is instantly felt everywhere — an idealization that is harmless for heat but must be repaired (relativistic diffusion, telegrapher's equation) when signal speed matters.

Worked examples

1 · How deep does the seasons' warmth reach? (fully worked)

Soil has thermal diffusivity $\kappa\approx5\times10^{-7}\ \mathrm{m^2\,s^{-1}}$ . Treating the surface temperature as an annual cosine, find the skin depth of the seasonal wave, the depth at which the swing is cut to about 4 percent, and the temperature a wine cellar at 3 m sees if the surface swings $\pm15\ \mathrm{K}$ .

Step 1. Set the drive frequency. One year is $\tau=3.156\times10^{7}\ \mathrm s$ , so $\omega=2\pi/\tau=1.99\times10^{-7}\ \mathrm{rad\,s^{-1}}$ . The half-space solution to $u_t=\kappa u_{xx}$ with a cosine boundary is $u=\Delta T\,e^{-x/d}\cos(\omega t-x/d)$ , obtained by seeking $u\propto\mathrm{Re}\,e^{i\omega t-qx}$ with $q^2=i\omega/\kappa$ , giving $q=(1+i)/d$ .

Step 2. Compute the skin depth $d=\sqrt{2\kappa/\omega}=\sqrt{2(5\times10^{-7})/(1.99\times10^{-7})}=\sqrt{5.02}=\mathbf{2.24\ m}$ . At this depth the annual swing is attenuated by $1/e\approx0.37$ . This is the natural yardstick for how far seasons soak into the ground.

Step 3. Find the half-year-lag depth. The phase lag is $x/d$ ; a lag of $\pi$ (half a year, seasons flipped) occurs at $x=\pi d=\mathbf{7.0\ m}$ , where the amplitude is $e^{-\pi}=0.043$ — about 4.3 percent of the surface swing. At 7 m the ground is coldest exactly when the surface is hottest.

Step 4. Evaluate the cellar at $x=3\ \mathrm m$ . The attenuation is $e^{-3/2.24}=e^{-1.34}=0.26$ , so a $\pm15\ \mathrm K$ surface swing becomes $\pm15\times0.26=\pm\mathbf{3.9\ K}$ , lagged by $3/2.24=1.34\ \mathrm{rad}\approx2.6$ months. A 3 m cellar therefore stays within about $\pm4\ \mathrm K$ of the annual mean, and its "summer" arrives in autumn — the everyday engineering payoff of a $\sqrt t$ law, and the reason water mains are buried below the frost line.

2 · Cooling a copper rod, and the price of simulating it

A copper rod of length $L=0.10\ \mathrm m$ ( $\kappa=1.16\times10^{-4}\ \mathrm{m^2\,s^{-1}}$ ) is heated, then its ends are clamped to 0 °C. Estimate how long it takes to cool, and find the largest stable time step for an explicit (FTCS) simulation on $N=100$ cells.

Step 1. Identify the slowest mode. From equation (5) the longest-lived term is $n=1$ , decaying with rate $\kappa\lambda_1=\kappa(\pi/L)^2$ . Its time constant is $\tau_1=\dfrac{L^2}{\pi^2\kappa}=\dfrac{(0.10)^2}{\pi^2(1.16\times10^{-4})}=\mathbf{8.7\ s}$ . After a few $\tau_1$ the rod is essentially cold; the $n=3$ mode ( $\tau_3=\tau_1/9\approx1.0\ \mathrm s$ ) and finer wiggles are long gone.

Step 2. Set up the explicit scheme. FTCS updates $u_j^{m+1}=u_j^m+r\,(u_{j+1}^m-2u_j^m+u_{j-1}^m)$ with diffusion number $r=\kappa\,\Delta t/(\Delta x)^2$ . Von Neumann analysis inserts $u_j^m=g^m e^{ikj\Delta x}$ and gives the amplification factor $g=1-4r\sin^2(k\Delta x/2)$ ; requiring $|g|\le1$ for the worst case $\sin^2=1$ demands $r\le\tfrac12$ .

Step 3. Put in the numbers. With $N=100$ cells, $\Delta x=L/N=10^{-3}\ \mathrm m$ , so the stability bound is

\Delta t\le\frac{(\Delta x)^2}{2\kappa}=\frac{(10^{-3})^2}{2(1.16\times10^{-4})}=\mathbf{4.3\times10^{-3}\ s}.

Halving the mesh spacing quarters the allowed step — the notorious cost of explicit diffusion solvers.

Step 4. Count the work. Reaching one time constant needs $\tau_1/\Delta t\approx8.7/4.3\times10^{-3}\approx\mathbf{2000}$ steps; a full cool-down of a few $\tau_1$ costs several thousand. The stiff $(\Delta x)^2$ scaling is exactly why implicit schemes (backward-Euler, Crank–Nicolson), which are unconditionally stable, are preferred for real diffusion problems.

Physics of the portrait

The heat equation is the birthplace of two subjects at once — the physics of conduction and the mathematics of Fourier analysis — and it carries the name of the man who wrote both into being while a prefect under Napoleon.

Origins · 1 minJoseph Fourier and how this was discovered

Joseph Fourier

1768–1830, French mathematician and physicist · the analytic theory of heat and the birth of harmonic analysis

Jean-Baptiste Joseph Fourier came to the theory of heat by an unlikely road: orphaned young, trained for the priesthood, swept into the French Revolution (and briefly imprisoned), then conscripted as a scientific adviser on Napoleon's Egyptian campaign. Back in France as prefect of Isère, he worked on the propagation of heat in solids and, in 1807, submitted a memoir to the Paris Academy claiming that any temperature profile could be written as a sum of sines and cosines and evolved term by term. The committee — Lagrange, Laplace, Legendre — balked at the audacity and the shaky convergence arguments, and withheld publication, though they awarded him the Academy's prize on the subject in 1811. Undeterred, Fourier expanded the work into his 1822 masterpiece, Théorie analytique de la chaleur, which not only set down $u_t=\kappa\nabla^2u$ and solved it by separation of variables, but in the process invented the trigonometric series and integral transform that now bear his name — tools that would outgrow heat entirely and reorganize much of physics and signal processing. He also, in passing, was the first to describe what we now call the greenhouse effect. The rigor he lacked was supplied later by Dirichlet and Riemann; the physical instinct was his from the start.

History

1807

— Fourier submits his memoir on the propagation of heat to the Paris Academy

1822

— Fourier publishes Théorie analytique de la chaleur, founding harmonic analysis

1855

— Fick casts the diffusion of matter in the very same equation as heat

1905

— Einstein links diffusion to Brownian motion, deriving mean-square displacement 2Dt

1949

— Von Neumann's stability analysis fixes the time-step limit for explicit schemes

Going further

Texts. Stone & Goldbart, Mathematics for Physics, Ch. 8–11 (PDEs, Green's functions) · Arfken & Weber, Mathematical Methods for Physicists, 7th ed., Ch. 9 · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed., Ch. 20–21 · Carslaw & Jaeger, Conduction of Heat in Solids, 2nd ed. · Crank, The Mathematics of Diffusion, 2nd ed. Key self-consistency checks for the results in this article: substituting $G$ into $u_t=\kappa u_{xx}$ confirms $\int G\,\mathrm dx=1$ and $\sigma^2=2\kappa t$ ; the rod orthogonality integral evaluates to $\int_0^L\sin^2(n\pi x/L)\,\mathrm dx=L/2$ ; the FTCS stability bound $r\le1/2$ follows from von Neumann analysis; and the soil skin depth $d=\sqrt{2\kappa/\omega}=2.24\ \mathrm m$ together with the copper diffusivity $\kappa=k/\rho c\approx1.16\times10^{-4}\ \mathrm{m^2\,s^{-1}}$ are consistent with standard tabulated values.

Prerequisites. Classification of Second-Order PDEs · Boundary Conditions and Well-Posedness · Separation of Variables in Cartesian Coordinates · Sturm–Liouville Eigenfunction Expansions. Leads to. Green's Functions for PDEs · The Helmholtz Equation · Laplace's and Poisson's Equations · The Wave Equation.

Ask the tutor: Ask me to work through the cooling of any specific rod, slab, or sphere — or to show exactly how the Gaussian kernel, the sine series, and the erfc solution are three views of the same physics.

Why does the heat equation smooth data instantly — is the infinite signal speed a real problem?

The Gaussian kernel $G(x,t)$ is strictly positive for every $x$ the instant $t>0$ , so heat deposited at one point is mathematically felt everywhere immediately, however faintly. That is genuine infinite propagation speed, and it follows from the parabolic structure. For heat and diffusion it is harmless: the amount reaching a far point in a short time is exponentially tiny, $\sim e^{-x^2/4\kappa t}$ , far below any measurable threshold. It only becomes a defect when you truly care about a front's speed limit — relativistic diffusion, or fast transients in nerve fibers — where the fix is a hyperbolic correction (the telegrapher's or Cattaneo equation) that restores a finite speed while reproducing the heat equation at long times.

If molecular dynamics is time-reversible, how can the heat equation be irreversible?

Nothing is added to the microscopic laws; irreversibility is emergent and statistical. The heat equation describes the coarse-grained density of an enormous number of particles, and while any single trajectory is reversible, the overwhelming majority of microstates consistent with a smooth profile lead forward to an even smoother one. Equation (6) makes this precise: the $L^2$ energy can only decrease under diffusion. Formally, the reversal is not forbidden by the microscopic dynamics but is infinitely improbable and ill-posed at the macroscopic level — reconstructing the exact past would require impossibly precise knowledge of every high-wavenumber mode. The arrow of time here is the arrow of ignorance made quantitative.

What is the difference between thermal conductivity k and thermal diffusivity κ?

They answer different questions. Conductivity $k$ (units $\mathrm{W\,m^{-1}\,K^{-1}}$ ) sets the steady flux for a given temperature gradient: how much power flows through a wall. Diffusivity $\kappa=k/(\rho c)$ (units $\mathrm{m^2\,s^{-1}}$ ) sets the speed at which a temperature pattern rearranges itself in time. A material can conduct a lot yet respond sluggishly if it also stores a lot of heat (large $\rho c$ ). That is why the timescale to equilibrate is $L^2/\kappa$ , not $L^2/k$ : what governs transient cooling is the ratio of how readily heat moves to how much heat a lump must shed to change its temperature.

Why does the diffusion length grow as √t rather than linearly?

Because diffusion is the sum of unbiased random steps, not directed motion. In a random walk the mean displacement is zero, but the variance — the sum of the individual step variances — grows linearly with the number of steps, hence linearly with time. The typical distance is the square root of the variance, so it grows as $\sqrt t$ . The equation says the same thing: the kernel's variance is $\sigma^2=2\kappa t$ , so $\sigma\propto\sqrt t$ . A practical corollary is that diffusion is fast over small distances and hopeless over large ones — doubling the distance quadruples the time — which is why cells are small and why stirring (advection, which is linear in time) beats waiting.

How is the heat equation related to Brownian motion and the Fokker–Planck equation?

They are the same equation wearing different clothes. If $u(x,t)$ is the probability density of a particle undergoing unbiased Brownian motion with no drift, its evolution is exactly $u_t=\kappa\nabla^2u$ — the driftless Fokker–Planck (Smoluchowski) equation, with $\kappa$ the diffusion coefficient $D$ . The heat kernel is then the transition probability, and its variance gives Einstein's mean-square displacement $\langle x^2\rangle=2Dt$ in one dimension. Add a force and you add a drift term $-\nabla\!\cdot(\mathbf{v}u)$ , turning it into the full Fokker–Planck equation. So Einstein's 1905 analysis of Brownian motion and Fourier's theory of heat are two readings of one parabolic operator.

Why is deblurring an image so much harder than blurring it?

Blurring an image is running the heat equation forward — a stable convolution with a Gaussian that damps high frequencies. Deblurring means inverting that, which is the ill-posed backward heat problem: the amplitudes that were multiplied by $e^{-\kappa k^2 t}$ must now be divided by it, multiplying the finest details — and any sensor noise sitting there — by an enormous factor. Naive inversion produces an explosion of high-frequency garbage. Real deconvolution therefore never inverts directly; it regularizes, trading exact recovery for stability (Wiener filtering, Tikhonov regularization, total-variation priors), or, as in modern diffusion models, learns a stochastic reverse process instead of computing the deterministic inverse.

Is the heat equation just the Schrödinger equation with imaginary time?

Formally yes, and it is a deep and useful fact. The free Schrödinger equation $i\hbar\,\psi_t=-\tfrac{\hbar^2}{2m}\nabla^2\psi$ becomes the heat equation under the Wick rotation $t\to-i\tau$ , with the diffusivity playing the role of $\hbar/2m$ . But the physics differs sharply: the heat equation's real exponential $e^{-\kappa k^2 t}$ damps modes, giving irreversible smoothing, whereas Schrödinger's imaginary exponential $e^{-i\hbar k^2 t/2m}$ rotates them, giving reversible, oscillatory, unitary evolution that conserves probability. The imaginary-time trick is nonetheless the backbone of the Feynman–Kac formula and of quantum Monte Carlo, where imaginary-time diffusion projects onto the ground state.