Quantum Solution of Classical Turbulence. Decaying Energy Spectrum

1. Prologue

Turbulence appears as an overwhelmingly complex phenomenon. As depicted in Figure ?? from a recent simulation [48], vortex lines of various shapes and sizes are entangled much like spaghetti. This visual complexity raises the question: How can such complexity be described analytically? Yet, it also sparks hope for a simplified statistical description.

With its myriad interacting particles, molecular dynamics similarly presents an intricate challenge. However, despite its complexity, a straightforward statistical description emerges that grows increasingly precise with the escalating complexity of the dynamical system. Maxwell, Boltzmann, and Gibbs demonstrated that Newton’s mechanics uniformly cover the energy surface over time, laying the groundwork for statistical mechanics—a robust theory, albeit sometimes computationally challenging, as in critical phenomena.

Why, then, should Navier-Stokes turbulence be any different?

Regrettably, to date, no known analog of the Gibbs distribution exists for turbulent flows. Therefore, a foundational element of turbulence theory must be to devise a substitute for the Gibbs distribution.

Hopf initiated this exploration in 1952 (see the recent review in [45]), formulating a functional equation that the probability distribution of the turbulent velocity field must satisfy. Through iterative application of this equation to the nonlinear term in the Navier-Stokes equation, one can generate an expansion in inverse powers of viscosity. The core challenge of turbulence theory is solving the Hopf equation in the opposite limit of low viscosity.

This beautiful equation is mathematically as intricate as the vortex spaghetti depicted earlier. Such complexity places turbulence high within the hierarchy of physics theories, nestled between critical phenomena and the quark confinement problem.

Our theory, initiated in 1993 (see Figure ?? for the historical outline), proposes a simpler variant of the Hopf equation—the loop equation—which suffices to define the statistics.

Figure 1. The path to the microscopic theory.

Figure 1. The path to the microscopic theory.

The loop equation corresponds to the Schrödinger equation in loop space. This profound analogy is not a poetic metaphor but a precise mathematical equivalence with significant implications, such as quantum interference effects affecting the scaling laws of classical turbulence.

Using the loop equation, we have identified a new instance of duality between the strong coupling phase of one theory (a fluctuating velocity field in three dimensions) and the weak coupling phase of another (a one-dimensional quantum ring of Fermi particles).

This weak coupling limit can be analytically solved, providing explicit formulas for observables in decaying turbulence, such as the energy dissipation decay index $n\left(t\right)=t\frac{\mathcal{E}\left(t\right)}{E\left(t\right)}$.

Experimental data from decaying grid turbulence (see Figure 2) corroborate our prediction that $n(\infty )=\frac{5}{4}$, within a 2% experimental error margin. We anticipate that more precise future measurements will validate this prediction more accurately.

Recent DNS and experiments [23,46] compute the energy spectrum and the velocity correlation function, matching our theory and challenging the Kolmogorov scaling. The accuracy is lower here due to large experimental errors. We found the way to reduce experimental errors by computing effective index for the second moment of velocity difference using numerical Fourier transform. This method dramatically improved the quality of the fit for the effective index in our theory, compared to traditional numerical differentiation of the energy spectrum. New large-scale experiments (both real and numerical) are welcome to verify our theory.

2. Definitions and notations

We use the units where the constant fluid density $\rho =1$. The 3D vectors and the dot products are denoted like this: $\overrightarrow{v}(\overrightarrow{r},t),\overrightarrow{k}·\overrightarrow{r}$. We also use the Einstein tensor notation with summation over repeated Greek indexes $v_{\alpha },r_{\beta },{\omega }_{\alpha \beta }=e_{\alpha \beta \gamma }{\partial }_{\beta }{v}_{\gamma }$, and so on. In these cases, these indexes run from 1 to d, where d is the dimension of space. In the remaining cases, we work only in 3D space. We always consider the Euclidean metric and Cartesian coordinates, so we make no distinction between upper and lower tensor and vector indexes. The parameterization of the loops is chosen to run from 0 to 1 with periodicity implied beyond these limits. The space loops $\overrightarrow{C}\left(\xi \right)$ are assumed to be continuous and periodic but not differentiable. In addition, the loop can have some extra periods, in which case this loop represents several separate loops. In particular, the same geometric loop $\overrightarrow{C}$ can be covered several times with ${\overrightarrow{C}}_n\left(\xi \right)=C\left(n\xi \right)$. The momentum loops $\overrightarrow{P}(\xi ,t)$ depend on time, whereas the spacial loops are static. The momentum loops can have discontinuities $\mathsf{\Delta }\overrightarrow{P}(\theta ,t)=\overrightarrow{P}(\theta +0,t)-\overrightarrow{P}(\theta -0,t)$. The momentum loops are independent of the space loops, being the vector functions of $\xi ,t$. The circulation $\mathsf{\Gamma }={\oint }_{C}d\overrightarrow{r}·\overrightarrow{v}(\overrightarrow{r},t)={\int }_0^{1}d\xi {\overrightarrow{C}}^{\prime }\left(\xi \right)·\overrightarrow{v}(\overrightarrow{C}\left(\xi \right),t)$. The whole theory, starting with circulation, is parametric invariant $\xi \Rightarrow f\left(\xi \right),f^{\prime }\left(\xi \right)>0$. The operators are denoted like $\widehat{a},{\widehat{a}}^{†},\widehat{X}$. We use three types of symbols for differentials. d is an ordinary differential of variable like $d\xi$ or vector variable like $d\overrightarrow{C}\left(\theta \right)=d\theta {\overrightarrow{C}}^{\prime }\left(\theta \right)$. D symbolizes a measure for a path integral $D\alpha$; it can be strictly defined by the limit of the multidimensional integral over the discrete values or as a limit of a multidimensional integral over all Fourier harmonics (this is cleaner). We do not need to specify the Gaussian measure, as it is well-known in physics and math. Finally, $\delta$ is a functional variation, and $\frac{\delta }{\delta }$ is a functional derivative like in $\frac{\delta }{\delta \varphi \left(\xi \right)}{\mathsf{\Lambda }}_N\left[\varphi \right]=\frac{\delta {\mathsf{\Lambda }}_N\left[\varphi \right]}{\delta \varphi \left(\xi \right)}$ with the $L_2$ norm in functional space. We also use notations ${\partial }_t=\frac{\partial }{\partial t},{\partial }_{\beta }=\frac{\partial }{\partial r_{\beta }}$ for the time derivative and for components of the spatial derivatives.

Figure 2. The experimental data [9,10] for decaying turbulence behind the oscillating grid. Reproduced with permission from "Comte-Bellot G, Corrsin S. The use of a contraction to improve the isotropy of grid-generated turbulence". Journal of Fluid Mechanics. 1966; 25(4):657-682, © 1966 Cambridge University Press. The two lines correspond to log-log plots of $1/v_{\perp }^2,1/v_{\Vert }^2$ for the flow behind the grid. The total inverse kinetic energy i$K=2/(v_{\perp }^2+v_{\Vert }^2)$, would have the the mean slope $(1.24+1.27)/2\approx 1.255\pm 0.02$ in agreement with our prediction $n(\infty )=\frac{5}{4}$. This plot was shown by K.R.Sreenivasan at the ICTS in Dec’23 in Bangalore, India [54] as the most reliable experimental data. Other measurements and simulations, at different conditions reviewed in [46], provide mismatching data influenced by initial and boundary conditions. The only large-scale DNS [40,46] covering three time decades also matches our predictions.

Figure 2. The experimental data [9,10] for decaying turbulence behind the oscillating grid. Reproduced with permission from "Comte-Bellot G, Corrsin S. The use of a contraction to improve the isotropy of grid-generated turbulence". Journal of Fluid Mechanics. 1966; 25(4):657-682, © 1966 Cambridge University Press. The two lines correspond to log-log plots of $1/v_{\perp }^2,1/v_{\Vert }^2$ for the flow behind the grid. The total inverse kinetic energy i$K=2/(v_{\perp }^2+v_{\Vert }^2)$, would have the the mean slope $(1.24+1.27)/2\approx 1.255\pm 0.02$ in agreement with our prediction $n(\infty )=\frac{5}{4}$. This plot was shown by K.R.Sreenivasan at the ICTS in Dec’23 in Bangalore, India [54] as the most reliable experimental data. Other measurements and simulations, at different conditions reviewed in [46], provide mismatching data influenced by initial and boundary conditions. The only large-scale DNS [40,46] covering three time decades also matches our predictions.

3. Summary

Here are the main results reported in this paper.

• We review the theory of Navier-Stokes loop equation, its relation to the Hopf functional equation, and the representation of the loop functional in terms of momentum loop.
• We present the solution of the loop equation in the inviscid limit of the three-dimensional Navier-Stokes theory in terms of the Euler ensemble. This ensemble consists of a one-dimensional ring of Ising spins in an external field related to random fractions of $\pi$.
• The continuum limit of this solution, $N\to \infty$, corresponds to the inviscid limit of the decaying turbulence in the Navier-Stokes equation. Effective turbulent viscosity is $\tilde{\nu }=\nu N^2\to \mathit{const}$.
• We derived an analytic formula for energy spectrum and dissipation in finite system (58a), (58c), (A188) and investigated it in Appendix K.
• The energy spectrum decays asymptotically as ${\left(\tilde{\nu }\right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{t}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\kappa }^{-\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ where $\kappa =k\sqrt{\tilde{\nu }t}$.
• The turbulent kinetic energy decays as $E\left(t\right)\propto t^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$
• Both effective indexes $n\left(t\right)=-\frac{dlogE}{dlogt},\mu \left(\kappa \right)=\frac{\partial logE(k,t)}{\partial logk}$ are nontrivial functions of the logarithm of scale and time, approaching $n(\infty )=\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.,\mu (\infty )=-\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, (see Figure 8, Figure 6).
• The 1966 experimental values [9,10], Figure 2 of $n(\infty )\approx 1.25\pm 0.02$ agree with our prediction.
• The results of DNS for the energy spectrum in decaying turbulence (Figure 4 of the review paper [46]) shows the energy spectrum decaying faster than $k^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$. The shape of a log-parabolic curve in DNS matches our prediction Figure 5.
• In the Section 8, we verify our prediction for the second velocity moment $〈\mathsf{\Delta }{v}^2〉\left(r\right)$ by numerical Fourier transform of the raw spectrum data from [40,46]. The K41 $\frac{2}{3}$ scaling law is ruled out by this DNS, but there is a match within the DNS errors of the scale-dependent effective index with our theory in a wide range of distances ( Figure 13 ).
• We computed the spectrum of singularity indexes $p_n$ in the product $\overrightarrow{v}(0,t)·\overrightarrow{v}(r,t)\sim \sum C_n\left(t\right)r^{p_n}$, similar but different from the OPE in CFT. Some of these singularity indexes come in complex conjugate pairs related to the zeros of the Riemann $\zeta$ function (84),(88),(93)).

Figure 3. Green curve is the inverse energy $1/E\left(t\right)$ as a function of time. It slowly approaches from above its asymptotic law $1/E\left(t\right)\propto t^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$ shown in red.

Figure 3. Green curve is the inverse energy $1/E\left(t\right)$ as a function of time. It slowly approaches from above its asymptotic law $1/E\left(t\right)\propto t^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$ shown in red.

Figure 4. The log-log plot of deviation from K41 energy spectrum in decaying turbulence, taken from [40,46], Panickacheril John John, Donzis Diego A. and Sreenivasan Katepalli R. 2022, Laws of turbulence decay from direct numerical simulationsPhil. Trans. R. Soc. A.38020210089, licensed under a Creative Commons Attribution (CC BY) license. The observed curve approximately matches our theoretical curve in Figure 5. The K41 spectrum would correspond to a horizontal line, a total mismatch.

Figure 4. The log-log plot of deviation from K41 energy spectrum in decaying turbulence, taken from [40,46], Panickacheril John John, Donzis Diego A. and Sreenivasan Katepalli R. 2022, Laws of turbulence decay from direct numerical simulationsPhil. Trans. R. Soc. A.38020210089, licensed under a Creative Commons Attribution (CC BY) license. The observed curve approximately matches our theoretical curve in Figure 5. The K41 spectrum would correspond to a horizontal line, a total mismatch.

Figure 5. Log-log plot of the universal function $H\left(\kappa \right){\kappa }^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$. It starts growing, reaches the maximum, then turns down, asymptotically decaying as ${\kappa }^{-\raisebox{1ex}{$11$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}$. This qualitatively matches the DNS plot at Figure 4. A larger DNS data at higher Reynolds numbers would be required to verify our theory with more precision.

Figure 5. Log-log plot of the universal function $H\left(\kappa \right){\kappa }^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$. It starts growing, reaches the maximum, then turns down, asymptotically decaying as ${\kappa }^{-\raisebox{1ex}{$11$}\!\left/ \!\raisebox{-1ex}{$6$}\right.}$. This qualitatively matches the DNS plot at Figure 4. A larger DNS data at higher Reynolds numbers would be required to verify our theory with more precision.

Figure 6. The theoretical plot of effective spectral index $\mu \left(T\right)=k{\partial }_{k}logE(k,t)$ as a function of $T=k^2\tilde{\nu }t$. It starts at 0 and goes down to the limit $-\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$. The DNS plot Figure 7 starts higher, at $\mu =2$, then goes down and stops at approximately $-1$ for accessible decay time t. The part from $-0.5$ to $-1$ qualitatively matches our curve, but clearly, more data at higher Reynolds numbers is required to verify this prediction of our theory.

Figure 6. The theoretical plot of effective spectral index $\mu \left(T\right)=k{\partial }_{k}logE(k,t)$ as a function of $T=k^2\tilde{\nu }t$. It starts at 0 and goes down to the limit $-\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$. The DNS plot Figure 7 starts higher, at $\mu =2$, then goes down and stops at approximately $-1$ for accessible decay time t. The part from $-0.5$ to $-1$ qualitatively matches our curve, but clearly, more data at higher Reynolds numbers is required to verify this prediction of our theory.

Figure 7. The DNS plot of spectral index $\mu \left(t\right)=k{\partial }_{k}logE(k,t)$ as a function of $logt$ at fixed values of k, taken from [46], Panickacheril John John, Donzis Diego A. and Sreenivasan Katepalli R. 2022, Laws of turbulence decay from direct numerical simulations. Phil. Trans. R. Soc. A.38020210089, licensed under a Creative Commons Attribution (CC BY) license. The curve goes from $\mu =2$ to $\mu \approx -1$, which partly overlaps with our theoretical curve in Figure 6. Larger dataset with higher Reynolds numbers is needed for a quantitative match with our theory.

Figure 7. The DNS plot of spectral index $\mu \left(t\right)=k{\partial }_{k}logE(k,t)$ as a function of $logt$ at fixed values of k, taken from [46], Panickacheril John John, Donzis Diego A. and Sreenivasan Katepalli R. 2022, Laws of turbulence decay from direct numerical simulations. Phil. Trans. R. Soc. A.38020210089, licensed under a Creative Commons Attribution (CC BY) license. The curve goes from $\mu =2$ to $\mu \approx -1$, which partly overlaps with our theoretical curve in Figure 6. Larger dataset with higher Reynolds numbers is needed for a quantitative match with our theory.

4. Introduction

The original version of this paper was overloaded with formulas, contradicting the preferred 21st-century style. Modern researchers prefer to follow the flow of ideas, with heavy computations hidden under the hood. Eventually, it will become a job for AI to verify computations using Mathematica^® and digest the results for busy readers.

Considering this, we moved all computations into a series of appendices, leaving only the general ideas, concepts, and results in the main body of the paper. This transformation provides distinct pathways for understanding and applying our study’s results based on the reader’s background and interests.

Portions of the upcoming discussion are borrowed from our first paper [30] in this series. We included them to clarify the big picture of the theory of turbulence but significantly modified these sections to reflect our deepened understanding of this theory.

(1)

Introduction for Physicists. The physics introduction discusses the potential correspondence between our theoretical developments and decaying turbulence as observed in real-world or numerical experiments. For physicists, this theory offers a solution to the Hopf functional equation for the statistical distribution of the velocity field in the unforced Navier-Stokes equation. This distribution represents a much-needed analog of the Gibbs distribution for decaying turbulence. There are strong indications that our theory is relevant to one of the two universality classes observed in these experiments.

(2)

Introduction for Mathematicians. This section summarizes the mathematical framework behind the loop equation [28] and its solution [30] in terms of the Euler ensemble. Addressed to mathematicians, this introduction allows those focusing on rigorous mathematical theory to bypass the more physically oriented discussions and delve directly into the Euler ensemble as a novel Number Theory set with conjectured connections to decaying turbulence. Pure mathematicians may want to prove, refine, or disprove the open mathematical problems and unproven conjectures left in this paper.

(3)

Guidance for Applied Mathematicians and Engineers. Applied mathematicians and engineers, primarily interested in practical applications rather than abstract theoretical constructs, are directed to this document’s Section 7. Here, they will find final formulas (58c), (58a), and accompanying Mathematica^®code [32,36] that facilitate the computation of both the energy spectrum and dissipation rates. These formulas are compared with real and numerical experiments in Section 8.

(4)

Notes for the Curious and Skeptical. The fourth category of readers—those curious yet skeptical about applying quantum mechanics to solve complex problems in classical physics—might still harbor doubts after reading the main text of this paper. For them, we have dedicated the last Section 10, which addresses some of their lingering questions and perhaps reassures their skepticism. These readers may want to dive into the Appendixes to learn the details of our theory after this discussion, hopefully eliminating their doubts.

In the following sections, we explore these themes in depth, aiming to provide clarity and actionable insights for all readers, regardless of their expertise or interest.

4.1. Physical Introduction: The Energy Flow and Random Vorticity Structures

Decaying turbulence is an old topic, traditionally examined within a weak turbulence framework—utilizing a truncated perturbative expansion in inverse powers of viscosity $\nu$ in the forced Navier-Stokes equation. Various phenomenological models have also been aligned with experimental observations, as discussed in the recent review [46].

However, the comprehensive turbulence theory requires solving the Navier-Stokes equations in the strong coupling limit $\nu \to 0$, the direct opposite of weak coupling. The universality of strong turbulence, with or without random forcing, poses the initial question in constructing such a theory.

Direct Numerical Simulation (DNS) data for energy decay in turbulent flows, detailed in [46], suggest a decay of the dissipation rate $E\sim t^{-n}$ with $n\approx 1.25$ or $1.60$ depending on the initial conditions (finite total momentum or zero total momentum but finite total angular momentum, see [46]). Thus, two universality classes of decaying turbulence have been identified.

It remains unclear which, if any, of the data in [40,46] reached the homogeneous isotropic turbulence limit corresponding to our regime. Moreover, stochastic forces added to the Navier-Stokes equation in simulations might contaminate the natural decay of turbulence. These forces are intended to initiate and enhance the spontaneous stochasticity of turbulent flow. However, in our theory [30], this inherent stochasticity is related to a dual quantum system and is discrete.

The Gaussian forcing can distort these quantum stochastic phenomena by stirring the flow ubiquitously and constantly. When the forcing is switched off, allowing the turbulence to decay towards a universal stage, energy dissipation should occur inside vorticity structures deeply embedded in the flow by the pure turbulent dynamics we study.

The forthcoming calculation supports the above relationship between singular vorticity distribution and energy flow. In the pure Navier-Stokes scenario, the energy balance reduces to energy dissipation by enstrophy within the bulk, counterbalanced by energy input from boundary forces (e.g., a large sphere encompassing the flow).

The general identity derived from the Navier-Stokes equations by multiplying both sides by $\overrightarrow{v}$ and averaging over an ensemble of stochastic solutions is:

$$\begin{array}{ccc}& & {\partial }_t〈E〉+{\int }_V{d}^{3}r〈\nu {\overrightarrow{\omega }}^2〉=\hfill \\ & & -{\int }_V{d}^{3}r{\partial }_{\beta }〈v_{\beta }\left(p+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.v_{\alpha }^2\right)+\nu v_{\alpha }({\partial }_{\beta }{v}_{\alpha }-{\partial }_{\alpha }{v}_{\beta })〉\hfill \end{array}$$

(1)

Applying the Stokes theorem, the right side reduces to the flow through the boundary $\partial V$ of the integration region V. The left side represents the dissipation within this volume:

$$\begin{array}{c}\hfill {\partial }_t〈E〉+{\mathcal{E}}_V=-{\int }_{\partial V}d\overrightarrow{\sigma }·〈\overrightarrow{v}\left(p+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.v_{\alpha }^2\right)+\nu \overrightarrow{\omega }\times \overrightarrow{v}〉\end{array}$$

(2)

This identity is valid for any volume, with the left side indicating viscous dissipation inside V and the right side representing the energy flow through the boundary $\partial V$.

Should a finite collection of vortex structures exist within the bulk, expanding this volume to infinite sphere results in the $\overrightarrow{\omega }\times \overrightarrow{v}$ term disappearing, as no vorticity persists at infinity.

Additionally, the velocity dictated by the Biot-Savart law diminishes as $|\overrightarrow{r}{|}^{-3}$ at infinity, so only the $\overrightarrow{v}p$ term remains significant:

$$\begin{array}{c}\hfill {\partial }_t〈E〉+〈{\mathcal{E}}_V〉\to -{\int }_{\partial V}d\overrightarrow{\sigma }·〈\overrightarrow{v}p〉\end{array}$$

(3)

This representation of energy flow will remain finite even as the sphere expands if the pressure scales as $p\to -\overrightarrow{F}·\overrightarrow{r}$, where $\overrightarrow{F}$ is the local force at any given point on a large sphere:

$$\begin{array}{c}\hfill {\partial }_t〈E〉+〈\mathcal{E}〉=f_{\alpha }{\displaystyle \underset{R\to \infty }{lim}}{R}^3{\int }_{S_2}{n}_{\alpha }{n}_{\beta }〈v_{\beta }\left(R\overrightarrow{n}\right)〉\end{array}$$

(4)

What about the Kolmogorov energy flow? It persists within any finite volume surrounding the set of vortexes:

$$\begin{array}{ccc}& & {\partial }_t〈E〉+〈\mathcal{E}〉=-{\int }_V{d}^{3}r〈v_{\beta }{\partial }_{\beta }p+v_{\alpha }{v}_{\beta }{\partial }_{\beta }{v}_{\alpha }〉=\hfill \\ & & -{\int }_V{d}^{3}r〈v_{\alpha }{v}_{\beta }{\partial }_{\beta }{v}_{\alpha }〉-{\int }_{\partial V}d\overrightarrow{\sigma }·〈\overrightarrow{v}p〉\hfill \end{array}$$

(5)

The triple velocity term in the last equation describes the Kolmogorov energy flow inside the volume V, and the second term represents the energy flow through the boundary.

Without a finite force $\overrightarrow{F}$ acting on the boundary, such as with periodic boundary conditions, the boundary integral would be absent, and the Kolmogorov relation would be fully applicable:

$$\begin{array}{c}\hfill {\partial }_t〈E〉+〈v_{\alpha }{v}_{\beta }{\partial }_{\beta }{v}_{\alpha }〉=-〈\mathcal{E}〉/V;\end{array}$$

(6)

This relation, alongside spatial symmetry properties in ${\mathbb{R}}_d$, leads to the Kolmogorov three-point correlation in a steady state ${\partial }_t〈E〉=0$:

$$\begin{array}{ccc}\hfill & & 〈v_{\alpha }\left({\overrightarrow{r}}_0\right)v_{\beta }\left({\overrightarrow{r}}_0\right)v_{\gamma }(\overrightarrow{r}+{\overrightarrow{r}}_0)〉=\hfill \\ & & \frac{\mathcal{E}}{(d-1)(d+2)V}\left({\delta }_{\alpha \gamma }{r}_{\beta }+{\delta }_{\beta \gamma }{r}_{\alpha }-\frac{2}{d}{\delta }_{\alpha \beta }{r}_{\gamma }\right);\hfill \end{array}$$

(7)

In the conventional approach to the turbulence problem, periodic Gaussian random forces $\overrightarrow{F}(\overrightarrow{r},t)$ are added to the Navier-Stokes equations in the conventional approach, based on time averaging:

$$\begin{array}{ccc}& & 〈{\mathcal{E}}_V〉=-{\int }_V{d}^{3}r〈v_{\beta }{\partial }_{\beta }p-v_{\beta }{f}_{\beta }(\overrightarrow{r},t)+v_{\alpha }{v}_{\beta }{\partial }_{\beta }{v}_{\alpha }〉=\hfill \\ & & {\int }_V{d}^{3}r〈v_{\beta }{f}_{\beta }(\overrightarrow{r},t)〉\hfill \end{array}$$

(8)

As the force becomes uniformly distributed across space, we derive another definition with $\mathcal{E}=\overrightarrow{F}·\overrightarrow{P}$, where $\overrightarrow{P}={\int }_V{d}^{3}r\overrightarrow{v}$ is the total momentum.

The phenomenon of turbulence we study exhibits a universal spontaneous stochasticity that does not depend on boundary conditions.

As long as energy flows from the boundaries, confined turbulence in the middle will dissipate this energy through singular vortex tubes. This spontaneous stochasticity results from the random distribution of these singular tubes within the volume of the velocity flow [31]. The dual picture from our recent theory [30] represents these by random gaps in the momentum curve $\overrightarrow{P}\left(\theta \right)$, as we shall discuss in the following sections.

The relation between the energy pumping at the boundary and the distribution of vortex blobs in the bulk follows from the Biot-Savart integral:

$$\begin{array}{c}\hfill \overrightarrow{v}\left(r\right)=-\overrightarrow{\nabla }\times \int d^3{r}^{\prime }\frac{\overrightarrow{\omega }\left(r^{\prime }\right)}{4\pi |r-r^{\prime }|}\end{array}$$

(9)

Generally, a gradient of harmonic potential is added to the Biot-Savart integral, dependent on the boundary conditions. We consider the velocity decaying at infinity, thus not adding such a term.

The net linear momentum $〈\int d^{3}r\overrightarrow{v}\left(0\right)·\overrightarrow{v}\left(\overrightarrow{r}\right)〉$, generally, is not zero in our theory, as we impose no such restriction. This nonvanishing linear moment places our theory in the most general $k^2$ (or Saffman) class.

On a large sphere $\partial V$ with radius $R\to \infty$,

$$\begin{array}{ccc}& & {\displaystyle \underset{R\to \infty }{lim}}{R}^3\overrightarrow{v}\left(R\overrightarrow{n}\right)\propto \hfill \\ & & \frac{1}{4\pi }{\displaystyle \sum _{\mathrm{blobs}}}{\displaystyle \underset{\mathrm{blob}}{\int }}{d}^3{r}^{\prime }\overrightarrow{\omega }\left(r^{\prime }\right)\times \left({\overrightarrow{r}}^{\prime }-<\overrightarrow{r}{>}_{\mathrm{blob}}\right);\hfill \end{array}$$

(10)

Here $<\overrightarrow{r}{>}_{\mathrm{blob}}$ is the geometric center of each blob. Substituting this into the identity (4), we directly relate the energy pumped by the forces at the boundary with the blob’s dipole moments of vorticity.

No forcing inside the flow is needed for this energy pumping; the energy flow starts at the boundary and propagates to numerous singular vorticity blobs, where it is finally dissipated. The distribution of these vorticity blobs is all we need for the turbulence theory. The forcing is required only as a boundary condition at infinity.

Another critical comment: with the velocity correlations growing with distance by the approximate K41 law, even the forcing at the remote boundary would influence the potential part of velocity in bulk. This boundary influence makes the energy cascade picture non-universal; it may depend upon the statistics of the random forcing.

Two asymptotic regimes manifesting this non-universality were observed for the energy spectrum $E(k,t)$: one for initial spectrum $E(k,0)\propto k^2$ and another for $E(k,0)\propto k^4$. The potential velocity part differs for these regimes, as the second one adds a constant velocity to the Biot-Savart integral to cancel the total momentum $\int d^{3}r\overrightarrow{v}$. In the general case, it will be a harmonic potential flow with certain boundary conditions at infinity, with explicit continuous dependence of the boundary forces. The most general case is the $k^2$ class, which does not require any restrictions.

Only the statistics of the rotational part of velocity, i.e., vorticity, could reach some universal regime independent of the boundary conditions at infinity. Certain discrete universality classes could exist as it is common in critical phenomena.

Unlike the potential part of velocity, vorticity is localized in singular regions—tubes and sheets, filling the space, as observed in numerical simulations. The potential part of velocity drops in the loop equations, and the remaining stochastic motion of the velocity circulation is equivalent to the vorticity statistics. Therefore, our solutions [30] of the loop equations [28,29] describe the internal stochastization of the decaying turbulence by a dual discrete system.

4.2. Mathematical Introduction. The loop equation and its solution

We derived a functional equation for the so-called loop average or Wilson loop in turbulence in the early nineties. All the references to our previous works can be found in a recent review paper [29].

The path to an exact solution by a dimensional reduction in this equation was proposed in the 1993 paper but has just been explored (see Figure ??). At the time, we could not compare a theory with anything but crude measurements in physical and numerical experiments at modest Reynolds numbers. All these experiments agreed with the K41 scaling, so the exotic equation based on unjustified methods of quantum field theory was premature. The specific prediction of the Loop equation, namely the Area law, could not be verified in DNS at the time with existing computer power.

The situation has changed over the last decades. No alternative microscopic theory based on the Navier-Stokes equation emerged, but our understanding of the strong turbulence phenomena grew significantly. On the other hand, the loop equations technology in the gauge theory also advanced over the last decades. The correspondence between the loop space functionals and the original vector fields was better understood, and various solutions to the gauge loop equations were found. In particular, the momentum loop equation was developed, similar to our momentum loop used below [26,27]. Recently, some numerical methods were found to solve loop equations beyond perturbation theory [2,20,21]. The loop dynamics was extended to quantum gravity, where it was used to study nonperturbative phenomena [4,50].

All these old and new developments made loop equations a major nonperturbative approach to gauge field theory. So, it is time to revive the hibernating theory of the loop equations in turbulence, where these equations are much simpler. The latest DNS [3,18,19,42] with Reynolds numbers of tens of thousands revealed and quantified violations of the K41 scaling laws. These numerical experiments are in agreement with so-called multifractal scaling laws [47].

Theoretically, we studied the loop equation in the confinement region (large circulation over large loop C) and justified the Area law suggested in ’93 on heuristic arguments. This law says that the tails of velocity circulation PDF in the confinement region are functions of the minimal area inside this loop. It was verified in DNS a few years ago [19], which triggered the further development of the geometric theory of turbulence[3,29,42]. In particular, the Area law was justified for flat and quadratic minimal surfaces, and an exact scaling law in confinement region $\mathsf{\Gamma }\propto \sqrt{Area}$ was derived [29]. The area law was verified with better precision in [18].

In the previous paper, [30], we have found a family of exact solutions of the loop equation for decaying turbulence [28,29]. This family describes a fixed trajectory of solutions with the universal time decay factor. The solutions are formulated in terms of the Wilson loop or loop average

$$\begin{array}{ccc}\hfill & & \mathsf{\Psi }(\gamma ,C)={〈exp\left(\frac{ı\gamma }{\nu }\oint d\overrightarrow{C}\left(\theta \right)·\overrightarrow{v}\left(\overrightarrow{C}\left(\theta \right)\right)\right)〉}_{init};\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill & & \mathsf{\Psi }(\gamma ,C)\Rightarrow {〈exp\left(\frac{ı\gamma }{\nu }\oint d\overrightarrow{C}\left(\theta \right)·\overrightarrow{P}\left(\theta \right)\right)〉}_{\mathbb{E}};\hfill \end{array}$$

(12)

In the first equation (the definition), the averaging $〈\cdots 〉$ goes over initial data for the solutions of the Navier-Stokes equation for velocity field $\overrightarrow{v}\left(\overrightarrow{r}\right)$. In the second one (the solution), the averaging goes over the space of solutions $\overrightarrow{P}\left(\theta \right)$ of the loop equation [30]. We choose in this paper the parametrization of the loop with $\xi =\frac{\theta }{2\pi }$ to match with the fermionic coordinates below (the parametrization is arbitrary, in virtue of parametric invariance of the loop dynamics).

The loop equation for the momentum loop $\overrightarrow{P}\left(\theta \right)$ follows from the Navier-Stokes equation for $\overrightarrow{v}$

$$\begin{array}{ccc}& & {\partial }_t{v}_{\alpha }=\nu {\partial }_{\beta }{\omega }_{\beta \alpha }-v_{\beta }{\omega }_{\beta \alpha }-{\partial }_{\alpha }\left(p+\frac{v_{\beta }^2}{2}\right);\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & {\partial }_{\alpha }{v}_{\alpha }=0;\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & {\omega }_{\beta \alpha }={\partial }_{\beta }{v}_{\alpha }-{\partial }_{\alpha }{v}_{\beta }\hfill \end{array}$$

(15)

After some transformations, replacing velocity and vorticity with the functional derivatives of the loop functional, we found the following momentum loop equation in [29,30]

$$\begin{array}{ccc}\hfill & & \nu {\partial }_t\overrightarrow{P}=-{\gamma }^2{\left(\mathsf{\Delta }\overrightarrow{P}\right)}^2\overrightarrow{P}+\hfill \end{array}$$

$$\begin{array}{ccc}& & \mathsf{\Delta }\overrightarrow{P}\left({\gamma }^2\overrightarrow{P}·\mathsf{\Delta }\overrightarrow{P}+ı\gamma \left(\frac{{(\overrightarrow{P}·\mathsf{\Delta }\overrightarrow{P})}^2}{\mathsf{\Delta }{\overrightarrow{P}}^2}-{\overrightarrow{P}}^2\right)\right);\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & \overrightarrow{P}\left(\theta \right)\equiv \frac{\overrightarrow{P}\left({\theta }^{+}\right)+\overrightarrow{P}\left({\theta }^{-}\right)}{2};\hfill \end{array}$$

(17)

$$\begin{array}{ccc}& & \mathsf{\Delta }\overrightarrow{P}\left(\theta \right)\equiv \overrightarrow{P}\left({\theta }^{+}\right)-\overrightarrow{P}\left({\theta }^{-}\right);\hfill \end{array}$$

(18)

The momentum loop has a discontinuity $\mathsf{\Delta }\overrightarrow{P}\left(\theta \right)$ at every parameter $0<\theta \le 1$, making it a fractal curve in complex space ${\mathbb{C}}_d$. The details can be found in [29,30]. We will skip the arguments $t,\theta$ in these loop equations, as there is no explicit dependence of these equations on either of these variables. This Anzatz (12) represents a plane wave in loop space, solving the loop equation for the Wilson loop due to the lack of direct dependence of the loop operator on the shape of the loop.

The superposition of these plane wave solutions would solve the Cauchy problem in loop space: find the stochastic function $\overrightarrow{P}\left(\theta \right)$ at $t=0$, providing the initial velocity field distribution. Formally, the initial distribution $W_0\left[P\right]$ of the momentum field $\overrightarrow{P}\left(\theta \right)$ is given by inverse functional Fourier transform.

$$\begin{array}{ccc}& & W_0\left[P\right]=\int DC{\delta }^3\left(\overrightarrow{C}\left[0\right]\right)\mathsf{\Psi }{(\gamma ,C)}_{t=0}\hfill \\ & & exp\left(-\frac{ı\gamma }{\nu }\int d\overrightarrow{C}\left(\theta \right)·\overrightarrow{P}\left(\theta \right)\right)\hfill \end{array}$$

(19)

In Appendix A, we solve the Cauchy problem for an exact stationary solution of the Navier-Stokes equation corresponding to the global rotation with the Gaussian random rotation matrix. The stochastic function $\overrightarrow{P}\left(\theta \right)$ in (A8), (A14) has a nontrivial Gaussian distribution with discontinuity $\mathsf{\Delta }P\left(\theta \right)$ related to slow $1/n$ decay of its Fourier expansion on the parametric circle. This is the simplest example of the fractal curve we study below in a decaying solution of the loop equation.

Rather than solving the Cauchy problem, we are looking for an attractor: the fixed trajectory for $\overrightarrow{P}(\theta ,t)$ with some universal probability distribution related to the decaying turbulence statistics.

The following transformation reveals the hidden scaling invariance of decaying turbulence

$$\begin{array}{c}\hfill \overrightarrow{P}=\sqrt{\frac{\nu }{2(t+t_0)}}\frac{\overrightarrow{F}}{\gamma }\end{array}$$

(20)

The new vector function $\overrightarrow{F}$ satisfies an equation

$$\begin{array}{ccc}\hfill & & 2{\partial }_{\tau }\overrightarrow{F}=\left(1-{\left(\mathsf{\Delta }\overrightarrow{F}\right)}^2\right)\overrightarrow{F}+\hfill \end{array}$$

$$\begin{array}{ccc}& & \mathsf{\Delta }\overrightarrow{F}\left({\gamma }^2\overrightarrow{F}·\mathsf{\Delta }\overrightarrow{F}+ı\gamma \left(\frac{{(\overrightarrow{F}·\mathsf{\Delta }\overrightarrow{F})}^2}{\mathsf{\Delta }{\overrightarrow{F}}^2}-{\overrightarrow{F}}^2\right)\right);\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & \tau =log(t+t_0)\hfill \end{array}$$

(22)

This equation is invariant under translations of the new variable $\tau =log(t+t_0)$, corresponding to the rescaling/translation of the original time.

$$\begin{array}{c}\hfill t\Rightarrow \lambda t+(\lambda -1)t_0\end{array}$$

(23)

There are two consequences of this invariance.

• There is a fixed point for $\overrightarrow{F}$.
• The approach to this fixed point is exponential in $\tau$, which is power-like in original time.

Both of these properties were used in [30]: the first one was used to find a fixed point, and the second one was used to derive the spectral equation for the anomalous dimensions ${\lambda }_i$ of decay $t^{-{\lambda }_i}$ of the small deviations from these fixed points. In this paper, we only consider the fixed point, leaving the exciting problem of the spectrum of anomalous dimensions for future research.

4.3. The big and small Euler ensembles

Let us remember the basic properties of the fixed point for $\overrightarrow{F}$ in [30]. It is defined as a limit $N\to \infty$ of the polygon ${\overrightarrow{F}}_0\cdots {\overrightarrow{F}}_N={\overrightarrow{F}}_0$ with the following vertices

$$\begin{array}{ccc}\hfill & & {\overrightarrow{F}}_k=\frac{\left\{cos\left({\alpha }_k\right),sin\left({\alpha }_k\right),icos\left(\frac{\beta }{2}\right)\right\}}{2sin\left(\frac{\beta }{2}\right)};\hfill \end{array}$$

(24)

$$\begin{array}{ccc}& & {\theta }_k=\frac{k}{N};\beta =\frac{2\pi p}{q};N\to \infty ;\hfill \end{array}$$

(25)

$$\begin{array}{ccc}& & {\alpha }_{k+1}={\alpha }_k+{\sigma }_k\beta ;{\sigma }_k=\pm 1,\beta \sum {\sigma }_k=2\pi pr;\hfill \end{array}$$

(26)

The parameters $\widehat{\mathsf{\Omega }},p,q,r,{\sigma }_0\cdots {\sigma }_N={\sigma }_0$ are random, making this solution for $\overrightarrow{F}\left(\theta \right)$ a fixed manifold rather than a fixed point. We suggested calling this manifold the big Euler ensemble of just the Euler ensemble.

It is a fixed point of (21) with the discrete version of discontinuity and principal value:

$$\begin{array}{ccc}& & \mathsf{\Delta }\overrightarrow{F}\equiv {\overrightarrow{F}}_{k+1}-{\overrightarrow{F}}_k;\hfill \end{array}$$

(27)

$$\begin{array}{ccc}& & \overrightarrow{F}\equiv \frac{{\overrightarrow{F}}_{k+1}+{\overrightarrow{F}}_k}{2}\hfill \end{array}$$

(28)

Both terms of the right side (21) vanish; the term proportional to $\mathsf{\Delta }\overrightarrow{F}$ and the term proportional to $\overrightarrow{F}$. Otherwise, we would have $\overrightarrow{F}\Vert \mathsf{\Delta }\overrightarrow{F}$, leading to zero vorticity [30]. The ensemble of all the different solutions is called the big Euler ensemble. The integer numbers ${\sigma }_k=\pm 1$ came as the solution of the loop equation, and the requirement of the rational $\frac{p}{q}$ came from the periodicity requirement.

We can use integration (summation) by parts to write the circulation as follows (in virtue of periodicity):

$$\begin{array}{ccc}\hfill & & \oint d\overrightarrow{C}\left(\theta \right)·\overrightarrow{P}\left(\theta \right)=-\oint d\overrightarrow{P}\left(\theta \right)·\overrightarrow{C}\left(\theta \right);\hfill \end{array}$$

(29)

$$\begin{array}{ccc}& & {\displaystyle \sum _k}\mathsf{\Delta }{C}_k{\overrightarrow{P}}_k=-\oint \mathsf{\Delta }{\overrightarrow{P}}_k·{\overrightarrow{C}}_k;\hfill \end{array}$$

(30)

A remarkable property of this solution $\overrightarrow{P}(\theta ,t)$ of the loop equation is that even though it satisfies the complex equation and has an imaginary part, the resulting circulation (29) is real! The imaginary part of the $\overrightarrow{P}(\theta ,t)$ does not depend on $\theta$ and thus drops from the integral $\oint d\overrightarrow{P}(\theta ,t)·\overrightarrow{C}\left(\theta \right)$.

There is, in general, a larger manifold of periodic solutions to the discrete loop equation, which has all three components of ${\overrightarrow{F}}_k$ complex and varying along the polygon.

We could not find a global parametrization of such a solution1. Instead, we generated it numerically by taking a planar closed polygon and evolving its vertices ${\overrightarrow{F}}_k$ by a stochastic process in the local tangent plane to the surface of the equations in multi-dimensional complex space.

We could not submit such a solution to an extra restriction $\mathbf{Im}{\overrightarrow{F}}_k=\mathit{const}$ needed to make circulation real. We cannot prove that such a general solution does not exist but rather take the Euler ensemble as a working hypothesis and investigate its properties.

This ensemble can be solved analytically in the statistical limit and has nice physical properties, matching the expected behavior of the decaying turbulence solution.

We assign equal weights to all elements of this set; we call this conjecture the ergodic hypothesis. This prescription is similar to assigning equal weights to each triangulation of curved space with the same topology in dynamically triangulated quantum gravity [1]. Mathematically, this is the most symmetric weight assignment, and there are general expectations that various discrete theories converge into the same symmetry classes of continuum theories in the statistical limit. This method works remarkably well in two dimensions [6,11,14], providing the same correlation functions as continuum gravity (Liouville theory [22]).

The fractions $\frac{p}{q}$ with fixed denominator are counted by Euler totient function $\phi \left(q\right)$ [16]

$$\begin{array}{ccc}\hfill & & \phi \left(q\right)={\displaystyle \sum _{\begin{array}{c}p=1\\ (p,q)\end{array}}^{q-1}}1=q{\displaystyle \prod _{p|q}}\left(1-\frac{1}{p}\right);\hfill \end{array}$$

(31)

For example $\phi \left(16\right)=8$ and $\phi \left(17\right)=16$.

In some cases, one can analytically average over spins $\sigma$ in the big Euler ensemble, reducing the problem to computations of averages over the small Euler ensemble $\mathcal{E}\left(N\right):N,p,q,r$ with the measure induced by averaging over the spins in the big Euler ensemble.

5. The Fermi ring and its continuum limit

In this paper, we perform this averaging over $\sigma$ analytically, without any approximations, reducing it to a partition function of a certain quantum mechanical system with Fermi particles. The Quantum Trace Theorem, establishing this connection (A29) is proven in Appendix B. This partition function is calculable using a WKB approximation in the statistical limit $N\to \infty$. As we shall shortly see, in the continuum limit $N\to \infty$, the accumulated numbers of Fermi particles ${\nu }_k=1$ and Dirac holes ${\nu }_k=0$ tend to some classical function ${\sum }_{l<k}(2{\nu }_l-1)\to \alpha \left(\xi \right)$ of "position" $\xi =\frac{k}{N}$, leading to the exact solution.

Specifically, as we prove in Appendix C, the loop functional in continuum limit $N\to \infty$ reduces to a quantum mechanical path integral (A67) over the Fermion density $\alpha \left(\xi \right)$. The effective Action for this path integral is given by circulation $ı\mathsf{\Gamma }$ expressed in terms of this density plus a quadratic functional corresponding to Brownian measure (A64) for $\alpha \left(\xi \right)$.

Thus, there are three sources of fluctuations: There is a phase factor related to the circulation, there is a Brownian positive distribution of the trajectory $\alpha \left(\xi \right)$(Gaussian measure), and finally, the circulation depends on the random fraction $\frac{p}{q}$ distributed according to the small Euler ensemble. The continuum limit of the latter distribution is derived in Appendix E, using new cotangent sums derived in the previous paper [30].

This is a new kind of quantum mechanical system with a complex Action, reflecting the irreversibility of turbulence. The square root of viscosity enters the denominator of the effective Action, like a coupling constant. As we argued in the previous paper, the turbulent limit of our theory corresponds to

$$\begin{array}{c}\hfill \nu \to \frac{\tilde{\nu }}{N^2}\to 0\end{array}$$

(32)

The parameter $\tilde{\nu }$ remains a free parameter of our theory, playing the role of turbulent viscosity. In particular, there is an anomalous energy dissipation in this limit

$$\begin{array}{ccc}& & \mathcal{E}\left(t\right)=\frac{\tilde{\nu }}{t^2}F\left(k_0^2\tilde{\nu }t\right);\hfill \end{array}$$

(33)

$$\begin{array}{ccc}& & F\left(0\right)=\frac{{\pi }^2}{864\zeta \left(3\right)};\hfill \end{array}$$

(34)

Here $k_0$ is the lower cutoff of the energy spectrum (to be discussed below). The decaying spectrum at small wave vectors $k<k_0$ is related to the energy pumping at the initial moment $t=t_0$ and is time-independent.

The dimensionless parameter $N\to \infty$ plays the role of the Reynolds number. We derive the turbulent limit $N=2m\to \infty$ without any $1/N$ corrections. This solution will then apply to the inertial range of the physical turbulence in a system with a finite but large Reynolds number.

6. Instanton in the path integral

This classical equation for our path integral reads (with $\mathsf{\Omega }\in O\left(3\right)$ being a random rotation matrix):

$$\begin{array}{ccc}\hfill & & {\alpha }^{\prime \prime }=-ı\kappa \left({\mathcal{C}}_{\mathsf{\Omega }}^{\prime }exp\left(ı\alpha \right)+{\left({\mathcal{C}}_{\mathsf{\Omega }}^{\prime }\right)}^{★}exp\left(-ı\alpha \right)\right);\hfill \end{array}$$

(35)

$$\begin{array}{ccc}& & \kappa =\frac{1}{2\pi y\sqrt{X}\sqrt{2\nu t}};\hfill \end{array}$$

(36)

$$\begin{array}{ccc}& & {\mathcal{C}}_{\mathsf{\Omega }}\left(\theta \right)=\overrightarrow{C}\left(\theta \right)·\widehat{\mathsf{\Omega }}·\{ı,1,0\};\hfill \end{array}$$

(37)

The parameter $\kappa$ is distributed according to the distributions (A94), (A100) of the variables $X,y$ in a small Euler ensemble in the statistical limit.

This complex equation leads to a complex classical solution (instanton). It simplifies for $z=exp\left(ı\alpha \right)$:

$$\begin{array}{ccc}& & z^{\prime \prime }=\frac{{\left(z^{\prime }\right)}^2}{z}+\kappa \left({\mathcal{C}}_{\mathsf{\Omega }}^{\prime }{z}^2+{\left({\mathcal{C}}_{\mathsf{\Omega }}^{\prime }\right)}^{★}\right);\hfill \end{array}$$

(38)

$$\begin{array}{ccc}& & z\left(0\right)=z\left(1\right)=1\hfill \end{array}$$

(39)

This equation cannot be analytically solved for arbitrary periodic function $C_{\mathsf{\Omega }}^{\prime }\left(\xi \right)$.

The weak and strong coupling expansions by $\kappa$ are straightforward. At small $\kappa$

$$\begin{array}{ccc}& & z\left(\xi \right)\to 1+2\kappa \left(-A\xi +{\int }_0^{\xi }\mathbf{Re}{\mathcal{C}}_{\mathsf{\Omega }}\left({\xi }^{\prime }\right)d{\xi }^{\prime }\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & +O\left({\kappa }^2\right);\hfill \end{array}$$

(40)

$$\begin{array}{ccc}& & A={\int }_0^1\mathbf{Re}{\mathcal{C}}_{\mathsf{\Omega }}\left({\xi }^{\prime }\right)d{\xi }^{\prime }\hfill \end{array}$$

(41)

At large $\kappa$

$$\begin{array}{ccc}& & z\left(\xi \right)\to ıexp\left(-ıarg{\mathcal{C}}_{\mathsf{\Omega }}^{\prime }\left(\xi \right)\right)=ı\frac{\left|{\mathcal{C}}_{\mathsf{\Omega }}^{\prime }\left(\xi \right)\right|}{{\mathcal{C}}_{\mathsf{\Omega }}^{\prime }\left(\xi \right)}\hfill \end{array}$$

(42)

This solution is valid at intermediate $\xi$, not too close to the boundaries $\xi =(0,1)$. In the region near the boundaries $\xi (1-\xi )\ll \frac{1}{\sqrt{\kappa }}$, the following asymptotic agrees with the classical equation

$$\begin{array}{ccc}\hfill & & z\to 1\pm ı\xi \sqrt{2\kappa \mathbf{Re}{C}_{\mathsf{\Omega }}^{\prime }\left(0\right)}+O\left({\xi }^2\right);\hfill \end{array}$$

(43)

$$\begin{array}{ccc}& & z\to 1\pm ı(1-\xi )\sqrt{2\kappa \mathbf{Re}{C}_{\mathsf{\Omega }}^{\prime }\left(1\right)}+O\left({(1-\xi )}^2\right);\hfill \end{array}$$

(44)

One can expand in small or large values of $\kappa$ and use the above distributions for $X,y$ term by term.

As it was noticed above, the viscosity $\nu =\frac{\tilde{\nu }}{N^2}\to 0$ in our theory. This limit makes $\kappa \sim N\to \infty$, justifying the strong coupling limit for the Wilson loop solution.

The classical limit of the circulation in exponential of (A67)

$$\begin{array}{c}\hfill {\int }_0^{1}d\xi \mathbf{Im}\left({\mathcal{C}}_{\mathsf{\Omega }}^{\prime }\left(\xi \right)exp\left(ı\alpha \left(\xi \right)\right)\right)\to {\int }_0^{1}d\xi \left|{\mathcal{C}}_{\mathsf{\Omega }}^{\prime }\left(\xi \right)\right|\end{array}$$

(45)

becomes a positive definite function of the rotation matrix $\mathsf{\Omega }$. At large $\kappa$ the leading contribution will come from the rotation matrix minimizing this functonal.

Let us think about the physical meaning of this finding. We have just found the density of our Fermi particles on a parametric circle

$$\begin{array}{c}\hfill \alpha \left(\xi \right)=\frac{\pi }{2}-arg{\mathcal{C}}_{\mathsf{\Omega }}^{\prime }\left(\xi \right)\end{array}$$

(46)

This density does not fluctuate in a turbulent limit, except near the endpoints $\xi \to 0,\xi \to 1$. In the vicinity of the endpoints, there is a different asymptotic solution (43) for $\alpha \to (z-1)/ı$.

Compute the Wilson loop for a specific loop, say, the circle, is an interesting problem, but there is a simpler quantity. In the next section we are considering an important calculable case of the vorticity correlation function, where the full solution in quadratures is available. This function has been directly observed in grid turbulence experiments [1,10] more than half a century ago and is being studied in modern large-scale real and numerical experiments[23,40,46].

This is the vorticity correlation function [29], corresponding to the loop C backtracking between two points in space ${\overrightarrow{r}}_1=0,{\overrightarrow{r}}_2=\overrightarrow{r}$, with the vorticity operators are inserted at these two points (see [30] for details and the justification). The Fourier transform of this function describes the decaying energy spectrum, also measured experimentally.

In Appendixes Appendix F, Appendix G, Appendix H, we express this correlation function as a particular case of the second variation of the loop functional. Then, in Appendixes Appendix I, Appendix J, we compute the path integral in the leading WKB approximation. This is a one-dimensional version of the instanton computations, familiar to the gauge theory experts. In the turbulent limit $N\to \infty ,\tilde{\nu }=\mathit{const}$, there are no higher order corrections to this WKB approximation of the path integral.

7. The decaying energy in finite system

The vorticity correlation in Fourier space doubles as an energy spectrum

$$\begin{array}{c}\hfill E(k,t)=4\pi {\overrightarrow{k}}^2〈\overrightarrow{v}·{\overrightarrow{v}}_k〉=4\pi 〈\overrightarrow{\omega }·{\overrightarrow{\omega }}_k〉\end{array}$$

(47)

The energy spectrum in a finite system with size L is bounded from below. At low $|\overrightarrow{k}|\le \pi /L$, the spectrum is no longer related to the turbulence but is given by the energy pumping by external forces at the boundaries.

This energy pumping [40,46] takes place at $t<t_0$, after which the pumping stops. At this moment, the energy spectrum is growing with wavevector by one of two possible laws (with P being the Birkhoff-Saffman invariant of the fluid and M being the Loitzansky invariant)

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}& E(k,t)\propto P{k}^2\hfill \\ & E(k,t)\propto M{k}^4\hfill \end{array}\right.\end{array}$$

(48)

The small k limit of the spectrum is time-independent as both P and M do not depend of time.

At $t>t_0$, without the forcing, the pumped energy dissipates at large k corresponding to smaller spatial scales of the hierarchy of vortex structures of all scales, ending with dissipative scales, or wavevectors $k\gg \pi /L$. After sufficient time, the universal regime kicks in, corresponding to the decaying turbulence. It is implied that a large amount of energy was pumped in, so it takes a long time to reach this decaying regime, corresponding to some fixed trajectory.

Our solution does not impose any restrictions on the SB invariant P and thus applies to the most general, first regime with $k^2$ spectrum at small k and some universal decay at large k, reflecting these distributed vortex structures.

The decaying energy $E_d\left(t\right)$, given by the part of the spectrum $k>k_0\sim 1/L$, has the following form

$$\begin{array}{ccc}\hfill & & E_d\left(t\right)={\int }_{k_0}^{\infty }dkE(k,t)\propto \frac{4\pi \tilde{\nu }}{t}{\int }_{k_0\sqrt{\tilde{\nu }t}}^{\infty }H\left(\kappa \right)d\kappa \hfill \end{array}$$

(49)

On top of the trivial decrease $\frac{\tilde{\nu }}{t}$, as prescribed by dimensional counting in an infinite system, there is some extra decrease related to the increase of the lower limit.

7.1. Computation of the energy dissipation

The energy in our theory does not have a finite statistical limit because of the contribution from the unknown potential part of velocity. This contribution could be infinite in the infinite system. Thus, we compute the energy by integrating its time derivative (i.e., the dissipation rate $-\mathcal{E}\left(t\right)$ ) given the zero energy left at infinite time.

$$\begin{array}{c}\hfill E\left(t\right)={\int }_t^{\infty }\mathcal{E}\left(t^{\prime }\right)d{t}^{\prime }\end{array}$$

(50)

This energy dissipation rate $\mathcal{E}\left(t^{\prime }\right)$ is calculable

$$\begin{array}{ccc}\hfill & & \mathcal{E}\left(t\right)\propto 4\pi \nu \int dk{k}^2〈\overrightarrow{\omega }\left(\overrightarrow{0}\right)·\overrightarrow{\omega }\left(\overrightarrow{k}\right)〉=\hfill \\ & & 4\pi \frac{\tilde{\nu }}{t^2}{\int }_{k_0\sqrt{\tilde{\nu }t}}^{\infty }{\kappa }^{2}H\left(\kappa \right)d\kappa \hfill \end{array}$$

(51)

In our theory, this integral has a finite limit in an infinite system ($k_0=0$).

This limit was computed in [30] in a slightly different grand canonical ensemble, where N was fluctuating with the weight $exp\left(-\mu N\right),\mu \to 0$.

With our current ensemble of fixed even $N\to \infty$ the results of [30] read:

$$\begin{array}{ccc}\hfill & & {\mathcal{E}}_{\infty }\left(t\right)=\frac{\tilde{\nu }}{t^2}\frac{{\pi }^2}{864\zeta \left(3\right)};\hfill \end{array}$$

(52)

In our present theory, the same quantity is given by the above integral at $k_0=0$

$$\begin{array}{c}\hfill {\mathcal{E}}_{\infty }\left(t\right)=4\pi \frac{\tilde{\nu }}{t^2}{\int }_0^{\infty }{\kappa }^{2}H\left(\kappa \right)d\kappa \end{array}$$

(53)

Comparing these two expressions, we get the normalization of $H\left(\kappa \right)$

$$\begin{array}{c}\hfill 4\pi {\int }_0^{\infty }{\kappa }^{2}H\left(\kappa \right)d\kappa =\frac{{\pi }^2}{864\zeta \left(3\right)}\end{array}$$

(54)

The integral on the left can be further reduced [37] to the following normalization condition:

$$\begin{array}{ccc}& & \mathcal{Z}=\frac{276480\zeta \left(3\right)\zeta \left(\frac{15}{2}\right)}{119\zeta \left(\frac{17}{2}\right)}{\int }_{{\mathsf{\Delta }}_1}^{{\mathsf{\Delta }}_2}d\mathsf{\Delta }\frac{P\left(\mathsf{\Delta }\right)}{Q\left(\mathsf{\Delta }\right)};\hfill \\ & & P\left(\mathsf{\Delta }\right)=(1-\mathsf{\Delta })\sqrt{r_0\left(\mathsf{\Delta }\right)}S\left(\mathsf{\Delta }\right)Q_{\alpha }(\mathsf{\Delta },1)\hfill \end{array}$$

(55)

$$\begin{array}{ccc}& & \left(\pi J\left(\mathsf{\Delta }\right)\left(r_0\left(\mathsf{\Delta }\right)+12\right)+L\left(\mathsf{\Delta }\right)\left(r_0\left(\mathsf{\Delta }\right)-6\right)\right);\hfill \end{array}$$

(56)

$$\begin{array}{ccc}& & Q\left(\mathsf{\Delta }\right)=L{\left(\mathsf{\Delta }\right)}^2\left(r_0\left(\mathsf{\Delta }\right)+12\right);\hfill \end{array}$$

(57)

This normalization constant $\mathcal{Z}$ can be used in equation (49) for the energy decay in a finite system. All the functions of $\mathsf{\Delta }$ were defined above.

As for the energy spectrum, this is not an independent function in our theory. Comparing the two expressions (47) and (51), we arrive at the following relation

$$\begin{array}{ccc}\hfill & & t^2\mathcal{E}\left(t\right)=4\pi \tilde{\nu }F\left(k_0\sqrt{\tilde{\nu }t}\right);\hfill \end{array}$$

(58a)

$$\begin{array}{ccc}& & F\left(\kappa \right)=\underset{\kappa }{\overset{\infty }{\int }}H\left(x\right)x^{2}dx,\hfill \end{array}$$

(58b)

$$\begin{array}{ccc}\hfill & & \sqrt{t}E(k,t)\propto 4\pi \tilde{\nu }\sqrt{\tilde{\nu }}H\left(k\sqrt{\tilde{\nu }t}\right);\hfill \end{array}$$

(58c)

Both the energy dissipation and the energy spectrum are related to the same function $H\left(\kappa \right)$, but the energy spectrum related to this function at large argument $\kappa =k\sqrt{\tilde{\nu }t}$, whereas the energy dissipation is related to integral of this function from small argument $\kappa =k_0\sqrt{\tilde{\nu }t}$ to infinity.

Our theory does not have the Saffmann part of the spectra; it only applies to an infinite system described by the region $k>k_0$. At the boundary of this region, we have the value of the spectrum (assuming $k_0\sqrt{\tilde{\nu }t}\ll 1$)

$$\begin{array}{c}\hfill E(k_0,t)\propto \frac{4\pi }{\sqrt{t}}\tilde{\nu }\sqrt{\tilde{\nu }}H\left(k_0\sqrt{\tilde{\nu }t}\right)\sim \frac{4\pi H\left(0\right)\tilde{\nu }\sqrt{\tilde{\nu }}}{\sqrt{t}}\end{array}$$

(59)

This will match the Saffman spectrum $P{k}^2$ at small k

$$\begin{array}{c}\hfill k\sim {\tilde{\nu }}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}{t}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}{P}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\end{array}$$

(60)

This boundary value decays as $t^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$; it is below the value $k_0\sim 1/L$, which is time-independent. We are not considering extremely large times such that $k_0\sqrt{\tilde{\nu }t}\gg 1$. At these times, the decay is over: there is not enough energy left for turbulence, and the whole pumped energy has dissipated.

We computed the universal function $H\left(\kappa \right)$ in the Appendix K and numerically integrated it to obtain $F\left(\kappa \right),\kappa =k\sqrt{\tilde{\nu }t}$. Asymptotically, at large $k\sqrt{\tilde{\nu }t}$

$$\begin{array}{ccc}& & E(k,t)\propto \frac{\tilde{\nu }\sqrt{\tilde{\nu }}}{\sqrt{t}}{\left(k\sqrt{\tilde{\nu }t}\right)}^{-\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$2$}\right.};\hfill \end{array}$$

(61)

$$\begin{array}{ccc}& & \mathcal{E}\left(t\right)\propto t^{-\raisebox{1ex}{$9$}\!\left/ \!\raisebox{-1ex}{$4$}\right.};\hfill \end{array}$$

(62)

$$\begin{array}{ccc}& & E\left(t\right)={\int }_t^{\infty }\mathcal{E}\left(t^{\prime }\right)d{t}^{\prime }\propto t^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}\hfill \end{array}$$

(63)

Figure 8. The effective index $n\left(t\right)=-\frac{t\mathcal{E}\left(t\right)}{E\left(t\right)}$ compared with asymptotic value $n(\infty )=\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.$.

Figure 8. The effective index $n\left(t\right)=-\frac{t\mathcal{E}\left(t\right)}{E\left(t\right)}$ compared with asymptotic value $n(\infty )=\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.$.

We computed the effective critical indexes

$$\begin{array}{ccc}& & n\left(t\right)=-t{\partial }_{t}logE=-\frac{t\mathcal{E}\left(t\right)}{E\left(t\right)};\hfill \end{array}$$

(64)

$$\begin{array}{ccc}& & s\left(\kappa \right)=-\kappa {\partial }_{\kappa }logH\left(\kappa \right)\hfill \end{array}$$

(65)

numerically in the Appendix K, using Mathematica^®. The accuracy is just 4-5 digits, but it can be easily improved by taking more CPU time once experimental data gets more precise. These curves are universal, and they change the regime before approaching their limits from below $n(\infty )=\frac{5}{4},s(\infty )=\frac{7}{2}$. These regime changes are due to the quantum effects (complex zeros of the zeta function contributing to the energy spectrum’s Mellin transform, as shown in Appendix K).

The experimental data [9,10] yields $n(\infty )\approx 1.25\pm 0.02$, which agrees with our theoretical prediction in Figure 8. Our universal curves for $n\left(t\right),s\left(\kappa \right)$ were computed directly from the analytic solution of the loop equation in the turbulent limit without any fitting parameters. It will be very interesting to compare these curves with more precise experiments (real or numerical).

7.2. The energy normalization problem

The above formulas do not specify the energy spectrum’s normalization, just the energy dissipation’s normalization. When one tries to recover the normalization of the energy spectrum, the following problem arises.

The normalization of the decaying energy (49) seems incompatible with its time derivative (51). In conventional turbulence models, the integral for dissipation diverges at large wavelengths, reflecting the singular vortex structures such as the Burgers vortex filaments [55] with viscous thickness $w\sim \sqrt{\nu }$. Integrating the square of vorticity in the Burgers vortex in the transverse plane, we get a large factor $1/\nu$; this compensation leads to anomalous dissipation $\mathcal{E}=\nu 〈{\overrightarrow{\omega }}^2〉$. independent of $\nu$.

In our dual theory, this factor of $\nu$ is compensated by a different mechanism: the vorticity is represented as a discontinuity at the curve $\overrightarrow{P}\left(\theta \right)$ in our solution: $\overrightarrow{\omega }\propto \overrightarrow{P}\times \mathsf{\Delta }\overrightarrow{P}$. Summing over a large number $N\to \infty$ of these discontinuities leads to the factor $N^2$ compensating small viscosity $\nu$ factor in front of the enstrophy. The enstrophy integral $\int dk{k}^2〈\overrightarrow{\omega }·{\overrightarrow{\omega }}_k〉$ converges at large k.

As a consequence, the vorticity correlation $〈\overrightarrow{\omega }·{\overrightarrow{\omega }}_k〉\sim 1/\nu$ is large at all k, not just at large k. Our limit $\nu \to 0,N\to \infty$ applies to the computation of integrated quantities such as total decaying energy but not to the energy spectrum as a function of wavelength.

The problem boils down to the following. The turbulent limit differs from the inviscid limit of the Navier-Stokes equation. In the turbulent limit, the average circulation $\mathsf{\Gamma }$ is much larger than viscosity, but the dimensional scales, determined by viscosity, stay finite.

To be more specific, the enstrophy in our system has the structure

$$\begin{array}{c}\hfill 〈{\overrightarrow{\omega }}^2〉=\frac{\mathcal{E}(t,|\mathsf{\Gamma }|/\nu )}{\nu }\end{array}$$

(66)

The dissipation $\mathcal{E}(t,|\mathsf{\Gamma }|/\nu )$ depends on time and the effective Reynolds number $\mathbf{Rey}=\left|\mathsf{\Gamma }\right|/\nu$ where $\mathsf{\Gamma }$ is a typical velocity circulation in our problem.

The turbulent limit corresponds to $\left|\mathsf{\Gamma }\right|\gg \nu$, rather than $\nu \to 0$. Mathematically, we are looking for a residue of the enstrophy at zero viscosity

$$\begin{array}{c}\hfill 〈{\overrightarrow{\omega }}^2〉\to \frac{\mathcal{E}(t,\infty )}{\nu }\end{array}$$

(67)

However, in the real physical world (or in the DNS), we shall use this formula for finite $\nu$ corresponding to the actual viscosity of water. We only use our dual theory to compute this residue ${\mathrm{Res}}_{\nu =0}〈{\overrightarrow{\omega }}^2〉=\mathcal{E}(t,\infty )$.

$$\begin{array}{c}\hfill 〈{\overrightarrow{\omega }}^2〉\to \frac{\mathcal{E}{(t,\infty )}_{dual}}{\nu }\end{array}$$

(68)

In particular, in the infinite system ($k_0=0$), we have the value (52) computed in the previous paper.

This comparison restores the normalization of the vorticity correlation and the energy spectrum

$$\begin{array}{ccc}\hfill & & 〈\overrightarrow{\omega }\left(\overrightarrow{0}\right)·\overrightarrow{\omega }\left(\overrightarrow{k}\right)〉=\frac{{\tilde{\nu }}^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}H\left(k\sqrt{\tilde{\nu }t}\right)}{\nu \sqrt{t}};\hfill \end{array}$$

(69)

$$\begin{array}{ccc}& & \sqrt{t}E(k,t)=\frac{4\pi {\tilde{\nu }}^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}H\left(k\sqrt{\tilde{\nu }t}\right)}{\nu }\hfill \end{array}$$

(70)

where $\nu$ is the physical viscosity of the fluid under consideration (say, water) and $\tilde{\nu }$ is the auxiliary scale parameter that comes from the dual theory in the turbulent limit.

Our theory has two unknown parameters: $\tilde{\nu }$ and $t_0$. The energy spectrum decreased faster than $k^{-3}$ so that the enstrophy integral $\int H\left(\kappa \right){\kappa }^{2}d\kappa$ converges at large $\kappa$ and so does the energy integral $\int H\left(\kappa \right)d\kappa$.

As we shall see in the next section, this fast decay of the energy spectrum also happens in DNS, in qualitative agreement with our decay ${\kappa }^{-\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$.

8. Comparing our theory with the DNS

As the first such test of our theory, we took the raw DNS data from [46], provided to us by the authors. This data is now available online [40] and can be downloaded without permission. We only compared the data corresponding to our $k^2$ case and restricted ourselves to four samples with the largest grid ${1024}^3$, labeled as sample $13,14,15,16$.

First, we verify the decay of the Reynolds number for each sample, as our theory only applies to the large Reynolds numbers. These plots are shown in Figure 9. As we see from these plots, all four Reynolds numbers are modest, but the sample 14 stands out as the closest to the strong turbulence we seek.

The next test is the effective length scale, which we define as

$$\begin{array}{ccc}\hfill & & L\left(t\right)=\frac{\int E(k,t)kdk}{\int E(k,t)k^{2}dk};\hfill \end{array}$$

(71)

The effective length $L\left(t\right)$ as a function of time is shown in Figure 10. The statistical equilibrium was not yet reached at $t<1000$, so we discarded this period. The late stages of decay where $L\left(t\right)\propto t$ correspond to low Reynolds number and do not agree with our theory: we interpret it as the non-turbulent stage of decay when the remaining energy is insufficient for the strong turbulence phase.

The next test is the energy decay curve $E\left(t\right)$ for the turbulent region of sample 14.

$$\begin{array}{ccc}\hfill & & E\left(t\right)=\int E(k,t)dk;\hfill \end{array}$$

(72)

$$\begin{array}{ccc}& & logE\left(t\right)\approx a+f(b+logL(t\left)\right)\hfill \end{array}$$

(73)

Our theory has two free parameters $t_0,\tilde{\nu }$. The fitting of $t_0$ is not trivial as we do not know which time range corresponds to the universal regime of the decay but still contains enough energy left for the strong turbulence. Following the suggestion of [51], we avoid fitting $t_0$ by investigating the decaying energy as a function of the decay length $L\left(t\right)$, which in our case scales as $\sqrt{\tilde{\nu }(t+t_0)}$. The precise definition was given in (71).

The parameters $a,b$ were fitted by nonlinear regression using "NonlinearModelFit" in Mathematica^®. The resulting log-log plot is shown in(Figure 11). The fit is perfect, with less than one percent of the standard deviation. This relation is approximately linear with the slope $-\frac{5}{2}$. Note also that the asymptotic index $-\frac{5}{2}$ comes as a ratio of the energy decay index $-\frac{5}{4}$ to the index $m=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ for $L\left(t\right)$ which we already tested. Let us now turn to the energy spectrum. The data is not as good here as the energy decay data for two reasons. First, each wave vector component has only $L=512$ independent values on a 1024 grid.

The energy spectrum is a function of the length of the wavevector, which is taking $O\left(L\right(L+1\left)\right(L+2)/6)$ different values between 0 and $L\frac{\sqrt{3}}{2}$. Unfortunately, the available data [40,46] aggregates the statistics at 512 equidistant bins in $|\overrightarrow{k}|$, reducing statistics. The large number of rank bins (by an equal number of data points in each bin) would give us much more information about the spectrum.

We have a scaling law $E(k,t)=\frac{H\left(k\sqrt{t+t_0}\right)}{\sqrt{t+t_0}}$, which means that the two-dimensional array of the data for $E(k,t)$ must collapse at one-dimensional subset.

We already saw the consequence of that collapse in the scaling law for $L\left(t\right)$. However, the low k part of the spectrum is discrete and corresponds to lattice artifacts.

We found the following method to avoid choosing the range of discrete wavelengths or fitting any scale parameters.

We consider the second moment of velocity, related to the energy spectrum by Fourier transform

$$\begin{array}{ccc}& & 〈\mathsf{\Delta }{\overrightarrow{v}}^2〉\left(r\right)=\int \frac{d^{3}k}{4{\pi }^3}\frac{\left(1-e^{ı\overrightarrow{k}·\overrightarrow{r}}\right)E(k,t)}{4\pi {\overrightarrow{k}}^2};\hfill \end{array}$$

(74)

$$\begin{array}{ccc}& & 〈\mathsf{\Delta }{\overrightarrow{v}}^2〉\left(r\right)=〈{\left(\overrightarrow{v}(\overrightarrow{r},t)-\overrightarrow{v}(0,t)\right)}^2〉\hfill \end{array}$$

(75)

There is a sharper test, namely the effective index ${\xi }_2(r,t)$ defined as a log-log derivative of this second moment

$$\begin{array}{ccc}& & {\xi }_2(r,t)=r{\partial }_{r}log〈\mathsf{\Delta }{v}^2〉\left(r\right)=\hfill \end{array}$$

$$\begin{array}{ccc}& & \frac{{\int }_0^{\infty }dk\left(-cos\left(kr\right)+\frac{sin\left(kr\right)}{kr}\right)E(k,t)}{{\int }_0^{\infty }dk\left(1-\frac{sin\left(kr\right)}{kr}\right)E(k,t)};\hfill \end{array}$$

(76a)

$$\begin{array}{ccc}& & {〈{\xi }_2(xL\left(t\right),t)〉}_t=f\left(x\right);\hfill \end{array}$$

(76b)

This universal function $f\left(x\right)$ is numerically well defined as the integration over the spectrum suppresses the noise related to lattice discretization unless the dimensionless coordinate x is too large.

In the K41 theory, this index is $\frac{2}{3}$, so one would expect at least a plateau of nearly K41 values in the inertial range of space scales. This index would have a higher constant value in multifractal models, around $0.7$. Figure 12 shows what we found instead for the DNS data [40,46].

This time, the theoretical curve (green line) goes outside the error bars, which means a lack of a global fit. Both curves are far from any constant value, so Kolmogorov and multifractal are out of the competition. However, the left parts of the theoretical and DNS curves, with $f>0.5$, are almost parallel in $logx$ scale. The larger values of x do not match so closely, and the error bars are bigger there, likely because of lattice artifacts in DNS. These left parts, however, can be moved on top of each other by properly choosing the length scale in our theory.

Figure 13. The plot of the effective index $f\left(x\right)$ in (Section 8) as a function of $logx$ (red dots with error bars) in the turbulent range. The K41 $\frac{2}{3}$ law (blue dashed line) is totally off the charts.

Figure 13. The plot of the effective index $f\left(x\right)$ in (Section 8) as a function of $logx$ (red dots with error bars) in the turbulent range. The K41 $\frac{2}{3}$ law (blue dashed line) is totally off the charts.

We tested this hypothesis by selecting the left side of this plot and adjusting the length scale in the theoretical line to get on top of the DNS data (shifting the green curve horizontally by a mean distance to the red curve). Here is the result of this selection/shifting (see Figure 13).

The theoretical curve (green line) is shifted left by some amount $\delta logx=s$ to minimize the mean square of the horizontal distance in the turbulent range. This shift corresponds to the adjustment of the length scale in our theory. The curves now match up to the standard deviation of the DNS. The time limit $T_0<t<T_1$ of decaying turbulence was chosen to ensure the diffusion law $L\left(t\right)\propto \sqrt{t}$. These DNS completely rule out the K41 scaling law $\frac{2}{3}$.

We recently encountered real experiments for compressible decaying turbulence in wind tunnels and atmosphere [23]. Our theory assumes incompressibility, so it doesn’t apply to these data in the air turbulence. Also, the magnitude of the errors in [23] is unclear.

Nevertheless, we still compared our prediction for the second moment with the air data from [23]. The shape of the experimental curve for $〈\mathsf{\Delta }{v}^2〉\left(r\right)$ (Figure 1 in [23] ) is similar to ours. It significantly deviates from the K41 prediction: instead of a straight line in the log-log scale, there is a curved line with the slope varying from 2 to 0 as r varies from zero to infinity, the same as our curves. There is no plateau at $\frac{2}{3}$ slope (Figure 1 in [23]; the slope linearly decreases with $logr$, similar to our decreasing slope in Figure 12.

However, the numerical values for the slope and curvature in [23] are quite different from those we derived above from the incompressible DNS [40,46], which fits our theory within the error bars. The origin of such a big discrepancy is unclear to us. Compressibility alone seems unlikely to explain such a large deviation from incompressible DNS.

It is difficult to estimate the slope of the measured data numerically [23] without increasing errors, which are already large enough in the data in [40,46]. Unlike numerical differentiation of the wind tunnel data, the Fourier integration of the DNS data suppresses the noise, which makes the DNS data for effective index ${\xi }_2\left(r\right)$ much more accurate.

We conclude that our theory passed its first test with flying colors, but a more detailed comparison with new large DNS or experiments is desirable.

9. The spectrum of scaling dimensions

In addition to numerically computing the energy spectrum and plotting effective critical indexes, we can go one step further in the mathematical analysis of the concept of the scaling laws in turbulence.

In a scale-invariant theory, the Mellin transform of the correlation function in coordinate or momentum space is a meromorphic function. The spectrum of critical indexes is given by the positions of poles of this function in a complex plane. The whole spectrum of critical indexes is real in the theory of critical phenomena, described by a Conformal Field Theory (CFT).

As we shall prove below, our theory is scale-invariant by this definition, but it is not a CFT. In particular, some critical indexes come in complex conjugate pairs, reflecting the dissipative nature of our turbulence theory.

Figure 14. Complex poles of the Mellin transform $h\left(p\right)$ of $H\left(\kappa \right)$ in (81) defining the critical indexes of the energy spectrum as a function of $\kappa =|\overrightarrow{k}|\sqrt{\tilde{\nu }t}$.

Figure 14. Complex poles of the Mellin transform $h\left(p\right)$ of $H\left(\kappa \right)$ in (81) defining the critical indexes of the energy spectrum as a function of $\kappa =|\overrightarrow{k}|\sqrt{\tilde{\nu }t}$.

9.1. The energy spectrum

The spectral density $H\left(\kappa \right)$ is in (A188). The Mellin transform is the following integral, assuming our function decreases at infinity.

$$\begin{array}{ccc}& & H\left(\kappa \right)={\int }_{-ϵ-ı\infty }^{-ϵ+ı\infty }\frac{dp}{2\pi ı}h\left(p\right){\kappa }^p;\hfill \end{array}$$

(77)

$$\begin{array}{ccc}& & h\left(p\right)={\int }_0^{\infty }\frac{d\kappa }{\kappa }{\kappa }^{-p}H\left(\kappa \right)\hfill \end{array}$$

(78)

We change the sign of p in conventional definition to better describe our functions.

The pure power law of decay would correspond to the Mellin transform $h\left(p\right)$ having a single pole in the left semi-plane. The position $p=-a$ of this pole becomes an index of the power law.

$$\begin{array}{ccc}& & H\left(\kappa \right)\propto {\kappa }^{-a}\hfill \end{array}$$

(79)

The next level of complexity would be a function, depending on an extra parameter n, such that the Mellin transform has a simple pole, moving with this parameter $p=a\left(n\right)$. In application to the moments of velocity difference, this parameter n is the degree of the moment $M_n=〈\mathsf{\Delta }{v}^n〉$. This pole will produce multifractal scaling laws

$$\begin{array}{c}\hfill M_n\sim r^{a\left(n\right)}\end{array}$$

(80)

Our energy spectrum has a more complex singularity structure in its Mellin transform

$$\begin{array}{ccc}\hfill & & h\left(p\right)=\frac{f\left(p\right)\zeta \left(p+\frac{15}{2}\right)\mathsf{\Gamma }(-p)}{(2p+7)(2p+17)\zeta \left(p+\frac{17}{2}\right)};\hfill \end{array}$$

(81)

$$\begin{array}{ccc}& & f\left(p\right)=20{\int }_{{\mathsf{\Delta }}_1}^{{\mathsf{\Delta }}_2}d\mathsf{\Delta }(1-\mathsf{\Delta })C^{p-1}(AC-Bp)\hfill \end{array}$$

(82)

where $A,B,C$ are some smooth positive functions of $\mathsf{\Delta }$ varying in finite limits (see Appendix K). Given these properties, it is simple to prove that the Taylor series of $f\left(p\right)$ at the origin converges as the expansion coefficients decrease as a factorial of the expansion order.

Figure 15. Complex poles of the Mellin transform $e\left(q\right)$ of $E\left(t\right)$ in (87) defining the spectrum of critical indexes of remaining energy as a function of time.

Figure 15. Complex poles of the Mellin transform $e\left(q\right)$ of $E\left(t\right)$ in (87) defining the spectrum of critical indexes of remaining energy as a function of time.

Figure 16. The spectrum of (complex) indexes of the power expansion of the velocity correlation function. The poles in the right semiplane determine the small-distance expansion of the correlation function $\overrightarrow{v}\left(0\right)·\overrightarrow{v}\left(r\right)\sim {\sum }_n{O}_n{\left|\overrightarrow{r}\right|}^{p_n}$.

Figure 16. The spectrum of (complex) indexes of the power expansion of the velocity correlation function. The poles in the right semiplane determine the small-distance expansion of the correlation function $\overrightarrow{v}\left(0\right)·\overrightarrow{v}\left(r\right)\sim {\sum }_n{O}_n{\left|\overrightarrow{r}\right|}^{p_n}$.

This convergence makes this function an entire function without any finite singularities. The values of $C\left(\mathsf{\Delta }\right)$ are bounded by two positive limits

$$\begin{array}{c}\hfill 0.0541984<C<0.0630755,\forall {\mathsf{\Delta }}_1<\mathsf{\Delta }<{\mathsf{\Delta }}_2\end{array}$$

(83)

Therefore, the entire function $f\left(p\right)$ decreases as $e^{-2.76342p}$ in the right semiplane and grows as $e^{-2.9151p}$ in the left semiplane. It oscillates along an imaginary axis, which is our integration path.

Theoretically, $f\left(p\right)$ could vanish at one or more positions of the poles of the remaining meromorphic function, eliminating these poles. We computed $f\left(p\right)$ at the lowest poles and ensured it was far from zero. However, canceling some higher poles by the root of $f\left(p\right)$ remains an open problem.

The singularities of the Mellin transform $h\left(p\right)$ for the energy spectrum $H\left(\kappa \right)$ are given by the following table of simple poles

(84)

Here $\pm t_n$ are imaginary parts of zeros of the $\zeta$ function on the critical line $z=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.+ıt$. About ${10}^{13}$ zeros are already known, though the Riemann hypothesis (no other complex zeros) still needs to be proven. These poles are shown in Figure 14.

Only poles with negative real parts contribute to the power expansion at $\kappa \to \infty$. The real negative poles yield decaying power terms, but the infinite series of complex poles at Riemann zeros adds the oscillations in log scale

$$\begin{array}{c}\hfill |W_n|{\kappa }^{-8}cos\left(t_{n}log\kappa +arg{W}_n\right)\end{array}$$

(85)

These slow oscillations are visible as regime change in the log-log plots of effective indexes in Figure 8, Figure 6.

As we already discussed above, the energy spectrum decays as $t^{-\raisebox{1ex}{$9$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}{k}^{-\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$. The wavelength decay at a fixed time is faster than K41 $k^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$. There is no theoretical reason (even at the level of heuristic) for K41 in decaying turbulence, as the dissipation $\mathcal{E}\left(t\right)$ is not a constant, so it cannot be used as a single scaling parameter.

9.2. The energy decay

The energy $E\left(t\right)$ is related to the same function

$$\begin{array}{ccc}& & E\left(t\right)={\int }_t^{\infty }dt\mathcal{E}\left(t\right)=\hfill \\ & & 4\pi \tilde{\nu }{\int }_t^{\infty }\frac{dT}{T^2}{\int }_{k_0\sqrt{\tilde{\nu }T}}^{\infty }H\left(x\right)x^{2}dx\hfill \end{array}$$

(86)

Substituting the Mellin transform for $H\left(x\right)$ and integrating twice, we get the Mellin transform for the energy

$$\begin{array}{ccc}& & E\left(t\right)={\int }_{-1-ϵ-ı\infty }^{-1-ϵ+ı\infty }\frac{dq}{2\pi ı}e\left(q\right){\left(k_0^2\tilde{\nu }t\right)}^q;\hfill \end{array}$$

(87a)

$$\begin{array}{ccc}& & e\left(q\right)=2\pi k_0^2{\tilde{\nu }}^2\frac{h(2q-1)}{q(q+1)}\hfill \end{array}$$

(87b)

The table of complex poles of this function is

(88)

These poles are shown in Figure 15. The leading pole is at $p=-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.$, corresponding to the asymptotic decay we compared to the grid turbulence decay data. Only poles with negative real parts contribute to the large t expansion.

9.3. The velocity correlation function in coordinate space

Let us transform the velocity correlation back to coordinate space from Fourier space

$$\begin{array}{ccc}& & 〈\mathsf{\Delta }{\overrightarrow{v}}^2〉\left(r\right)=\hfill \\ & & 2\int \frac{d^{3}k}{{\left(2\pi \right)}^3}\left(1-exp\left(ı\overrightarrow{k}·\left(\overrightarrow{r}\right)\right)\right)\frac{E(k,t)}{4\pi {\overrightarrow{k}}^2}=\hfill \\ & & \frac{{\tilde{\nu }}^2}{{\pi }^2\nu t}{\int }_{ϵ-ı\infty }^{ϵ+ı\infty }\frac{dp}{2\pi ı}{\left(\frac{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}{\sqrt{\tilde{\nu }t}}\right)}^p\hfill \end{array}$$

$$\begin{array}{ccc}& & cos\left(\frac{\pi p}{2}\right)\mathsf{\Gamma }(-p-1)h(-1-p);\hfill \end{array}$$

(89)

$$\begin{array}{ccc}& & 〈\mathsf{\Delta }{\overrightarrow{v}}^2〉\left(r\right)=〈{\left(\overrightarrow{v}(\overrightarrow{r},t)-\overrightarrow{v}(0,t)\right)}^2〉\hfill \end{array}$$

(90)

The Mellin transform simplifies to

$$\begin{array}{ccc}\hfill & & 〈\mathsf{\Delta }{\overrightarrow{v}}^2〉\left(r\right)=\hfill \end{array}$$

$$\begin{array}{ccc}& & \frac{{\tilde{\nu }}^2}{\nu t}{\int }_{ϵ-ı\infty }^{ϵ+ı\infty }\frac{dp}{2\pi ı}V\left(p\right){\left(\frac{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}{\sqrt{\tilde{\nu }t}}\right)}^p;\hfill \\ & & V\left(p\right)=\hfill \end{array}$$

(91)

$$\begin{array}{ccc}& & -\frac{f(-1-p)\zeta \left(\frac{13}{2}-p\right)csc\left(\frac{\pi p}{2}\right)}{16{\pi }^2(p+1)(2p-15)(2p-5)\zeta \left(\frac{15}{2}-p\right)}\hfill \end{array}$$

(92)

The poles of this function in the right semiplane represent the indexes of the power singularities of the velocity correlation function at coinciding points $(0,\overrightarrow{r})$. Should our theory be a CFT (which it is not), this spectrum would be related to the spectrum of anomalous dimensions ${\mathsf{\Delta }}_n$ in the OPE:

$$\mathrm{CFT}:\overrightarrow{v}\left(0\right)·\overrightarrow{v}\left(r\right)\sim \sum _n{O}_n{\left|r\right|}^{p_n};p_n=2{\mathsf{\Delta }}_v-{\mathsf{\Delta }}_n;$$

We have such an expansion with $\frac{|\overrightarrow{r}|}{\sqrt{\tilde{\nu }t}}$ in place of $\left|r\right|$ and a factor $\frac{{\tilde{\nu }}^2}{\nu t}$ in front. The spectrum of these scaling indexes $p_n$ (unrelated to a dilatation operator as far as we know) is given in the following table:

(93)

Figure 17. The universal function (91) as a function of $\rho =|\overrightarrow{r}|/\sqrt{\tilde{\nu }t}$. The turnover is caused by subleading terms in the power expansion starting with ${\rho }^2$. The next terms involve quantum oscillations manifesting as a turnover from power growth to power decay. Asymptotic at large $\rho$ is $\mathit{const}$ as it follows from the table (93) of the decay indexes.

Figure 17. The universal function (91) as a function of $\rho =|\overrightarrow{r}|/\sqrt{\tilde{\nu }t}$. The turnover is caused by subleading terms in the power expansion starting with ${\rho }^2$. The next terms involve quantum oscillations manifesting as a turnover from power growth to power decay. Asymptotic at large $\rho$ is $\mathit{const}$ as it follows from the table (93) of the decay indexes.

Figure 18. Oscillations of the effective index ${\xi }_2\left(r\right)$ at large ${log}_{10}r$. This is a theoretical curve corresponding to the zoom into a region of large separations, currently inaccessible by DNS with required accuracy.

Figure 18. Oscillations of the effective index ${\xi }_2\left(r\right)$ at large ${log}_{10}r$. This is a theoretical curve corresponding to the zoom into a region of large separations, currently inaccessible by DNS with required accuracy.

We have no CFT but a calculable spectrum of scaling dimensions. Unlike the CFT in three dimensions, this spectrum is complex.

Only poles with positive/(non-positive) real parts contribute to the power expansion at $|\overrightarrow{r}|\to 0/\infty$. The leading term at $r\to 0$ is $r^2$, which is calculable in general form from its definition after expanding the exponential and averaging over directions of $\overrightarrow{k}$

$$\begin{array}{ccc}& & 〈\mathsf{\Delta }{\overrightarrow{v}}^2〉\left(r\right)\to \frac{r^2{\mathcal{E}}_{\infty }\left(t\right)}{24{\pi }^3\nu }=\frac{r^2\tilde{\nu }}{20736\pi \zeta \left(3\right)\nu t^2}\hfill \end{array}$$

(94)

The next term is $r^{\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ as it follows from the table in (93). As mentioned about the energy spectrum, there is no K41 scaling index $p=\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$. This omission is not a contradiction, as K41 does not apply to decaying turbulence. Instead of pure scaling laws with single decay indexes, we found an infinite spectrum of scaling indexes, some of which come as complex conjugate pairs, which leads to quantum oscillations: see Figure 18 for the oscillation of the index ${\xi }_2$ at the latest stage. This region is inaccessible with modern scale of the DNS.

The theoretical curve in Figure 17 agrees with the Fourier-transformed data of [40,46] but deviates from [23]. The probable reason is the compressibility of the air in real experiments in [23].

The imaginary parts of these complex scaling dimensions coincide with those of the famous Riemann zeta zeros, establishing an intriguing relation between Turbulence and Number Theory.

10. Discussion

This section will try to reconcile traditional perspectives on turbulence phenomena, including enduring beliefs and myths, with our new theory.

10.1. Myth and Reality of Turbulent Scaling Laws

More than eighty years ago, Kolmogorov and Obukhov made a breakthrough in turbulence theory by establishing the relation (7) for the three-point correlation function of the velocity field in turbulent flow. This formula only selects the potential part of the triple velocity correlation function by taking two coincident points. When taking the curl, we get zero:

$$〈v_{\alpha }\left({\overrightarrow{r}}_0\right)v_{\beta }\left({\overrightarrow{r}}_0\right){\omega }_{\gamma }(\overrightarrow{r}+{\overrightarrow{r}}_0)〉=0$$

This relation indicates no constraints on the rotational part concerning triple vorticity correlations and sheds no light on the scale invariance of turbulence theory.

Moreover, the linear term in coordinates does not have any support in the Fourier spectrum: this is an example of the harmonic terms added to the Biot-Savart integral for the velocity field.

The K41 scaling law was introduced as a phenomenological model, not intended to replace the missing microscopic theory. It was based on the assumption that the local dissipation density does not fluctuate—a limitation its creators were aware of, prompting them to propose a log-normal distribution for this variable later on. However, even this modified model lacked a microscopic justification and failed to fully correspond with empirical observations.

Subsequent experiments and DNS [53,56,57] have invalidated the K41 scaling laws (including the log-normal model) over the past thirty years. Regarding decaying turbulence, the experimental data [46] have diverged even further from Kolmogorov scaling laws despite all attempts to stretch this data or discard the non-fitting region as "erosion." We highlight significant deviations—six orders of magnitude—from the $k^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$ scaling in Figure 4. There are also recent measurements [23] with significant deviations of the log-log derivative of the second velocity moment $〈\mathsf{\Delta }{v}^2〉\left(r\right)$ from the $\frac{2}{3}$ predicted by K41.

A broader assumption posited that power laws with anomalous dimensions might exist in the inertial range. The assumed analogy to critical phenomena led to the proposal of multifractal scaling laws [47], which, as a phenomenological model, successfully described observed deviations from the K41 laws in forced turbulence[53,56].

However, there are no theoretical grounds for conformal symmetry in turbulence; the ’current conservation’ conditions ${\partial }_{\alpha }{v}_{\alpha }=0,{\partial }_{\gamma }{\omega }_{\gamma }=0$ in the CFT would prescribe both velocity and vorticity dimensions of $d-1=2$, contradicting the fact that vorticity is a curl of velocity.

Moreover, the anomalous dimension would not explain decaying turbulence, as the log-log plots would remain straight lines, though the slopes would become irrational numbers. We must allow nonlinear correlation functions on a log-log scale, as indicated by the data in Figure 2 (top) in [23]. Here, energy spectra for various parameters converge into a universal curve on the decreasing part of the spectrum, which is curved on the log-log scale, indicating that a simple power law cannot describe it. Instead, it is a nontrivial universal function of $logk$, spanning several decades.

Both the DNS and experimental papers [23,46] noted significant deviations from scaling laws (whether K41 or multifractal). The conclusion of [46] was cautiously negative: "it is somewhat disappointing that the results are not more closely aligned with theoretical arguments." The most recent paper [23] made a stronger negative claim: "Our results point to a Reynolds number-independent logarithmic correction to the classical power law for decaying turbulence that calls for theoretical understanding."

Our recent paper [31] presented the theoretical argument for the breaking of scaling laws due to logarithmic divergences in a dilute gas of vortex filaments. In this approximation, there were logarithmic terms in the effective energy for the filament, leading to violations of scaling laws akin to asymptotic freedom in QCD. This approximation does not apply to decaying turbulence with a large density of vortex structures, but at least it identifies a dynamical mechanism for the deviations from the scaling laws.

In the present paper, we used raw data from the DNS [40,46] to compute the effective index of the velocity correlation by numerical Fourier transform of their energy spectrum (see Section 8). Our effective index is plotted in Figure 13. The K41 scaling law $\mathsf{\Delta }{v}^2\sim r^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$ is very far from reality, as it is clear from these plots. Our theory is much closer, and by fitting our arbitrary length scale, we obtained a very good fit in the turbulent range within experimental errors.

The microscopic theory developed here is not conformally invariant but retains a critical aspect of CFT. The Mellin transform of the vorticity field’s correlation function in coordinate space is a meromorphic function of the Mellin parameter p. This characteristic implies some underlying scale invariance with an infinite discrete spectrum of complex anomalous dimensions (93).

In this way, our theory extends the multifractal scaling laws by accommodating an infinite discrete spectrum of scaling dimensions. According to recent DNS results [54], our predictions align with observed slopes, unlike conventional scaling models such as those proposed by Kolmogorov and Saffman.

This result suggests that experiments and DNS in decaying turbulence should be conducted at larger scales and higher Reynolds numbers, fitting the data for logarithms of the spectrum and energy dissipation decay as nonlinear functions of the logarithm of the product of wavevector and the square root of time, as we (successfully) did in Section 8.

As part of this reevaluation, one should magnify and study the decaying part of the spectrum way beyond its middle part, roughly described by $-\frac{5}{3}$ law with logarithmic corrections. This decaying part, the "dissipative subrange," was discarded as an unfitting puzzle piece in conventional data fitting, but it fits well in our theory.

Heisenberg [17] and Chandrasekhar [8] proposed in the middle of the last century for the "dissipative subrange," the spectrum decay $k^{-7}$, based on a model equation by Heisenberg. At that time, there were no mathematical tools to solve the turbulence problem exactly, so the model equations like that one passed as theories. The fame of two Nobel laureates involved added weight to this model assumption, so it stays alive to this day.

K.R. Sreenivasan dispelled this die-hard myth in his paper [52].

Chandra’s initial enthusiasm for Heisenberg’s work was moderated when he learned from J. von Neumann, in a colloquium that Chandra gave at Princeton in the spring of 1949, that the $k^{-7}$ power law in the far-dissipation range did not have experimental support.

Our theory also contradicts the $k^{-7}$ law: no pole exists between $-13/2$ and $-8\pm ıt_n$ in our Mellin transform spectrum (84). Instead, we have nontrivial dynamics at this "dissipative subrange," not just in the "inertial range" between energy pumping and dissipation. The full plot of the effective index for the spectrum is shown in Figure 6, asymptotically approaching $-\frac{7}{2}$. As long as enough energy is left for a turbulent flow, our theory has a universal decaying spectrum spanning several decades and a strongly curved second moment $〈\mathsf{\Delta }{v}^2〉\left(r\right)$ with effective index (log derivative) shown in Figure 12. We call the corresponding range of scales "turbulent range", combining old inertial and dissipation ranges. Our theory perfectly matches DNS/experiment in the whole turbulent range without any dimensionless fitting parameters.

In conclusion, single-power scaling laws cannot describe the observed critical phenomena in decaying turbulence. Instead, compare these phenomena with the microscopic theory, which goes beyond empirical laws, replacing them with universal nonlinear functions for the energy spectrum, energy decay, and velocity correlation.

10.2. Stochastic solution of the Navier-Stokes equation and ergodic hypothesis

Richard Feynman wrote about turbulence in his Lectures in Physics [12]:

Nobody in physics has really been able to analyze it mathematically satisfactorily in spite of its importance to the sister sciences. It is the analysis of circulating or turbulent fluids.... What we really cannot do is deal with actual, wet water running through a pipe. That is the central problem which we ought to solve some day, and we have not.

These words were written over sixty years ago, but the problem remains unsolved.

We address this problem by seeking a stochastic solution to the unforced Navier-Stokes equation, covering a universal manifold over an infinite time. Our solution reveals power-like singularities in correlation functions, which emerge after averaging across this manifold in the statistical limit as its dimension approaches infinity.

We already have a partial answer to the question posed by Feynman about water flowing through a pipe: its local kinetic energy density decays with time (and also with distance from the grid at the entrance) as $t^{-\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$. A more detailed answer for the pressure as a function of the total amount of water pushed through the pipe would require some future investigation of our solution.

These singularities originate from the infinite time required to cover this manifold uniformly.

We identify this manifold (the Euler ensemble) by solving the loop equation—a subset of the Hopf functional equation for the generating functional of velocity field probabilities. Notably, none of the solutions within this manifold experiences finite-time blow-ups. Instead, we encounter singularities from the fixed trajectory of the loop equation, not from its finite-time solutions.

We adopt the most natural invariant measure from the perspective of number theory: each element of the Euler ensemble is weighted equally, an assumption we term the quantum ergodic hypothesis.

With this invariant measure, the Euler ensemble stands out because the loop functional is equal to the trace of the evolution operator in a quantum system—the Fermi particles on a ring interacting with a quantum field made of fractions of $\pi$. Every distinct state, including every distinct fraction, contributes equally to the quantum trace in this discrete system. For our purposes, this means treating every element in the Euler ensemble with equal weight.

Our quantum ergodic hypothesis thus stipulates an exact equivalence between the loop functional and the quantum trace of an evolution operator for the one-dimensional ring of Fermi particles. This quantum analogy has paved the way for an analytical solution in the turbulent limit. This limit corresponds to the quasiclassical limit of this Fermi system, where viscosity acts like Planck’s constant.

The quantum ergodic hypothesis results from a more general relation between fluid dynamics and quantum mechanics. The loop equation, in the general case, with finite viscosity and external stochastic forces in the Navier-Stokes equation, represents the Schrödinger equation in loop space [28,29].

The time evolution of the wave function, which is the loop functional, is given by the sum over "classical" histories, corresponding to this loop space Hamiltonian. Dirac and Feynman established that the weight of each history $exp\left(ıS/\hslash \right)$ for any quantum system with the Action S.

Comparing the Dirac-Feynman rule with the definition of the loop functional, we conclude that the velocity circulation ${\mathsf{\Gamma }}_C\left[v\right]={\oint }_C\overrightarrow{v}·d\overrightarrow{r}$ plays the role of the Action in the loop quantum mechanics, and viscosity plays the role of Planck’s constant. The sum goes over the classical solutions of the Navier-Stokes equation with various initial data, with equal weight for each solution.

We cannot describe all these solutions for the velocity field, but surprisingly, we can compute the weighted sum of all these solutions, i.e., the loop functional. Let us stress that this is not an asymptotic solution of the loop equation, with some terms neglected at large times. Our solution (12), (20), (24) exactly satisfies the Navier-Stokes equation at a finite time for the loop functional in the turbulent limit $N\to \infty ,\nu \to 0$.

The wave functional is not localized in the weak turbulence phase (small circulations compared to viscosity), so states are not quantized. This quantization occurs only in the strong turbulent phase (large circulations); the Euler ensemble or the Fermi ring describes it.

In the same way, as one-dimensional quantum mechanical motion in external potential becomes finite and quantized when potential well becomes deep enough, our loop functional at large time transforms from continuous distribution in loop space to the quantized finite motion characterized by the momentum loop $\overrightarrow{P}\left(\theta \right)$. The continuous quantum mechanical integral over phase space with equal weight $1/\left(2\pi \right)$ per DOF becomes the discrete sum over all distinct quantum levels with unit weight.

We hope our quantum ergodic hypothesis can be proven from the Navier-Stokes equation, starting with the quantum representation of the loop equation as a Schrödinger equation in loop space. If confirmed, this hidden quantum mechanics of classical turbulence may become a law of Nature rather than a computational method.

Like the classical ergodic hypothesis, this may take another hundred years. Theoretical physics does not wait for rigorous proof but rather explores the consequences of the conjectured theory and compares them with physical and numerical experiments.

10.3. The physical meaning of the loop equation and dimensional reduction

The long-term evolution of Newton’s dynamical system with many particles eventually covers the energy surface (microcanonical ensemble). The ergodic hypothesis, accepted in Physics but still not proven mathematically, states that this energy surface is covered uniformly. Turbulence theory aims to find a replacement for the microcanonical ensemble for the Navier-Stokes equation. This surface would also participate in the decay in the pure Navier-Stokes equation without artificial forcing.

In both cases, Newton and Navier-Stokes, the probability distribution must satisfy the Hopf equation, which follows from the dynamics without specifying the mechanism of the stochastization. Indeed, the Gibbs and microcanonical distributions in Newton’s dynamics satisfy the Hopf equation in a rather trivial way: it reduces to the conservation of the probability measure (Liouville theorem), which suggests the energy surface as the only additive integral of motion to use in the exponent of the fixed point distribution.

The loop technology has been thoroughly discussed in the last few decades in gauge theories, including QCD [2,20,21,26,39], where the loop equations were first derived [24,25].

In the case of decaying turbulence, the loop equations represent a closed subset of the Hopf equations, which is still sufficient to generate the statistics of vorticity. In this case, the exact solution we have found for the loop functional also follows from the integrals of motion, this time, the conservation laws in the loop space.

The loop space Hamiltonian we derived from the unforced Navier-Stokes equation does not have any potential terms (those with explicit dependence upon the shape of the loop). The Schrödinger equation with only kinetic energy in the Hamiltonian conserves the momentum. The corresponding wave function is a superposition of plane waves $exp\left(ı\overrightarrow{p}·\overrightarrow{x}\right)$. This superposition is the solution we have found, except the dot product $\overrightarrow{p}·\overrightarrow{x}$ becomes a symplectic form $\oint \overrightarrow{P}\left(\theta \right)·d\overrightarrow{C}\left(\theta \right)$ in the loop space.

Our momentum $\overrightarrow{P}(\theta ,t)$ is not an integral of motion, but simple scaling properties of the pure Navier-Stokes equation lead to the solution with $\overrightarrow{P}(\theta ,t)\propto \overrightarrow{F}\left(\theta \right)/\sqrt{t}$, with $\overrightarrow{F}\left(\theta \right)$ being the integral of motion at large time (i.e., a fixed point). The rest is a purely technical task: substituting this scaling solution into the Navier-Stokes equation and solving the resulting universal equation for a fixed point $\overrightarrow{F}\left(\theta \right)$.

This equation led us to the Fermi ring in the quasiclassical limit. The solution of the Fermi ring in this limit resulted in the energy spectrum and dissipation in a finite system found in Section 7.

10.4. Classical Flow and Quantum Mechanics

Our computations rely significantly on number theory, particularly Jordan’s multitotients, ${\phi }_l\left(q\right)$, which extend the Euler totient function [41]. What could number theory share with turbulent flow? The quantization of parameters in the fixed manifold of decaying turbulence originates from the quantum correspondence identified in the nineties [28]. The statistical distribution of a nonlinear classical Navier-Stokes (NS) PDE is related to the wave functional of a linear Schrödinger equation in loop space, as detailed in the previous section.

Is quantum mechanics at work in a water faucet with a grid filter? Yes and no.

This relationship is indirect: the loop functional, a Fourier transform of the classical probability distribution for circulation, equals the complex quantum amplitude of the loop space theory represented by a Fermi ring. Probability is real and positive, while its Fourier transform is complex, reflecting the irreversibility of the Navier-Stokes dynamics.

A probability distribution for circulation satisfies another loop equation [28,29], with all coefficients being real. The time evolution of this probability spans alternative histories, as it typically does in statistical mechanics, but the weights of each history remain real and positive.

On the other hand, the complex loop functional adheres to the quantum mechanical evolution equation, resulting in quantum interference of alternative histories. The quantum interference is quite significant here, with the dominant complex trajectory—the instanton—describing notable quantum effects such as exponential cancellation of contributions from alternative histories and penetration into classically forbidden regions within loop space.

The quantization mechanism of the parameters in the plane wave solution mirrors that of conventional quantum mechanics: the solution’s periodicity $\overrightarrow{P}(\xi +1)=\overrightarrow{P}\left(\xi \right)$.

From the conventional perspective, the fractal curve in complex momentum space, $\overrightarrow{P}\left(\theta \right)\in {\mathbb{C}}_3$, or Fermi particles on a circle, may seem unrelated to turbulent flow. One could ask how turbulence can be addressed without directly studying the velocity field.

The well-established duality phenomenon, known as ADS/CFT duality, equates the strong coupling phase of a conformal field theory to the weak coupling phase of quantum geometry in another dimension. This relationship is more than just a method for calculating correlation functions of a strongly fluctuating vector field; it reveals a second identity of the original theory.

The quantum Fermi ring’s particle density fluctuations disappear in the turbulent limit. In contrast, the original theory’s fluctuations are so intense that the vorticity field ceases to exist. This Fermi ring arguably reveals the true identity of decaying turbulence, as a classical function describes a smooth Fermion density—the instanton solution we identified.

Coming back to decaying turbulence in a water faucet, the probability distribution of velocity circulation in the water stream you wash your hands with is classical, of course; however, this classical distribution decays by a complex law based on quantum interference for its Fourier transform. This Fourier transform (loop functional) adds up from alternative histories with quantum mechanical complex weights.

Mathematical physics sometimes has dual representations for the same phenomena, such as the duality between particles and waves in quantum mechanics, or between matrix models [6,11,15] and Liouville theory [22,49] in 2D quantum gravity.

Additional complexities arise in gauge theory due to short-distance singularities involving the infinite fluctuating degrees of freedom in quantum field theory. Wilson loop functionals in coordinate space are singular in gauge field theory and cannot be multiplicatively renormalized.

There are no short-distance divergences in the Navier-Stokes equations and NS loop dynamics. The Euler equations represent a singular limit that, as argued, should be resolved through singular topological solitons regularized by the Burgers vortex.

In the dual theory of this paper, the singularities exist in the dual space ${\mathbb{C}}_3$. Anomalous dissipation is achieved through numerous finite discontinuities of the fractal curve $\overrightarrow{P}\left(\xi \right)\in {\mathbb{C}}_3$.

However, these singularities only occur in the inviscid limit, $\nu \propto 1/N^2\to 0$, representing Euler singularities like line vortices that are regularized by finite viscosity, just like our singularities.

10.5. Renormalizability of the Inviscid Limit of the Loop Equation

Let us examine the relationship $\nu N^2=\mathit{const}$ between diminishing viscosity and the increasing number N of discontinuities on the momentum loop $\overrightarrow{P}\left(\theta \right)$.

The Navier-Stokes (Navier-Stokes) equation is essentially an idealization of molecular dynamics, approximating nonlocal theory by a truncated expansion in powers of gradients.

In the case of laminar flow, this truncated expansion poses no issues. However, in our solution for the loop equation, the velocity field becomes singular in the local limit.

Mathematically, velocity and vorticity are not ordinary functions in ${\mathbb{R}}_3$ but stochastic variables, with $\mathsf{\Delta }\overrightarrow{v}\sim {\left(\mathsf{\Delta }\overrightarrow{r}\right)}^{\alpha }$.

Fractal Calculus [13] was introduced to generally describe such fields. Yet, this alone does not account for turbulence, particularly since the more general power laws with multifractal dimensions $〈{\left(\mathsf{\Delta }\overrightarrow{v}\right)}^n〉\sim {\left(\mathsf{\Delta }\overrightarrow{r}\right)}^{{\zeta }_n}$ cannot be explained in this manner.

In our theory, there is a universal function (91) for the velocity correlation function as a function of the separation $\overrightarrow{r}$. However, the $\overrightarrow{r}$ dependence is not defined by a single power: the Mellin transform of this function reveals infinitely many poles (93) in the complex plane. The absence of branch cuts supports the idea of scale invariance, but the analogy with CFT ends here.

Our solution characterizes a fractal vorticity field, yet this field does not conform to known fractal types. The fractal curve $\overrightarrow{P}\left(\theta \right)$ is not merely a random walk on a circle: dynamic restrictions, such as periodicity, lead to nontrivial critical behavior not describable by any finite set of fractal power laws.

We define the turbulent limit of the velocity field distribution by discretizing the loop equation (replacing the continuum loop with a polygon with an increasing number of vertices).

The relation between vanishing viscosity and an increasing number of degrees of freedom parallels the renormalization group (RG) relation between the bare coupling constant $g_0$ and the lattice spacing a in QCD: $a\propto g_0^{\alpha }exp\left(-\beta /g_0^2\right)$.

The naive local limit $a\to 0,g_0=\mathit{const}$ does not exist in QCD, but the RG limit effectively describes the strong interaction of hadrons.

Similarly, in our approach, renormalizability is present. The dissipation rate remains finite as $\tilde{\nu }=\nu N^2\to \mathit{const}$. Furthermore, the energy spectrum, expressed as a function of $k\sqrt{\tilde{\nu }t}$, remains finite in the local limit.

Thus, akin to renormalizable quantum field theory (QFT), a "dimensional transmutation" occurs where infinities are absorbed into the dimensional parameter $\tilde{\nu }$, defining the time scale.

The energy scale is related to the fluid’s physical viscosity $\nu$ in the denominator. Formally setting $\nu =0$ would tend the physical energy scale to infinity. The fictitious limit $\nu \to 0$ was only taken to compute the turbulent limit of the energy spectrum as a function of Reynolds number $\mathsf{\Gamma }/\nu$.

In the real world, this number is large because of the large circulation $\mathsf{\Gamma }$ compared to finite viscosity $\nu$. Still, we can compute these dimensionless functions by taking the inviscid limit, as described by equation (69) in Section 7.2.

In particular, the energy spectrum is proportional to $1/\nu$ times the function of the effective Reynolds number, time, and wave vector. We tend the Reynolds number to infinity in the turbulent limit but keep a finite physical viscosity factor in the energy spectrum.

10.6. Relation of Our Solution to Weak Turbulence

The solution of the loop equation with finite area derivative, satisfying the Bianchi constraint, belongs to the category of Stokes-type functionals [24], similar to the Wilson loop for gauge theory and fluid dynamics.

The Navier-Stokes (NS) Wilson loop represents a case of the Abelian loop functional, characterized by commuting components of the vector field $\overrightarrow{v}$. As extensively discussed in [24,25,29], any Stokes-type functional $\mathsf{\Psi }(\gamma ,C)$ that satisfies the boundary condition at a shrunk loop $\mathsf{\Psi }[\gamma ,0]=1$ and solves the loop equation can be iterated in the nonlinear term in the NS equations. This iteration is applicable in conditions of high viscosity.

The resulting expansion in inverse powers of viscosity (representing weak turbulence) coincides with the standard perturbation expansion of the NS equations for the velocity field, averaged over the initial data distribution.

We have shown in [28,29] (and also here, in Section 4.2) how the velocity distribution for random uniform vorticity in the fluid was reproduced by a singular momentum loop $\overrightarrow{P}\left(\theta \right)$.

The solution for $\overrightarrow{P}\left(\theta \right)$ at this specific fixed point of the loop equation was complex and had slowly decreasing Fourier coefficients, leading to a discontinuity $sign(\theta -{\theta }^{\prime })$ in a pair correlation function (A14). The corresponding Wilson loop is equated to the Stokes-type functional (A6).

Using this example, we demonstrated how a discontinuous random momentum loop describes the vorticity distribution in the stochastic NS flow. Here, the vorticity acts as a global random variable corresponding to a random uniform fluid rotation—a well-known exact solution of the NS equation.

This example corresponds to a special fixed point for the loop equation. Although not general enough to describe turbulent flow, it is an ideal mathematical model for loop technology. It illustrates how the momentum loop solution aggregates all terms of the $1/\nu$ expansion in the NS equation.

In our general solution, with the Euler ensemble, the summation of a divergent perturbation expansion occurs at an extreme level, leading to a universal fixed point in the limit of vanishing viscosity.

Given initial conditions, after a finite time, the solution will still depend on viscosity and the initial conditions. We expect no singularities for a smooth initial velocity field.

11. Remaining problems

• The loop functional for the circular loop is the simplest object in this theory. It can be computed using the methods developed in this paper, yielding even simpler results. In this case, the classical equation is trivial, so the computations reduce to the functional determinant and the resolvent. On the other hand, this is an observable quantity and could be measured in DNS. It would be an interesting problem to solve and compare with the DNS.
• The higher moments of circulation or velocity differences are calculable from this general WKB approximation for the path integral at $\nu \to 0,N\to \infty$. The nth moment of $\mathsf{\Delta }\overrightarrow{v}=\int d^{3}k\frac{ı\overrightarrow{k}\times {\overrightarrow{\omega }}_k}{{\overrightarrow{k}}^2}\left(exp\left(ı\overrightarrow{k}·\overrightarrow{r}\right)-exp\left(ı\overrightarrow{k}·{\overrightarrow{r}}^{\prime }\right)\right)$ reduces to the loop functional for the same backtracking "hairpin" traversed n times, with vorticity inserted n times at the ends. This computation will produce analytical results for the multifractal scaling laws for velocity moments.
• The spectrum of indexes for deviations from our fixed trajectory [30] can be evaluated to compute vorticity correlation functions in the Navier-Stokes equation with an infinitesimal random force.

12. Conclusion

We have established an exact duality between decaying classical turbulence in 3+1 dimensions and a solvable one-dimensional quantum theory of Fermi particles on a ring. In this framework, Fermi particles are confined within the vorticity field, where strong vorticity fluctuations correspond to weak fluctuations in the Fermion density. Elaborating on this theory, we present an analytical solution for decaying turbulence in quadrature.

Our theory replaces the old scaling laws (including multifractal versions) with universal functions, nonlinear in a log-log scale. This microscopic theory challenges the prevailing paradigms of the past eighty years, but it perfectly fits the available experimental data [9,10] and recent DNS data [40,46,51]. It offers a far more intriguing perspective, revealing unexpected connections between nonlinear classical physics and quantum mechanics.