On the blow-up of some complex solutions of the 3-d Navier-Stokes Equations: Theoretical Predictions and Computer simulations

TL;DR

This paper investigates complex solutions to the 3D Navier-Stokes equations that blow up in finite time, combining theoretical analysis and computer simulations to reveal detailed blow-up behaviors and energy concentration phenomena.

Contribution

It identifies two distinct types of blow-up solutions and provides detailed simulation results showing energy concentration and new features not previously predicted.

Findings

Two types of blow-up solutions with different divergence rates

Energy and enstrophy concentrate around singular points

Simulations reveal features beyond current theoretical predictions

Abstract

We consider some complex-valued solutions of the Navier-Stokes equations in $R^{3}$ for which Li and Sinai proved a finite time blow-up. We show that there are two types of solutions, with different divergence rates, and report results of computer simulations, which give a detailed picture of the blow-up for both types. They reveal in particular important features not, as yet, predicted by the theory, such as a concentration of the energy and the enstrophy around a few singular points, while elsewhere the "fluid" remains quiet.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\sum _{j=1}^3{u}_j\frac{\partial }{\partial x_j}\mathbf{u}=\Delta \mathbf{u}-\nabla p,\mathbf{x}=(x_1,x_2,x_3)\in {\mathbb{R}}^3.$$

$$\nabla \cdot \mathbf{u}=0,\mathbf{u}(\cdot ,0)={\mathbf{u}}_0.$$

$$\mathbf{v}(\mathbf{k},t)=\frac{i}{(2\pi {)}^3}{\int }_{{\mathbb{R}}^3}\mathbf{u}(\mathbf{x},t)e^{i⟨\mathbf{k},\mathbf{x}⟩}d\mathbf{x},\mathbf{k}=(k_1,k_2,k_3)\in {\mathbb{R}}^3,$$

$$\mathbf{v}(\mathbf{k},t)=e^{-t{\mathbf{k}}^2}{\mathbf{v}}_0(\mathbf{k})+$$

$$+{\int }_0^t{e}^{-(t-s){\mathbf{k}}^2}ds{\int }_{{\mathbb{R}}^3}⟨\mathbf{v}(\mathbf{k}-\mathbf{k}{}^{\prime },s),\mathbf{k}⟩P_{\mathbf{k}}\mathbf{v}(\mathbf{k}{}^{\prime },s)d\mathbf{k}{}^{\prime },$$

$$P_{\mathbf{k}}\mathbf{v}=\mathbf{v}-\frac{⟨\mathbf{v},\mathbf{k}⟩}{{\mathbf{k}}^2}\mathbf{k}.$$

$$E(t)=\frac{1}{2}{\int }_{{\mathbb{R}}^3}|\mathbf{u}(\mathbf{x},t){|}^{2}d\mathbf{x},S(t)={\int }_{{\mathbb{R}}^3}|\nabla \mathbf{u}(\mathbf{x},t){|}^{2}d\mathbf{x}={\int }_{{\mathbb{R}}^3}|\omega (\mathbf{x},t){|}^{2}d\mathbf{x}.$$

$$E(t)+{\int }_0^{t}S(\tau )d\tau =E(0),$$

$${\mathbf{v}}_A(\mathbf{k},t)=A{e}^{-t{\mathbf{k}}^2}{\mathbf{v}}_0(\mathbf{k})+{\int }_0^t{e}^{-{\mathbf{k}}^2(t-s)}\sum _{p=2}^{\infty }{A}^p{\mathbf{g}}^{(p)}(\mathbf{k},s)ds.$$

$${\mathbf{g}}^{(2)}(\mathbf{k},s)={\int }_{{\mathbb{R}}^3}⟨{\mathbf{v}}_0(\mathbf{k}-\mathbf{k}{}^{\prime }),\mathbf{k}⟩P_{\mathbf{k}}{\mathbf{v}}_0(\mathbf{k}{}^{\prime })e^{-s(\mathbf{k}-\mathbf{k}{}^{\prime }{)}^2-s(\mathbf{k}{}^{\prime }{)}^2}d\mathbf{k}{}^{\prime },$$

$${\mathbf{g}}^{(p)}(\mathbf{k},s)=$$

$$={\int }_0^{s}d{s}_2{\int }_{{\mathbb{R}}^3}⟨{\mathbf{v}}_0(\mathbf{k}-\mathbf{k}{}^{\prime }),\mathbf{k}⟩P_{\mathbf{k}}{\mathbf{g}}^{(p-1)}({\mathbf{k}}^{{}^{\prime }},s_2)e^{-s(\mathbf{k}-\mathbf{k}{}^{\prime }{)}^2-(s-s_2){(\mathbf{k}}^{{}^{\prime }}{)}^2}d{\mathbf{k}}^{{}^{\prime }}+$$

$$+\sum _{\genfrac{}{}{0ex}{}{p_1+p_2=p}{p_1,p_2>1}}{\int }_0^{s}d{s}_1{\int }_0^{s}d{s}_2{\int }_{{\mathbb{R}}^3}⟨{\mathbf{g}}^{(p_1)}(\mathbf{k}-\mathbf{k}{}^{\prime },s_1),\mathbf{k}⟩\cdot$$

$$\cdot P_{\mathbf{k}}{\mathbf{g}}^{(p_2)}(\mathbf{k}{}^{\prime },s_2)e^{-(s-s_1)(\mathbf{k}-\mathbf{k}{}^{\prime }{)}^2-(s-s_2)(\mathbf{k}{}^{\prime }{)}^2}d\mathbf{k}{}^{\prime }+$$

$$+{\int }_0^{s}d{s}_1{\int }_{{\mathbb{R}}^3}⟨{\mathbf{g}}^{(p-1)}(\mathbf{k}-\mathbf{k}{}^{\prime },s_1),\mathbf{k}⟩P_{\mathbf{k}}{\mathbf{v}}_0(\mathbf{k}{}^{\prime })e^{-(s-s_1)(\mathbf{k}-\mathbf{k}{}^{\prime }{)}^2-s(\mathbf{k}{}^{\prime }{)}^2}d\mathbf{k}{}^{\prime }.$$

$$\mathbf{k}=p{\mathbf{k}}^{(0)}+\sqrt{p}\mathbf{Y},{\mathbf{h}}^{(p)}(\mathbf{Y},s)={\mathbf{g}}^{(p)}(p{\mathbf{k}}^{(0)}+\sqrt{p}\mathbf{Y},s),s_j=s\left(1-\frac{{\theta }_j}{p_j^2}\right),j=1,2.$$

$${\mathbf{h}}^{(p)}(\mathbf{Y},s)=p^{\frac{5}{2}}\sum _{\genfrac{}{}{0ex}{}{p_1+p_2=p}{p_1,p_2>\sqrt{p}}}\frac{1}{p_1^2{p}_2^2}{\int }_{{\mathbb{R}}^3}{P}_{{\mathbf{e}}_3+\frac{\mathbf{Y}}{\sqrt{p}}}{\mathbf{h}}^{(p_2)}\left(\frac{\mathbf{Y}{}^{\prime }}{\sqrt{1-\gamma }},s\right)\cdot$$

$$\cdot ⟨{\mathbf{h}}^{(p_1)}\left(\frac{\mathbf{Y}-\mathbf{Y}{}^{\prime }}{\sqrt{\gamma }},s\right),{\mathbf{e}}_3+\frac{\mathbf{Y}}{\sqrt{p}}⟩d\mathbf{Y}{}^{\prime }\left(1+o(1)\right),$$

$${\mathbf{h}}^{(p)}(\mathbf{Y},s)=\left(H_1^{(p)}(\mathbf{Y},s),H_2^{(p)}(\mathbf{Y},s),\frac{F^{(p)}(\mathbf{Y},s)}{\sqrt{p}a}\right),$$

$$Y_1{H}_1^{(p)}(\mathbf{Y},s)+Y_2{H}_2^{(p)}(\mathbf{Y},s)+F^{(p)}(\mathbf{Y},s)=\mathcal{O}(p^{-\frac{1}{2}}{a}^{-1}),$$

$${\mathbf{h}}^{(p)}(\mathbf{Y},s)=p(\Lambda (s){)}^p\prod _{j=1}^3{g}_{{\sigma }_j}(Y_j)\left(\mathbf{H}(\mathbf{Y})+{\delta }^{(p)}(\mathbf{Y},s)\right).$$

$$\mathbf{H}(\mathbf{Y})=\left(H_1(\mathbf{Y}),H_2(\mathbf{Y}),0\right),$$

$${\delta }^{(p)}(\mathbf{Y},s)=\left({\delta }_1^{(p)}(\mathbf{Y},s),{\delta }_2^{(p)}(\mathbf{Y},s),{\delta }_3^{(p)}(\mathbf{Y},s)\right)$$

$$g_1^{(2)}(\mathbf{Y})\mathbf{H}(\mathbf{Y})={\int }_0^{1}d\gamma {\int }_{{\mathbb{R}}^2}{g}_{\gamma }^{(2)}(\mathbf{Y}-\mathbf{Y}{}^{\prime })g_{1-\gamma }^{(2)}(\mathbf{Y}{}^{\prime })\mathcal{L}(\mathbf{H};\gamma ,\mathbf{Y},\mathbf{Y}{}^{\prime })\mathbf{H}\left(\frac{\mathbf{Y}{}^{\prime }}{\sqrt{1-\gamma }}\right)d\mathbf{Y}{}^{\prime }$$

$$\mathcal{L}(\mathbf{H};\gamma ,\mathbf{Y},\mathbf{Y}{}^{\prime })=-(1-\gamma {)}^{\frac{3}{2}}⟨\frac{\mathbf{Y}-\mathbf{Y}{}^{\prime }}{\sqrt{\gamma }},\mathbf{H}\left(\frac{\mathbf{Y}-\mathbf{Y}{}^{\prime }}{\sqrt{\gamma }}\right)⟩+$$

$$+{\gamma }^{\frac{1}{2}}(1-\gamma )⟨\frac{\mathbf{Y}{}^{\prime }}{\sqrt{1-\gamma }},\mathbf{H}\left(\frac{\mathbf{Y}{}^{\prime }}{\sqrt{1-\gamma }}\right)⟩.$$

$$H_j(\mathbf{Y})=\sum _{m_1,m_2=0}^{\infty }{\ell }_{m_1{m}_2}^{(j)}{\mathrm{He}}_{m_1}(Y_1){\mathrm{He}}_{m_2}(Y_2),j=1,2.$$

$${\mathbf{H}}^{(0)}=c{\mathbf{Y}}_{\perp },{\mathbf{Y}}_{\perp }=(Y_1,Y_2,0),$$

$${\mathbf{v}}^{(p_0,\alpha )}(\mathbf{k},t)=C\sum _{p=p_0}^{\infty }p{\sigma }_{\alpha }(p)g^{(3)}({\mathbf{Y}}^{(p)}){\mathbf{Y}}_{\perp }^{(p)}{\int }_0^t{e}^{-{\mathbf{k}}^2(t-s)}{\left(\frac{\Lambda (s)}{\Lambda (\tau )}\right)}^{p}ds,\alpha =I,II.$$

$$\ln \frac{\Lambda (s)}{\Lambda (\tau )}=-\kappa (\tau -s)(1+r(t-s)),\kappa =\frac{{\Lambda }^{\prime }(\tau )}{\Lambda (\tau )}$$

Full text

Dedicated to Yakov Grigorievich Sinai on the occasion of his 80-th anniversary 3-d Navier Stokes equations. Blow-up. Global regularity problem

1 Introduction

Since the beginning of the modern mathematical theory of the Navier-Stokes (NS) equations, with Jean Leray in 1934 [6], one of the main open questions is whether the solutions of the initial value problem in ${\mathbb{R}}^{3}$ with smooth initial data and in absence of external forces can become singular in a finite time (“global regularity problem”).

Leray believed that the singular solutions do exist, and are related to turbulence. Although we now have a developed theory of turbulence which has no connection to the NS singularities, it is clear that the existence of singular solutions could provide a deeper understanding of physical phenomena such as the fast concentration of energy in a small space region, as we see in hurricanes. We know in fact [12] that a loss of smoothness implies the divergence of the solution at some point of the physical space (and this is also shown by the Li-Sinai solutions). There is at present no effective model for such phenomena.

Much work has been devoted to the global regularity problem, both theoretical and by computer simulations. An important result was recently obtained by T. Tao [13], who proved a finite-time blow-up for modified NS equations which satisfy the energy conservation. The model is obtained by replacing the quadratic term in the NS equations by a suitable average, and is related to the so-called “dyadic” model of Katz and Pavlovic, for which a finite-time blow-up can also be proved [5]. Other theoretical results have been obtained for suitable modifications of the equations (see [3], [13] and references therein).

Li and Sinai [8] consider a class of complex-valued solutions of the classical initial value problem in the whole space ${\mathbb{R}}^{3}$, in absence of boundary conditions and external forces

[TABLE]

[TABLE]

where ${\mathbf{u}}=(u_{1},u_{2},u_{3}):{\mathbb{R}}^{3}\times[0,\infty)\to{\mathbb{C}}^{3}$ is the complex “velocity” field, $p$ is the (complex) pressure, $\mathbf{u}_{0}$ is the initial data, and the viscosity is $\nu=1$ (which is always possible by rescaling). If we introduce the modified Fourier transform

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the scalar product in ${\mathbb{R}}^{3}$, we are led, by a Duhamel formula, to the integral equation

[TABLE]

[TABLE]

where ${\mathbf{v}}_{0}$ is the transform of ${\mathbf{u}}_{0}$, and $P_{{\mathbf{k}}}$ is the solenoidal projector expressing incompressibility

[TABLE]

Li and Sinai [8] proved that there is a class of real solutions of the integral equation (3), corresponding to complex solutions of the NS equations (1), which become singular at a finite time (blow-up). Observe that we get a real NS solutions for odd initial data $\mathbf{v}_{0}$.

As we discuss below, the main feature of such solutions, which is of great help both for the theoretical analysis and for computer simulations, is that they have a rather simple structure in ${\mathbf{k}}$-space: their support is concentrated in a thin cone around a fixed direction.

The total energy $E(t)$ and the totale enstrophy $S(t)$ (which measures the intensity of the vorticity field $\omega(\mathbf{x},t)=\nabla\times\mathbf{u}(\mathbf{x},t)$)

[TABLE]

are related by the law of conservation of energy

[TABLE]

which however for complex solutions is not coercive. In fact, for the Li-Sinai solutions both $E(t)$ and $S(t)$ diverge as we approach the critical time.

We recall that for real solutions (see, e.g., [7]), if $\mathbf{u}_{0}$ is such that $E(0)$ and $S(0)$ are finite, there is a unique regular solution in a maximal interval $t\in[0,T_{c})$, where the critical time $T_{c}$ depends on $\mathbf{u}_{0}$. Moreover global regularity (i.e., $T_{c}=+\infty$) holds if $S(t)$ is finite for all $t>0$, a result which follows from an “enstrophy inequality” if $E(0)$ and $S(0)$ are small enough. Therefore in the real case, if $T_{c}$ is finite, $S(t)$ diverges as $t\uparrow T_{c}$, while $E(t)$ decreases by energy conservation, implying a transfer of energy to the high $|{\mathbf{k}}|$ modes as $t\uparrow T_{c}$.

The main purpose of the present paper is to obtain a better understanding of the relevant features of the complex Li-Sinai solutions, in particular of the peculiar “mechanism ” by which the high ${\mathbf{k}}$-modes are enhanced, leading to the blowup. As we discuss in the next section, such mechanism is inherited by the real solutions obtained by antisymmetrizing the initial data $\mathbf{v}_{0}$ of Li and Sinai, although it is unclear whether it is “strong enough” for a blow-up. In any case the results of the present paper can be applied to the study of those solutions, a program that is already under way, shedding light on a new class of (real) fluid flows.

The computer simulations are of great help, as they reveal important properties which are not, so far, predicted by the theory. They indicate, in particular, point-wise convergence as $t\uparrow T_{c}$ of the solutions in the physical ${\mathbf{x}}$-space, everywhere except at a few points (one point for solutions of type $I$ and two points for type $II$) around which the energy concentrates.

The computer simulations reported in the present paper required the implementation of a computational scheme for the integral equation (3) which can be used for a class of complex and real solutions of the NS equation, including solutions that blow up.

The plan of the paper is as follows. In §2 we briefly report the Li-Sinai theory, and we identify two different types of solutions with different behavior near the critical time. We show that the solutions converge point-wise in ${\mathbf{k}}$-space as $t\uparrow T_{c}$, and derive the divergence rates for $E(t)$ and $S(t)$ for the two types. In the following sections, §3, and §4, we report the results of computer simulations. Finally, §5 gives some technical details on the computation, and §7 is devoted to concluding remarks.

We would like to remark that we simulate the solution of the integral equation (3) for a $3$-d vector in ${\mathbb{R}}^{3}$, with support going away to infinity, a task which requires a grid of the order of $10^{9}$ points, an unusual challenge even for modern supercomputers. Reasonable results could be obtained only thanks to the relatively simple structure of the Li-Sinai solutions in ${\mathbf{k}}$-space.

Computer simulations of a blowup for the integral equation (3) were first reported in [1]. However, due to computational limitations, their results only give a qualitative description of the divergence of energy and enstrophy.

Complex-valued solutions which blow up in a finite time with a similar behavior in $k$-space have been found also for the Burgers equations [9] and other models [10]. Computer simulations of blow-up solutions for the two-dimensional Burgers equations, which are much easier to handle, are reported in [2].

2 The Li-Sinai complex solutions and their divergence rates at the blow-up

We begin by a briefly description of the Li-Sinai theory, which is an excellent guideline to the understanding of the main features of the solutions. We refer the reader to the paper [8] for the proofs and further details.

As we mentioned above, we consider real solutions of the equation (3), which in general correspond to complex solutions ${\mathbf{u}}({\mathbf{x}},t)$. The initial data ${\mathbf{v}}_{0}$ are localized around a point ${\mathbf{k}}^{{(0)}}$, at a certain distance from the origin. In our simulations we took ${\mathbf{k}}^{{(0)}}=(0,0,a)$ with $5\leq a\leq 25$, and the support of $\mathbf{v}_{0}$ in a circle with center ${\mathbf{k}}^{{(0)}}$ and radius $r<a$.

Multiplying the initial data ${\mathbf{v}}_{0}$ by a positive parameter $A$ and iterating the Duhamel formula we can write the solution as a power series

[TABLE]

Substituting into the equation, we see that the functions $\mathbf{g}^{(p)}(\mathbf{k},s)$ satisfy a recursive relation which reminds the BBGKY hierarchy of Statistical Physics. Setting $\mathbf{g}^{(1)}(\mathbf{k},s)=e^{-s{{\mathbf{k}}}^{2}}{\mathbf{v}}_{0}({\mathbf{k}})$ and

[TABLE]

we find for $p>2$ the recursive relation

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

If $C=\rm{supp}\;{\mathbf{v}}_{0}$, then, by iteration of the convolution, the support of $\mathbf{g}^{(p)}$ will be $\underbrace{C+\ldots+C}_{p\;times}$. As $C$ is around ${\mathbf{k}}^{{(0)}}=(0,0,a)$, the support of the solution extends along the $k_{3}$-axis.

By analogy with the theory of probability, where convolution corresponds to the distribution of a sum of random variables, we know that for large $p$ the support of $\mathbf{g}^{(p)}$ is around $p{\mathbf{k}}^{{(0)}}$, in a region with transversal dimensions of the order $\sqrt{p}$. Moreover if $p$ is large the terms of the sum for which $\max\{p_{1},p_{2}\}\leq p^{1\over 2}$ can be neglected, and the Gaussian densities give a significant contribution to the integrals only for $s_{1},s_{2}$ near the endpoint $s$. Therefore we introduce the new variables and functions

[TABLE]

Integrating over $\theta_{j},j=1,2$ and setting $\gamma={p_{1}\over p}$ we get

[TABLE]

[TABLE]

where $\mathbf{e}_{3}=(0,0,1)$. We write the components of ${\mathbf{h}}^{(p)}({\mathbf{Y}},s)$ in the form

[TABLE]

and as, by incompressibility, ${\mathbf{h}}^{(p)}$ is orthogonal to ${\mathbf{k}}=(\sqrt{p}Y_{1},\sqrt{p}Y_{2},pa+\sqrt{p}Y_{3})$, we get

[TABLE]

which shows that $F^{(p)}({\mathbf{Y}},s)$ is of finite order. Therefore $\mathbf{h}^{(p)}({\mathbf{Y}},s)$, is essentially transversal to the $k_{3}$-axis, and as $p\to\infty$, $P_{\mathbf{e}_{3}+{{\mathbf{Y}}\over\sqrt{p}}}\mathbf{h}^{(p_{2})}\to\mathbf{h}^{(p_{2})}$, i.e., the solenoidal projector in (9) tends to the identity.

The fundamental Ansatz is that for some set of initial data $\mathbf{v}_{0}$, when $p$ is large and $s$ in some interval of time, the recursive relation (9) has an approximate solution which is asymptotically of the form

[TABLE]

Here ${\Lambda}(s)$ is a positive function, which will be discussed below, $g_{\sigma}(x)={e^{-{x^{2}\over 2\sigma}}\over\sqrt{2\pi\sigma}}$ denotes the centered Gaussian density on ${\mathbb{R}}$, $\sigma_{1},\sigma_{2},\sigma_{3}$ are positive constants, $\mathbf{H}$ is a vector function independent of time, orthogonal to $\mathbf{e}_{3}$, and depending only on $Y_{1},Y_{2}$,

[TABLE]

and the remainder

[TABLE]

is such that $\delta^{(p)}({\mathbf{Y}},s)\to 0$ as $p\to\infty$.

It is important to observe that by the Ansatz (11), the function $\mathbf{h}^{(p)}({\mathbf{Y}},s)$ is proportional to a product of Gaussian functions, and the time dependence of its leading term is determined by the function $\Lambda(s)$.

In view of possible rescalings it is not restrictive to set $\sigma_{i}=1,i=1,2,3$. Inserting (11) into (9), treating $\gamma$ as a continuous variable, neglecting the remainders, and integrating over $Y_{3}$, one can see that $\mathbf{H}({\mathbf{Y}})$ is a solution of the integral fixed point equation

[TABLE]

where, by abuse of notation, $\mathbf{H}=(H_{1},H_{2})$ and ${\mathbf{Y}}=(Y_{1},Y_{2})$ are in ${\mathbb{R}}^{2}$, $g^{(2)}_{\sigma}({\mathbf{Y}})={e^{-{Y_{1}^{2}+Y_{2}^{2}\over 2\sigma}}\over 2\pi\sigma}$, and

[TABLE]

[TABLE]

The solutions, or “fixed points”, of the functional equation (13) are found by expanding $\mathbf{H}$ in Hermite polynomials ${\rm He}_{k},k=0,1,\ldots$, which are orthogonal with respect to the standard Gaussian,

[TABLE]

There are infinitely many fixed points, and there is a class ${\cal C}$ of them such that for $\mathbf{H}\in{\cal C}$ there is an open set of initial data ${\mathbf{v}}_{0}$ for which the solution satisfies the Ansatz (11) for $s\in S=[s_{-},s_{+}]$, where $S$ is a non-empty time interval. The open set is constructed by linearized stability analysis, and the proofs are based on the renormalization group method.

Following Li and Sinai we consider the behavior of the solution near the blow-up time for a fixed point which is proportional to the radial vector

[TABLE]

where $c>0$ is a suitable constant determined by the fixed point analysis, and from now on we use the pedex $\perp$ to denote vector components perpendicular to the $k_{3}$-axis.

A delicate analysis (see [8] and references therein) shows that the function $\Lambda(s)$ is differentiable and strictly increasing.

Observe that the recursive relation (9) is unchanged if we replace $\mathbf{h}^{(p)}$ with $(-1)^{p}\mathbf{h}^{(p)}$, and the fixed point equation (13) is also unchanged. Therefore in the fundamental Ansatz (11) we can replace the term $(\Lambda(s))^{p}$ by $(-1)^{p}(\Lambda(s))^{p}$.

We get two types of solutions, those for which $h^{(p)}$ is given by (11), which will be called solutions of type $I$, and those for which $h^{(p)}$ is replaced by $(-1)^{p}h^{(p)}$, which will be called solutions of type $II$. Observe that if the initial data ${\mathbf{v}}_{0}$ is chosen according to the prescriptions in [8] and leads to a solution of type $I$ the fixed point $\mathbf{H}^{{(0)}}$, then the initial data $-{\mathbf{v}}_{0}$ lead to a solution of type $II$ with the same fixed point. As we shall see below, the two types behave rather differently.

We now set $A={1\over\Lambda(\tau)}$, for $\tau\in S$, and consider the solution (6) for $t<\tau$. Taking into account (15), setting $\sigma_{I}(p)=1$, $\sigma_{II}(p)=(-1)^{p}$, $p=1,2,\ldots$, the tail of the series in equation (6) can be replaced for both types by its asymptotics

[TABLE]

Here we write for clarity ${\mathbf{Y}}^{(p)}={{\mathbf{k}}-p{\mathbf{k}}^{{(0)}}\over\sqrt{p}}$, $p_{0}$ is large enough for the asymptotic behavior to set in, $C$ is a constant, and $g^{(3)}$ is the three-dimensional standard Gaussian density. The explicit asymptotics (16) shows that the main support of the solution extends along the direction $\mathbf{k}^{{(0)}}$, i.e., along the positive $k_{3}$-axis in a thin cone of transversal diameter proportional to $\sqrt{k_{3}}$.

As $\Lambda$ is strictly increasing, ${\Lambda(s)\over\Lambda(\tau)}<{\Lambda(t)\over\Lambda(\tau)}<1$, and the terms of the series (16) fall off exponentially fast in $p$, with an exponential rate which vanishes as $t\uparrow\tau$. Therefore the critical time is $T_{c}=\tau$, and, as the single terms of the series (6) tend to a finite limit as $t\uparrow\tau$, any divergence is due to the tail.

We now pass to a discussion of the main properties of the blow-up, based on the explicit asymptotic expression of the tail series (16). Although the results that follow are new, we will only give a sketch of the proofs. Full proofs are straighforward, but would require taking care of the corrections, which at some points is a rather cumbersome procedure. In what follows we may denote by ${\rm{const}}$ different constants.

As a first step we give a more convenient representation of the tail series (16), which is valid when $t$ is close to $\tau$. (Observe that it is not restrictive to assume that the initial time is close to $\tau$.) Since $\Lambda(s)$ is differentiable, we have, as $s\uparrow\tau$

[TABLE]

where $\kappa>0$ and $r(s)\to 0$ as $s\to 0$. Neglecting the remainder, and integrating over $s$, if $p_{0}$ is large, we get as $t\uparrow\tau$, for ${\alpha}=I,II$,

[TABLE]

Let now ${\mathbf{v}}^{({\alpha})}_{A(\tau)}({\mathbf{k}},t)$, ${\alpha}=I,II$, be the solution expressed by the series (6), where $A=A(\tau)={1\over{\Lambda}(\tau)}$, with $\tau>0$ and the initial data $\pm{\mathbf{v}}_{0}$ leading to the fixed point $\mathbf{H}^{{(0)}}$ chosen as explained above (and in more detail in [8]). We first show that, as $t\uparrow\tau$, there is a point-wise limit, i.e., for any fixed ${\mathbf{k}}$

[TABLE]

In fact the series (18) is bounded in absolute value by

[TABLE]

and the terms for $p>{\mathbf{k}}^{2}$ are bounded by $c\;p^{-{1\over 2}}\;e^{-{pa^{2}\over 4}}$, where the constant $c$ depends on ${\mathbf{k}}$. Hence the series (18) converges absolutes, and so does the original series (6).

Another point which was not considered in the previous literature is that of the rate of divergence of the totale energy $E(t)$ and the totale enstrophy $S(t)$. We will now show that the divergence rates are different for solutions of the two types. More precisely, as $t\uparrow\tau$, we have

[TABLE]

where $\beta_{I}=1$, $\beta_{II}={1\over 2}$ and $C^{({\alpha})}_{E},C^{({\alpha})}_{S}$ are constants depending on the initial data.

By what we said above, the divergence rate is determined by the asymptotic tail $\tilde{\mathbf{v}}^{(p_{0},{\alpha}}$ in (18). For solutions of type $I$, as $\langle{\mathbf{Y}}^{(p_{1})}_{\perp},{\mathbf{Y}}^{(p_{2})}_{\perp}\rangle={k_{1}^{2}+k_{2}^{2}\over\sqrt{p_{1}p_{2}}}$, and setting $R_{p}({\mathbf{k}})=\sqrt{p}{g^{(3)}({\mathbf{Y}}^{(p)})\over|{\mathbf{k}}|^{2}+\kappa p}$ we have

[TABLE]

The product $R_{p}({\mathbf{k}})R_{p+j}({\mathbf{k}})$ has a factor $\exp\{-{1\over 2}(|{\mathbf{Y}}^{(p+j)}|^{2}+|{\mathbf{Y}}^{(p)}|^{2})\}$, and writing ${\mathbf{Y}}^{(p+j)}$, $j=0,1,\ldots$, in terms of ${\mathbf{Y}}:={\mathbf{Y}}^{(p)}$ we find

[TABLE]

so that for the integrals of the terms of the series (21) we get

[TABLE]

Let $j_{0}(p):=[p^{{1\over 2}+{\epsilon}}]$, ${\epsilon}\in(0,{1\over 2})$, where $[\cdot]$ is the integer part. Due to the factor $\exp\{-{1\over 2}{j^{2}a^{2}\over 2p+j}\}$ we have

[TABLE]

Observe that due to the Gaussian factor, the integral (23) in the region $|\tilde{\mathbf{Y}}^{p,j}|>p^{{\epsilon}\over 2}$ falls off faster than any inverse power of $p$, and if $j\leq j_{0}(p)$ and $|\tilde{\mathbf{Y}}^{p,j}|\leq p^{{\epsilon}\over 2}$, setting $s_{j}={j\over\sqrt{p}}$, we have

[TABLE]

For $p$ large $\tilde{\mathbf{Y}}^{p,j}\sim\tilde{{\mathbf{Y}}}^{(s_{j})}:={\mathbf{Y}}-{s_{j}\over 2}{\mathbf{k}}_{0}$, and observing that $|{\mathbf{k}}_{\perp}|^{2}=p|\tilde{\mathbf{Y}}_{\perp}^{(s_{j})}|^{2}$, and changing the integration variable to $\tilde{\mathbf{Y}}^{(s_{j})}$, we get the following asymptotics, uniformly in $j=0,1,\ldots,j_{0}(q)$,

[TABLE]

Finally, as $s_{j}\leq s_{j_{0}(q)}\leq p^{{\epsilon}}$, we have, as $p\to\infty$,

[TABLE]

Therefore the sum $\sum_{j=1}^{j_{0}(p)}e^{-\kappa j(\tau-t)}I_{p,j}$ tends to a constant for large $p$, and gives the main contribution to the energy. In conclusion for $t\uparrow\tau$ we have

[TABLE]

[TABLE]

which proves the first asymptotics (20) for solutions of type I.

For $\tilde{\mathbf{v}}^{(p_{0},II)}$ the series in (18) has alternating signs, and (21) is replaced by

[TABLE]

[TABLE]

The contribution of the diagonal terms to $E(t)$ is the same as for the solutions of type $I$, and, taking into account the asymptotics (26), we have as $t\uparrow\tau$,

[TABLE]

For the off-diagonal terms, observe that, as for type $I$, the sum for $j>j_{0}(p)$ gives a negligible contribution for large $p$. Furthermore, we have

[TABLE]

[TABLE]

Writing ${\mathbf{Y}}^{p+\ell-1},{\mathbf{Y}}^{p+\ell}$ in terms of ${\mathbf{Y}}^{p}=:{\mathbf{Y}}$ we see that

[TABLE]

For $j\leq j_{0}(q)$, taking into account (31), where we neglect the term of order ${\cal O}(\tau-t)$, we need to compute the integrals

[TABLE]

where the exponential in the definition of $F_{p,2j}$ is written in the form (32). We then change the integration variable to $\tilde{\mathbf{Y}}^{(s_{2j-1})}$, and observe that in the region $|\tilde{\mathbf{Y}}^{(s_{2j-1})}|\leq p^{{\epsilon}\over 2}$, neglecting terms of the order ${\cal O}(p^{-1+2{\epsilon}})$, we have

[TABLE]

The term $\tilde{Y}_{3}^{(s_{2j-1})}$ given no contribution to the integral by parity, so that for large $p$ $J_{p,j}\sim{a^{2}\over 2\sqrt{p}}s_{2j-1}I_{p,2j-1}$, where $I_{p,2j-1}$ is given by (26). Summing over $j$, as in (27), we find

[TABLE]

where the constants are positive. Proceeding as for the diagonal part (30) we see that the off-diagonal terms give a contribution which is asymptotically of order $(\tau-t)^{-{1\over 2}}$.

The first relation (20) is proved.

For the second relation (20) it is enough to observe that the enstrophy density in ${\mathbf{k}}$-space is proportional to the energy density multiplied by a factor $|{\mathbf{k}}|^{2}$, which, taking into account the asymptotics (25) produces an additional factor $p^{2}$ in the analogues of the integrals $I_{p,j}$.

As a final remark for this section, we would like to point out that, as shown by the series representation (6), the Li-Sinai solutions describe a peculiar mechanism of enhancing the high ${\mathbf{k}}$-modes, which is due to convolutions of modulated Gaussian terms with centers on the $k_{3}$-axis at the points $k_{3}=pa$, $p\in{\mathbb{Z}}_{+}$.

It is easy to see that the a similar mechanism is going to work also for the real solutions obtained by antisymmetrizing the Li-Sinai initial data, although it is not clear how effective it will be, because the Gaussian centers will be a sum of positive and negative terms and tend to be closer to the origin.

3 Computer Simulations: the Li-Sinai solutions in ${\mathbf{k}}$-space

We simulate the integral equation (3) on a regular mesh in ${\mathbf{k}}$-space contained in the region $R=[-127,127]\times[-127,127]\times[-19,L]$. The parameter $L$, as we describe below, is of critical importance, and for the simulations reported below it takes the values $2028,2528,3028$. In what follows we will only indicate the simulation range of the longitudinal variable $k_{3}$. More details are given in §6.

The initial data ${\mathbf{v}}_{0}$ are chosen according to the prescriptions in [8] (§7, formula (39)), with support concentrated around the point $\mathbf{k}^{(0)}=(0,0,a)$, $a>0$, and leading to the fixed point (15). A full screening for the “best” cases with a large random choice of the initial parameters, as in our paper on the 2-d Burgers equations [2], was not possible because for the 3-d NS equations it takes too much computer time. We only considered about a hundred cases and followed up the most promising ones.

For all results reported in this paper the initial data are of the following form

[TABLE]

where $a=20$, ${\mathbb{I}}_{D}$ is the indicator function of the support $D=\{{\mathbf{k}}:|{\mathbf{k}}|\leq 17\}$, and the positive constant $C$ controls the initial energy and enstrophy. The initial data ${\mathbf{v}}_{0}^{+}$ lead to solutions of type $I$, and the initial data ${\mathbf{v}}_{0}^{-}$ to solutions with oscillating sign of type $II$.

We studied in more detail the solutions of type II, as their behavior is more similar to that of the real solutions obtained by antisymmetrizing the initial data, which is an object of our present research.

In describing the blow-up an important role is played, due to the structure of the solutions, by the marginal densities for energy and enstrophy along the $k_{3}$ axis in ${\mathbf{k}}$-space

[TABLE]

[TABLE]

By our definition of the transform (2) the total energy and enstrophy in (4) are given by

[TABLE]

The marginal densities along the third axis in ${\mathbf{x}}$-space are

[TABLE]

[TABLE]

We also consider the transverse marginals $E_{j}(k_{j},t)$, $\tilde{E}_{j}(x_{j},t)$, $S_{j}(k_{j},t)$, $\tilde{S}_{j}(x_{j},t)$, $j=1,2$. which are defined in an obvious way.

For the initial data (35), if the constant $C$ is large enough, beyond some critical value corresponding to the initial energy to $E_{0}>25,000$, the solution blows up after a time of the order $10^{-4}$ time units. However, as in our screening we never went beyond a time of the order $10^{-2}$, it may well be that some of the cases with lower initial energy do in fact blow up at a later time. The question whether the critical value is real remains open.

The qualitative behavior of the solution does not change when we increase the constant $C$ in (35) beyond the (apparent) critical value. Therefore we only report results with initial data (35) where $C$ corresponds to the initial energy is $E_{0}=200\times(2\pi)^{3}\approx 49,500$. That is, we report results of only two solutions, of type $I$ and $II$.

The rapid growth of the energy and enstrophy takes place in a very short time with respect to the total run time, as shown in Fig. 1 for solutions of type $I$ and Fig. 2 for solutions of type $II$. Observe also that the enstrophy starts growing earlier than the energy, and the critical time is much smaller for the solution of type $I$, although the initial energy and enstrophy are the same.

The concentration of the support of the solutions around the $k_{3}$-axis and the role of the fixed point $\mathbf{H}^{{(0)}}$ are already visible at the beginning of the blow-up and for relatively small values of $k_{3}$. As shown in Fig. 3 for the solution of type $II$, the direction of the vector field on a plane orthogonal to the $k_{3}$-axis is approximately radial, and the absolute value $|{\mathbf{v}}({\mathbf{k}},t)|$ falls off sharply as we move away from the $k_{3}$-axis.

The different behavior of the solutions of type $I$ and $II$ is illustrated in Fig. 4, which shows a local plot of the marginal enstrophy density $S_{3}(k_{3},t)$. For the solution of type $II$ we see sharp peaks, approximately at the points $k_{3}=pa$, $p=1,2,\ldots$, corresponding to the centers of the Gaussian factors in the expansion (18), with zeroes between them, due to interference of neighboring peaks with different sign. The solution of type $I$ shows instead very mild peaks, and only for low values of $p$.

As we shall see in the next section, the different behavior in ${\mathbf{k}}$-space implies a different location of the singular point in physical ${\mathbf{x}}$- space.

The main difficulty in following the blow-up by computer simulations is that, as we get close to the critical time, the essential support of the solutions moves away along the $k_{3}$-axis. The divergence of the energy could not be followed in a proper way, because even in the initial times of the blow-up its growth is due to modes with $k_{3}$ exceeding the maximal value available to us $L=3025$.

The divergence of the enstrophy is easier to follow, because, as we mentioned above, it starts growing much earlier than the energy, when the range of the values of $k_{3}$ which contribute to the growth is contained within the simulation region $[-19,L]$, for $L\geq 2028$.

Fig. 5 shows the behavior in time of the enstrophy marginal density $S_{3}(k_{3},t)$ along the $k_{3}$-axis. It can be seen that in the range of times under consideration there is a significant increase of the enstrophy due to the modes with $k_{3}$ within the simulation range. The transversal distribution of the enstrophy given in Fig. 6, for type $II$ and the same times as in Fig. 5, shows that the relevant support of the enstrophy grows very slowly in the transversal direction and is well contained inside the computation region $[-127,127]\times[-127,127]$.

Fig. 5 also shows that, as we approach the critical time, the main support of $S_{3}(k_{3},t)$ moves away to the high $k_{3}$ region, while the function for low $k_{3}$ values, to the left of the growing maxima, changes very little. This behavior is in accordance with the point-wise convergence in ${\mathbf{k}}$-space, which was deduced in Section 2.

The rate of convergence of the solution for low $k_{3}$ values is given for type $II$ in Fig. 7, where a transversal component of ${\mathbf{v}}({\mathbf{k}},t)$ is reported as a function of $k_{3}$, for $k_{1},k_{2}$ fixed, and several values of $t$. Observe that the behavior is a kind of damped oscillation with period $2a$, and zeroes at the points $k_{3}\approx{1\over 2}(2j+1)\;a$, $j=1,\ldots$. The other transversal component behaves in the same way.

The problem of estimating the critical time $\tau$ is rather challenging because the divergence is due to the large $k_{3}$ modes, which act in a different way on the enstrophy and energy. The most reliable way of estimating $\tau$ is that based on formula (18), which indicates that, except for power-law corrections, the absolute value of the solution $|{\mathbf{v}}_{A(\tau)}({\mathbf{k}},t)|$ for ${\mathbf{k}}_{\perp}$ fixed, falls off exponentially fast in $k_{3}$ with a rate proportional to $(\tau-t)$.

The same property extends to the marginal energy density $E_{3}(k_{3},t)$, and indeed, as shown by Fig. 8 for type $II$, it is verified for $k_{3}\geq 400$ already at times relatively far from the critical time. The points align on a straight line with great accuracy, and the slope decreases with time. (The dots represent the local maxima, as the marginal density $E_{3}$, in analogy with the enstrophy density $S_{3}$ shown in Figg. 4, is wildly oscillating.) Moreover the exponential decay rate is remarkably stable with respect to the longitudinal simulation range for $L\geq 2028$, for solutions of both types.

If we now plot the exponential decay rate vs time, as shown on Fig.9 for both types, we get an estimate of the critical time by looking at the intercept with the horizontal axis. The estimates are $\tau\approx 1110\times 10^{-7}$ for type $I$ and $\tau\approx 1630\times 10^{-7}$ for type $II$. Observe that all estimates are to be considered as overestimates, because the simulations do not take into account the contribution of the modes outside the simulation region, which enhances the divergence.

We conclude this section with a check of the rates of divergence (20) predicted in Section 2. We restrict our considerations to the divergence of the total enstrophy, because, as we explained above, the data on the growth of the energy that we could collect are not good enough. In Fig.10 we report a plot of $\log S(t)$ versus $\log{1\over\tau_{*}-10^{7}\;t}$, where $\tau_{*}$ is the estimated critical time in magnified units ($\tau_{*}=1110$ for type $I$, and $\tau_{*}=1630$ for type $II$) in the range of times which are sufficiently large and within the reliable range. The results should be compared with the prediction of (20), which is a power ${\alpha}=3$ for type $I$ and ${\alpha}=2.5$ for type $II$. Taking into account that that the slopes in Fig. 10 are very sensitive to the value of the predicted critical time, we conclude that the results, namely ${\alpha}\approx 3,12$ for type $I$ and ${\alpha}\approx 2,6$ for type $II$, are compatible with the predictions in (20).

More computer resources are needed in order to get a deeper understanding of the role of the high $k_{3}$ modes on the divergence of the total energy and the total enstrophy as $t\uparrow\tau$.

4 Computer simulations: the Li-Sinai solutions in $\mathbf{x}$-space

In ${\mathbf{x}}$-space the main support of the solution is contained within a small volume around the origin, and as $t\uparrow\tau$, it concentrates in sharp “spikes”, which become singular in the limit.

As shown in Fig. 7, the solution of type $II$ oscillates in $k_{3}$ with a period $T\approx 2a$ (with $a=20$), so that we expect large values of the velocity ${\mathbf{u}}({\mathbf{x}},t)$ around the points $\mathbf{x}^{{(0)}}_{\pm}=(0,0,\pm x_{3}^{{(0)}})$ with $x^{{(0)}}_{3}\approx 0.16\approx{\pi\over a}$. The behavior of the marginal energy densities $\tilde{E}_{3}(x_{3},t)$ and $\tilde{E}_{1}(x_{1},t)$ for two different times given in Fig. 11 and Fig. 12. It can be seen that $\tilde{E}_{3}(x_{3},t)$ shows two spikes at the points $\pm x_{3}^{{(0)}}$, while the transversal marginal $\tilde{E}_{1}(x_{1},t)$ has a single spike at the origin. Both spikes increase indefinitely as $t\uparrow\tau$, while everywhere else both marginals seem to converge to a finite limit.

A possible conclusion, which however would require further investigation, is that for all ${\mathbf{x}}\neq\mathbf{x}^{{(0)}}_{\pm}$ the function ${\mathbf{u}}({\mathbf{x}},t)$ itself also converges. A result of this kind was in fact proved for the complex Burgers equations (see [9] and [2] for computer simulations).

The marginal $\tilde{E}_{3}(x_{3},t)$ for the solution of type $I$, as given in Fig. 13, also for two different times, differs from the previous one, in that it has a single spike, suggesting that the only singularity of ${\mathbf{u}}({\mathbf{x}},t)$ as $t\uparrow\tau$ is at the origin. Here also the marginals indicate divergence at ${\mathbf{x}}=0$ and convergence elsewhere.

An important role with respect to possible singularities of the NS equations is played by the vorticity stretching vector $\mathbf{w}({\mathbf{x}},t)=\omega({\mathbf{x}},t)\cdot\nabla\mathbf{u}(\mathbf{x},t)$, where $\omega({\mathbf{x}},t)=\nabla\times{\mathbf{u}}({\mathbf{x}},t)$ is the vorticity (see, e.g, [11]). In Fig. 14 we show, for the solution of type $II$, a joint logarithmic plot of the enstrophy marginal $\tilde{S}_{3}(x_{3},t)$ and of the corresponding marginal for vorticity stretching

[TABLE]

($W_{3}$, for dimensional homogeneity, is divided by $E_{0}$). Observe that the spikes of the enstrophy marginal are more enhanced than for the energy, and those of $W_{3}$ even more so. The transverse marginals $\tilde{S}_{1}(x_{1},t)$, and the corresponding one for the vorticity stretching have a spike of similar magnitude at the origin.

5 Some information on the computer simulations

As we said above, we simulate the integral equation (3) in Fourier ${\mathbf{k}}$-space, and it is of great help the fact that the support of the solutions in ${\mathbf{k}}$-space is concentrated along the $k_{3}$-axis.

The discretization is implemented by a regular mesh of points containing the origin, and such that it contains the region where the solution is significantly non-zero. For all simulations reported in this paper the mesh in ${\mathbf{k}}$-space is taken uniform with step $1$, and is a set of the type $R=[-127,127]\times[-127,127]\times[-19,L]\subset{\mathbb{Z}}^{3}$ (the brackets $[\ldots]$ denote intervals in ${\mathbb{Z}}$), with $L=2028,2528,3028$. Control simulations with a refined mesh were performed to check stability, which showed that the results are remarkably stable with respect to refinements of the mesh, in accordance with the fact that the solution ${\mathbf{u}}({\mathbf{x}},t)$ in ${\mathbf{x}}$-space is essentially concentrated in a small region around the origin (see §5).

However, as we approach the critical time, the simulations are very sensitive to the time step, and to the strip length $L$. We always checked that the time step $\delta_{t}=10^{-7}$ and the given value of $L$ are such as to ensure stability in the range of times under consideration.

The discrete computation is obtained by the so-called Nyström method with respect to the ${\mathbf{k}}$ variables and a predictor-corrector scheme with respect to the time variable, which is based on the Euler method and the Trapezoidal method.

The predictor-corrector scheme iterates a recursive procedure to compute successive approximations $\mathbf{V}^{(j)}({\mathbf{k}},n\delta_{t})$, $j=1,2,\ldots$, of the function ${\mathbf{v}}({\mathbf{k}},n\delta_{t})$, until a convergence criterion is satisfied: $|\mathbf{V}^{(j+1)}({\mathbf{k}},n\delta_{t})-\mathbf{V}^{(j+1)}({\mathbf{k}},n\delta_{t})|\leq{\rm tol}$, where the tolerance ${\rm tol}$ is set at $10^{-8}$.

The procedure is computationally challenging, and in fact for each ${\mathbf{k}}$ in the mesh it requires the evaluation of a three-dimensional integral. The integral is a convolution, and by using the fast Fourier transform (FFT) we can reduce the computational cost. Note that, as we are simulating the NS equations in ${\mathbb{R}}^{3}$, the convolution is not periodic, and has to be implemented on a computational grid which is doubled in size.

The accuracy of the approximated solution corresponding to the chosen initial data is evaluated on an experimental basis by comparing the results obtained for different discretization parameters.

Our computer simulations were performed at CINECA of Bologna (Italy) on the FERMI Supercomputer (Model: IBM-BlueGene/Q; Architecture: $10$ BGQ Frame; Processor Type: IBM PowerA2, $1.6$ GHz; Computing Cores $163840$; Computing Nodes $10240$; RAM: 1GByt/core)

The computation method was implemented in Fortran 90 (IBM Fortran compiler) with MPI library for parallel computations, and 2Decomp&FFT for the parallel computation of the fast Fourier transform.

6 Concluding Remarks

We deduced in Section 2 some important consequences of the work of Li and Sinai [8] which clarify the behavior of the solutions near the blow-up. With the help of computer simulation it was then possible to check the predictions and to estimate the critical time. Moreover, what is perhaps more relevant, the computer simulations give a detailed picture of the behavior in ${\mathbf{x}}$-space, indicating important properties of the solutions, such as the point-wise convergence as $t\uparrow\tau$ in ${\mathbf{x}}$-space, except for the singular points, which are not so far predicted by the theory, but can hopefully be proved rigorously in the near future.

We would like to remark that simulations of the solutions of the 3-d NS equations are usually computationally onerous and sometimes unreliable, especially for flows with large values of the enstrophy (see [4] for a review). It is a remarkable fact that the singular complex solutions proposed by Li and Sinai, due to their simple structure in ${\mathbf{k}}$-space, are relatively easy to follow by computer simulations on the supercomputers of the last generation.

It is of particular importance in this respect that the solution can be represented as a power series (6), where the parameter $A$ governs the blow-up time. A great help also comes from the stability of the computation with respect to the discretization step in ${\mathbf{k}}$-space, which is due to the confinement of the energy in a small region of ${\mathbf{x}}$-space. In fact, in the original paper [8] the infinite extension of the domain in ${\mathbf{x}}$-space does not seem to be essential, except for the absence of boundary conditions, and the proofs can possibly be adapted to the periodic case on the torus $T^{3}$.

The general picture that comes out is that of a motion in which the “fluid” moves very fast in a small region around the origin for solutions of type $I$ or around two symmetric points close to the origin for solutions of type $II$, along flow lines with high curvature, as indicated by the behavior of the vorticity and the vorticity stretching described in §5.

By antisymmetrizing the initial data one gets real-valued solutions which share some basic properties of the complex solutions, such as the restriction of the support to a thin cone in $k$ space. The study of such solutions, which for some interval of time show strong similarities to solutions of type $II$, is in progress.

7 Acknowledgements

We thank Prof. Ya. G. Sinai and Dr. D. Li for their constant interest on our work, and for many discussions and suggestions. We acknowledge the CINECA award under the ISCRA initiative IsB10_3DNS (2014) and IsC23_3DNS (2015), for the availability of high performance computing resources and support.

Funding To S.F.: Ministero dell’Istruzione, dell’Università e della Ricerca, COFIN Mat/07 2012.