The Millennium Prize Problem

Learning Objectives

By the end of this section, you will be able to:

1. State the Clay Millennium problem for the 3D incompressible Navier-Stokes equations in plain English and in formal mathematics.
2. Explain the “race” between nonlinear vortex stretching and viscous dissipation that sits at the heart of the question.
3. Derive the basic energy estimate $\tfrac{d}{dt}\,\|\mathbf{u}\|_{L^2}^2 = -2\nu\,\|\nabla\mathbf{u}\|_{L^2}^2$ and recognise why it is not enough.
4. Reproduce the Beale-Kato-Majda criterion: smoothness up to time $T$ persists iff $\int_0^T \|\boldsymbol\omega(t)\|_{L^\infty}\,dt < \infty$.
5. Distinguish classical, weak (Leray-Hopf), and suitable weak solutions, and quote the Caffarelli-Kohn-Nirenberg partial-regularity result.
6. Implement a Python integrator for the toy 1D blow-up ODE and a PyTorch 1D viscous Burgers solver, and use them to feel why 1D and 2D are tractable while 3D resists every known technique.

The Million-Dollar Question

In May 2000 the Clay Mathematics Institute announced seven open problems, each carrying a US$1 000 000 prize. One of them — the only one drawn directly from applied mathematics — is the Navier-Stokes existence and smoothness problem. Twenty-five years later it is still open.

Put as simply as possible: start with a smooth, gentle 3D velocity field and let it evolve under the Navier-Stokes equations. Will the flow remain smooth forever, or can it develop an infinite-velocity gradient — a “singularity” — at some finite future time?

The official problem (Clay, 2000): Prove that for any smooth, divergence-free initial velocity field with finite energy on $\mathbb{R}^3$, the 3D incompressible Navier-Stokes equations have a global, smooth solution — OR exhibit a smooth initial datum for which the solution breaks down in finite time.

The asymmetry is striking. Existence of a solution for short time is elementary. Uniqueness for the short-time solution is elementary. Existence of a weak solution for all time was proved by Jean Leray in 1934. The unsolved question lives in a narrow gap: does the smooth solution that exists for a moment continue to exist for all time?

A Quick Recap: The 3D Navier-Stokes Equations

For an incompressible Newtonian fluid filling a domain $\Omega \subset \mathbb{R}^3$, the velocity $\mathbf{u}(\mathbf{x},t) \in \mathbb{R}^3$ and pressure $p(\mathbf{x},t) \in \mathbb{R}$ satisfy

The five symbols, decoded:

The Race: Stretching vs Dissipation

Strip the Navier-Stokes equations down to the two terms that fight each other. The convection term $(\mathbf{u}\cdot\nabla)\mathbf{u}$ is quadratic in $\mathbf{u}$ — it amplifies what is already there. The viscous term $\nu\,\Delta\mathbf{u}$ is linear — it smooths what is already there. The Millennium question, deeply, is about which one wins as features get sharper.

Here is the heart of the difficulty, written as a back-of-the-envelope scaling argument. Suppose $\mathbf{u}$ develops a feature of size $U$ on a length scale $\ell$. Then derivatives scale like $U/\ell$, and:

The ratio is the local Reynolds number $R_\ell = U\ell/\nu$. When $R_\ell \gg 1$ convection dominates and the feature tries to grow. When $R_\ell \ll 1$ diffusion dominates and the feature relaxes. The terrifying observation is that as a vortex stretches, ℓ shrinks faster than U grows, and $R_\ell$ can — at least a priori — keep climbing.

Before tangling with vector PDEs, let's capture this race in the simplest possible scalar ODE. It is the 1D cartoon of the Millennium problem and is small enough to play with on a slider.

Interactive: The Toy Blow-up ODE

Imagine a single “vorticity unit” $\omega(t)$. Two forces act on it:

• Self-stretching: $+\omega^2$ — a tube of vorticity stretched by the flow concentrates its rotation, which makes the stretching even stronger. Pure positive feedback.
• Viscous damping: $-\nu\,\omega$ — the standard linear decay you would write for any diffusing scalar.

The ODE is

Separation of variables gives the exact solution, valid as long as the denominator stays nonzero:

When $\omega_0 > \nu$ the bracket is positive, the exponential drives the denominator to zero, and the solution explodes at the finite blow-up time

When $\omega_0 \le \nu$ the bracket is non-positive, the denominator stays bounded below by $1$, and $\omega(t) \to 0$. Two regimes, one knob — exactly the shape of the Millennium question.

What “Smoothness” Really Means

When mathematicians say a solution is smooth, they mean it has continuous partial derivatives of every order — i.e. it lies in $C^\infty$. In practice we work in function spaces that quantify how many derivatives are well-behaved and in what averaged sense.

The reason these distinctions matter: in 3D a function can have bounded energy ($L^2$) and yet be unbounded pointwise, the way $f(\mathbf{x}) = |\mathbf{x}|^{-1/2}$ on the unit ball is in $L^2$ but has a spike at the origin. The Millennium problem is whether the Navier-Stokes equations can drive an initially smooth solution into exactly that kind of singularity.

The Formal Clay Statement

Clay phrased four sub-problems. Solving any one of them — A or B — wins the prize. C and D are the corresponding statements for the periodic version.

The community is split, but the consensus tilts gently toward (A) being true — based on numerical evidence, physical intuition, and the partial-regularity theorems below. Tao's 2016 averaged-equation blow-up result is the strongest hint that (B) might be the right answer, because it shows that a close cousin of Navier-Stokes does blow up.

Energy Estimates: What We CAN Prove

Take the inner product of the momentum equation with $\mathbf{u}$ itself and integrate over the whole domain. Three of the four terms simplify beautifully:

• $\int \mathbf{u}\cdot\partial_t\mathbf{u}\,dx = \tfrac{1}{2}\tfrac{d}{dt}\int|\mathbf{u}|^2\,dx$ — the kinetic-energy time derivative.
• $\int \mathbf{u}\cdot(\mathbf{u}\cdot\nabla)\mathbf{u}\,dx = \tfrac{1}{2}\int \mathbf{u}\cdot\nabla|\mathbf{u}|^2\,dx = 0$ — the nonlinear term integrates to zero by $\nabla\cdot\mathbf{u}=0$ and integration by parts. The bad term vanishes in the energy budget!
• $\int \mathbf{u}\cdot\nabla p\,dx = 0$ — same reason, pressure drops out.
• $\int \mathbf{u}\cdot\nu\Delta\mathbf{u}\,dx = -\nu\int|\nabla\mathbf{u}|^2\,dx$ — integration by parts on the Laplacian.

The result is the famous basic energy identity:

Integrating in time gives the cornerstone a-priori bound

What this says, in words: energy can only decrease, and the integral of $\|\nabla\mathbf{u}\|_{L^2}^2$ over all time is finite. Both are wonderful — and both are spectacularly not enough.

The energy bound controls $\|\nabla\mathbf{u}\|_{L^2}$ in an integral sense over time. To rule out a singularity, we would need an instantaneous pointwise bound on $\nabla\mathbf{u}$. There is a giant gap between “finite on average” and “never infinite”.

The Vorticity Equation: Where 3D Bites

Define vorticity $\boldsymbol\omega = \nabla\times\mathbf{u}$. Take the curl of the momentum equation. Pressure drops out (curl of a gradient is zero), and you get

The new term $(\boldsymbol\omega\cdot\nabla)\mathbf{u}$ is the vortex-stretching term. It is the soul of the Millennium problem.

In 2D the vortex-stretching term $(\boldsymbol\omega\cdot\nabla)\mathbf{u}$ is identically zero (vorticity is perpendicular to the plane of the flow, gradient lies in the plane, dot product vanishes). Vorticity is then transported by the flow and dissipated by viscosity. It can never amplify. This is why 2D is solved.

In 3D the stretching term is generically nonzero, and the formal argument $d|\boldsymbol\omega|/dt \sim |\boldsymbol\omega|^2$ is exactly the toy blow-up ODE you just played with — minus the depletion effects coming from geometry. The Millennium question is whether those geometric effects are powerful enough to globally beat the quadratic feedback.

Leray's Weak Solutions (1934)

Faced with an apparent dead-end on the classical problem, Jean Leray invented a now-standard trick: drop the requirement that the equation be satisfied pointwise. Multiply by a smooth test function $\boldsymbol\varphi$ and integrate, moving derivatives onto $\boldsymbol\varphi$. Any function $\mathbf{u}$ satisfying the resulting integral identity for all test functions is called a weak solution.

Leray proved: For any $\mathbf{u}_0 \in L^2$ divergence-free, there exists a weak solution defined for all $t \in [0,\infty)$ that satisfies the energy inequality. That is a remarkable theorem — global existence of some kind of solution from any finite-energy start. Two giant caveats:

1. Regularity unknown. The Leray solution might be discontinuous, might have unbounded gradients — it just needs to satisfy the integral identity.
2. Uniqueness unknown. We do not know that the Leray solution is the only weak solution with the given initial data. (Recent work by Buckmaster-Vicol shows that, for the related super-critical class of weak solutions, uniqueness fails.)

The Millennium gap, restated: we have a global weak solution. We have a short-time classical solution. We do not know that they agree, and we do not know that the classical solution can be extended for all time.

The Beale-Kato-Majda Criterion

Beale, Kato, and Majda (1984) gave the cleanest characterisation of blow-up:

Suppose the classical solution $\mathbf{u}$ exists smoothly on $[0,T)$. Then the solution can be extended past $T$ if and only if

$\displaystyle\int_0^T \|\boldsymbol\omega(t)\|_{L^\infty}\,dt < \infty.$

In words: the only way a smooth Navier-Stokes solution can die at time T is for the supremum norm of vorticity to be non-integrable on the approach to T. This is enormously useful because it converts the abstract question “does smoothness survive?” into a concrete diagnostic that can be tracked in any DNS simulation: just plot $\|\boldsymbol\omega\|_{L^\infty}$ against time and see whether the time-integral is starting to diverge.

Caffarelli-Kohn-Nirenberg Partial Regularity

If you cannot rule out singularities everywhere, the next best thing is to bound the size of the set where they could live. That is the Caffarelli-Kohn-Nirenberg theorem (1982), the deepest known result on 3D Navier-Stokes:

CKN, 1982: The set $S \subset \mathbb{R}^3 \times (0,\infty)$ of possible singularities of a suitable weak solution has parabolic Hausdorff dimension at most $1$.

Translation: the singular set, if non-empty, is at most a curve in space-time. It cannot be a 2D sheet, cannot fill any open region. Whatever exotic object a 3D Navier-Stokes singularity is, it is extremely thin.

CKN does not say $S$ is empty. That is the prize. It does say that solving the prize is equivalent to ruling out a one-dimensional pencil of nasty points.

Interactive: The Kolmogorov Energy Cascade

The physical picture motivating the math is Andrei Kolmogorov's 1941 theory of turbulence. Energy is injected at a large length scale $L$, cascades through an “inertial range” of progressively smaller eddies without losing energy along the way, and is finally dissipated by viscosity at the Kolmogorov scale $\eta = (\nu^3/\varepsilon)^{1/4}$.

Between $L$ and $\eta$ the energy spectrum follows the celebrated $-5/3$ law:

The range $L/\eta$ grows like $\mathrm{Re}^{3/4}$. For atmospheric flows $\mathrm{Re}\sim 10^9$, so the inertial range stretches over six decades. That is the regime where a direct numerical simulation would need $\sim \mathrm{Re}^{9/4} \sim 10^{20}$ grid points — utterly out of reach. It is also the regime where the Millennium question lives.

Why 2D Is Solved (and 3D Is Not)

In 2D Navier-Stokes, vorticity $\omega = \partial_x v - \partial_y u$ is a scalar (the perpendicular component of the 3D vorticity vector). Its evolution equation is

Notice what is not there: no vortex-stretching term. The scalar equation has the structure of a convection-diffusion equation, and admits a maximum principle: the supremum of $|\omega|$ is non-increasing along the flow. So $\|\omega\|_{L^\infty}$ is bounded by its initial value forever, the BKM-style integral is automatically finite, and global smoothness follows without breaking a sweat. Leray-Ladyzhenskaya-Lions (1969) proved global regularity in 2D once and for all.

In 3D the stretching term $(\boldsymbol\omega\cdot\nabla)\mathbf{u}$ re-enters. The maximum principle for vorticity collapses. The toy ODE you played with becomes an honest possibility. Everything in 3D regularity theory boils down to trying to control this one term that 2D does not have.

Worked Example: Computing the Blow-up Time

Let's solve the toy ODE by hand for $\omega_0 = 2$, $\nu = 1$, then sanity-check against the interactive panel.

Python: Integrating the Blow-up ODE

Now let's do the same calculation numerically. We'll build a 4th-order Runge-Kutta integrator with an automatic blow-up detector and confirm the closed-form answer to four decimal places.

Read the annotation cards on the right — every numerical value printed by this code is hand-computed in the explanations, so you can verify each step yourself.

1import numpy as np
2import matplotlib.pyplot as plt
3
4
5def vorticity_rhs(w, nu):
6    """
7    Right-hand side of the toy vorticity model.
8
9        dω/dt = ω² − ν · ω
10
11    The ω² term is the 1D cartoon of *vortex stretching* — a tube
12    of vorticity stretched by the flow concentrates its rotation,
13    amplifying ω. The −ν·ω term is *viscous diffusion* spread out
14    over the unit cell — it relaxes ω back toward zero.
15    """
16    return w * w - nu * w
17
18
19def rk4_step(w, dt, nu):
20    """One classical 4th-order Runge-Kutta step on the toy ODE."""
21    k1 = vorticity_rhs(w, nu)
22    k2 = vorticity_rhs(w + 0.5 * dt * k1, nu)
23    k3 = vorticity_rhs(w + 0.5 * dt * k2, nu)
24    k4 = vorticity_rhs(w + dt * k3, nu)
25    return w + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
26
27
28def integrate(w0, nu, t_end, n_steps, cap=1e8):
29    """
30    Integrate the toy ODE and stop once ω exceeds 'cap'
31    (so we can SEE the finite-time blow-up).
32    """
33    dt = t_end / n_steps
34    ts = [0.0]
35    ws = [w0]
36    w = w0
37    for k in range(n_steps):
38        w = rk4_step(w, dt, nu)
39        t = (k + 1) * dt
40        if not np.isfinite(w) or w > cap:
41            ts.append(t)
42            ws.append(cap)
43            print(f"BLOW-UP detected at t ≈ {t:.4f}, ω ≈ {w:.3e}")
44            return np.array(ts), np.array(ws)
45        ts.append(t)
46        ws.append(w)
47    return np.array(ts), np.array(ws)
48
49
50def exact_blowup_time(w0, nu):
51    """
52    Closed-form blow-up time for dω/dt = ω² − ν·ω, ω(0) = ω0 > ν.
53
54        t*  =  (1/ν) · ln( ω0 / (ω0 − ν) )
55
56    If ω0 ≤ ν the solution is global (no blow-up); we return ∞.
57    """
58    if w0 <= nu:
59        return np.inf
60    return (1.0 / nu) * np.log(w0 / (w0 - nu))
61
62
63if __name__ == "__main__":
64    # Case A: stretching wins (ω0 = 2, ν = 1) → blow-up at t = ln 2.
65    ts_a, ws_a = integrate(w0=2.0, nu=1.0, t_end=0.69, n_steps=2000)
66    t_star_a   = exact_blowup_time(2.0, 1.0)
67    print(f"Case A   t* (formula) = {t_star_a:.4f}")
68
69    # Case B: viscosity wins (ω0 = 0.5, ν = 1) → global decay.
70    ts_b, ws_b = integrate(w0=0.5, nu=1.0, t_end=4.0, n_steps=2000)
71    print(f"Case B   ω(4) = {ws_b[-1]:.3e}  (should be ~0)")

BLOW-UP detected at t ≈ 0.6930, ω ≈ 1.20e+08
Case A   t* (formula) = 0.6931
Case B   ω(4) = 9.158e-03  (should be ~0)

PyTorch: A Differentiable Viscous Burgers Equation

The toy ODE is just one variable. Let's climb one rung of the ladder: the 1D viscous Burgers equation $\partial_t u + u\,\partial_x u = \nu\,\partial_x^2 u$ is a genuine PDE with the same competition between a quadratic nonlinearity and linear dissipation — but without the divergence-free constraint and without 3D vortex-stretching. 1D Burgers IS proven to have global smooth solutions for any $\nu > 0$.

We'll implement it in PyTorch — same code runs on CPU and GPU — and use it to watch what happens as we drop viscosity through four orders of magnitude. The maximum velocity gradient grows like $\sqrt{1/\nu}$, big but bounded; in 3D Navier-Stokes nobody knows.

1import torch
2
3
4def viscous_burgers_step(u, nu, dx, dt):
5    """
6    One explicit step of the 1D viscous Burgers equation
7
8        ∂u/∂t  +  u · ∂u/∂x   =   ν · ∂²u/∂x²
9
10    on a periodic grid. This PDE is the cleanest 1D analogue of the
11    3D Navier-Stokes equations:
12
13        the term u·∂u/∂x is the 1D inertial (nonlinear) term;
14        the term ν·∂²u/∂x² is the viscous diffusion.
15
16    Unlike 3D NS, the 1D Burgers equation IS proven to have a
17    global smooth solution for any ν > 0 and any smooth initial
18    condition. The reason is exactly the maximum principle:
19    diffusion in 1D smooths everything out before steepening
20    can run away. The Millennium question asks whether the
21    analogous statement is true in 3D.
22
23    Inputs
24        u  : tensor of shape (N,), the velocity on the grid
25        nu : float, kinematic viscosity
26        dx : grid spacing
27        dt : time step
28    Returns
29        u_new : tensor of shape (N,), one step forward in time.
30    """
31    # Periodic finite differences
32    u_left  = torch.roll(u,  1)
33    u_right = torch.roll(u, -1)
34
35    # Convective term: upwinded central difference  u · ∂u/∂x
36    # We use a central scheme because the question is about
37    # smoothness, not shocks; central is 2nd-order accurate.
38    conv = u * (u_right - u_left) / (2.0 * dx)
39
40    # Diffusive term: standard 3-point Laplacian  ν · ∂²u/∂x²
41    diff = nu * (u_right - 2.0 * u + u_left) / (dx * dx)
42
43    # Forward-Euler update.
44    return u + dt * (-conv + diff)
45
46
47def run(nu, n_x=256, t_end=0.5, cfl=0.4):
48    """
49    Drive a smooth initial bump u0(x) = sin(2πx) and watch what
50    steepens / smooths / blows up as we shrink ν.
51    """
52    L  = 1.0
53    dx = L / n_x
54    # CFL: dt ≤ cfl · dx² / (2ν)  for stability of the diffusion.
55    dt = cfl * dx * dx / max(2.0 * nu, 1e-8)
56    n_steps = int(t_end / dt)
57
58    x = torch.linspace(0.0, L, n_x + 1)[:-1]
59    u = torch.sin(2.0 * torch.pi * x)
60
61    for k in range(n_steps):
62        u = viscous_burgers_step(u, nu, dx, dt)
63
64    grad_max = ((torch.roll(u, -1) - torch.roll(u, 1)) / (2.0 * dx)).abs().max()
65    return u, grad_max.item()
66
67
68if __name__ == "__main__":
69    for nu in [1e-1, 1e-2, 1e-3, 1e-4]:
70        u, g = run(nu)
71        print(f"ν = {nu:.0e}   max |∂u/∂x| = {g:8.2f}")

Recent Progress and Open Avenues

The problem has resisted every attempt for 90+ years, but the field is far from stagnant. Some of the most influential recent directions:

Why This Matters: From Forecasts to Aircraft

It is tempting to dismiss the problem as “just a math question”. Engineers and physicists routinely simulate Navier-Stokes flows that look perfectly fine, even at Reynolds numbers far beyond anything provable. So why care?

Summary

The Millennium Prize problem for Navier-Stokes is not exotic mathematics. It is the question of whether the most successful equation in continuum mechanics is, internally, consistent with the smoothness it appears to produce.

Key Takeaways

1. The Clay statement asks for global smoothness from smooth, divergence-free initial data on $\mathbb{R}^3$ — or a finite-time singularity from such data.
2. The heart of the difficulty is the vortex-stretching term $(\boldsymbol\omega\cdot\nabla)\mathbf{u}$ in 3D, which can a priori amplify vorticity quadratically — the 1D cartoon is exactly $d\omega/dt = \omega^2 - \nu\omega$.
3. We have a basic energy identity $\tfrac{d}{dt}\|\mathbf{u}\|_{L^2}^2 + 2\nu\|\nabla\mathbf{u}\|_{L^2}^2 = 0$ but it controls only integral norms, not pointwise gradients.
4. Leray (1934) gives global weak solutions of unknown regularity and uniqueness.
5. Beale-Kato-Majda (1984): smoothness extends past $T$ iff $\int_0^T \|\boldsymbol\omega\|_{L^\infty}\,dt < \infty$.
6. CKN (1982): the singular set of a suitable weak solution has parabolic Hausdorff dimension at most $1$.
7. 2D is solved by Leray-Ladyzhenskaya-Lions; the vortex-stretching term identically vanishes and a vorticity maximum principle holds.
8. Numerical experiments — including our 1D Burgers solver — show gradients staying finite even as $\nu \to 0$. That is suggestive, not conclusive.

Coming Next: In Section 8 we leave the question of whether smooth solutions exist and learn how engineers compute the ones they need anyway — the world of computational fluid dynamics.