From Jean Leray to the millennium problem: the Navier–Stokes equations

Abstract

One of the Millennium problems of Clay Mathematics Institute from 2000 concerns the regularity of solutions to the instationary Navier–Stokes equations in three dimensions. For an official problem description, we refer to an article by Charles Fefferman (Existence and smoothness of the Navier–Stokes equation. http://www.claymath.org/sites/default/files/navierstokes.pdf) on the whole space problem as well as the periodic setting in $\mathbb {R}^3$. In this article, we will introduce the Navier–Stokes equations, describe their main mathematical problems, discuss several of the most important results, starting from 1934 with the seminal work by Jean Leray, and proceeding to very recent results on non-uniqueness and examples involving singularities. On this tour de force we will explain the open problem of regularity.

1 Introduction

1.1 The model of Navier–Stokes

The Navier–Stokes equations are the standard model to describe the flow of a viscous, incompressible fluid such as water; however, polymeric fluids, blood and oil need more sophisticated models involving further nonlinear terms. Dating back to Newton, Daniel Bernoulli and Euler the laws of conservation of mass and of linear momentum yield for the velocity field $u=(u_1, u_2, u_3)$ and the hydrostatic pressure p of a viscous and incompressible fluid the set of partial differential equations

$$\begin{aligned} u_t + u \cdot \nabla u - \nu \varDelta u + \nabla p= & {} f \end{aligned}$$

(1)

$$\begin{aligned} \text {div }u= & {} 0. \end{aligned}$$

(2)

The only nonlinear term $u \cdot \nabla u$ stands for $\sum \nolimits _{j=1}^3 u_j \frac{\partial u}{\partial x_j}$; this term naturally appears when computing the total time derivative of the velocity u(x(t), t) along the trajectory, x(t), of a particle, i.e., the acceleration

$$\begin{aligned} a(x(t),t) = \frac{\mathrm{d}}{\mathrm{d}t}\, u (x(t),t)= \sum \limits _{j=1}^3 {\dot{x}}_j(t) \frac{\partial u}{\partial x_j}+ \frac{\partial u}{\partial t}. \end{aligned}$$

Based on Newton’s law of viscosity and the so-called Stokes’ hypothesis (1845) which was more or less already formulated by A.J.C. Barré de Saint-Venant in 1843 the easiest way to model viscosity is given by the term $-\nu \varDelta u + \nabla p$ where $\nu >0$ is the constant coefficient of viscosity. For simplicity, we put the constant density $\rho$ to 1. In more complicated models of viscous fluids, the term $-\nu \varDelta u + \nabla p$ is replaced by nonlinear terms leading to non-Newtonian fluids or by either higher order or nonlocal terms with respect to u. From this point of view, the Navier–Stokes equations are the simplest model for viscous, incompressible fluids. The system (1), (2) is to be completed by an initial value $u_0$ for the velocity u at time $t=0$, say, and boundary conditions for u on the boundary $\partial \varOmega$ of a domain $\varOmega \subset \mathbb {R}^3$ unless $\varOmega =\mathbb {R}^3$ or $\varOmega$ equals the torus ${\mathbb {T}}^3 = (\mathbb {R}/(2\pi \mathbb {Z}))^3$ with periodic boundary conditions. The most common boundary condition is given by the adhesion or no-slip condition $u=0$ on $\partial \varOmega$, i.e., particles on $\partial \varOmega$ at time $t=0$ stay on the boundary for all times.

1.2 Basic observations

The Navier–Stokes system (1)–(2) together with initial-boundary value conditions define a coupled system of nonlinear partial differential equations so that—as for a nonlinear ordinary differential equation—the question on global-in-time existence and uniqueness of solutions might be nontrivial. Unfortunately, compared with the Poisson problem $-\varDelta u=f$ or the heat equation $u_t -\varDelta u=f$, there exists no maximum principle for the velocity field u. Given boundary data $u \not = 0$ on $\partial \varOmega$, we cannot expect to get the a priori estimate $\underset{x \in \varOmega }{\max }\, |u(x)| \le \underset{y \in \partial \varOmega }{\max }\, |u(y)|$ in the stationary case. Actually, considering the flow in a finite channel with a bottleneck between the inlet and outlet the velocity vector will attain its maximum somewhere in the bottleneck due to the condition $\text {div }u=0$, the incompressibility of the fluid.

Moreover, the pressure plays a hindering but nevertheless crucial role. For the stationary Stokes problem, i.e., for the system

$$\begin{aligned} -\nu \varDelta u + \nabla p=0, \quad \text {div }u=0 \text { in } \varOmega , \quad u=0 \text { on } \partial \varOmega , \end{aligned}$$

(3)

the solenoidality condition $\text {div }u =0$ can be considered as a linear constraint requiring a Lagrange multiplier so that (3) can be solved. This multiplier is given by the pressure. Furthermore, in (1)–(2), there is no time derivative for p so that with a solution p(x, t) (and u(x, t)) also $p(x,t)+q(t)$ for any function q(t) which may be even non-measurable defines a solution.

Usually the pressure will be eliminated by the so-called Helmholtz projection P which is given by the solution to a weak Poisson problem with Neumann boundary condition. The solution operator, a pseudodifferential operator of order 0, acts nonlocally in the domain $\varOmega$ so that in the Navier–Stokes system the pressure depends nonlocally on the term $u\cdot \nabla u$. The Laplacian $-\varDelta$ will be replaced by the Stokes operator $A=-P \varDelta$ which partly has properties as the Laplacian $-\varDelta$, in particular in the whole space $\varOmega =\mathbb {R}^3$ or the periodic setting $\mathbb {T}^3$; but it is much more difficult to handle in domains with boundary since the operator P maps vector fields u to solenoidal vector fields with vanishing normal component $(P u) \cdot n=0$ on $\partial \varOmega$ so that in general $-P\varDelta \ne -\varDelta P$. In other words, the lack of the maximum principle for the Stokes problem (3) compared to the Laplace equation is related to the nonlocality of the Stokes operator A.

1.3 The state of the art before 1934

There are only very few explicit solutions to the Navier–Stokes system. These examples are solutions in special geometries like an infinite tube (Hagen–Poiseuille flow) or in an infinite rotating cylinder with inner and outer cylinder (Taylor–Couette flow). The common ingredient is the fact that for these solutions the nonlinear term $u \cdot \nabla u$ vanishes. The so-called Jeffery–Hamel flow in an infinite two-dimensional wedge with in- or outflow through the vertex of the wedge admits a radially symmetric solution in the sense that ru(x), $r=|x|$, solves a nonlinear second-order ordinary differential equation. Solutions to the stationary Stokes system (3) in more general domains, either bounded or exterior domains, were found by the method of single- and double-layer hydrodynamic potentials leading to boundary integral equations and exploiting Fredholm theory. The instationary Navier–Stokes system in the whole space $\mathbb {R}^3$ was solved by Oseen [47, §7.4] using an ansatz by power series. However, the power series will converge only on a small time interval $[0,T_0]$, and Oseen remarks [47, §7.5] that a further iterative process on a sequence of consecutive intervals $[T_0, T_1]$, $[T_1, T_2]$, ...will probably stop in finite time. The main problem in his and other approaches is the fact that the authors were working in spaces of classical smooth functions for which they could not (and we cannot yet even today) find uniform a priori estimates globally in time. To be more precise, these spaces were not adapted to the underlying physics of the Navier–Stokes system.

2 Weak solutions

2.1 Construction of weak solutions

Adequate a priori estimates are given by the energy (in-)equality obtained by testing (1) with the solution u itself, i.e., taking the $L^2(\varOmega )$-scalar product of (1) with u and applying integration by parts. Since $u{ |{_{\partial \varOmega }}}=0$ and $\text {div }u=0$ in $\varOmega$, (1) yields the equation

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t} \frac{1}{2} \int \nolimits _\varOmega |u|^2\,\mathrm{d}x + \nu \int \nolimits _\varOmega |\nabla u|^2\,\mathrm{d}x + \int \nolimits _\varOmega \nabla p \cdot u\,\mathrm{d}x + \int \nolimits _\varOmega (u \cdot \nabla u) \cdot u\,\mathrm{d}x=0 \end{aligned}$$

where $\int \nolimits _\varOmega \nabla p \cdot u\,\mathrm{d}x=0$ and $\int \nolimits _\varOmega (u \cdot \nabla u) \cdot u\,\mathrm{d}x= \frac{1}{2} \int \nolimits _\varOmega u \cdot \nabla |u|^2 \,\mathrm{d}x=0$. Hence, an integration in time on [s, t] yields the energy equality

$$\begin{aligned} \frac{1}{2} \Vert u(t)\Vert _2^2 + \nu \int \nolimits _s^t \Vert \nabla u\Vert _2^2 \,\mathrm{d}\tau = \frac{1}{2} \Vert u(s)\Vert _2^2; \end{aligned}$$

(4)

here, $\Vert \cdot \Vert _2$ denotes the $L^2(\varOmega )$-norm of functions, vector fields, etc. With $s=0$ and $u(0)= u_0 \in L^2(\varOmega )$, we conclude that

$$\begin{aligned} u \in L^\infty \big (0,T;L^2(\varOmega )\big ), \quad \nabla u \in L^2 \big (0,T;L^2(\varOmega )\big ), \end{aligned}$$

i.e., the flow field u has a kinetic energy $\frac{1}{2} \Vert u(t)\Vert _2^2$ which is uniformly bounded for all times; the same holds for the dissipation energy, $\nu \int \nolimits _0^t \Vert \nabla u\Vert _2^2\,\mathrm{d}\tau$, which is caused by friction of particles close to each other, but with different velocities (strain rate tensor $(\nabla u + (\nabla u)^\top )(x,t) \not =0$). The idea to use energy estimates dates back to Leray [41]. He was the first to construct global-in-time weak solutions called solutions turbulentes.

This procedure motivates the definition of weak solutions which—for simplicity—will be stated for bounded domains $\varOmega \subset \mathbb {R}^3$ only. In this definition, $L^2_\sigma (\varOmega )=PL^2(\varOmega )$ is the set of all weakly solenoidal vector fields v in $L^2(\varOmega )$ such that the normal component, $v \cdot n$, vanishes on $\partial \varOmega$ in a generalized sense, and $C^\infty _{0,\sigma }(\varOmega )$ denotes the set of all solenoidal test vector fields in $C_0^\infty (\varOmega )$. Finally, $H_0^1(\varOmega )$ is the Sobolev space of $L^2$-vector fields v such that $v{ |{_{\partial \varOmega }}}=0$ and the weak gradient, $\nabla v$, belongs to $L^2(\varOmega )$.

Definition 1

A solenoidal vector field u in the Leray–Hopf class

$$\begin{aligned} {\mathcal {L}} {\mathcal {H}}_T = L^\infty \big (0,T;L^2_\sigma (\varOmega )\big ) \cap L^2 \big (0,T; H_0^1(\varOmega )\big ) \end{aligned}$$

(5)

is called a weak solution in the sense of Leray and Hopf of the Navier–Stokes system (1), (2) with Dirichlet boundary data $u{ |{_{\partial \varOmega }}}=0$ and initial value $u_0 \in L^2_\sigma (\varOmega )$ if for all smooth compactly supported and solenoidal (with respect to x) vector fields $\varphi$ ($\varphi \in C^\infty ([0,T); C^\infty _{0,\sigma }(\varOmega ))$)

$$\begin{aligned} \begin{aligned} -\int \nolimits _0^T \int \nolimits _\varOmega&u \varphi _t \,\mathrm{d}x \,\mathrm{d}\tau + \nu \int \nolimits _0^T \int \nolimits _\varOmega \nabla u \cdot \nabla \varphi \,\mathrm{d}x \,\mathrm{d}\tau \\&-\int \nolimits _0^T \int \nolimits _\varOmega u \otimes u \cdot \nabla \varphi \,\mathrm{d}x \,\mathrm{d}\tau = \int \nolimits _\varOmega u_0 (x) \varphi (0,x)\,\mathrm{d}x \end{aligned} \end{aligned}$$

(6)

and u(t) satisfies for all $t\in [0,T)$ the energy inequality

$$\begin{aligned} \frac{1}{2} \Vert u(t)\Vert _2^2 + \nu \int \nolimits _0^t \Vert \nabla u\Vert _2^2 \,\mathrm{d}\tau \le \frac{1}{2} \Vert u_0\Vert _2^2. \end{aligned}$$

(7)

Weak solutions of the Navier–Stokes system in the sense of Leray and Hopf can be constructed for bounded domains and all kinds of unbounded domains, even with rough boundaries. The classical procedure is by the Galerkin approximation in which the nonlinear PDE (1)–(2) is replaced by a sequence of nonlinear $(N \times N)$-ODE systems of first order in time yielding a sequence of approximate solutions $(u_N)_{N \in \mathbb {N}}$. Using the energy equality (4) for each $u_N$, the sequence $(u_N)$ is shown to be bounded in the Leray–Hopf class ${\mathcal {L}}{\mathcal {H}}_T$. Then, by functional analytic principles, $(u_N)$ possesses a weakly/weakly-$(*)$-convergent subsequence, called $(u_N)$ again, such that

$$\begin{aligned} u_N \overset{*}{\rightharpoonup } u \text { in } L^\infty \big (0,T;L^2_\sigma (\varOmega )\big ), \quad u_N \rightharpoonup u \text { in } L^2 \big (0,T;H_0^1(\varOmega )\big ). \end{aligned}$$

Unfortunately, weak convergence does not allow to pass to the limit in the energy equality (4), but, due its lower semi-continuity with respect to norms, guarantees the energy inequality (7) for the weak limit u. Actually, it is still an open problem whether any weak solution $u \in {\mathcal {L}} {\mathcal {H}}_T$ of the variational problem (6) satisfies the energy equality (4) or the energy inequality (7).

The crucial step in the passage to the limit $u_N \rightarrow u$ is the nonlinear term $u_N \cdot \nabla u_N$; for that step a compactness argument is needed to conclude from further properties of the sequence of approximate solutions $(u_N)$ that $u_N \rightarrow u$ in $L^2(0,T; L^2(\varOmega ))$. The Galerkin approximation method was exploited by Hopf [35] to construct global weak solutions of the Navier–Stokes equations in an arbitrary domain $\varOmega \subset \mathbb {R}^3$ the shape of which may even depend on time. Leray [41] accomplished this analysis for the whole space case $\varOmega =\mathbb {R}^3$ in which approximate solutions and global a priori estimates are easily obtained by the use of the heat kernel; however, his compactness argument on $\mathbb {R}^3$ was quite involved.

2.2 Fundamental properties of weak solutions in $\mathbb {R}^2$ and $\mathbb {R}^3$

Mathematically, it makes sense to discuss the Navier–Stokes system (1)–(2) also in two dimensions. Indeed, a flow field $u=(u_1, u_2, u_3)$ in a three-dimensional domain $\varOmega =\varOmega ' \times \mathbb {R}\subset \mathbb {R}^3$ such that $u_3 \equiv 0$ and $u_1, u_2,\, p$ are independent of the third variable $x_3$, is a solution of (1)–(2) in the two-dimensional domain $\varOmega ' \subset \mathbb {R}^2$.

In the following, let us compare properties of weak solutions in $\mathbb {R}^d$, $d=2$ or 3. To this aim, we introduce the Serrin number S(s, q) for the Bochner space $L^s(L^q) := L^s(0,T;L^q(\varOmega ))$:

$$\begin{aligned} u\in L^s(L^q)\quad \Rightarrow \quad S_u := S(s,q) := \frac{2}{s}+ \frac{d}{q}. \end{aligned}$$

Then, for $u\in L^\infty (L^2)$ we have $S_u=\frac{d}{2}$. Moreover, by Sobolev’s embedding theorems, $u\in L^2(H_0^1) \subset L^2(L^{2^*})$ where $2^*= \frac{2d}{d-2}$ so that $S_u=\frac{d}{2}$ as above (actually, when $d=2$, there holds only $H^1 \subset L^q$ for all $2 \le q < \infty$). Thus, the Serrin number for weak solutions equals $S_u=\frac{d}{2}$.

Returning to the energy (in-)equality let us test the Navier–Stokes solution with the weak solution u itself. The crucial term is the nonlinear one, $u \cdot \nabla u$, leading to the integral

$$\begin{aligned} \int \nolimits _0^T \int \nolimits _\varOmega (u \cdot \nabla u) \cdot u \,\mathrm{d}x \,\mathrm{d}t. \end{aligned}$$

(8)

Although it is easily seen by an integration by parts that for almost all $t \in (0,T)$ the integral $\int \nolimits _\varOmega (u \cdot \nabla u) \cdot u\,\mathrm{d}x$ is well-defined and vanishes, the double integral (8) requires in view of $\nabla u \in L^2(L^2)$ more or less that $u \in L^4(L^4)$. For $d=2$, $S(4,4)=1=\frac{d}{2}$ so that the “necessary” condition $u \in L^4(L^4)$ is satisfied and consequently the integral in (8) vanishes. The final conclusion is that weak solutions in the 2d-setting satisfy the energy equality. However, in $\mathbb {R}^3$, $S(4,4)=\frac{5}{4} < \frac{3}{2}=\frac{d}{2}$, i.e., the condition $u \in L^4(L^4)$ cannot be guaranteed for weak solutions and the validity of the energy equality remains open. On the other hand, Galdi [24] proved that the condition $u \in L^4(L^4)$ yields the validity of the energy equality even for weak solutions in the sense of distributions, i.e., u satisfies (6) in Definition 1, but a priori is not necessarily an element of the Leray–Hopf class ${\mathcal {L}}\mathcal H_T$; in this case, the condition $u\in {\mathcal {L}}{\mathcal {H}}_T$ follows automatically.

There is an intrinsic argument why Serrin’s condition on u is important. The Navier–Stokes equations have the remarkable property that with a solution u and an associated pressure p also the pair $(u_\lambda , p_\lambda )$,

$$\begin{aligned} u_\lambda (x,t)= \lambda u(\lambda x, \lambda ^2 t), \quad p_\lambda (x,t)= \lambda ^2 p(\lambda x, \lambda ^2 t), \end{aligned}$$

(9)

for any $\lambda >0$ is a solution of (1)–(2) for the same coefficient of viscosity $\nu$. This property guarantees that a “solution” controlled in a wind tunnel experiment for a scaled model is related—together with all its measurands—to a real world solution via the scaling (9). Mathematically, this property would open up the opportunity to work with $u_\lambda$ for sufficiently large or small $\lambda$ rather than u itself and exploit the smallness of some adequate norm of $u_\lambda$. However, it is easy to check that

$$\begin{aligned} \Vert u_\lambda \Vert _{L^s(0,\infty ; L^q(\mathbb {R}^3))} = \lambda ^{1-\big (\frac{2}{s}+\frac{3}{q}\big )}\Vert u\Vert _{L^s(0,\infty ; L^q(\mathbb {R}^3))} \end{aligned}$$

so that $\Vert u_\lambda \Vert _{L^s(0,\infty ; L^q(\mathbb {R}^3))} = \Vert u\Vert _{L^s(0,\infty ; L^q(\mathbb {R}^3))}$ if and only if $\frac{2}{s}+ \frac{3}{q}=1$, i.e., $S_u=1$. Therefore, the space $L^s(L^q)$ with $S(s,q)=1$ is called scale invariant. We conclude that the analysis in scale invariant function spaces yields in some sense optimal results which cannot be improved by replacing u with $u_\lambda$.

Looking at uniqueness and regularity of weak solutions, we put the nonlinear term to the right-hand side as a given external force. Then, a rigorous analysis based on Sobolev embeddings and Gagliardo–Nirenberg inequalities shows that $u \cdot \nabla u$ can be absorbed by the term $u_t -\nu \varDelta u$ on the left-hand side in the two-dimensional case, but not in the three-dimensional one. A formal check is based on Serrin numbers: Since $\nabla u \in L^2(L^2)$, we have $S_{\nabla u}= \frac{d}{2}+1$ so that with $S_u=\frac{d}{2}$ and Hölder’s inequality $S_{u\cdot \nabla u} = \frac{d}{2}+ (\frac{d}{2}+1)=d+1$. On the other hand, formally $S_{u_t}=S_{\varDelta u} = \frac{d}{2}+2$. In the two-dimensional case where $d+1=\frac{d}{2}+2$, the nonlinear term can be controlled by linear ones. Actually, a more careful analysis proves that weak solutions in $\mathbb {R}^2$ are unique and smooth in the classical sense for all times $t >0$. However, in the 3d case, where $d+1 >\frac{d}{2}+2$, the nonlinear term cannot be absorbed by $u_t-\nu \varDelta u$ and uniqueness as well as regularity remain open questions—unless the weak solution u satisfies further conditions. This topic will be discussed in the next section.

2.3 Basic results on regularity of weak solutions

Assume in the 3d case that a weak solution u lies in the so-called Serrin class, i.e.,

$$\begin{aligned} u \in L^s(0,T; L^q(\varOmega )), \quad S_u=1, \quad 2<s<\infty . \end{aligned}$$

(10)

Then, u is unique among all weak solutions with the same data satisfying the energy inequality (weak–strong uniqueness, see [51]). The same condition, (10), implies that u is a regular solution. For instance, if $\varOmega \subset \mathbb {R}^3$ is a “uniform $C^2$-domain” and $u_0$ is regular, see (20), (21) in Sect. 2.5 for details, then for each $0<T'<T$

$$\begin{aligned} \begin{aligned} u&\in L^\infty \big ([0,T']; H^1(\varOmega )\big ) \cap L^2\big ([0,T']; H^2(\varOmega )\big ),\\ u_t&\in L^2 \big ([0,T']; L^2(\varOmega )\big ); \end{aligned} \end{aligned}$$

(11)

in particular, all derivatives in the Navier–Stokes system exist in the weak sense and the equations are satisfied in the pointwise sense almost everywhere in time and space. Even $C^\infty$-regularity up to the boundary $\partial \varOmega$ can be obtained assuming that $\partial \varOmega$ is of class $C^\infty$. Precise results and proofs are found in the monograph by H. Sohr [53, Theorems V. 1.8.1 and V. 1.8.2], whereas the result on uniqueness dates back to Leray [41] and Serrin [50, 51].

On the other hand, if u is only a solution in the sense of distributions lying in Serrin’s class (10) for each interval $(\delta ,T)$, $0<\delta <T$, and converges to $u_0$ weakly in $L^2$ as $t\rightarrow 0+$, then u is actually a Leray–Hopf solution satisfying (5) in Definition 1 and even $u(t)\rightarrow u_0$ in $L^2$ as $t\rightarrow 0$, see Galdi [25]. A similar result holds in the limit case when $L^s(L^q)$ is replaced by $C^0(L^3)$.

The limit case $s=\infty$, $q=3$ in (10) posed further problems concerning uniqueness as well as regularity since $L^\infty$ is neither reflexive nor separable. The corresponding regularity problem was solved not before 2003 by Escauriaza et al. [18] for the whole space $\varOmega =\mathbb {R}^3$; the main idea is a backward uniqueness property for the heat equation (inequality) for the vorticity $\omega = \mathrm{rot}\, u$: roughly spoken, if $|\omega _t - \nu \varDelta \omega | \le M(|\omega |+|\nabla \omega |)$ on $\mathbb {R}^3 \times (-T,0)$, $\omega (t=0)=0$ on $\mathbb {R}^3$ and $\omega$ is restricted by some rough growth condition, then $\omega$ vanishes on $\mathbb {R}^3 \times (-T,0)$. Subsequent articles deal with the case of domains.

In addition to Serrin’s condition $S_u=1$, there are various other conditions on u, $\nabla u$, on specific components of either u or $\nabla u$, on the vorticity $\omega =\mathrm{curl}\,u$, the pressure p and further physical or mathematical quantities which yield regularity and classical smoothness of weak solutions. Numerous results of this type can be found in the literature under the name conditional regularity; see Sect. 3.1 below for a discussion of selected results.

In particular, the vorticity,

$$\begin{aligned} \omega =\mathrm{curl}\,u= (\partial _2 u_3-\partial _3 u_2, \partial _3 u_1 -\partial _1 u_3, \partial _1 u_2 -\partial _2 u_1) \in \mathbb {R}^3, \end{aligned}$$

has interesting mathematical and physical properties which are different in $\mathbb {R}^2$ and $\mathbb {R}^3$. For a planar flow, $d=2$, the derivatives $\partial _3$ and the component $u_3$ vanish so that only the third component $\omega _3$ of the vorticity $\omega$ may be nonzero. In this latter case, the scalar vorticity $\omega (:=\omega _3)$ solves the scalar vorticity transport equation

$$\begin{aligned} \omega _t -\nu \varDelta \omega + u \cdot \nabla \omega =0, \end{aligned}$$

(12)

a parabolic equation for $\omega$ which admits control of $\omega$ by the maximum and minimum principle since in each term $\omega$ appears under a differential operator. However, note that in a domain $\varOmega \subset \mathbb {R}^2$ the Dirichlet boundary condition $u=0$ on $\partial \varOmega$ does not yield enough information for $\omega$ on the boundary to allow for a rigorous application of the maximum principle. In the three-dimensional case, the vorticity transport equation reads

$$\begin{aligned} \omega _t -\nu \varDelta \omega +u \cdot \nabla \omega -\omega \cdot \nabla u=0, \end{aligned}$$

(13)

a vectorial equation with the additional term $\omega \cdot \nabla u$, called the vortex stretching term. For instance, if locally, $\partial _3 u_3 >0$ (e.g., $u_3 (x_1,x_2,x_3)\equiv x_3$) and $\omega _3 >0$, but $|\omega _1|+|\omega _2|\ll 1$, then $(\omega \cdot \nabla u)_3 \approx \omega _3 \partial _3 u_3 >0$ so that by (13) $\omega _3$ is likely to increase with respect to t and $x_3$ and could tend to infinity. For further details on regularity related to vorticity, we refer to Sect. 3.1.

2.4 Construction of regular solutions

Although there are several well-established methods to construct global-in-time weak solutions, the problem to prove that they are globally regular is still unsolved—it is one of the Millennium Prize problem.

Another idea is to construct regular solution from the very beginning. Concerning Galerkin’s method, we can use smooth ansatz functions $u_j$ given by an orthonormal basis of eigenfunctions of the Stokes operator $A=-P\varDelta$ in a bounded domain. Actually, similarly to the Laplacian on a bounded smooth domain, A has a compact resolvent, and, consequently, there exists a complete set of eigenfunctions $u_j$ in the Sobolev space $H^2(\varOmega )\cap H^1_0(\varOmega )$. The crucial point to construct a regular solution u on a given interval (0, T) with the property (11) is to get adequate higher order a priori estimates of u. To this aim, the Navier–Stokes system (1) is tested with Au, rather than with u, and yields, amongst others, the integral

$$\begin{aligned} I=\int \nolimits _\varOmega (u \cdot \nabla u) \cdot A u \,\mathrm{d}x. \end{aligned}$$

This integral is estimated by Hölder’s inequality, the interpolation estimate

$$\begin{aligned} \Vert \nabla u\Vert _3 \le c \Vert \nabla u\Vert _2^{1/2} \, \Vert \nabla ^2 u\Vert _2^{1/2} \le c \Vert \nabla u\Vert _2^{1/2} \, \Vert A u\Vert _2^{1/2} \end{aligned}$$

(14)

and the Sobolev embedding $H_0^1(\varOmega ) \subset L^6(\varOmega )$ as follows:

$$\begin{aligned} |I|&\le \Vert u\Vert _6\, \Vert \nabla u\Vert _3 \, \Vert A u\Vert _2 \le c \Vert \nabla u\Vert _2\Vert \nabla u\Vert _3\, \Vert A u\Vert _2\\&\le c \Vert \nabla u\Vert ^{3/2}\, \Vert A u\Vert _2^{3/2} \le \frac{\nu }{2} \Vert A u\Vert _2^2 + \frac{c}{\nu ^3} \Vert \nabla u\Vert _2^6. \end{aligned}$$

The term $\frac{\nu }{2} \Vert A u\Vert _2^2$ can be absorbed by the term $\nu \Vert A u\Vert _2^2$ which appears on the left-hand side of the Navier–Stokes system. However, the function

$$\begin{aligned} y(t)= \Vert \nabla u(t)\Vert _2^2 \end{aligned}$$

can only be shown to satisfy a differential inequality of the form $y' \le C y^3$ and stays bounded for

$$\begin{aligned} t \le T_0 := \frac{\nu ^3/c}{\Vert \nabla u_0\Vert _2^4}. \end{aligned}$$

(15)

These formal arguments can be made rigorous and prove that the Galerkin approximation method yields a strong solution at least on the interval $[0,T_0)$, the length of which depends on norms of the initial value. By the way, in the planar case, the interpolation estimate (14) can be replaced by the much stronger estimate $\Vert \nabla u\Vert _4 \le c \Vert \nabla u\Vert _2^{1/2} \, \Vert \nabla ^2 u\Vert _2^{1/2}$.

On the other hand, Leray [41] has shown that a weak solution (in the whole space case) getting singular as $t \rightarrow T_0-$ in the sense that $\Vert \nabla u(t)\Vert _2 \rightarrow \infty$ blows up at least with the rate

$$\begin{aligned} \Vert \nabla u(t)\Vert _2 \ge \frac{c \nu ^{3/4}}{(T_0-t)^{1/4}}, \quad t \rightarrow T_0-. \end{aligned}$$

Moreover, since friction (dissipation) leads to a loss of kinetic energy, for any global weak solution u, there exists $T_\infty > T_0$ such that u is regular on $(T_\infty ,\infty )$. These two facts lead to the famous Structure Theorem of Leray. Similar results hold for weak solutions on bounded and exterior domains satisfying a so-called strong energy inequality. To be more precise, for a given global weak solution on $[0,\infty )$ there exist finitely or at most countably many open intervals of regularity, $I_j$, in $(0,T_\infty )$, such that u is regular on each interval $I_j$ as well as on $(T_\infty ,\infty )$, and the complement satisfies—in terms of the Hausdorff measure ${\mathcal {H}}^{1/2}$—the condition

$$\begin{aligned} {\mathcal {H}}^{1/2} \Big ([0,\infty \setminus \bigcup \nolimits _j I_j \setminus (T_\infty ,\infty ) \Big ) =0, \end{aligned}$$

(16)

see, e.g. [23]. In particular, the set of temporal singularities is a Lebesgue set of measure 0, but could contain uncountably many epochs.

2.5 Construction of regular solutions by semigroups techniques

Another important method for the construction of regular solutions is based on analytic semigroup theory. It is well-known that the Stokes operator $A=-P \varDelta$ generates a bounded analytic semigroup $e^{-tA}$, $t \ge 0$, for bounded and exterior domains as well as for layers, tubes and half spaces (and variants of them). Then, Duhamel’s formula yields the nonlinear integral equation

$$\begin{aligned} u(t)= e^{-\nu tA} u_0 -\int \nolimits _0^t e^{-\nu (t-\tau )A} P(u \cdot \nabla u) (\tau )\,\mathrm{d}\tau . \end{aligned}$$

(17)

For the analysis of this equation, we need the maximal regularity property of the Stokes operator, i.e., the a priori estimate

$$\begin{aligned} \Vert u_t\Vert _{L^s(0,T,L^q)} + \Vert A u\Vert _{L^s(0,T;L^q)} \le c \Vert f\Vert _{L^s(0,T;L^q)} \end{aligned}$$

(18)

for solutions of the instationary Stokes system

$$\begin{aligned} u_t + \nu A u = f \text { in } \varOmega \times (0,T), \quad u(0)=0. \end{aligned}$$

(19)

Under suitable assumptions on $\varOmega$, this estimate holds for all $1<s, q<\infty$. The first result of this type is due to Solonnikov [54] for $s=q$; for a more functional analytic proof see Geissert et al. [28]. Moreover, Noll et al. [46] proved that the Stokes operator has a bounded ${\mathcal {H}}^\infty$-calculus which, among others, allows to characterize the domains of fractional powers of the Stokes operator; for the latter result, we refer also to Giga [29]. Based on (18) and fractional powers of A, a solution u of (17) can be found by Banach’s fixed point theorem. However, in order to get a contraction and self-map on a closed ball in some Banach space smallness assumptions either on the initial value $u_0$ or on the length of the interval of existence [0, T) must be posed. For a modern review of the semigroup approach, we refer to a handbook article by Hieber and Saal [34].

Since this method is strongly based on $L^q$-theory, $q \not =2$, a careful analysis of the Helmholtz projection and the Stokes operator on $L^q$-spaces is indispensable. However, due to a famous counter-example by Bogovskiǐ [60] there exist smooth unbounded domains on which the Helmholtz projection fails to satisfy the usual mapping properties of existence and uniqueness. In such a case, the Stokes operator $A=-P\varDelta$ is not well-defined. On the other hand, Hilbert space methods, e.g., the Lemma of Lax–Milgram, allow for any bounded and unbounded domain, even with rough boundary, to define the operators P and A on $L^2$-spaces, see [53]. A combination of local $L^q$-theory and global-in-space $L^2$-theory was the key to the analysis of P and A in general smooth unbounded domains in spaces of type

$$\begin{aligned} {\tilde{L}}^q (\varOmega )= {\left\{ \begin{array}{ll} L^2(\varOmega ) + L^q(\varOmega ), &{} 1<q<2\\ L^2(\varOmega ) \cap L^q(\varOmega ), &{} 2 \le q <\infty \end{array}\right. }, \end{aligned}$$

see Kozono, Sohr and the author [19, 20].

Another interesting observation called “Weak Neumann implies Stokes” is described by Geissert et al. in [27]: For a domain with possibly non-compact smooth boundary, assume that the Helmholtz projection P exists on some $L^q$ space, $1< q < \infty$. Then, the Stokes operator generates an analytic semigroup on $L^q$ satisfying maximal regularity; hence, a unique local mild solution to the Navier–Stokes equations can be constructed provided $q > d$.

Looking for strong solutions u in Serrin’s class $L^s(0,T;L^q)$ where $\frac{2}{s}+ \frac{3}{q}=1$, $3<q<\infty$, also the initial value $u_0$ must have more regularity than $L^2_\sigma (\varOmega )$. For the Stokes system (19) with $f=0$, but $u(0)= u_0 \not = 0$, it is obvious that Serrin’s condition for its solution $u(t)=e^{-t A} u_0$ is equivalent to the property

$$\begin{aligned} \int \nolimits _0^\infty \Vert e^{-t A} u_0\Vert _q^s \,\mathrm{d}t < \infty \end{aligned}$$

(20)

which in terms of real interpolation theory can be rewritten in the form

$$\begin{aligned} A^{-1}u_0 \in \big (L^q_\sigma (\varOmega ), {\mathcal {D}}(A_q)\big )_{1-1/s,s} \end{aligned}$$

(21)

where $L^q_\sigma (\varOmega )=\{u \in L^q(\varOmega ): \text {div }u=0 \text { in } \varOmega ,\, u \cdot n=0 \text { on } \partial \varOmega \}$ and ${\mathcal {D}}(A_q)$ is the domain of the Stokes operator A on $L^q_\sigma (\varOmega )$. Moreover, the space defined by (20), (21) equals a “solenoidal” subspace of the Besov space $B_{q,s}^{-1+3/q} (\varOmega )$. Without going into details of Besov spaces on domains, we mention that in view of Serrin’s condition the order of regularity $-1+ \frac{3}{q}$ of the space $B_{q,s}^{-1+3/q}(\varOmega )$ is negative. It was shown by Sohr et al. [21] that the optimal (largest) space of initial values for the nonlinear Navier–Stokes system to yield an at least local-in-time weak and also strong solution in Serrin’s class is given by the “Besov space” $B_{q,s}^{-1+3/q}(\varOmega ) \cap L^2_\sigma (\varOmega )$. Moreover, there exists an absolute constant $\varepsilon _*= \varepsilon _*(q,\varOmega ) >0$ such that any $T >0$ satisfying the smallness condition

$$\begin{aligned} \int \nolimits _0^T \Vert e^{-\tau A} u_0\Vert _q^s \,\mathrm{d}\tau \le \varepsilon _*\end{aligned}$$

defines an initial interval of regularity [0, T), see [22].

3 Advanced analysis of weak and regular solutions

3.1 Conditional regularity

Besides the Serrin condition (10) for the velocity field u, there are numerous related conditions on components of u or $\nabla u$, etc. A classical generalization by Beirão da Veiga [1] concerns $\nabla u$ such that $S_{\nabla u}=2$. More technical proofs are necessary when only one or two components of u or $\nabla u$ are taken into account, and it is not immediate to decide whether the same Serrin condition applies to $u_3$ as for u, say. Most of the following results are proven for the case $\varOmega =\mathbb {R}^3$, but we are not touching on problems and details for domains $\varOmega \ne \mathbb {R}^3$.

In [44], Neustupa et al. proved that if one component of a (suitable) weak solution locally satisfies $S_{u_3}=\frac{1}{2}$, say, then the solution is locally regular. Recently, an optimal result for one column, $\partial _3 u$, of $\nabla u$ was found by Z. Skalák [52]: If $S_{\partial _3 u}=2$, then the solution is regular. Looking for the row $\nabla u_3$ of $\nabla u$, the $L^2$-case $\nabla u_3\in L^4(L^2)$ where $S_{\nabla u_3}=2$ was solved by Han et al. [33], together with further results on fractional derivatives. Considering only one component of the tensor $\nabla u$ results are less complete in the sense that $S_{\partial _iu_j}$ is assumed to be less than 2. Here, for diagonal element $\partial _3u_3$ results are better than for off-diagonal elements. We refer to Cao and Titi [9] for both types of results which are too technical for this short review. If mixed norms are used, i.e., with different integration exponents in $x_1$-, $x_2$- and $x_3$-direction, then the critical Serrin number $S=2$ can be approached for diagonal elements $\partial _iu_i$, but not for off-diagonal elements so far; see Guo et al. [32].

Concerning the vorticity $\omega$, it is not the size of $\omega$ itself which is responsible for the possible existence of singularities as proposed by the vortex stretching mechanism. Indeed, Constantin and Fefferman [15] found that the crucial term is the angle $\theta (x,y,t)= \sphericalangle (\omega (x,t), \omega (y,t))$ which should not vary too fast for large $|\omega |$; if

$$\begin{aligned} |\sin \theta (x,y,t)| \le K |x-y|^{1/2} \end{aligned}$$

for large vorticities $|\omega (x,t)|$, $|\omega (y,t)| \ge M$, then the solution will not blow up. A substantial improvement was found by Beirão da Veiga and Berselli [2] in the sense that it suffices to get for large $\omega (x,t), \omega (y,t)$ a bound $|\sin \theta (x,y,t)| \le g(t,x) |x-y|^\alpha$ where $\alpha \in [\frac{1}{2},1]$ and $g\in L^s(L^q)$ with $S_g=\alpha -\frac{1}{2}$. Returning to the 2d case in which $\omega \widehat{=} (0,0,\omega )^T$, the angle $\theta (x,y,t)$ is constant to either 0 or $\pm \pi$ so that the above condition is satisfied everywhere. Actually, if two components of the vorticity in the 3d case can be controlled, then the 3d-situation is close to the 2d one yielding regularity; however, it is an open problem whether the control of only one component of $\omega$ is sufficient to get regularity.

For the pressure p, there are less criteria available in the literature. Recall that p is unique only up to functions depending on t so that assumptions like $p\ge 0$ or $p\ge -C$ are meaningless. However, Seregin and Šverák [49] proved that a weak solution is regular provided there exists a function $g\ge 0$ such that $p(x,t)\ge -g(x,t)$ and g satisfies the condition

$$\begin{aligned} \sup _{x_0\in \mathbb {R}^3} \sup _{t_0-r_0^2\le t<t_0} \int \nolimits _{B_{r_0}(x_0)} \frac{g(x,t)}{|x-x_0|}\,\mathrm{d}x <\infty \end{aligned}$$

for each $t_0>0$; moreover, the function $t\mapsto \int \nolimits _{B_{r_0}(x_0)} \frac{g(x,t)}{|x-x_0|}\,\mathrm{d}x$ is assumed to be continuous at $t_0$ from the left. The same result holds when the total head pressure $\frac{1}{2} |v|^2+p$ is bounded above by g. The condition above generalizes the pointwise condition $p\ge -C$ and prevents the fluid from developing the phenomenon of cavitation.

A completely different criterion is based on the geometry of the rate of deformation tensor, the symmetrized velocity gradient $D(u) = \frac{1}{2}\big (\nabla u + (\nabla u)^\top \big )$. Since D is symmetric with vanishing trace due to the solenoidality of u, D has three real eigenvalues $\lambda _1\le \lambda _2\le \lambda _3$ summing up to 0. Obviously, $\lambda _1 \le 0$ and $\lambda _3 \ge 0$ without any a priori knowledge about the sign of $\lambda _2$. Given an orthonormal basis of $\mathbb {R}^3$ of eigenvectors of D, $\{e_1,e_2,e_3\}$, it can be shown from (2) that an infinitely small test volume of the fluid is compressed or stretched in the direction $e_i$ if $\lambda _i<0$ or $\lambda _i>0$, respectively. Hence, in the case when $\lambda _2\le 0$ a test cube will be compressed to a thin but long tube, whereas in the second case ($\lambda _2>0$) the test cube will be pressed flat to a thin sheet. Then, if the positive part of $\lambda _2$ as a function of (x, t) can be controlled locally by some integrability condition, the solution is locally regular. This result is due to Neustupa and Penel [45].

3.2 Analysis in scale invariant borderline spaces

The classical borderline case $C^0[0,T]; L^3(\varOmega )) \subset L^\infty (0,T; L^3(\varOmega ))$ for strong solutions can neither be reached by the Galerkin method nor the semigroup approach in Sect. 2.5 directly; one reason is the elementary fact that the norm in $L^\infty (0,T)$ does not necessarily converge to 0 for a bounded function as $T \rightarrow 0$. This problem was solved for the first time by Kato [37] by a more involved iteration scheme than the classical fixed point iteration of Banach. The idea is to start a first iteration for $(u_j)_{j \in \mathbb {N}}$ based on (17) such that

$$\begin{aligned} t^{\frac{1}{2}(1-\frac{3}{q})} u_j(t)&\in C \big ([0,T]; L^q(\varOmega )\big ),\\ t^{\frac{1}{2}} \nabla u_j(t)&\in C [0,T]; L^3(\varOmega )\big ). \end{aligned}$$

Here, $q >3$ is any exponent. Proving the convergence of this scheme, in a second step the convergence of $(u_j)$ in $C ([0,T; L^3(\varOmega ))$ and the property

$$\begin{aligned} u \in C \big ([0,T]; L^3(\varOmega )\big ) \end{aligned}$$

are proved. The proof by Kato was designed for the whole space $\mathbb {R}^3$, but can be adapted to many other domains. Moreover, if $u_0 \in L^3(\varOmega )$ is sufficiently small, then $T=\infty$ is possible. A variant of Kato’s method was described by Giga [30] where the gradient term in the iteration is no longer needed.

The largest space of initial values allowing for local and/or global in time solutions in scale invariant function spaces was found so far by Koch and Tataru [39]. Assume that $u_0$ lies in the scaling invariant space $BMO^{-1}(\mathbb {R}^3)$, i.e., $u_0$ can be written in the form $u_0=\text {div }f$ with $f\in BMO(\mathbb {R}^3)$ and has small norm in this space. Then, there exists a unique local/global solution $u\in L_{\mathrm{loc}}^2(\mathbb {R}^3\times [0,\infty ))$ such that the scaling invariant norm

$$\begin{aligned} \Bigg (\sup _{x\in \mathbb {R}^3,R>0} \frac{1}{R^3} \int \nolimits _0^{R^2}\! \int \nolimits _{\{y:\, |y-x|<R\}} |u|^2\, \mathrm{d}y\,\mathrm{d} \tau \Bigg )^{1/2} \end{aligned}$$

related to the $BMO^{-1}$-norm is finite.

However, the largest scale invariant space for the Navier–Stokes system is the homogeneous Besov space $\dot{B}^{-1}_{\infty ,\infty }$ (for a proof see Cannone [8, Proposition 7]) containing the space $BMO^{-1}$; hence, conditions formulated in this space are among the most interesting. Since $B^{-1}_{\infty ,\infty }$ is larger than $\dot{B}^{-1}_{\infty ,\infty }$, initial values in $B^{-1}_{\infty ,\infty }$ will not be considered, but rather regularity question. E.g., Cheskidov and Shvydkoy [13] proved that a weak Leray–Hopf solution $u\in C^0((0,T]; B^{-1}_{\infty ,\infty })$ or with sufficiently small jumps in the $B^{-1}_{\infty ,\infty }$-norm is regular on (0, T]. On the other hand, weak solutions of Leray–Hopf type may be discontinuous at $t=0$ in the topology of $B^{-1}_{\infty ,\infty }$. Actually, for the torus ${\mathbb {T}}^3$, Cheskidov and Shvydkoy [12] construct an initial value $u_0\in \bigcap _{1<q\le \infty } B^{3/q-1}_{q,\infty } \cap L^2_\sigma$ such that each Leray–Hopf type weak solution u with $u(0)=u_0$ satisfies the discontinuity estimate

$$\begin{aligned} \limsup _{t\rightarrow 0+} \Vert u(t)-u_0\Vert _{B^{-1}_{\infty ,\infty }} \ge \delta ; \end{aligned}$$

here, $\delta >0$ is independent of u.

The Navier–Stokes system is ill-posed in the space $\dot{B}^{-1}_{\infty ,\infty }$, see Bourgain and Pavlović [3]. Actually, the authors proved that for any $\delta >0$ there exists an initial value $u_0$ in the Schwartz class $\mathcal S(\mathbb {R}^3)$, a $T\in (0,\delta )$ and a corresponding solution (u, p) to the Navier–Stokes system such that

$$\begin{aligned} \Vert u_0\Vert _{\dot{B}^{-1}_{\infty ,\infty }} \le \delta \quad \text {and}\;\; \Vert u(T)\Vert _{\dot{B}^{-1}_{\infty ,\infty }} \ge \frac{1}{\delta }. \end{aligned}$$

This result was extended by Yoneda [59] to the spaces $\dot{B}^{-1}_{\infty ,r}$ for $r>2$, and by Wang [57] to the homogeneous Besov spaces $\dot{B}^{-1}_{\infty ,r}$ for any $1\le r\le 2$.

On the other hand, there exist initial values of arbitrarily large norms in $\dot{B}^{-1}_{\infty ,\infty }$ allowing for global smooth solutions of the Navier–Stokes system. In [11], Chemin et al. construct initial values of the form

$$\begin{aligned} u_{0,\varepsilon } = \Big ( v_0^h(x_h,\varepsilon x_3), \frac{1}{\varepsilon } v_0^3(x_h,\varepsilon x_3)\Big ), \quad 0<\varepsilon < \varepsilon _0, \end{aligned}$$

on the domain $(0,2\pi )^2\times \mathbb {R}$ where $x_h=(x_1,x_2)$ and $x_3\in \mathbb {R}$. Note that due to the factor $\frac{1}{\varepsilon }$ the norm $\Vert u_{0,\varepsilon } \Vert _{\dot{B}^{-1}_{\infty ,\infty }}$ increases like $\frac{1}{\varepsilon }$ as $\varepsilon \rightarrow 0$, see [11, p. 989]. Further examples of highly oscillating initial values $u_\varepsilon \in \dot{B}^{-1}_{\infty ,\infty }(\mathbb {R}^3)$ with $\Vert u_\varepsilon \Vert _{B^{-1}_{\infty ,\infty }} \nearrow \infty$ as $\varepsilon \rightarrow 0$ are found in [10].

For a recent survey on the analysis in critical function spaces, we refer to Gallager [26].

3.3 Local regularity and possible singularities

In contrast to global-in-space regularity results for weak solutions u looking for singular epochs as in Leray’s Structure Theorem, see (16), local regularity theory is interested in space-time points $z_0=(x_0,t_0)\in \mathbb {R}^4$ such that $u\notin L^\infty (U_\delta (z_0))$ for any neighborhood $U_\delta \subset \mathbb {R}^4$ of $z_0$. Due to the local character, the weak solution has to satisfy a localized version of the usual energy inequality and is called a suitable weak solution. The proof of existence of suitable weak solutions needs a more careful construction, in particular, the use of cut-off functions leads to pressure terms which must be estimated. Suitable weak solutions are known to exist for bounded, exterior and more general unbounded smooth domains, cf. [7, 19]. The celebrated result of Caffarelli et al. [7] states that for a suitable weak solution u the set of singular space-time points, ${\mathcal {S}} \subset \mathbb {R}^4$, satisfies

$$\begin{aligned} {\mathcal {H}}^1({\mathcal {S}}) = 0. \end{aligned}$$

(22)

Due to the different scaling in time and space, the Hausdorff measure ${\mathcal {H}}^1$ can be replaced by an even stronger parabolic variant, ${\mathcal {P}}^1$, of ${\mathcal {H}}^1$. This result implies that the set ${\mathcal {S}}$ cannot contain a line or any smooth path; thus, for axially symmetric flows, singularities can occur only on the axis of rotation. The result (22) is an integrated version of a local regularity analysis on space-time cylinders $Q_r(z_0) = B_r(x_0) \times (t_0,t_0-r^2)$ using either the criterion

$$\begin{aligned} \mathop {\text {lim sup}}\limits _{r\rightarrow 0+} \frac{1}{r} \iint _{Q_r(z_0)} |\nabla u|^2 \,\mathrm{d}x\, \mathrm{d}\tau <\varepsilon _0 \end{aligned}$$

or

$$\begin{aligned} \frac{1}{r^2} \iint _{Q_r(z_0)} |u|^3 \,\mathrm{d}x\, \mathrm{d}\tau <\varepsilon _0, \end{aligned}$$

where the smallness parameter $\varepsilon _0$ is an absolute constant. We note that the above integrals are scale-invariant quantities with respect to the scaling (9). Moreover, the term $\nabla u$ can be replaced by $\omega =\mathrm{rot}\, u$ and the $L^3(L^3)-$norm of u can be replaced by any scale-invariant local, weighted $L^s(L^q)$-norm ($3\le s,q\le \infty$) of u. For proofs and further results, we refer to [7, 40, 42, 58].

A classical idea of Leray [41] to construct solutions with singularities uses an ansatz of self-similar solutions of the form

$$\begin{aligned} u(x,t) = \frac{1}{\sqrt{-t}} V\Big (\frac{x}{\sqrt{-t}}\Big ) \end{aligned}$$

(23)

with $V\in L^3(\mathbb {R}^3)$, getting singular as $t\rightarrow 0-$. Since u is a weak solution of the Navier–Stokes equations, the vector field V satisfies a t-independent modified Navier–Stokes equations with an additional y-dependent term of the form $y\cdot \nabla _y V$. However, the problem to find a nonzero field V leading to a singular weak solution was left open and answered not before the end of the last century by Nečas et al. [43] in the negative ($V\equiv 0$). The same negative answer was given by Tsai [56] in case that either $V\in L^q(\mathbb {R}^3)$ for any $q>3$ or the self-similar solution u has bounded energy on the space-time cylinder $(-1,0)\times B_1(0)$ in the sense that

$$\begin{aligned} \mathop {\text {ess sup}}\limits _{-1<t<0} \int \nolimits _{B_1(0)} |u|^2\,\mathrm{d}x + \int \nolimits _{-1}^0 \int \nolimits _{B_1(0)} |\nabla u|^2\,\mathrm{d}x\, \mathrm{d}\tau < \infty . \end{aligned}$$

In the pointwise sense, the singularities discussed above are of the form $|u(x,t)|\sim \frac{c}{|x|+\sqrt{-t}}$ and called singularities of Type I. They can be excluded in the ansatz of self-similar solutions as above and for axisymmetric flows, see [14, 38]. All other singularities are called of Type II.

4 Recent approaches

4.1 The counter-example by Terence Tao

A solution with smooth initial value and singularity in finite time was constructed by Tao [55] for a modified averaged Navier–Stokes system on $\mathbb {R}^3$. The counter-example is based on a blow-up solution for the Euler system so that the dissipation term $-\varDelta u$ does act only as a pertubation and finally allows for a singularity of Type II. The first idea is a discrete (dyadic) dynamical system of scalar functions $(X_n(t))_{n\in \mathbb {Z}}$ of the form

$$\begin{aligned} \partial _t X_n(t) = -\lambda ^{2n\alpha }X_n + \lambda ^{n-1} X_{n-1}^2 - \lambda ^n X_n X_{n+1} \end{aligned}$$

(24)

with parameter $\lambda =(1+\varepsilon _0)^{5/2}$, $0<\varepsilon _0\ll 1$. Actually, in the final analysis $(X_n)$ is a family of vector fields in $\mathbb {R}^4$ with pairwise disjoint supports in Fourier space such that energy is either pumped or drained or rotated between different components, either slowly or abruptly. Then, the solution u blowing up in finite time has the form $u(x,t) \sim \sum _n X_n(t) \lambda ^{3n/5} \psi (\lambda ^{2n/5}x)$ and is based on a Schwartz function $\psi$ with Fourier transform supported in an annulus of radii $\lambda ^{-5/2}$ and $\lambda ^{5/2}$. To control the nonlinear transport of energy in (24) it would be helpful if at a given time only one pair $(X_n,X_{n+1})$ interacts. This can be achieved by a introducing a piecewise constant function $n=n(t)$ in (24), but would lead to a t-dependent nonlinear term in the Navier–Stokes system. Here and for the transformation to a continuous model a sophisticated averaging process plays a crucial role. To this aim, we return to the bilinear operator $B(u,v) := -\frac{1}{2}(u\cdot \nabla v + v\cdot \nabla u)$ and the trilinear form

$$\begin{aligned} \langle B(u,v),w \rangle _{\mathbb {R}^3} = -\frac{1}{2}\int \nolimits _{\mathbb {R}^3} \big (u\cdot \nabla v + v\cdot \nabla u\big )w\,\mathrm{d}x \end{aligned}$$

which vanishes for solenoidal vector fields $u=v=w$, see (8). Then, Tao introduces the operators ${\mathcal {R}} u(x):= R_\xi u(R_\xi ^{-1} x)$ defining the rotation around the axis $\xi \in \mathbb {R}^3$ which map solenoidal vector fields to solenoidal ones. Moreover, classical multiplier operators $\mathcal Mu = \mathcal F^{-1}m(\cdot ){\mathcal {F}} u$ on $L^p(\mathbb {R}^3)$ are needed. With these definitions, the trilinear form $\langle B(u,v),w \rangle _{\mathbb {R}^3}$ is replaced by the averaged form

$$\begin{aligned} \langle {{\tilde{B}}}(u,v),w \rangle _{\mathbb {R}^3} = \int \nolimits _\varOmega \big \langle B({\mathcal {M}}^{(1)}_\omega {\mathcal {R}}^{(1)}_\omega u,{\mathcal {M}}^{(2)}_\omega {\mathcal {R}}^{(2)}_\omega v),{\mathcal {M}}^{(3)}_\omega {\mathcal {R}}^{(3)}_\omega w \big \rangle _{\mathbb {R}^3} \,\mathrm{d}\mu \end{aligned}$$

where $\omega$ is a random variable defined on a subset $\varOmega$ of the space $SO(\mathbb {R}^3)\times SO(\mathbb {R}^3)\times SO(\mathbb {R}^3)\times \mathbb {R}^3 \times \mathbb {R}^3 \times \mathbb {R}^3$. The crucial properties of ${{\tilde{B}}}$ are the cancellation property $\langle {{\tilde{B}}}(u,u),u\rangle _{\mathbb {R}^3} = 0$ and energy estimates as satisfied by B. The analysis is performed in Fourier space in which the energy cascade is mainly controlled by balls tending to infinity as time t approaches the finite blow-up time $T_*$ more or less by jumps on a sequence of epochs $(t_n)$ tending to $T_*$ exponentially fast. Tao’s example implies that an affirmative answer to the regularity problem must either be based on more sophisticated estimates and hidden properties of B which do not hold for ${{\tilde{B}}}$ or on a further balance among the terms $u_t+u\cdot \nabla u$, $-\varDelta u$ and the pressure $\nabla p$.

4.2 Recent approaches to non-uniqueness

In [36], Jia and Šverák construct a situation with initial values admitting multiple solutions. Starting point is a self-similar initial value $u_0$ on $\mathbb {R}^3$ such that the Navier–Stokes system possesses for each $\sigma u_0$, $0\le \sigma \ll 1$, a unique global self-similar solution $u_\sigma (x,t) = t^{-1/2} U_\sigma (x t^{-1/2})$, cf. (23). Linearizing the modified Navier–Stokes equation for $U_\sigma$ around $U_\sigma$ the linear operator ${\mathcal {L}}_\sigma$,

$$\begin{aligned} {\mathcal {L}}_\sigma \varphi = \varDelta \varphi + \frac{x}{2}\cdot \nabla \varphi +\frac{1}{2}\varphi - U_\sigma \cdot \nabla \varphi - \varphi \cdot \nabla U_\sigma +\nabla p, \end{aligned}$$

is the crucial operator for a further bifurcation analysis. It leads to certain assumptions on eigenvalues of ${\mathcal {L}}_\sigma$ such that for $\sigma>\sigma _0>0$ a pitchfork bifurcation will occur. The theoretical assumptions could not yet be verified rigorously in concrete cases. However, numerical tests by Guillod and Šverák [31] indicate that the above scenario will develop.

The method of convex integration developed about 70 and 50 years ago by Nash and Gromov is the basic tool of De Lellis und Székelyhidi Jr. to prove the existence of weak solutions to the 3D Euler equations in the class $C^\beta$, $0<\beta <\tfrac{1}{3}$, in space-time with prescribed kinetic energy $E_{\mathrm{kin}}(t)\ge 0$, s [6, 16, 17]. In the last years, these results were generalized to the Navier–Stokes equations where the smoothing property by dissipation (friction) was the crucial problem: there exist “weak” solutions u in the class $C^0([0,T];L^2(\varOmega ))$ such that u will neither necessarily satisfy the energy inequality nor have finite dissipation energy $\int \nolimits _0^t \Vert \nabla u(\tau )\Vert _2^2 \,{\text {d}}\tau$. Buckmaster and Vicol [4] prescribe a smooth function $E_{\mathrm{kin}}(t)\ge 0$ and construct a weak solution $u\in C^0([0,T];L^2({\mathbb {T}}^3)) \cap W^{1,1+\beta } ({\mathbb {T}}^3))$ with $\beta > 2^{-16}$ and $\frac{1}{2}\Vert u(t)\Vert _2^2 = E_{\mathrm{kin}}(t)$. Since $E_{\mathrm{kin}}$ is not assumed to be monotonically decreasing, u will not be a physically reasonable solution. This result was extended by those authors and Colombo [5] as follows: given two smooth solutions $u^{(1)}$, $u^{(2)}$ of the Navier–Stokes system on ${\mathbb {T}}^3 \times [0,T]$ with $f\equiv 0$ there exists a weak solution in the above sense with $u(0) = u^{(1)}(0)$ such that

$$\begin{aligned} u\equiv u^{(1)}\; \text { in }\;[0,T/3] \quad \text {und }u\equiv u^{(2)}\quad \text {in }\;[2T/3,T]. \end{aligned}$$

In particular, choosing $u^{(1)}\equiv 0$ and $u^{(2)} \ne 0$ there exists a solution which is at rest on some initial time interval, but bifurcates to a nonzero smooth solution. On the other hand, u may start with $u\equiv u^{(1)}\ne 0$ on [0, T/3] but come to rest at $t=2T/3$. Nevertheless, the set ${\mathcal {S}}$ of singular instants of u is small in the sense ${\mathcal {H}}^{1-\beta }({\mathcal {S}}) =0$.

Since these solutions cannot be shown to satisfy the energy inequality on [T/3, 2T/3], the question arises whether there exist weak solutions in the sense of Leray–Hopf with given physically reasonable kinetic energy function. A partial answer is given by Ożański [48] who prescribes an arbitrary continuous, non-increasing energy function $E_{\mathrm{kin}}(t)$ and constructs a suitable weak solution u satisfying the localized energy inequality such that $\big |\frac{1}{2} \Vert u(t)\Vert _2^2 - E_\mathrm{kin}(t)\big |$ is bounded by any given positive constant. However, u satisfies a Navier–Stokes inequality, i.e., u is a solution in the Leray–Hopf class with the property that

$$\begin{aligned} f:=u_t-\nu \varDelta u +u\cdot \nabla u +\nabla p,\quad f\cdot u\le 0. \end{aligned}$$

In other words, the force f is not given, but is defined a posteriori and satisfies $f\cdot u\le 0$; i.e., formally f acts against the direction of the flow field.

5 Conclusion

Even after more than 85 years when Leray constructed global weak solutions, the regularity problem for global-in-time weak solutions of the Navier–Stokes equations or the problem of the direct construction of global-in-time regular solutions is still open. After an analysis in classical $L^2$, $L^p-$spaces, Lorentz and Sobolev spaces, with and without weights, modern analysis also works in Besov spaces, but there are also results on regularity in spaces of type BMO, Hardy ${\mathcal {H}}^1$, Triebel–Lizorkin, Morrey, and others. However, the gap between weak solutions with Serrin number $S=\frac{3}{2}$ and strong solutions with $S=1$ could not be closed. Further work on the construction of solutions with finite-in-time blow up either failed (e.g., self-similar solutions by Leray) or lead to solutions of a modified Navier–Stokes system (Tao) where it is not clear whether and how the construction can be transferred to the original PDE system. A recent trend is to prove non-uniqueness in a new class of weak solutions. However, these solutions which have striking and unexpected behavior do not satisfy the physical properties (energy inequality, finite dissipation energy) as required from the physical point of view.