# Initial evolution of three-dimensional turbulence en route to the Kolmogorov state: Emergence and transformations of coherent structures, self-similarity, and instabilities

Giorgio Krstulovic, Sergey Nazarenko

► To cite this version:

Giorgio Krstulovic, Sergey Nazarenko. Initial evolution of three-dimensional turbulence en route to the Kolmogorov state: Emergence and transformations of coherent structures, self-similarity, and instabilities. Physical Review Fluids, 2024, 9 (3), pp.034610. 10.1103/PhysRevFluids.9.034610. hal-04800934

HAL Id: hal-04800934

https://hal.science/hal-04800934v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 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

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³ 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 ≈ -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 an instability should lead to formation of another type of quasi-singular vortex structures–thin intense vortex filaments. Later in [11], pseudo-spectral 256³ 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 × 1024 × 2048 in [15] and 972 × 2048 × 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 3D decaying solutions $\mathbf{u}(\mathbf{x}, t)$ of the incompressible hyper-viscous Navier-Stokes equations:

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = \nabla p - \nu(-\Delta)^{\alpha} \mathbf{u} \qquad (1) $$

$$ \nabla \cdot \mathbf{u} = 0. \qquad (2) $$

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

$$ E(t) = \frac{1}{2V} \int |\mathbf{u}|^2 d\mathbf{x}, \quad \text{and} \quad \epsilon(t) = \frac{\nu}{V} \int \mathbf{u}(-\Delta)^{\alpha} \mathbf{u} d\mathbf{x} = \nu \sum_{\mathbf{k}} |\mathbf{k}|^{2\alpha} |\hat{\mathbf{u}}(\mathbf{k})|^2, \qquad (3) $$

where $\hat{\mathbf{u}}(\mathbf{k})$ is the Fourier transform of $\mathbf{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

$$ \hat{u}_i(\mathbf{k}, 0) = \sum_{|\mathbf{k}|<k_I} \frac{\sqrt{C}}{k} e^{i\phi_k}, \quad i = 1, 2, 3. \qquad (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_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. We observe three stages which will be

<figure><img src="image_2.png" /><figcaption>FIG. 1. Temporal evolution of the energy and energy dissipation. Run at resolution $1024^3$.</figcaption></figure>

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 see later, this phase ends when the propagating spectrum reaches 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). During the evolution, vorticity blobs shrink, creating thinner and thinner

<figure><img src="image_3.png" /><figcaption>FIG. 2. Visualizations of the entrophy 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$.</figcaption></figure>

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), \qquad (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 asymptotically 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:

$$E_{\text{fit 0}}(k, t) = E_0(t)k^{-m_0(t)}e^{-2\delta_0(t)k} \qquad (6)$$

$$E_{\text{fit 1}}(k, t) = E_1(t)k^{-m_1(t)}e^{-2\delta_1(t)k-(\eta(t)k)^{n_\alpha}} \qquad (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 = \text{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-similarity 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 = \text{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³.

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

<figure><img src="image_4.png" /><figcaption>FIG. 3. Temporal evolution of energy spectra and fits. Fit 6 is plotted in magenta and 7 in cyan. Run at resolution 2048³.</figcaption></figure>

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{d\eta}{dt}$, we infer that the sudden jump of the derivative indicates that

$$t_* \approx 1.24.$$

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

$$k_*^{\text{raw}}(t) = \max[k|E(k, t) > \epsilon] \qquad (8)$$

$$\text{or } k_*^\delta(t) = \lambda/\delta(t). \qquad (9)$$
<figure><img src="image_5.png" /><figcaption>FIG. 4. Evolution of the fitting parameters in Eqs. 6 and 7. Run at resolution 2048³.</figcaption></figure>

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_*^{\text{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_{\text{max}}$. We observe a clear exponential growth of $k^*$ in the range $k^*(t) < k_{\text{max}}$, which, at the current resolution,

<figure><img src="image_6.png" /><figcaption>FIG. 5. Temporal evolution of $k^*t$. The maximum wavenumber of the simulation is $k_{\text{max}} = 2048/3$.</figcaption></figure>

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 streamlines 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,

$$ error(x) = \int_{T_0}^{T_1} \int_{T_0}^{T_1} \left[ \int_{\eta_i}^{\eta_f} \left( \log \tilde{F}(\eta, t) - \log \tilde{F}(\eta, s) \right)^2 d\eta \right] ds dt \quad (10) $$

$$ \eta = \frac{k}{k_*(t)}, \quad \tilde{F}(\eta, t) = k_*(t)^{-x} E(k, t). \quad (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 within 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). We use now the obtained value of $x$ to show the

<figure><img src="image_7.png" /><figcaption>FIG. 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_i = 0.0126$, $\eta_f = 1.12$ (green dashed vertical lines) and $T_0 = .7475$, $T_1 = 1.2275$. Run at resolution $2048^3$.</figcaption></figure>

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. 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
<figure><img src="image_8.png" /><figcaption>FIG. 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$.</figcaption></figure>

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><img src="image_9.png" /><figcaption>FIG. 8. Energy spectra at intermediate times. Run at resolution 1024³</figcaption></figure>

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 behaviour as apparent in Fig. 9, typical of turbulent flows. Indeed, turbulence is fully developed and

<figure><img src="image_10.png" /><figcaption>FIG. 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$.</figcaption></figure>

the energy spectrum displays a Kolmogorov $-5/3$ scaling as shown in Fig. 10. 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 $\epsilon(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 [17, 22–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
<figure><img src="image_11.png" /><figcaption>FIG. 10. Temporal evolution of compensated energy spectrum at larger times. Run at resolution 1024³.</figcaption></figure>

$\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><img src="image_12.png" /><figcaption>FIG. 11. Temporal evolution of compensated energy spectrum at larger times. Run at resolution 1024³.</figcaption></figure>

## 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 were 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 reveal 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.

[1] U. Frisch, *Turbulence: The Legacy of A. N. Kolmogorov* (Cambridge University Press, 1995).

[2] C. Leith, *The Physics of Fluids* **10**, 1409 (1967).

[3] C. Connaughton and S. Nazarenko, *Physical review letters* **92**, 044501 (2004).

[4] S. A. Orszag, *Journal of Fluid Mechanics* **41**, 363-386 (1970).

[5] V. Grebenev, S. Nazarenko, S. Medvedev, Y. Chirkunov, and I. Schwab, *Journal of Physics A: Mathematical and Theoretical* **47**, 025401 (2014).

[6] W. J. Bos, C. Connaughton, and F. Godeferd, *Physica D: Nonlinear Phenomena* **241**, 232 (2012), special Issue on Small Scale Turbulence.

[7] S. GALTIER, S. V. NAZARENKO, A. C. NEWELL, and A. POUQUET, *Journal of Plasma Physics* **63**, 447-488 (2000).

[8] C. Sulem, P.-L. Sulem, and H. Frisch, *Journal of Computational Physics* **50**, 138 (1983).

[9] S. Thalabard, S. Nazarenko, S. Galtier, and S. Medvedev, *48*, 285501 (2015).

[10] 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.

[11] T. Passot, H. Politano, P. Sulem, J. Angilella, and M. Meneguzzi, *Journal of Fluid Mechanics* **282**, 313-338 (1995).

[12] P. Saffman, *Vortex Dynamics*, Cambridge Monographs on Mechanics (Cambridge University Press, 1995).

[13] S. J. Lin and G. M. Corcos, *Journal of Fluid Mechanics* **141**, 139-178 (1984).

[14] K. Beronov and S. Kida, *Physics of Fluids* **8**, 1024 (1996).

[15] 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.

[16] D. S. Agafontsev, E. A. Kuznetsov, and A. A. Mailybaev, *Journal of Fluid Mechanics* **813**, R1 (2017).

[17] S. V. Nazarenko and V. N. Grebenev, *Journal of Physics A: Mathematical and Theoretical* **50**, 035501 (2017).

[18] S. Thalabard, S. Nazarenko, S. Galtier, and S. Medvedev, *Journal of Physics A: Mathematical and Theoretical* **48**, 285501 (2015).

[19] C. Sulem, P.-L. Sulem, and H. Frisch, *Journal of Computational Physics* **50**, 138 (1983).

[20] G. K. BATCHELOR and I. PROUDMAN, *The Quarterly Journal of Mechanics and Applied Mathematics* **7**, 83 (1954), https://academic.oup.com/qj mam/article-pdf/7/1/83/5300208/7-1-83.pdf.

[21] S. NAZARENKO, N. K.-R. KEVLAHAN, and B. DUBRULLE, *Journal of Fluid Mechanics* **390**, 325-348 (1999).

[22] C. C. LIN, *Proceedings of the National Academy of Sciences of the United States of America* **34**, 540 (1948).

[23] 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.

[24] J. Panickacheril John, D. A. Donzis, and K. R. Sreenivasan, *Philosophical Transactions of the Royal Society A* **380**, 20210089 (2022).