# Initial evolution of 3D turbulence en route to the Kolmogorov state:
emergence and transformations of coherent structures, self-similarity and
instabilities

Giorgio Krstulovic  Université Côte d’Azur, Observatoire de la Côte d’Azur,
CNRS, Laboratoire Lagrange, Boulevard de l’Observatoire CS 34229 – F 06304
Nice Cedex 4, France  Sergey Nazarenko  Université Côte d’Azur, CNRS, Institut
de Physique de Nice - INPHYNI, Parc Valrose, 06108 Nice, France

###### Abstract

In this work, we study numerically the temporal evolution of an initially
random large-scale velocity field under governed by the hyperviscous
incompressible Navier-Stoke equations. Three stages are clearly observed
during the evolution. First, the initial condition development is
characterized by a spectrum evolving in a self-similar way with a wave number
front $k_{*}(t)$ propagating toward high values exponentially in time. This
evolution corresponds to the formation and shrinking of think vortex pancakes
exponentially in time, as it has been previously reported in simulations of
the incompressible Euler equations. At the second stage, pancakes become
unstable, rolling up on the edges and breaking up in the middle, leading to
the emergence of _vortex ribs_ – quasi-periodic arrangements of vortex
filaments. Those filaments then twist creating structures akin to ropes. At
the last stage, a fully developed turbulent state is observed, characterized
by a Kolmogorov energy spectrum and exhibiting a decay law compatible with
standard turbulent predictions in the case of an integral scale saturated at
the scale of the box.

##  I Introduction

Traditionally, studies of hydrodynamic turbulence focus on statistically
steady states. These include testing the degree of validity of the Kolmogorov
theory and investigating deviations from this theory [1]. This is related to
the view that turbulence properties are most universal, and therefore most
interesting, when a stationary state is reached. However, non-stationary
turbulence is also interesting and important, and it deserves greater
attention. Particularly interesting is the temporal evolution en route to the
steady state. Here, some unexpected nontrivial behavior was discovered for
phenomenological descriptions of turbulence: the Leith and the eddy-damped
quasi-normal Markovian (EDQNM) models. The Leith model is a nonlinear singular
diffusion equation for the energy spectrum [2, 3], whereas the EDQNM model is
a nonlinear integro-differential equation [4, 5]. It was found, first for the
Leith model and later for the EDQNM model, that turbulence initially excited
in a finite range of scales develops a power-law spectrum with an exponent
$x_{*}$ which is different from Kolmogorov’s $-5/3$ behind the front
propagating toward higher wave numbers. Namely, for the Leith model $x_{*}$
was found to be $\approx-1.85$ [3, 5] and for the EDQNM closure $\approx-2$
[6]. It is interesting that such an exponent is model-dependent but
independent of the shape of the initial spectrum. The Kolmogorov spectrum in
both of these models forms after time $t_{*}$ when the initial spectral front
reaches the dissipative scale: it propagates as a wave in the $k$-space
propagating from the dissipative scale toward the initial scale, i.e. from
large to small wave numbers. This evolution scenario is typical for all
finite-capacity turbulent systems, i.e. systems whose constant-flux
(Kolmogorov-type) spectra correspond to a finite physical-space density of the
cascading invariant (e.g. energy) when the dissipation coefficient (e.g.
viscosity) tends to zero. Besides the Leith model and EDQNM, this behavior was
also observed for the Magneto-Hydrodynamic (MHD) turbulence [7] and turbulence
described by the Burgers equation [8], and later for a great variety of Wave
Turbulence systems [9].

Analytical treatment is naturally easier for the differential (e.g. Leith)
models than for the integro-differential ones (e.g. EDQNM). In the Leith model
case, the pre-$t_{*}$ and the post-$t_{*}$ stages to the self-similar
solutions of the second and the third kinds respectively. In the self-
similarity of the second kind the similarity coefficients are found by solving
a nonlinear eigenvalue problem, whereas in the third kind self-similarity they
are fixed by the previous self-similar stage. It was shown that the anomalous
power-law exponent c for the pre-$t_{*}$ spectrum is related is a heteroclinic
bifurcation of a dynamical system associated with the self-similar class of
solutions. The post-$t_{*}$ stage (with the Kolmogorov scaling invading the
$k$-range from right to left) has the similarity coefficients inherited from
the pre-$t_{*}$. Similar behavior is typical for the integro-differential
models; for example it was demonstrated for the kinetic equation of the weak
wave MHD turbulence [7].

A revealing picture arises when analyzing a similar scenario in the Burgers
equation model. Since the Burgers equation is integrable, its analysis can be
done fully rigorously. It was found [8] that the pre-$t_{*}$ scaling is
associated with formation of the so-called pre-shock: a cubic-root singularity
arising in the physical at the wave breaking moment. The post-$t_{*}$ scaling
is related to the shock singularity giving the Burgers power-law spectrum
which is an analog of the Kolmogorov spectrum. It is tempting to extend such a
two-singularities scenario to the hydrodynamic turbulence. However, the
closure models which only involve the spectrum do not allow to study
singularities because they do not contain information about the phases of the
Fourier amplitudes of velocity. On the other hand, a two-singularity scenario
was indeed discussed in [10], who used pseudo-spectral direct numerical
simulations (DNS) of the Euler equations at resolution $256^{3}$ and observed
that pancake-like vorticity structures form at the initial (inviscid) stage of
the turbulent evolution. Appearance of the vortex pancakes, and accompanying
it power-law with exponent $\approx-4.5$, turned out to be a rather robust
process–it was observed in runs with initial conditions in the form of the
Taylor-Green vortex as well as in runs with random initial conditions. The
pancakes were observed to be catastrophically shrinking in thickness–a process
which should eventually be arrested by viscosity. They further suggested (but
not demonstrated numerically) that the pancakes should become unstable at the
moment when their thickness becomes comparable to the dissipative scale. Such
a instability should lead to formation of another type of quasi-singular
vortex structures–thin intense vortex filaments. Later in [11], pseudo-
spectral $256^{3}$ DNS of the Navier-Stokes equations was performed using the
same random-mode initial conditions as in [10]. Presence of viscosity and use
of a state of the art (for that time) visualization technique allowed them to
observe instability of the vortex pancakes leading to their breakup into
filamentary vortices. Also, the authors of [11] have suggested that the
prototype for the observed vortex pancakes are the so-called Burgers vortex
layers in which an interplay of an external strain and the viscosity results
in a stationary sheet-like vortex structure [12]. They performed a linear
analysis of instability of the Burgers vortex layer in a way which was
previously used to demonstrate instability of the mixing layers in [13]. Later
such a linear instability analysis was further developed in [14].

The process of formation of shrinking vortex pancakes in evolving turbulence
was recently reproduced numerically in high-resolution DNS of the Euler
equations: at resolution $486\times 1024\times 2048$ in [15] and $972\times
2048\times 4096$ – in [16]. Self-similarity in evolving pancakes was observed
and characterized. However, no break-up of the vortex pancakes into
filamentary was observed. Later, the same group of authors claimed that no
Kelvin-Helmholtz type instability is possible for the rapidly shrinking vortex
pancakes. These observations are in agreement with the suggestion that
including dissipation is essential for the vortex pancake instability because
it is dissipation that slows down the pancake collapse and makes them behave
like quasi-stationary Burgers vortex layers. It was also claimed in [15, 16]
that the shrinking vortex pancakes give rise to appearance of the Kolmogorov
-5/3 spectrum.

The present paper is devoted to a study of the initial evolution of
hydrodynamic turbulence by DNS of the Navier-Stokes equations with hyper-
viscosity using a pseudo-spectral method and random wave initial conditions
similar to the ones used in [10]. Our motivation is twofold: we would like to
check if the Fourier-space self-similarity (previously observed in the Leith
and the EDQNM closures) is also realized in DNS of hydrodynamic turbulence, as
well as to examine, at high resolution, realizability of the previously
conjectured two-singularity scenario in where the initial evolution first
leads to formation of shrinking vorticity pancakes which subsequently become
unstable and break up into thin quasi-singular vortex filaments.

We find that the initial evolution stage is characterized by a spectrum
evolving in a self-similar way with a front $k_{*}(t)$ propagating toward high
$k$’s exponentially in time. Note that the exponential shrinking of the
minimal size $\delta(t)\sim 1/k_{*}(t)$ was also reported in [10] and later in
[15] and [16]. We observe that the spectral front reaches the maximum
wavenumber in a finite time $t_{*}\approx 1.24$. On the other hand, at the
present available resolution we were not able to detect any front acceleration
which could be interpreted as a change from the exponentially shrinking regime
to a finite-time blowup behaviour. Thus, our tentative observation is that the
value of $t_{*}$ is resolution dependent and would increase if the inertial
range of wavenumbers was widened, and that no evidence for the finite time
singularity formation could be seen.

In the physical space, like previously reported in [10, 15] and [16], at the
early stage we see formation of vorticity pancakes whose thickness is
shrinking. After the minimal thickness have reached the dissipative scale, the
pancakes go unstable: they roll on the edges and break up in the middle parts
into ”vortex ribs”–quasi-periodic sets of intense vortex filaments. Thus, the
scenario of two consecutive (quasi-)singularities is indeed observed. Note
that the pancakes form at before $t_{*}$ and persist/grow well beyond $t_{*}$.
Their instability and breakup into quasi-periodic arrays of co-directed vortex
ribs starts much later, at $t\gtrsim 3.5t_{*}$. After forming, the vortex ribs
pull together and twist forming vortex “ropes”. The Kolmogorov spectrum does
indeed form at the later stage, $t\gtrsim 4.6t_{*}$, after the vortex ribs and
ropes become numerous and entangled, forming a rather chaotic vorticity field.
Note that, in contrast to what was claimed in [15, 16] the Kolmogorov spectrum
formation occurs much later than the shrinking pancaked form and disappear. In
fact, no spectrum numerical result for the energy spectrum was presented in
[15, 16] to support their claim about the Kolmogorov spectrum forming at the
shrinking pancake stage.

We believe that the mechanism of the pancake instability in the system with
hyperviscosity is similar to the Burgers layer instability in the system with
regular kinematic viscosity previously considered in [11]. Respectively, we
believe that the Navier-Stokes system with normal viscosity should behave in a
similar way to what we observe in the Navier-Stokes system with hyperviscosity
considered in our paper. Namely, the route to fully developed turbulence could
be roughly divided into three stages: (i) vortex pancake formation and
exponential shrinking (ii) pancake roll-ups at the edges and instability in
the middle parts leading to the vortex ribs, and (iii) twisting of vortex ribs
and forming vortex ropes with further entanglement of the rib and rope
structures and formation of a fully developed chaos in the vorticity field.
This latter stage is characterised by formation of a Kolmogorov scaling in the
turbulent energy spectrum.

##  II Mathematical model and numerical setup

We consider $3$D decaying solutions ${\bf u}({\bf x},t)$ of the incompressible
hyper-viscous Navier-Stokes equations:

$\displaystyle\frac{\partial{\bf u}}{\partial t}+{\bf u}\cdot\nabla{\bf u}$
$\displaystyle=$ $\displaystyle\nabla p-\nu(-\Delta)^{\alpha}{\bf u}$ (1)
$\displaystyle\nabla\cdot{\bf u}$ $\displaystyle=$ $\displaystyle 0.$ (2)

The total energy and total dissipation of the system are given by

$E(t)=\frac{1}{2V}\int|{\bf u}|^{2}\mathrm{d}{{\bf x}},\quad{\rm
and}\quad\epsilon(t)=\frac{\nu}{V}\int{\bf u}(-\Delta)^{\alpha}{\bf
u}\mathrm{d}{{\bf x}}=\nu\sum_{\textbf{k}}|\textbf{k}|^{2\alpha}|\widehat{{\bf
u}}(\textbf{k})|^{2},$ (3)

where $\widehat{{\bf u}}(\textbf{k})$ is the Fourier transform of ${\bf u}$.
They naturally obey the balance equation $\frac{dE}{dt}=-\epsilon(t)$.

We consider the evolution of a random large-scale initial condition, defined
as

$\widehat{u}_{i}(\textbf{k},0)=\sum_{|\textbf{k}|<k_{\rm
I}}\frac{\sqrt{C}}{k}e^{i\phi_{k}},\hskip 28.45274pti=1,2,3.$ (4)

The phases $\phi_{k}$ are chosen randomly and $C$ is a normalization constant
such that $E_{0}=E(t=0)$.

In the numerical simulations presented in this work we have set $\alpha=4$,
$k_{\rm I}=1.75$ and $E_{0}=1$. The cubic domain is of size $2\pi$. We adapted
the standard pseudo-spectral code GHOST [https://github.com/pmininni/GHOST]]
and used resolutions of $1024^{3}$ and $2048^{3}$ colocation points with the
hyper-viscous dissipation coefficient set to $\nu=4\times 10^{-19}$ and
$\nu=1.52\times 10^{-21}$, respectively. Both simulations are started from
exactly the same initial conditions (same large-scale Fourier modes). Run at
resolution $2048^{3}$ is evolved up to $t=2$ and it used to perform a fine
analysis of the energy spectrum at early stages. Run at resolution $1024^{3}$
goes up to $t=10$ and allows us to perform visualization of the enstrophy
field and to study later stages of the evolution.

The temporal evolution of energy and dissipation are shown in Fig. 1.

Figure 1: Temporal evolution of the energy and energy dissipation. Run at
resolution $1024^{3}$.

We observe three stages which will be reported in the following. In early
evolution (stydied in section III) dissipation plays no role, and the dynamics
is effectively driven by the Euler equation ($\nu=0$). As we will se later,
this phase ends when the propagating spectrum reach the viscous scale at
$t^{*}\approx 1.225$. During the second stage, between $t^{*}$ and the maximum
of dissipation around $t=7$ dissipation comes in, and different structures
appear in the flow (see section IV). Finally, closely after the maximum of
dissipation turbulence is fully developed exhibiting Kolmogorov scaling (see
section V).

##  III Early stages: inviscid evolution before $t_{*}$

The initial condition consists on large scale structures, oriented randomly in
the box. Visualizations of the initial condition are displayed in Fig.2 (top
row).

Figure 2: Visualizations of the enstrophy field at $t=0$ (top row) and $t=1.2$
(bottom row). The left panels show 3D visualizations of the field, whereas the
center and right columns displays cuts of the $xy$\- and $xz$-planes. Run at
resolution $1024^{3}$.

During the evolution, vorticity blobs shrink, creating thinner and thinner
vortex sheets, following a process similar to the one reported in earlier
papers [10, 11, 15, 16] (see Fig.2 (bottom row)). A natural question arises,
are those vortex sheets shrinking in a completely self-similar manner? If so,
is such a self-similarity of the blowup (second kind) type as it is in the
simplified turbulence models, e.g. the Leith model and the EDQNM model
discussed in the introduction. This type of behaviour would correspond to a
precursor of an eventual finite-time singularity of the inviscid dynamics. In
the following, we will address these questions and study the self-similarity
of the early stage evolution of our simulations.

A self-similar solution for the energy spectrum has the following form, of the
second kind (keeping notation of Nazarenko and Grebenev [17]):

$E(k,t)=g(t)F(\eta),\quad\eta=k/k_{*}(t),$ (5)

where $k_{*}(t)$ is the characteristic wave number of the propagating front of
the spectrum. For the second-kind self-similarity, $g(t)=k_{*}(t)^{x}$ and
$k_{*}(t)=(t_{*}-t)^{b}$, where $x$ and $b$ are some similarity constants
which are to be found from solving a ”nonlinear eigenvalue problem” of
simultaneous matching the two boundary conditions: $F(\eta)\to\eta^{-x}$ for
$\eta\to 0$ and $F(\eta)\to 0$ (sufficiently fast) at $\eta\to\infty$. (Note
that in such as solution, the spectrum is asympotically stationary at $\eta\to
0$). Such kind of self-similarity was found for the Leith model [3, 5, 18] and
the EDQNM model [6], both being simplified models for 3D hydrodynamic
turbulence. This leads to a natural question if the second-kind kind self-
similarity could also be observed in the DNS of the 3D Navier-Stokes
turbulence. As we will see later, this form of self-similar solution is not
observed in our simulations which, however, does not rule out a possibility
that it could appear in simulations at much higher resolution at times very
close to $t_{*}$. Our present results indicate that the solution is also
approximately self-similar with $k_{*}(t)$ function of time increasing
exponentially. In this case, $t_{*}$ has a meaning of the time at which the
self-similar stage ends (when the dissipative and smallest resolved scales are
reached) rather than a ”blowup” time.

In order to study the pre-$t^{*}$ self-similar evolution and beyond, we fit
the spectra with the following trial functions:

$\displaystyle E_{\rm fit\,0}(k,t)$ $\displaystyle=$ $\displaystyle
E_{0}(t)k^{-m_{0}(t)}e^{-2\delta_{0}(t)k}$ (6) $\displaystyle E_{\rm
fit\,1}(k,t)$ $\displaystyle=$ $\displaystyle
E_{1}(t)k^{-m_{1}(t)}e^{-2\delta_{1}(t)k-(\eta(t)k)^{n_{\alpha}}}$ (7)

The interpretation is the following. At early times, the dissipation is
negligible and in absence of finite time singularities, the velocity field is
analytic and (at least) an exponential decay is expected in the high-$k$
range. In Euler flows, singularities appear as poles in the complex plane of
the analytical continuation of the velocity flow [19]. The scale $\delta(t)$
is usually called the analyticity strip and corresponds to the distance from
the real plane to such singularities. If $\delta$ vanishes, then a singularity
has crossed the real axes corresponding to an actual singularity of the flow.
At a later stage, when dissipation start acting, a stronger decay is expected.
Then $\eta(t)$ would correspond to the Kolmogorov scale in fully developed
turbulent states. We expect that at early times $\eta(t)$ will be very small,
followed by a sudden increase at $t^{*}$. The proper value of $n_{\alpha}$
certainly depends on the kind of dissipation. It turns out that in our hyper-
viscous setup the best fitting value is $n_{\alpha}=4$. Note that self-
similarity implies $m_{0}=m_{1}=$ const, $\delta_{0}(t)=\delta_{1}(t)\propto
1/k_{*}(t)$. Also, as $\eta$ represents a dissipative scale, it can not
decrease with time, and therefore self-simlarity imposes $\eta=0$, which is
expected for $t<t^{*}$. We will see that, like in the second-type self-
similarity, our spectrum appears to be stationary at the IR end, which implies
$E_{0}=E_{1}=$ const. We will not fix the time dependence of $k_{*}$, and will
find it from fitting the data. In this section we analyze the data of the run
at resolution $2048^{3}$.

The energy spectra at different moments of time, together with the fits (6)
and (7) are shown in Fig. 3.

Figure 3: Temporal evolution of energy spectra and fits. Fit 6 is plotted in
magenta and 7 in cyan. Run at resolution $2048^{3}$.

We see that, somehow, the fits could be used as an extrapolation of the
numerical results beyond numerical precision. The temporal evolution of the
fitting parameters $E_{0,1},\delta_{0,1},m_{0,1}$ and $\eta$ are displayed in
Fig. 4. From the inset of Fig. 4.(d), which shows $\frac{{\rm d}\eta}{{\rm
d}t}$, we infer that the sudden jump of the derivative indicates that

$t_{*}\approx 1.24.$

Figure 4: Evolution of the fitting parameters in Eqs.6 and 7. Run at
resolution $2048^{3}$.

Again, the suggested fits indicate approximate self-similarity only when the
resulting parameters $E_{0},E_{1}$ and $m_{0},m_{1}$ are relatively constant
in time. From the evolution of the fitting parameters shown in Fig. 4, we see
that the period of time for which this condition is approximately satisfied is
$0.6\lesssim t\lesssim t_{*}$. From the lower-left frame of Fig. 4, we see
that at the pre-$t_{*}$ stage the exponent of the power law behind the
propagating front of the spectrum is $m\sim 3.25$, which is consistent with
the previous numerical results obtained at low resolutions in [10].

Let us seek self-similar solutions by determining the temporal evolution of
$k_{*}(t)$ numerically. We can define $k_{*}(t)$ by using the fits Eqs. 6 and
7 or by measuring it directly from the spectrum, namely

$\displaystyle k_{*}^{\rm raw}(t)$ $\displaystyle=$
$\displaystyle\max{[k\,|E(k,t)>\epsilon}]$ (8) $\displaystyle\hbox{or}\quad
k_{*}^{\delta}(t)$ $\displaystyle=$ $\displaystyle\lambda/\delta(t).$ (9)

We expect that, up to a constant, all definitions present the same trend. In
numerics we have set $\epsilon=10^{-30}$. The pre-factor $\lambda$ is set in
order to match the values of $k_{*}^{\rm raw}$ and $k_{*}^{\delta}(t)$ during
the self-similar evolution.

In Fig. 5 we present the temporal evolution of $k_{*}(t)$, normalized by the
maximum wavenumber of the simulation $k_{\rm max}$.

Figure 5: Temporal evolution of $k^{*}t$. The maximum wavenumber of the
simulation is $k_{\rm max}=2048/3$.

We observe a clear exponential growth of $k^{*}$ in the range $k^{*}(t)<k_{\rm
max}$, which, at the current resolution, excludes the second-kind self-
similarity with a finite-time singularity scenario. Again, we do not rule out
a theoretical possibility that at much higher resolution one might see a
transition to an accelerated (second-kind) self-similar regime very close to
$t_{*}$.

On the other hand, the exponential growth of $k^{*}(t)$, and the respective
exponential shrinking of characteristic physical scale, has a clear physical
interpretation as a process of shrinking of the pancake-like vorticity
structures embedded to a strong strain field produced by the large scales. The
fact that the large-scale strain is indeed very strong is explained by a very
steep initial energy spectrum. For a very strong strain, the action of the
pancakes onto themselves can be neglected compared to the action of the large-
scale strain on them. Such a nonlocal (in scale) dynamics can be described the
Rapid Distortion Theory (RDT) [20, 21], and the exponential shrinking of the
scale in this approach is simply due to the motion of the vorticity along the
exponentially converging steamlines of the large-scale straining motion.

Let us now embark on a more formal search of the self-similar solution that
provides the best fit to our numerical data before $t_{*}$. We look for a
self-similar solution given by equation (5) with $g(t)=k_{*}(t)^{x}$. We use
$k_{*}(t)$ found from fitting the data for finding the scaling parameter $x$
via minimizing the following error,

$\displaystyle error(x)$ $\displaystyle=$
$\displaystyle\int_{T_{0}}^{T_{1}}\int_{T_{0}}^{T_{1}}\left[\int_{\eta_{\rm
i}}^{\eta_{\rm
f}}\left(\log{\tilde{F}(\eta,t)}-\log{\tilde{F}(\eta,s)}\right)^{2}{\rm
d}\eta\right]{\rm d}s{\rm d}t$ (10) $\displaystyle\eta=\frac{k}{k_{*}(t)}$ ,
$\displaystyle\tilde{F}(\eta,t)=k_{*}(t)^{-x}E(k,t).$ (11)

Note that for a perfect self-similar evolution, $\tilde{F}(\eta,t)=F(\eta)$,
and then for the correct value of the exponent $x$ we thus have $error(x)=0$.
The resulting error of minimizing (10) is displayed on the inset of figure 6,
and gives the optimal value $x=-3.49$. The error is computed in the interval
$\eta\in(0.0126,1.1220)$ (represented by dashed green vertical lines). The
value of the exponent $-3.49$ is consistent withithin $7\%$ with the result
for the power-law exponent behind the propagating front $-m\approx 3.25$, as
seen in the inset of Fig. 4 (b).

Figure 6:  Self-similar form Eq. (5) for different times before $t*$. Error
function (arbitrary normalization) (10) as a function of the scaling parameter
$\gamma$. Bounds of the integrals are set to $\eta_{\rm i}=0.0126$, $\eta_{\rm
f}=1.12$ (green dashed vertical lines) and $T_{0}=.7475$ , $T_{1}=1.2275$. Run
at resolution $2048^{3}$.

We use now the obtained value of $x$ to show the self-similar behavior of the
energy spectra. The function $F(\eta)$ at different times is displayed in
Figure 6. The collapse onto a single curve is excellent.

##  IV Intermediate times, rolling pancakes and vortex ribs

After $t^{*}$, dissipation start acting, taking the system beyond the purely
self-similar shrinking of pancakes. Visualizations of the enstrophy field at
$t=3.6$ and $t=4.8$ are displayed in Fig. 7.

Figure 7: Visualizations of the enstrophy field at $t=3.6$ (top row) and
$t=4.8$ (middle row). The left panels of the top and the middle rows show 3D
visualizations of the field, whereas the center and right columns display cuts
of the $xy$\- and $xz$-planes. The bottom row show several consecutive frames
zoomed into a region with a typical vortex sheet roll-up and formation of the
vortex ribs. Run at resolution $1024^{3}$.

Early post-$t_{*}$ evolution is characterised by persistence of the vortex
pancakes which continue to shrink and whose vorticity continues to grow for a
time of about three to four $t_{*}$. In the top row of Fig. 7 we see that at
$t=3.6$ some of the vortex pancakes are rolling up at the edges, quickly
exciting small scales which are purged by hyperviscosity. Such a purging
allows the pancakes to keep shrinking in other (more central) areas until the
hyper-viscosity counteracts with the shrinking by (hyper-) diffusing it out of
the thin layers thereby leading to a slowdown and eventual arrest of the
shrinking process. After that some pancakes become unstable which leads to
their breakup into sets of quasi-periodic quasi-colinear thin vortex tubes,
“vortex ribs”, see the middle row of Fig. 7 for the enstrophy field at
$t=4.8$. Transverse cross-section of such ribs shows a vorticity pattern
typical for a vortex street arising after the Kelvin-Helmholtz type
instability of a vortex layer. Similar behaviour was previously observed in
$256^{3}$ simulations of Navier-Stokes equations, and the rollup of the vortex
pancakes into the array of vortex filaments was explained as a result of
instability which could be understood as an instability of the Burgers vortex
layer [11]. The dissipation mechanism in our simulations is hyperviscosity
rather than a regular viscosity, but clearly similar quasi-stationary vortex
layers in which the vortex stretching by external large-scale strain is
balanced by a layer (hyper-)diffusion must appear and become unstable in our
system too. After forming, vortex ribs start approaching to each other and
form twisted vortex bundles reminiscent of ropes, as also seen in the middle
row of Fig. 7 (left panel). The bottom panel of Fig. 7 shows a zoom of the box
where the shrinking of pancakes, the rolling and the formation of vortex ribs
is clearly visible. Further in time, the number of vortex ribs and ropes
proliferate and become abundant throughout the fluid, and their mutual
interactions lead to an increasingly chaotic vortex tangle at $t\gtrsim
4.5t_{*}$, see e.g. the transition in the enstrophy field character from
$t=4.8$, the bottom row of Fig. 7, to $t=7.05$ in Fig. 9.

In terms of the energy spectrum, the evolution for $t_{*}<t<4.5t_{*}$ is not
too eventful: it is characterised by a gradual increase of the power-law
exponent from $\sim-3.4$ to $\sim-5/3$.

Figure 8: Energy spectra at intermediate times. Run at resolution $1024^{3}$

##  V Long times

We now address the later stage, $t\gtrsim 7$, which corresponds to the
dynamics after the maximum of dissipation is reached (see Fig. 1). At this
stage, vortex ribs are not longer present and filaments and vortex structures
display a more chaotic behviour as apparent in Fig. 9, tipical of turbulent
flows.

Figure 9: Visualizations of the enstrophy field at $t=7.05$. The left panel
shows a 3D visualization of the field, whereas the center and right columns
display cuts of the $xy$\- and $xz$-planes. Run at resolution $1024^{3}$.

Indeed, turbulence is fully developped and the energy spectrum displaus a
Kolmogorov $-5/3$ scaling as shwon in Fig. 10.

Figure 10: Temporal evolution of compensated energy spectrum at larger times.
Run at resolution $1024^{3}$.

We notice the emergence of a clear Kolmogorov scaling for $5<k<50$. This is
apparent in the left panel of Fig. 10 where $k^{5/3}$-compensated spectra are
displayed. For $k>50$ we see a bump which is typical in simulations with
hyperviscosity: it corresponds to a tendency to the turbulence thermalization
due to presence of the cutoff wave number $k_{max}$. The right panel of Fig.
10 displays the spectra compensated by $\varepsilon(t)^{2/3}k^{5/3}$. Here we
can see that the whole spectrum, including both the Kolmogorov and the thermal
bump parts, follow approximately the same law of decay at the late stages. In
the other words, the late stage of the spectrum evolution is also self-
similar.

The decay of energy and dissipation observed in Fig. 1 is typical of turbulent
flows. Several power laws have been proposed and reported depending on the
nature of the large-scale configuration of the turbulent flow [22, 23, 17,
24]. Turbulent decay laws are constructed under the assumption of existence of
an inertial range. Assuming that the energy dissipation is equal to the energy
input at large scales, a simple balance equation can be easily derived
$\epsilon=-\frac{dE}{dt}\sim E^{3/2}/\ell$, where $\ell$ is the integral scale
of the flow. In general, the difficulty lies in determining the temporal
evolution of the integral scale $\ell$. Such phenomenological arguments also
apply to our hyperviscous simulation as they only require a well-developed
turbulent flow obeying Kolmogorov phenomenology. In our particular case, the
situation is simpler as the integral scale is saturated at the size of the
box, and therefore $\ell$ remains constant in this stage. This assumption
leads to a decay low $E(t)\sim t^{-2}$. In our simulations, the decaying range
is too narrow to check this exponent directly, but our data is compatible with
the law $\epsilon(t)\sim E(t)^{3/2}$, as shown in Fig. 11.

Figure 11: Temporal evolution of compensated energy spectrum at larger times.
Run at resolution $1024^{3}$.

##  VI conclusions

In this paper, we have performed DNS of the hyper-viscous Navier-Stokes system
with no forcing and with only few randomly selected low-wavenumber Fourier
modes excited initially. The goal was to study transition to the fully
turbulent state, with specific attention to self-similar behaviour and to
different coherent quasi-singular vortex structures appearing at different
evolution stages. We we motivated by the previous results for the differential
(Leith) and integro-differential closures (EDQNM) for the hydrodynamic
turbulence as well as for the Burgers equation. For these systems, the second-
kind self-similarity was observed, where the propagating spectral front
reached the infinite wavenumber in a finite time $t^{*}$ after which a
reflected wave propagates toward the smaller wavenumbers bringing the
Kolmogorov scaling in the wake (Burgers scaling in the Burgers equation
model). In the Burgers case this evolution is related to the successive
appearance of singularities (quasi-singularities in presence of small
viscosity) of two types: the pre-shock and the shock respectively. Thus we
wanted to know if there is a self-similar evolution in the Navier-Stokes
system, if this self-similarity is of the second (blowup) type (with a
consequence of a blowup occuring in the Euler equations), and if we could
identify quasi-singular structures of different types during the evolution en
route to the Kolmogorov state.

Our analysis of the DNS data shows that at the inviscid stage the evolution is
indeed self-similar but it not of the second kind: instead of the blowup
propagation of the spectral front to infinite $k$ in a finite time we observe
an exponential in time motion. Second, we do observe quasi-singularities: (1)
shrinking vortex pancakes are forming at the inviscid stage via an RDT-type
mechanism of vortex stretching/shrinking by a large-scale strain and (2)
breakup of pancakes into quasi-periodic sets of quasi-singular vortex
filaments, vortex ribs, by a Burgers vortex sheath instability. The last stage
of transition to the Kolmogorov state is characterised by ribs twisting into
ropes followed by entanglement of the ribs and ropes into a visually
unstructured vortex tangle.

There is a natural question: why do the turbulence closures, Leith model and
EDQNM, both exhibit the second-kind self-similarity with a blowup spectral
front propagation whereas our DNS of the hyper-viscous Navier-Stokes does not?
The clue can be found in the RDT-type mechanism of the vortex pancake
stretching/squashing. Indeed, this mechanism is based on a nonlocal (in the
scale space) interaction: the strong large-scale strain affects the small
vorticity scales associated with the thin pancakes much stronger than the
small scales affect each other. But the Leith and EDQNM models do not capture
such effects: locality of interactions is ”built into” these two models. Thus
we arrive at a conjecture that it is the nonlocal interactions of scales that
upset the blowup behaviour. However, at this point it could not be ruled out
that much higher resolution simulations, or simulations with some special
initial conditions, might reveil a blowup which could supersede the
exponential scale shrinking.

##  VII Acknowledgements

This work was funded by the Simons Foundation Collaboration grant Wave
Turbulence (Award ID 651471). Computations were carried out at the Mésocentre
SIGAMM hosted at the Observatoire de la Coôte d’Azur.

## References

  * Frisch [1995] U. Frisch,  _Turbulence: The Legacy of A. N. Kolmogorov_ (Cambridge University Press, 1995).
  * Leith [1967] C. Leith, The Physics of Fluids 10, 1409 (1967).
  * Connaughton and Nazarenko [2004] C. Connaughton and S. Nazarenko, Physical review letters 92, 044501 (2004).
  * Orszag [1970] S. A. Orszag, Journal of Fluid Mechanics 41, 363–386 (1970).
  * Grebenev  _et al._ [2014] V. Grebenev, S. Nazarenko, S. Medvedev, Y. Chirkunov, and I. Schwab, Journal of Physics A: Mathematical and Theoretical 47, 025401 (2014).
  * Bos  _et al._ [2012] W. J. Bos, C. Connaughton, and F. Godeferd, Physica D: Nonlinear Phenomena 241, 232 (2012), special Issue on Small Scale Turbulence.
  * GALTIER  _et al._ [2000] S. GALTIER, S. V. NAZARENKO, A. C. NEWELL, and A. POUQUET, Journal of Plasma Physics 63, 447–488 (2000).
  * Sulem  _et al._ [1983a] C. Sulem, P.-L. Sulem, and H. Frisch, Journal of Computational Physics 50, 138 (1983a).
  * Thalabard  _et al._ [2015a] S. Thalabard, S. Nazarenko, S. Galtier, and S. Medvedev,  48, 285501 (2015a).
  * Brachet  _et al._ [1992] M. E. Brachet, M. Meneguzzi, A. Vincent, H. Politano, and P. L. Sulem, Physics of Fluids A: Fluid Dynamics 4, 2845 (1992), https://doi.org/10.1063/1.858513 .
  * Passot  _et al._ [1995] T. Passot, H. Politano, P. Sulem, J. Angilella, and M. Meneguzzi, Journal of Fluid Mechanics 282, 313–338 (1995).
  * Saffman [1995] P. Saffman,  _Vortex Dynamics_, Cambridge Monographs on Mechanics (Cambridge University Press, 1995).
  * Lin and Corcos [1984] S. J. Lin and G. M. Corcos, Journal of Fluid Mechanics 141, 139–178 (1984).
  * Beronov and Kida [1996] K. Beronov and S. Kida, Physics of Fluids 8, 1024 (1996).
  * Agafontsev  _et al._ [2015] D. S. Agafontsev, E. A. Kuznetsov, and A. A. Mailybaev, Physics of Fluids 27, 085102 (2015), https://aip.scitation.org/doi/pdf/10.1063/1.4927680 .
  * Agafontsev  _et al._ [2017] D. S. Agafontsev, E. A. Kuznetsov, and A. A. Mailybaev, Journal of Fluid Mechanics 813, R1 (2017).
  * Nazarenko and Grebenev [2017] S. V. Nazarenko and V. N. Grebenev, Journal of Physics A: Mathematical and Theoretical 50, 035501 (2017).
  * Thalabard  _et al._ [2015b] S. Thalabard, S. Nazarenko, S. Galtier, and S. Medvedev, Journal of Physics A: Mathematical and Theoretical 48, 285501 (2015b).
  * Sulem  _et al._ [1983b] C. Sulem, P.-L. Sulem, and H. Frisch, Journal of Computational Physics 50, 138 (1983b).
  * BATCHELOR and PROUDMAN [1954] G. K. BATCHELOR and I. PROUDMAN, The Quarterly Journal of Mechanics and Applied Mathematics 7, 83 (1954), https://academic.oup.com/qjmam/article-pdf/7/1/83/5300208/7-1-83.pdf .
  * NAZARENKO  _et al._ [1999] S. NAZARENKO, N. K.-R. KEVLAHAN, and B. DUBRULLE, Journal of Fluid Mechanics 390, 325–348 (1999).
  * LIN [1948] C. C. LIN, Proceedings of the National Academy of Sciences of the United States of America 34, 540 (1948).
  * Saffman [1967] P. G. Saffman, The Physics of Fluids 10, 1349 (1967), https://pubs.aip.org/aip/pfl/article-pdf/10/6/1349/12566054/1349_1_online.pdf .
  * Panickacheril John  _et al._ [2022] J. Panickacheril John, D. A. Donzis, and K. R. Sreenivasan, Philosophical Transactions of the Royal Society A 380, 20210089 (2022).