A gradient left to itself

What a parcel does in the first instant splits the velocity gradient into a stretch and a spin and stops there, one instant in. That is enough to say what a flow is doing now. It is not enough to say what the gradient itself will do next, because the gradient is not a fixed property of the place — it is carried by the parcel, and it changes, and the thing it changes in response to is mostly itself.

Differentiate the Euler equations once in space and follow a parcel. The gradient $A_{ij} = \partial u_i/\partial x_j$ obeys

$$\frac{D\mathbf A}{Dt} = -\mathbf A^2 - \mathbf H,$$

where $\mathbf H$ is the matrix of second derivatives of the pressure. The first term is local: the gradient multiplying itself, a thing any single parcel could compute from its own neighbourhood. The second is not local at all. The pressure is fixed by a Poisson equation over the whole of the fluid at once, which is the reason it has no speed, and its Hessian at one point depends on the gradient everywhere else.

This essay throws the non-local part away and asks what is left. Keep only as much of the pressure as incompressibility forces on one parcel by itself — its trace, spread evenly in every direction — and the equation closes on nine numbers:

$$\frac{d\mathbf A}{dt} = -\mathbf A^2 + \tfrac13 \operatorname{tr}\!\left(\mathbf A^2\right)\mathbf I.$$

It is called the restricted Euler equation. It is not the Euler equation, and the whole value of it is seeing exactly how it fails, because it fails in a way that says which part of the physics was doing the work.

Two numbers carry the whole trajectory #

A trace-free three-by-three matrix has two invariants, the ones the classification of a critical point in space is built on:

$$Q = -\tfrac12 \operatorname{tr}\mathbf A^2, \qquad R = -\tfrac13 \operatorname{tr}\mathbf A^3 .$$

$Q$ is positive where the spin outweighs the strain and negative where the strain outweighs the spin — it is a quarter of the squared vorticity minus half the squared strain rate, which is why it is also the quantity called the Q-criterion. $R$ carries the sign of the stretching: positive when two directions are being pulled apart and one squeezed, negative when one is pulled and two squeezed.

Take the trace of the restricted equation multiplied by $\mathbf A$, and again by $\mathbf A^2$, and the nine equations collapse onto two that mention nothing else:

$$\frac{dQ}{dt} = -3R, \qquad \frac{dR}{dt} = \tfrac23 Q^2 .$$

That closure is the first surprise, and it is exact: the invariants evolve without reference to the orientation of the matrix or the direction of its axes. It was checked here by differencing $Q$ and $R$ along the nine-component right-hand side for twenty random gradients, and the two relations hold to $9\times10^{-11}$.

The second surprise follows in one line. Multiply the first equation by $Q^2$ and the second by $\tfrac{27}{2}R$ and add:

$$\frac{d}{dt}\left(\tfrac{27}{4}R^2 + Q^3\right) = 0 .$$

So $D = \tfrac{27}{4}R^2 + Q^3$ never changes, and every trajectory is the explicit curve $Q = \sqrt[3]{D - \tfrac{27}{4}R^2}$. There is no time in it and no integration in it; the figure above is simply those curves drawn for several values of $D$. And because $dR/dt = \tfrac23 Q^2$ is never negative, $R$ only ever increases. Every curve is travelled from left to right, and as $R$ grows without bound the cube root sends $Q$ to minus infinity along the curve $D = 0$.

$D = 0$ is Vieillefosse’s line, $Q = -3(|R|/2)^{2/3}$, and it is the curve on which the matrix has a repeated eigenvalue. It is also the boundary a critical point in space crosses when it changes from a spiralling kind to a purely straining one. Here it has a second job: its right-hand branch is where every trajectory in the plane ends.

The sign that decides #

The cleanest way to see the end is the one case with a closed form. Let the gradient be an axisymmetric strain, $\mathbf A = a\,\mathrm{diag}(1, 1, -2)$. The matrix stays diagonal, and the restricted equation reduces to the one line

$$\frac{da}{dt} = a^2, \qquad a(t) = \frac{a_0}{1 - a_0 t}.$$

This is the same Riccati equation that describes a compression wave steepening into a shock: in one dimension, with no pressure at all, a velocity gradient obeys $d(u_x)/dt = -u_x^2$ and a compressive one reaches minus infinity at $t = 1/|u_x|$. The three-dimensional incompressible equation carries the same quadratic, and the isotropic pressure removes only the part that would change the volume.

Now the sign. With $a_0 > 0$, two directions are stretched and one is squashed. A small sphere of fluid becomes a pancake, and $a$ reaches infinity at $t^* = 1/a_0$. With $a_0 < 0$, one direction is stretched and two are squashed, the sphere becomes a needle — and $a$ does not blow up at all. It decays, as $1/t$, for ever.

The principal stretches of a material sphere are the exponentials of the integrated rates, and here they come out in closed form: $1/(1-a_0t)$ twice and $(1-a_0t)^2$ once. For the sheet that last factor goes to zero at a finite time, so the thickness of the pancake vanishes. For the tube the same formula with $a_0$ negative gives a length growing as $(1+|a_0|t)^2$ — polynomial, never singular. Both strains were normalised to the same size, $|\mathbf A| = 1$, so the sheet’s singular time is $\sqrt6 = 2.449$; the integration reached it to fourteen figures.

That asymmetry is the most useful thing the restricted equation says, and it is the opposite of the picture the word stretching suggests. A vortex tube being stretched is the canonical image of what turbulence does to vorticity, and the local, self-driven dynamics of the gradient do not favour it. Left to itself, a gradient runs away towards sheets.

How the end arrives #

A general gradient has spin in it, off-diagonal parts, no symmetry at all, and its trajectory has no closed form. It still has the conserved $D$, and that pins down how it approaches the singularity.

Near the end the trajectory hugs the right-hand branch of $D = 0$, on which $Q = -3(R/2)^{2/3}$. Substituting into $dR/dt = \tfrac23 Q^2$ gives an equation for $R$ alone whose solution is

$$R = \frac{2}{(t^*-t)^3}, \qquad Q = -\frac{3}{(t^*-t)^2}, \qquad |\mathbf A| \sim \frac{1}{t^*-t},$$

and the time still to go, read from where the trajectory is, is $t^* - t = 2^{1/3}R^{-1/3}$ exactly. Each invariant grows as the power of the gradient it is built from, and the gradient grows as one over the time left, as it does in the one-dimensional steepening.

The gradient in that figure was chosen to start in the upper half of the plane, where vorticity dominates: its $Q$ is positive and its $R$ is small. It circles over the top of the phase plane, crosses $Q = 0$, and falls onto the attracting branch. The regression of $\log|\mathbf A|$, $\log(-Q)$ and $\log R$ against $\log(t^*-t)$, over the stretch where the gradient is between a hundred and a hundred thousand times its starting size, gives $-1.0000$, $-2.0000$ and $-3.0000$.

The time was found two ways, and they had to agree. Once by integrating all nine components of the matrix with a fourth-order Runge–Kutta step of $h/|\mathbf A|$ — so the step shrinks exactly as fast as the solution grows — until the gradient had grown a million-fold, then adding the closed-form remainder. Once by integrating only the pair $(Q, R)$ in the same way. For twelve random starting gradients the two singular times agree to $1.2\times10^{-12}$. And $D$, which the nine-component integration knows nothing about, stays constant along it to $7\times10^{-14}$ of the size of its two terms.

What the singular matrix looks like #

Scaled by the time left, the gradient settles on a fixed matrix: $\mathbf A(t^*-t) \to \mathbf B$. Putting $\mathbf A = \mathbf B/(t^*-t)$ into the restricted equation gives the matrix equation

$$\mathbf B^2 + \mathbf B - 2\mathbf I = 0,$$

whose roots are $1$ and $-2$. With the trace zero, $\mathbf B$’s eigenvalues must be $1, 1, -2$: every singularity of the restricted equation is, in the end, a sheet — two directions stretching at the same rate and one collapsing at twice it. At the end of the integration the residual of that matrix equation is $4\times10^{-10}$.

What the equation does not fix is whether $\mathbf B$ is symmetric. A matrix with eigenvalues $1, 1, -2$ can have an off-diagonal part, and when it does the vorticity survives all the way to the singularity. For the gradient followed here, $|\boldsymbol\omega|$ ends at 0.905 of the strain rate’s magnitude. The spin does not die; it comes into line.

The alignment is exact, and the reason is short. $\mathbf B$ acts as the identity on a plane and as $-2$ on a line not in it. Take the one direction in that plane which is perpendicular to the line: $\mathbf B$ leaves it alone, and so does $\mathbf B^{\mathsf T}$, because nothing else $\mathbf B$ does has a component along it. A direction that a matrix and its transpose both leave alone is an eigenvector of the symmetric part with eigenvalue 1 and is annihilated by the skew part — which means the vorticity points along it. So the middle strain rate is the gradient’s own double eigenvalue, 1/(t − t), to seven figures, and the vorticity lies along its axis*. The other two strain rates depend on how far $\mathbf B$ is from symmetric, which depends on where the trajectory started; for this start they end at $1.581$ and $-2.581$.

This is where the calculation meets measurement. The best-known statistic of real turbulent velocity gradients, first reported from direct simulations by Ashurst, Kerstein, Kerr and Gibson in 1987, is exactly this alignment: the vorticity lies preferentially along the intermediate eigenvector of the strain, not the most-stretching one as the tube picture predicts. And the intermediate strain rate is preferentially positive, with average rates in something like the ratio $3 : 1 : -4$. The restricted equation predicts both qualitatively from nothing but the local quadratic term. It predicts them too strongly, and it predicts them at a singularity that real flows do not reach, which is the next section.

Every start, not a chosen one #

A single trajectory is a demonstration. The claim is about all of them.

The ensemble is isotropic and Gaussian, which is a choice rather than a claim about turbulence; any ensemble without a special relation between its nine entries will do the same, because the only starts that avoid the singularity sit at a point and on a single curve of a two-dimensional plane. The spread in the times is the spread in how close each start lands to them.

The point is the origin, and it holds more than the zero matrix. Simple shear is there. Its gradient has one entry, $A_{12}$, and its square is zero, so $Q = R = 0$ and the right-hand side of the restricted equation vanishes identically: a pure shear is an exact steady state of it. That is the flow the first instant found to be half spin and half strain, and here it is the one gradient of any size that the local dynamics leave alone. Every nearby gradient drifts off it, slowly at first because both invariants are small, which is why the slowest member of the ensemble was nearly a shear. How slowly is exact: a shear with a sheet-forming strain a thousandth of its size added to it reaches the singularity in a thousand time units, the inverse of the nudge, because the shear’s square contributes nothing and the added strain runs its own Riccati clock.

That curve is the left-hand branch of $D = 0$. A gradient started exactly on it — an axisymmetric tube, or anything with the same two invariants — slides along it towards the origin, slowing as it goes, and never arrives: its speed along the curve is $\tfrac23 Q^2$, which vanishes at the origin. It is the only way in, and it is a knife edge.

The second trajectory is worth dwelling on because it is not a failure of the calculation. In exact arithmetic it would slide in for ever. In floating point it carries an error of order $10^{-16}$ in $D$, and a non-zero $D$ is a different trajectory — one of the curves that swing away to the right. The error grows because the branch is unstable to exactly that departure, and the time it takes to matter is the time a perturbation of $10^{-16}$ takes to become of order one. The only gradients that do not blow up are a set of measure zero, and they are unstable even to rounding.

What was thrown away, and what it must be doing #

The restricted equation is a statement about what the gradient would do if the rest of the fluid had no say. The rest of the fluid does have a say, through the part of the pressure Hessian that was discarded, and through viscosity, which was never in the Euler equations at all.

Neither is small. The anisotropic pressure Hessian is the whole of the non-local coupling; it is what tells a parcel that its neighbours are being squeezed too, and that the fluid cannot collapse onto a sheet at one point without something giving way around it. Measured in simulations of real turbulence — not computed here — it is comparable in size to the $\mathbf A^2$ term and acts against it along the attracting branch, which is how the same local tendency produces bounded gradients and a characteristic pattern in the $(R, Q)$ plane rather than a singularity.

That pattern is the famous teardrop: the joint distribution of $Q$ and $R$ in a turbulent flow is concentrated in two quadrants — spinning with $R < 0$, straining with $R > 0$ — and has a long tail stretched out along exactly the right-hand branch of Vieillefosse’s line. The tail is the restricted equation’s attractor, visible in data; its finite length is the pressure and the viscosity cutting the runaway off. Nothing here computes the teardrop, which needs a turbulent field, and it is described rather than drawn for that reason.

So the calculation earns its place by being wrong in a specific way. It gets the direction of the local runaway right — sheets, vorticity on the middle axis, a positive intermediate strain — and it gets the outcome completely wrong, because the outcome is decided by a term it does not have. That is a sharper statement than “the pressure matters”: it says which way the pressure has to push and where in the plane it has to push hardest.

What the picture cannot show #

The phase-plane figure shows every trajectory at once and hides their speeds. Two points on the same curve can be seconds or ages apart, and the thin curves near the origin are traversed so slowly that a trajectory can spend most of its life on a few millimetres of the drawing. The ensemble figure is the correction: it is the same plane measured in time rather than in position.

Nor does any figure here show a flow. Each point is a single parcel’s gradient, with no neighbours and no place, and the singularity is a statement about that parcel’s matrix rather than about a velocity field. A velocity field whose gradient did this at one point would have to do something drastic everywhere around it, and the restricted equation cannot say what, because saying so is precisely the non-local pressure it has discarded. Whether the full Euler equations ever produce a singularity from smooth data is open — the stretching argument states the criterion such a singularity would have to meet — and nothing here bears on it.

The convention the result depends on #

Three choices sit under every number above, and each is named rather than assumed.

The gradient is $A_{ij} = \partial u_i/\partial x_j$, row index for the velocity component. The transpose convention exists in the literature and flips the sign of the vorticity read off the skew part; it does not change $Q$ or $R$.

Time is measured in units of $1/|\mathbf A_0|$, with $|\mathbf A_0|$ the Frobenius norm of the starting gradient. The equation has no other scale — replacing $\mathbf A$ by $\lambda\mathbf A$ and $t$ by $t/\lambda$ leaves it unchanged — so every time quoted is a pure number, and a real gradient of $1000\ \mathrm{s^{-1}}$ would reach its singularity a thousand times sooner than one of $1\ \mathrm{s^{-1}}$.

And the pressure kept is the isotropic part and only that. A model that kept some of the anisotropic part — there are several, and they are the question this essay ends on — is a different equation with a different answer.

Who found it, and when #

The closed pair for $Q$ and $R$ and the conserved discriminant are Vieillefosse’s, in two papers of 1982 and 1984, written to ask whether the local dynamics of an ideal fluid contained a runaway. Cantwell wrote the full matrix solution in 1992 and drew the phase plane in the form used here. The observation that turbulent vorticity aligns with the intermediate strain axis is Ashurst, Kerstein, Kerr and Gibson’s, from direct simulation in 1987, and it was made before the restricted equation was widely connected to it.

The $(R, Q)$ plane as a way of classifying local flow patterns belongs to Chong, Perry and Cantwell, from 1990, and is the same classification the critical points in space use. The observation that a turbulent field fills a teardrop in that plane came from the simulations of the 1990s, and it is the reason this equation, which is wrong about the outcome, stayed in use: it is the only model in which the teardrop’s tail has an explanation that can be written down.

Still open: what the discarded pressure has to supply #

The restricted equation fails by blowing up, and the natural question is how little of the missing pressure is needed to stop it. The models that answer it keep the gradient’s local closure and add back a stated approximation to the anisotropic Hessian — a pressure that remembers how the parcel’s neighbourhood has been deformed over a recent time window, or one built from a small tetrahedron of parcels rather than one — together with a viscous term expressed the same way.

Each such model is a closure, in exactly the sense the Reynolds-averaged equations need one, and each can be tested on whether it reproduces the teardrop. The calculation that would follow from here is the simplest of them: the deformation-history closure, which replaces the missing pressure by what an initially isotropic Hessian becomes after a fixed decorrelation time of deformation. The question it answers is whether a single remembered time is enough to turn the runaway of the figures above into a bounded, stationary distribution with a tail along the right branch — and if so, how long the memory has to be.

Beside it is the question the sign asymmetry raises about the tubes that are seen. The stretching rates of a real flow are not one number, and the tubes in simulated turbulence are the most intense structures in it, while the local dynamics computed here favour sheets. One account has the sheets rolling up into tubes by the instability that turns a shear layer into vortices, which is not a local process at all. How much of a real flow’s most intense vorticity is made that way, rather than by direct stretching, is a question about the non-local part, and it is where this equation stops being able to help.