Quantum solution of classical turbulence: Decaying energy spectrum

Migdal, Alexander

doi:10.1063/5.0228660

This paper presents a recent advancement that transforms the problem of decaying turbulence in the Navier–Stokes equations in 3 + 1 dimensions into a Number Theory challenge: finding the statistical limit of the Euler ensemble. We redefine this ensemble as a Markov chain, establishing its equivalence to the quantum statistical theory of N Fermions on a ring, interacting with an external field associated with random fractions of π. Analyzing this theory in the turbulent limit, where $N \to \infty$ and $ν \to 0$ ⁠, we discover the solution as a complex trajectory (instanton) that acts as a saddle point in the path integral over the density of these Fermions. By computing the contribution of this instanton to the vorticity correlation function, we obtain an analytic formula for the observable energy spectrum—a complete solution of decaying turbulence derived entirely from first principles without the need for approximations or fitted dimensionless parameters. Our analysis reveals the full spectrum of critical indices in the velocity correlation function in coordinate space, determined by the poles of the Mellin transform, which we prove to be a meromorphic function. Real and complex poles are identified, with the complex poles reflecting dissipation and uniquely determined by the famous complex zeros of the Riemann zeta function. Universal functions of the scaling variables supersede the traditional turbulent scaling laws (K41, Heisenberg, and multifractal). These functions for the energy spectrum, energy decay rate, and velocity correlation significantly deviate from power laws but closely match the results from grid turbulence experiments [Comte-Bellot and Corrsin, J. Fluid Mech. 48(2), 273–337 (1971); Comte-Bellot and Corrsin, J. Fluid Mech. 25(4), 657–682 (1966)] and recent direct numerical simulation data [Panickacheril John, Donzis, and Sreenivasan, Philos. Trans. A Math. Phys. Eng. Sci. 380(2218), 20210089 (2022)] within experimental error margins.

I. PROLOGUE

Turbulence appears as an overwhelmingly complex phenomenon. As depicted in Fig. 1 from a recent simulation,⁵⁰ 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.

FIG. 1.

View large Download slide

The three-dimensional vortices in quantum turbulent flow. Reproduced from Polanco, Müller, and Krstulovic, Nat. Commun. 12, 7090 (2021). Licensed under a Creative Commons Attribution (CC BY) license.

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 the Navier–Stokes (NS) 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 Ref. 47) 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 Fig. 2 for the historical outline), proposes a simpler variant of the Hopf equation—the loop equation—which suffices to define the statistics.

FIG. 2.

View large Download slide

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 (1D) 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 (t) = t \frac{E (t)}{E (t)}$ ⁠.

Experimental data from decaying grid turbulence (see Fig. 3) corroborate our prediction (Fig. 4) 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.

FIG. 3.

View large Download slide

The experimental data^9,10 for decaying turbulence behind the oscillating grid. Reproduced with permission from Comte-Bellot and Corrsin, J. Fluid Mech. 25(4), 657–682 (1966). Copyright 1966 Cambridge University Press. The two lines correspond to log-log plots of $1 / v_{⊥}^{2}, 1 / v_{∥}^{2}$ for the flow behind the grid. The total inverse kinetic energy i $K = 2 / (v_{⊥}^{2} + v_{∥}^{2})$ ⁠, would have 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,⁵⁶ as the most reliable experimental data. Other measurements and simulations, at different conditions reviewed in Ref. 48, provide mismatching data influenced by initial and boundary conditions. The only large-scale DNS^42,48 covering three time decades also matches our predictions.

FIG. 4.

Green curve is the inverse energy 1 / E ( t ) as a function of time. It slowly approaches from above its asymptotic law 1 / E ( t ) ∝ t 5 / 4 shown in red.

View large Download slide

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

Recent direct numerical simulation (DNS) and experiments^24,48 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 a way to reduce experimental errors by computing an effective index for the second moment of velocity difference using the 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. Indeed, new experimental data were published recently for decaying turbulence.¹⁵ Their measurements of time decay of the turbulent kinetic energy obey our asymptotic law $t^{- 5 / 4}$ ⁠, as shown in Fig. 4 of their paper. The large tank data relevant to our isotropic homogeneous turbulence are shown on the inlet. The authors of this paper are trying to explain their small tank data in terms of the Saffman–Kolmogorov model, dismissing the large tank asymptotic decay as “viscous effects.” We disagree with such dismissal, as our turbulence theory explains these data. The only problem with these data are that the Reynolds numbers are relatively small. We need confirmation with higher Reynolds numbers.

II. DEFINITIONS AND NOTATIONS

DNS = direct numerical simulation; OPE = operator product expansion; Wilson loop = loop average = average of phase factor of circulation as a functional of the shape of the loop. ν is the viscosity, $\vec{r}$ is the coordinate, $\vec{k}$ is the momentum (or wavevector), $ξ, η, ρ$ are dimensionless integration variables, α, β, γ, and other Greek indices used for vector components with implied summation for repeated Greek indexes. We use the units where the constant fluid density ρ = 1. The 3D vectors and the dot products are denoted like this: $\vec{v} (\vec{r}, t), \vec{k} \cdot \vec{r}$ ⁠. We also use the Einstein tensor notation with summation over repeated Greek indexes $v_{α}, r_{β}, ω_{α β} = e_{α β γ} \partial_{β} v_{γ}$ ⁠, 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 the 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 $\vec{C} (ξ)$ 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 $\vec{C}$ can be covered several times with ${\vec{C}}_{n} (ξ) = C (n ξ)$ ⁠. The momentum loops $\vec{P} (ξ, t)$ depend on time, whereas the spacial loops are static. The momentum loops can have discontinuities $Δ \vec{P} (θ, t) = \vec{P} (θ + 0, t) - \vec{P} (θ - 0, t)$ ⁠. The momentum loops are independent of the space loops, being the vector functions of $ξ, t$ ⁠. The circulation $Γ = \oint_{C} d \vec{r} \cdot \vec{v} (\vec{r}, t) = \int_{0}^{1} d ξ \vec{C}' (ξ) \cdot \vec{v} (\vec{C} (ξ), t)$ ⁠. The whole theory, starting with circulation, is parametric invariant $ξ \Rightarrow f (ξ), f' (ξ) > 0$ ⁠. The operators are denoted like $\hat{a}, {\hat{a}}^{†}, \hat{X}$ ⁠. We use three types of symbols for differentials. d is an ordinary differential of variable like $d ξ$ or vector variable like $d \vec{C} (θ) = d θ \vec{C}' (θ)$ ⁠. D symbolizes a measure for a path integral $D α$ ⁠; 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, δ is a functional variation, and $\frac{δ}{δ}$ is a functional derivative like in $\frac{δ}{δ ϕ (ξ)} Λ_{N} [ϕ] = \frac{δ Λ_{N} [ϕ]}{δ ϕ (ξ)}$ ⁠, with the L₂ norm in functional space. We also use notations $\partial_{t} = \frac{\partial}{\partial t}, \partial_{β} = \frac{\partial}{\partial r_{β}}$ for the time derivative and for components of the spatial derivatives.

III. SUMMARY

The main results reported in this paper are as follows:

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 π.
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{ν} = ν N^{2} \to const$ ⁠.
We derived an analytic formula for energy spectrum and dissipation in finite system (58a), (58c), and (J9) and investigated it in Appendix K.
The energy spectrum decays asymptotically as ${(\tilde{ν})}^{3 / 2} t^{- 1 / 2} κ^{- 7 / 2}$ ⁠, where $κ = k \sqrt{\tilde{ν} t}$ ⁠.
The turbulent kinetic energy decays as $E (t) \propto t^{- 5 / 4}$ ⁠.
Both effective indexes $n (t) = - \frac{d log E}{d log t}, μ (κ) = \frac{\partial log E (k, t)}{\partial log k}$ are nontrivial functions of the logarithm of scale and time, approaching $n (\infty) = 5 / 4, μ (\infty) = - 7 / 2$ (see Figs. 5 and 6).
The 1966 experimental values,^9,10 Fig. 3 of $n (\infty) \approx 1.25 \pm 0.02$ ⁠, agree with our prediction.
The results of DNS for the energy spectrum in decaying turbulence (Fig. 5 of the review paper⁴⁸) show the energy spectrum decaying faster than $k^{- 5 / 3}$ ⁠. The shape of a log-parabolic curve in DNS matches our prediction (Fig. 7).
In Sec. VIII, we verify our prediction for the second velocity moment $⟨ Δ v^{2} ⟩ (r)$ by numerical Fourier transform of the raw spectrum data from Refs. 42 and 48. 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 (Figs. 8 and 9).
We computed the spectrum of singularity indexes p_n in the product $\vec{v} (0, t) \cdot \vec{v} (r, t) \sim \sum C_{n} (t) 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 ζ function (88), (92), and (97).

FIG. 5.

The theoretical plot of effective spectral index μ ( T ) = k ∂ k log E ( k , t ) as a function of T = k 2 ν ̃ t . It starts at 0 and goes down to the limit − 7 / 2 . The DNS plot (Fig. 8) starts higher, at μ = 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 are required to verify this prediction of our theory.

View large Download slide

The theoretical plot of effective spectral index $μ (T) = k \partial_{k} log E (k, t)$ as a function of $T = k^{2} \tilde{ν} t$ ⁠. It starts at 0 and goes down to the limit $- 7 / 2$ ⁠. The DNS plot (Fig. 8) starts higher, at μ = 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 are required to verify this prediction of our theory.

FIG. 6.

The effective index n ( t ) = − t E ( t ) E ( t ) compared with asymptotic value n ( ∞ ) = 5 / 4 .

View large Download slide

The effective index $n (t) = - \frac{t E (t)}{E (t)}$ compared with asymptotic value $n (\infty) = 5 / 4$ ⁠.

FIG. 7.

View large Download slide

Log-log plot of the universal function $H (κ) κ^{5 / 3}$ ⁠. It starts growing, reaches the maximum, then turns down, asymptotically decaying as $κ^{- 11 / 6}$ ⁠. This qualitatively matches the DNS plot at Fig. 20. A larger DNS data at higher Reynolds numbers would be required to verify our theory with more precision.

FIG. 8.

The DNS plot of spectral index μ ( t ) = k ∂ k log E ( k , t ) as a function of log t at fixed values of k. Reproduced from Panickacheril John, Donzis, and Sreenivasan, Philos. Trans. R. Soc. A. 380(2218), 20210089 (2022), licensed under a Creative Commons Attribution (CC BY) license. The curve goes from μ = 2 to μ ≈ − 1 , which partly overlaps with our theoretical curve in Fig. 5. Larger dataset with higher Reynolds numbers is needed for a quantitative match with our theory.

View large Download slide

The DNS plot of spectral index $μ (t) = k \partial_{k} log E (k, t)$ as a function of $log t$ at fixed values of k. Reproduced from Panickacheril John, Donzis, and Sreenivasan, Philos. Trans. R. Soc. A. 380(2218), 20210089 (2022), licensed under a Creative Commons Attribution (CC BY) license. The curve goes from μ = 2 to $μ \approx - 1$ ⁠, which partly overlaps with our theoretical curve in Fig. 5. Larger dataset with higher Reynolds numbers is needed for a quantitative match with our theory.

FIG. 9.

The plot of the effective index f(x) in (80) as a function of log x (red dots with error bars) in the turbulent range. The K41 2 3 law (blue dashed line) is totally off the charts.

View large Download slide

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

IV. 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 Appendixes, leaving only the general ideas, concepts, and results in the main body of the paper. Furthermore, all the heavy computations, requiring tedious analytical transformations and numerical/analytical integrations, are performed in Mathematica^® notebooks. The resulting formulas and plots are then quoted in the Appendixes. Thus, we have four levels of hierarchy in this paper:

Summary
Sections of the main text
Appendixes
Mathematica^® notebooks

This transformation provides distinct pathways for understanding and applying our study's results based on the reader's background and interest.

We illustrate this hierarchical structure in this online document.³⁷ Browsing this document would give a reader a big picture with an interactive ability to zoom into various topics.

Portions of the upcoming discussion are borrowed from our first paper³¹ 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.

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.
Introduction for mathematicians. This section summarizes the mathematical framework behind the loop equation²⁹ and its solution³¹ 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.
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 Sec. VII. Here, they will find final formulas (58a) and (58c) and accompanying Mathematica^® code^33,38 that facilitate the computation of both the energy spectrum and dissipation rates. These formulas are compared with real and numerical experiments in Sec. VIII.
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 Sec. X, 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 Secs. IV A and IV B, we explore these themes in depth, aiming to provide clarity and actionable insights for all readers, regardless of their expertise or interest.

A. 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 ν in the forced Navier–Stokes equation. Various phenomenological models have also been aligned with experimental observations, as discussed in the recent review.⁴⁸

However, the comprehensive turbulence theory requires solving the Navier–Stokes equations in the strong coupling limit $ν \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 Ref. 48, suggest a decay of the dissipation rate $E \sim t^{- n}$ with $n \approx 1.2$ or 1.4 depending on the initial conditions (finite total momentum or zero total momentum but finite total angular momentum, see Ref. 48). Thus, two universality classes of decaying turbulence have been identified.

It remains unclear which, if any, of the data in Refs. 42 and 48 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,³¹ 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 toward 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

\vec{v}

and averaging over an ensemble of stochastic solutions is

\begin{matrix} \partial_{t} ⟨ E ⟩ + \int_{V} d^{3} r ⟨ ν {\vec{ω}}^{2} ⟩ \\ = - \int_{V} d^{3} r \partial_{β} ⟨ v_{β} (p + 1 / 2 v_{α}^{2}) + ν v_{α} (\partial_{β} v_{α} - \partial_{α} v_{β}) ⟩ . \end{matrix}

(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:

\partial_{t} ⟨ E ⟩ + E_{V} = - \int_{\partial V} d \vec{σ} \cdot ⟨ \vec{v} (p + 1 / 2 v_{α}^{2}) + ν \vec{ω} \times \vec{v} ⟩ .

(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 $\vec{ω} \times \vec{v}$ term disappearing, as no vorticity persists at infinity.

Additionally, the velocity dictated by the Biot–Savart law diminishes as

{| \vec{r} |}^{- 3}

at infinity, so only the

\vec{v} p

term remains significant:

\partial_{t} ⟨ E ⟩ + ⟨ E_{V} ⟩ \to - \int_{\partial V} d \vec{σ} \cdot ⟨ \vec{v} p ⟩ .

(3)

This representation of energy flow will remain finite even as the sphere expands if the pressure scales as

p \to - \vec{F} \cdot \vec{r}

⁠, where

\vec{F}

is the local force at any given point on a large sphere:

\partial_{t} ⟨ E ⟩ + ⟨ E ⟩ = f_{α} lim_{R \to \infty} R^{3} \int_{S_{2}} n_{α} n_{β} ⟨ v_{β} (R \vec{n}) ⟩ .

(4)

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

\begin{matrix} \partial_{t} ⟨ E ⟩ + ⟨ E ⟩ = - \int_{V} d^{3} r ⟨ v_{β} \partial_{β} p + v_{α} v_{β} \partial_{β} v_{α} ⟩ \\ = - \int_{V} d^{3} r ⟨ v_{α} v_{β} \partial_{β} v_{α} ⟩ - \int_{\partial V} d \vec{σ} \cdot ⟨ \vec{v} p ⟩ . \end{matrix}

(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

\vec{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:

\partial_{t} ⟨ E ⟩ + ⟨ v_{α} v_{β} \partial_{β} v_{α} ⟩ = - ⟨ E ⟩ / V .

(6)

This relation, alongside spatial symmetry properties in

ℝ_{d}

⁠, leads to the Kolmogorov three-point correlation in a steady state

\partial_{t} ⟨ E ⟩ = 0

⁠:

\begin{matrix} ⟨ v_{α} ({\vec{r}}_{0}) v_{β} ({\vec{r}}_{0}) v_{γ} (\vec{r} + {\vec{r}}_{0}) ⟩ \\ = \frac{E}{(d - 1) (d + 2) V} (δ_{α γ} r_{β} + δ_{β γ} r_{α} - \frac{2}{d} δ_{α β} r_{γ}) . \end{matrix}

(7)

In the conventional approach to the turbulence problem, periodic Gaussian random forces

\vec{F} (\vec{r}, t)

are added to the Navier–Stokes equations:

\begin{matrix} ⟨ E_{V} ⟩ = - \int_{V} d^{3} r ⟨ v_{β} \partial_{β} p - v_{β} f_{β} (\vec{r}, t) + v_{α} v_{β} \partial_{β} v_{α} ⟩ \\ = \int_{V} d^{3} r ⟨ v_{β} f_{β} (\vec{r}, t) ⟩ . \end{matrix}

(8)

As the force becomes uniformly distributed across space, we derive another definition with $E = \vec{F} \cdot \vec{P}$ ⁠, where $\vec{P} = \int_{V} d^{3} r \vec{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.³² The dual picture from our recent theory³¹ represents these by random gaps in the momentum curve $\vec{P} (θ)$ ⁠, as we shall discuss in Sec. IV B.

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

\vec{v} (r) = - \vec{\nabla} \times \int d^{3} r' \frac{\vec{ω} (r')}{4 π | r - r' |} .

(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 \vec{v} (0) \cdot \vec{v} (\vec{r}) ⟩$ ⁠, 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² (or Saffman) class.

On a large sphere

\partial V

with radius

R \to \infty

⁠,

\begin{matrix} {lim}_{R \to \infty}^{} R^{3} \vec{v} (R \vec{n}) \propto \frac{1}{4 π} \sum_{blobs} \int_{blob} d^{3} r' \vec{ω} (r') \times (\vec{r}' - ⟨ \vec{r} ⟩_{blob}) . \end{matrix}

(10)

Here

⟨ \vec{r} ⟩_{blob}

is the geometric center of each blob. Substituting this into the identity , 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 \vec{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² 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³¹ of the loop equations^29,30 describe the internal stochastization of the decaying turbulence by a dual discrete system.

B. 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.³⁰

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 Fig. 2). 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 past 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 past 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.^27,28 Recently, some numerical methods were found to solve loop equations beyond perturbation theory.^2,21,22 The loop dynamics was extended to quantum gravity, where it was used to study nonperturbative phenomena.^4,52

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,19,20,44 with Reynolds numbers of tens of thousands revealed and quantified violations of the K41 scaling laws. These numerical experiments are in agreement with the so-called multifractal scaling laws.⁴⁹

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,²⁰ which triggered the further development of the geometric theory of turbulence.^3,30,44 In particular, the Area law was justified for flat and quadratic minimal surfaces, and an exact scaling law in the confinement region $Γ \propto \sqrt{Area}$ was derived.³⁰ The area law was verified with better precision in Ref. 19.

In the previous paper,³¹ we have found a family of exact solutions of the loop equation for decaying turbulence.^29,30 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 as

Ψ (γ, C) = {⟨ exp (\frac{ı γ}{ν} \oint d \vec{C} (θ) \cdot \vec{v} (\vec{C} (θ))) ⟩}_{init},

(11)

Ψ (γ, C) \Rightarrow {⟨ exp (\frac{ı γ}{ν} \oint d \vec{C} (θ) \cdot \vec{P} (θ)) ⟩}_{E} .

(12)

In the first equation (the definition), the averaging

⟨ \dots ⟩

goes over initial data for the solutions of the Navier–Stokes equation for velocity field

\vec{v} (\vec{r})

⁠. In the second one (the solution), the averaging goes over the space of solutions

\vec{P} (θ)

of the loop equation.³¹ We choose in this paper the parametrization of the loop with

ξ = \frac{θ}{2 π}

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

\vec{P} (θ)

follows from the Navier–Stokes equation for

\vec{v}

\partial_{t} v_{α} = ν \partial_{β} ω_{β α} - v_{β} ω_{β α} - \partial_{α} (p + \frac{v_{β}^{2}}{2}),

(13)

\partial_{α} v_{α} = 0,

(14)

ω_{β α} = \partial_{β} v_{α} - \partial_{α} v_{β} .

(15)

After some transformations, replacing velocity and vorticity with the functional derivatives of the loop functional, we found the following momentum loop equations:^30,31

\begin{matrix} ν \partial_{t} \vec{P} = - γ^{2} {(Δ \vec{P})}^{2} \vec{P} + Δ \vec{P} (γ^{2} \vec{P} \cdot Δ \vec{P} + ı γ (\frac{{(\vec{P} \cdot Δ \vec{P})}^{2}}{Δ {\vec{P}}^{2}} - {\vec{P}}^{2})), \end{matrix}

(16)

\vec{P} (θ) \equiv \frac{\vec{P} (θ^{+}) + \vec{P} (θ^{-})}{2},

(17)

Δ \vec{P} (θ) \equiv \vec{P} (θ^{+}) - \vec{P} (θ^{-}) .

(18)

The momentum loop has a discontinuity $Δ \vec{P} (θ)$ at every parameter $0 < θ \leq 1$ ⁠, making it a fractal curve in complex space $ℂ_{d}$ ⁠. The details can be found in Refs. 30 and 31. We will skip the arguments $t, θ$ 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

\vec{P} (θ)

at t = 0, providing the initial velocity field distribution. Formally, the initial distribution

W_{0} [P]

of the momentum field

\vec{P} (θ)

is given by inverse functional Fourier transform as

\begin{matrix} W_{0} [P] = \int D C δ^{3} (\vec{C} [0]) Ψ {(γ, C)}_{t = 0} exp (- \frac{ı γ}{ν} \int d \vec{C} (θ) \cdot \vec{P} (θ)) . \end{matrix}

(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 $\vec{P} (θ)$ in (A8) and (A14) has a nontrivial Gaussian distribution with discontinuity $Δ P (θ)$ 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 $\vec{P} (θ, t)$ with some universal probability distribution related to the decaying turbulence statistics.

The following transformation reveals the hidden scaling invariance of decaying turbulence:

\vec{P} = \sqrt{\frac{ν}{2 (t + t_{0})}} \frac{\vec{F}}{γ} .

(20)

The new vector function

\vec{F}

satisfies an equation

\begin{matrix} 2 \partial_{τ} \vec{F} = (1 - {(Δ \vec{F})}^{2}) \vec{F} + Δ \vec{F} (γ^{2} \vec{F} \cdot Δ \vec{F} + ı γ (\frac{{(\vec{F} \cdot Δ \vec{F})}^{2}}{Δ {\vec{F}}^{2}} - {\vec{F}}^{2})), \end{matrix}

(21)

τ = log (t + t_{0}) .

(22)

This equation is invariant under translations of the new variable

τ = log (t + t_{0})

⁠, corresponding to the rescaling/translation of the original time given as

t \Rightarrow λ t + (λ - 1) t_{0} .

(23)

There are two consequences of this invariance:

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

Both of these properties were used in Ref. 31: 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 λ_i of decay $t^{- λ_{i}}$ of the small deviations from this fixed point. In this paper, we only consider the fixed point, leaving the exciting problem of the spectrum of anomalous dimensions for future research.

C. The big and small Euler ensembles

Let us remember the basic properties of the fixed point for

\vec{F}

in Ref. 31. It is defined as a limit

N \to \infty

of the polygon

{\vec{F}}_{0} \dots {\vec{F}}_{N} = {\vec{F}}_{0}

with the following vertices:

{\vec{F}}_{k} = \frac{{cos (α_{k}), sin (α_{k}), i cos (\frac{β}{2})}}{2 sin (\frac{β}{2})},

(24)

θ_{k} = \frac{k}{N}, β = \frac{2 π p}{q}, N \to \infty,

(25)

α_{k + 1} = α_{k} + σ_{k} β; σ_{k} = \pm 1, β \sum σ_{k} = 2 π p r .

(26)

The parameters

\hat{Ω}, p, q, r, σ_{0}, \dots, σ_{N} = σ_{0}

are random, making this solution for

\vec{F} (θ)

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 with the discrete version of discontinuity and principal value:

Δ \vec{F} \equiv {\vec{F}}_{k + 1} - {\vec{F}}_{k},

(27)

\vec{F} \equiv \frac{{\vec{F}}_{k + 1} + {\vec{F}}_{k}}{2} .

(28)

Both terms of the right side vanish; the term proportional to

Δ \vec{F}

and the term proportional to

\vec{F}

⁠. Otherwise, we would have

\vec{F} ∥ Δ \vec{F}

⁠, leading to zero vorticity.³¹ The ensemble of all the different solutions is called the big Euler ensemble. The integer numbers

σ_{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):

\oint d \vec{C} (θ) \cdot \vec{P} (θ) = - \oint d \vec{P} (θ) \cdot \vec{C} (θ),

(29)

\sum_{k} Δ C_{k} {\vec{P}}_{k} = - \sum_{k} Δ {\vec{P}}_{k} \cdot {\vec{C}}_{k} .

(30)

A remarkable property of this solution $\vec{P} (θ, 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 $\vec{P} (θ, t)$ does not depend on θ and thus drops from the integral $\oint d \vec{P} (θ, t) \cdot \vec{C} (θ)$ ⁠.

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

We could not find a global parametrization of such a solution.⁴⁶ Instead, we generated it numerically by taking a planar closed polygon and evolving its vertices ${\vec{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 $Im {\vec{F}}_{k} = 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.¹ 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,16 providing the same correlation functions as continuum gravity (Liouville theory²³).

The fractions

\frac{p}{q}

with fixed denominator are counted by Euler totient function

φ (q)

⁠, which is calculated as¹⁷

φ (q) = \sum_{\begin{matrix} p = 1 \\ (p, q) \end{matrix}}^{q - 1} 1 = q \prod_{p | q} (1 - \frac{1}{p}) .

(31)

For example,

φ (16) = 8

and

φ (17) = 16

⁠.

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

V. THE FERMI RING AND ITS CONTINUUM LIMIT

In this paper, we perform this averaging over σ 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 (B14) 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 $ν_{k} = 1$ and Dirac holes $ν_{k} = 0$ tend to some classical function $\sum_{l < k} (2 ν_{l} - 1) \to α (ξ)$ of “position” $ξ = \frac{k}{N}$ ⁠, leading to the exact solution.

Specifically, as we prove in Appendix C, the loop functional in the continuum limit $N \to \infty$ reduces to a quantum mechanical path integral (C29) over the Fermion density $α (ξ)$ ⁠. The effective Action for this path integral is given by circulation $ı Γ$ expressed in terms of this density plus a quadratic functional corresponding to Brownian measure (C26) for $α (ξ)$ ⁠.

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 $α (ξ)$ (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.³¹

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

ν \to \frac{\tilde{ν}}{N^{2}} \to 0 .

(32)

The parameter

\tilde{ν}

remains a free parameter of our theory, playing the role of turbulent viscosity. In particular, there is an anomalous energy dissipation in this limit

E (t) = \frac{\tilde{ν}}{t^{2}} F (k_{0}^{2} \tilde{ν} t),

(33)

F (0) = \frac{π^{2}}{864 ζ (3)} .

(34)

Here k₀ is the lower cutoff of the energy spectrum (to be discussed below). The decaying spectrum at small wavevectors

k < k_{0}

is related to the energy pumping at the initial moment t = t₀ and is time-independent.

The dimensionless parameter $N \to \infty$ plays the role of the Reynolds number. We derive the turbulent limit $N = 2 m \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.

VI. INSTANTON IN THE PATH INTEGRAL

The classical equation for our path integral reads [with

Ω \in O (3)

being a random rotation matrix]

α ″ = - ı κ (C_{Ω}^{'} exp (ı α) + {(C_{Ω}^{'})}^{⋆} exp (- ı α)),

(35)

κ = \frac{1}{2 π y \sqrt{X} \sqrt{2 ν t}},

(36)

C_{Ω} (θ) = \vec{C} (θ) \cdot \hat{Ω} \cdot {ı, 1, 0} .

(37)

The parameter κ is distributed according to the distributions and 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 (ı α)

⁠:

z ″ = \frac{{(z^{'})}^{2}}{z} + κ (C_{Ω}^{'} z^{2} + {(C_{Ω}^{'})}^{⋆}),

(38)

z (0) = z (1) = 1 .

(39)

This equation cannot be analytically solved for arbitrary periodic function

C'_{Ω} (ξ)

⁠.

The weak and strong coupling expansions by κ are straightforward. At small κ,

\begin{matrix} z (ξ) \to 1 + 2 κ (- A ξ + \int_{0}^{ξ} Re C_{Ω} (ξ') d ξ') + O (κ^{2}), \end{matrix}

(40)

A = \int_{0}^{1} Re C_{Ω} (ξ') d ξ' .

(41)

At large κ,

z (ξ) \to ı exp (- ı arg C_{Ω}^{'} (ξ)) = ı \frac{| C_{Ω}^{'} (ξ) |}{C_{Ω}^{'} (ξ)} .

(42)

This solution is valid at intermediate ξ, not too close to the boundaries

ξ = (0, 1)

⁠. In the region near the boundaries

ξ (1 - ξ) ≪ \frac{1}{\sqrt{κ}}

⁠, the following asymptotic agrees with the classical equation:

z \to 1 \pm ı ξ \sqrt{2 κ Re C'_{Ω} (0)} + O (ξ^{2}),

(43)

z \to 1 \pm ı (1 - ξ) \sqrt{2 κ Re C'_{Ω} (1)} + O ({(1 - ξ)}^{2}) .

(44)

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

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

The classical limit of the circulation in exponential of

\int_{0}^{1} d ξ Im (C_{Ω}^{'} (ξ) exp (ı α (ξ))) \to \int_{0}^{1} d ξ | C_{Ω}^{'} (ξ) |

(45)

becomes a positive definite function of the rotation matrix Ω. At large κ, the leading contribution will come from the rotation matrix minimizing this functional.

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

α (ξ) = \frac{π}{2} - arg C_{Ω}^{'} (ξ) .

(46)

This density does not fluctuate in a turbulent limit, except near the end points

ξ \to 0, ξ \to 1

⁠. In the vicinity of the end points, there is a different asymptotic solution for

α \to (z - 1) / ı

⁠.

Computing the Wilson loop for a specific loop, say, the circle, is an interesting problem, but there is a simpler quantity. In Sec. VII, 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.^24,42,48

This is the vorticity correlation function,³⁰ corresponding to the loop C backtracking between two points in space ${\vec{r}}_{1} = 0, {\vec{r}}_{2} = \vec{r}$ ⁠, with the vorticity operators are inserted at these two points (see Ref. 31 for details and the justification). The Fourier transform of this function describes the decaying energy spectrum, also measured experimentally.

In Appendixes F–H, we express this correlation function as a particular case of the second variation of the loop functional. Then, in Appendixes I and 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{ν} = const$ ⁠, there are no higher order corrections to this WKB approximation of the path integral.

VII. THE DECAYING ENERGY IN FINITE SYSTEM

The vorticity correlation in Fourier space doubles as an energy spectrum

E (k, t) = 4 π {\vec{k}}^{2} ⟨ \vec{v} \cdot {\vec{v}}_{k} ⟩ = 4 π ⟨ \vec{ω} \cdot {\vec{ω}}_{k} ⟩ .

(47)

The energy spectrum in a finite system with size L is bounded from below. At low $| \vec{k} | \leq π / 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^42,48 takes place at

t < t_{0}

⁠, after which the pumping stops. At this moment, the energy spectrum is growing with wavevector by one of the following two possible laws (with P being the Birkhoff–Saffman invariant of the fluid and M being the Loitzansky invariant):

{\begin{array}{l} E (k, t) \propto P k^{2}, \\ E (k, t) \propto M k^{4} . \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 ≫ π / 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² spectrum at small k and some universal decay at large k, reflecting these distributed vortex structures.

The decaying energy

E_{d} (t)

⁠, given by the part of the spectrum

k > k_{0} \sim 1 / L

⁠, has the following form:

E_{d} (t) = \int_{k_{0}}^{\infty} d k E (k, t) \propto \frac{4 π \tilde{ν}}{t} \int_{k_{0} \sqrt{\tilde{ν} t}}^{\infty} H (κ) d κ .

(49)

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

A. 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

- E (t)

] given the zero energy left at infinite time.

E (t) = \int_{t}^{\infty} E (t') d t' .

(50)

This energy dissipation rate

E (t')

is calculated as

\begin{matrix} E (t) \propto 4 π ν \int d k k^{2} ⟨ \vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k}) ⟩ \\ = 4 π \frac{\tilde{ν}}{t^{2}} \int_{k_{0} \sqrt{\tilde{ν} t}}^{\infty} κ^{2} H (κ) d κ . \end{matrix}

(51)

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

This limit was computed in Ref. 31 in a slightly different grand canonical ensemble, where N was fluctuating with the weight $exp (- μ N), μ \to 0$ ⁠.

With our current ensemble of fixed even

N \to \infty

⁠, the results of Ref. 31 read

E_{\infty} (t) = \frac{\tilde{ν}}{t^{2}} \frac{π^{2}}{864 ζ (3)} .

(52)

In our present theory, the same quantity is given by the above integral at

k_{0} = 0

E_{\infty} (t) = 4 π \frac{\tilde{ν}}{t^{2}} \int_{0}^{\infty} κ^{2} H (κ) d κ .

(53)

Comparing these two expressions, we get the normalization of

H (κ)

as

4 π \int_{0}^{\infty} κ^{2} H (κ) d κ = \frac{π^{2}}{864 ζ (3)} .

(54)

The integral on the left can be further reduced³⁹ to the following normalization condition:

Z = \frac{276 480 ζ (3) ζ (\frac{15}{2})}{119 ζ (\frac{17}{2})} \int_{Δ_{1}}^{Δ_{2}} d Δ \frac{P (Δ)}{Q (Δ)},

(55)

\begin{matrix} P (Δ) = (1 - Δ) \sqrt{r_{0} (Δ)} S (Δ) Q_{α} (Δ, 1) \\ (π J (Δ) (r_{0} (Δ) + 12) + L (Δ) (r_{0} (Δ) - 6)), \end{matrix}

(56)

Q (Δ) = L {(Δ)}^{2} (r_{0} (Δ) + 12) .

(57)

This normalization constant

Z

can be used in Eq. for the energy decay in a finite system. All the functions of Δ were defined above.

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

t^{2} E (t) = 4 π \tilde{ν} F (k_{0} \sqrt{\tilde{ν} t}),

(58a)

F (κ) = \int_{κ}^{\infty} H (x) x^{2} d x,

(58b)

\sqrt{t} E (k, t) \propto 4 π \tilde{ν} \sqrt{\tilde{ν}} H (k \sqrt{\tilde{ν} t}) .

(58c)

Both the energy dissipation and the energy spectrum are related to the same function

H (κ)

⁠, but the energy spectrum related to this function at large argument

κ = k \sqrt{\tilde{ν} t}

⁠, whereas the energy dissipation is related to integral of this function from small argument

κ = k_{0} \sqrt{\tilde{ν} t}

to infinity.

Our theory does not have the Saffman 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{ν} t} ≪ 1

⁠)

E (k_{0}, t) \propto \frac{4 π}{\sqrt{t}} \tilde{ν} \sqrt{\tilde{ν}} H (k_{0} \sqrt{\tilde{ν} t}) \sim \frac{4 π H (0) \tilde{ν} \sqrt{\tilde{ν}}}{\sqrt{t}} .

(59)

This will match the Saffman spectrum

P k^{2}

at small k with

k \sim {\tilde{ν}}^{3 / 4} t^{- 1 / 4} P^{- 1 / 2} .

(60)

This boundary value decays as

t^{- 1 / 4}

⁠; 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{ν} t} ≫ 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 (κ)

in Appendix K and numerically integrated it to obtain

F (κ), κ = k \sqrt{\tilde{ν} t}

⁠. Asymptotically, at large

k \sqrt{\tilde{ν} t}

⁠,

E (k, t) \propto \frac{\tilde{ν} \sqrt{\tilde{ν}}}{\sqrt{t}} {(k \sqrt{\tilde{ν} t})}^{- 7 / 2},

(61)

E (t) \propto t^{- 9 / 4},

(62)

E (t) = \int_{t}^{\infty} E (t') d t' \propto t^{- 5 / 4} .

(63)

We computed the effective critical indexes

n (t) = - t \partial_{t} log E = - \frac{t E (t)}{E (t)},

(64)

s (κ) = - κ \partial_{κ} log H (κ)

(65)

numerically in 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 yield $n (\infty) \approx 1.25 \pm 0.02$ ⁠, which agrees with our theoretical prediction in Fig. 6.

Our universal curves for $n (t), s (κ)$ 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).

B. 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⁵⁸ with viscous thickness $w \sim \sqrt{ν}$ ⁠. Integrating the square of vorticity in the Burgers vortex in the transverse plane, we get a large factor $1 / ν$ ⁠; this compensation leads to anomalous dissipation $E = ν ⟨ {\vec{ω}}^{2} ⟩$ ⁠, independent of ν.

In our dual theory, this factor of ν is compensated by a different mechanism: the vorticity is represented as a discontinuity at the curve $\vec{P} (θ)$ in our solution: $\vec{ω} \propto \vec{P} \times Δ \vec{P}$ ⁠. Summing over a large number $N \to \infty$ of these discontinuities leads to the factor N² compensating small viscosity ν factor in front of the enstrophy. The enstrophy integral $\int d k k^{2} ⟨ \vec{ω} \cdot {\vec{ω}}_{k} ⟩$ converges at large k.

As a consequence, the vorticity correlation $⟨ \vec{ω} \cdot {\vec{ω}}_{k} ⟩ \sim 1 / ν$ is large at all k, not just at large k. Our limit $ν \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 as follows: The turbulent limit differs from the inviscid limit of the Navier–Stokes equation. In the turbulent limit, the average circulation Γ 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

⟨ {\vec{ω}}^{2} ⟩ = \frac{E (t, | Γ | / ν)}{ν} .

(66)

The dissipation

E (t, | Γ | / ν)

depends on time and the effective Reynolds number

Rey = | Γ | / ν

⁠, where Γ is a typical velocity circulation in our problem.

The turbulent limit corresponds to

| Γ | ≫ ν

⁠, rather than

ν \to 0

⁠. Mathematically, we are looking for a residue of the enstrophy at zero viscosity

⟨ {\vec{ω}}^{2} ⟩ \to \frac{E (t, \infty)}{ν} .

(67)

However, in the real physical world (or in the DNS), we shall use this formula for finite ν corresponding to the actual viscosity of water. We only use our dual theory to compute this residue

\underset{ν = 0}{Res} ⟨ {\vec{ω}}^{2} ⟩ = E (t, \infty)

⁠.

⟨ {\vec{ω}}^{2} ⟩ \to \frac{E {(t, \infty)}_{dual}}{ν} .

(68)

In particular, in the infinite system (⁠

k_{0} = 0

⁠), we have the value computed in the previous paper.

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

⟨ \vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k}) ⟩ = \frac{{\tilde{ν}}^{5 / 2} H (k \sqrt{\tilde{ν} t})}{ν \sqrt{t}},

(69)

\sqrt{t} E (k, t) = \frac{4 π {\tilde{ν}}^{5 / 2} H (k \sqrt{\tilde{ν} t})}{ν},

(70)

where ν is the physical viscosity of the fluid under consideration (say, water) and

\tilde{ν}

is the auxiliary scale parameter that comes from the dual theory in the turbulent limit.

Our theory has two unknown parameters: $\tilde{ν}$ and t₀. The energy spectrum decreased faster than $k^{- 3}$ so that the enstrophy integral $\int H (κ) κ^{2} d κ$ converges at large κ and so does the energy integral $\int H (κ) d κ$ ⁠.

As we shall see in Sec. VIII, this fast decay of the energy spectrum also happens in DNS, in qualitative agreement with our decay $κ^{- 7 / 2}$ ⁠.

The presence of a hidden second scale $\tilde{ν}$ in our theory invalidates the formal scale invariance of Navier–Stokes equation used to determine the dependence of the energy spectrum from viscosity in the Onsager-like scaling theories.¹² This phenomenon of the hidden finite scale in the field theory is the same as dimensional transmutation in QCD, where the UV divergences combined with vanishing bare coupling constant create a finite mass scale. The physical origins of the hidden scale in Navier–Stokes equation are the thermal fluctuations of the initial velocity field at molecular scales, amplified in a turbulent flow to spontaneous stochasticity.

VIII. COMPARING OUR THEORY WITH THE DNS

As the first such test of our theory, we took the raw DNS data from Ref. 48, provided to us by the authors. These data are now available online⁴² and can be downloaded without permission. We only compared the data corresponding to our k² case and restricted ourselves to four samples with the largest grid 1024³, 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 Fig. 10. 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.

FIG. 10.

View large Download slide

The time decay of Reynolds numbers for each of the four samples from Refs. 42 and 48.

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

L (t) = \frac{\int E (k, t) kdk}{\int E (k, t) k^{2} d k} .

(71)

The effective length L(t) as a function of time is shown in Fig. 11.

FIG. 11.

The log-log plot of the effective length L(t) for the sample 14 from Refs. 42 and 48. The turbulent part, 1000 < t < 8000 , closely fits our theory L ( t ) ∝ t .

View large Download slide

The log-log plot of the effective length L(t) for the sample 14 from Refs. 42 and 48. The turbulent part, $1000 < t < 8000$ ⁠, closely fits our theory $L (t) \propto \sqrt{t}$ ⁠.

The statistical equilibrium was not yet reached at $t < 1000$ ⁠, so we discarded this period. The late stages of decay where $L (t) \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(t) for the turbulent region of sample 14, defined as

E (t) = \int E (k, t) d k,

(72)

log E (t) \approx a + f (b + log L (t)) .

(73)

Our theory has two free parameters

t_{0}, \tilde{ν}

⁠. The fitting of t₀ 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 Ref. 53, we avoid fitting t₀ by investigating the decaying energy as a function of the decay length L(t), which in our case scales as

\sqrt{\tilde{ν} (t + t_{0})}

⁠. The precise definition was given in .

The parameters a, b were fitted by nonlinear regression using “NonlinearModelFit” in Mathematica^®. The resulting log-log plot is shown in Fig. 12. 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 = 1 / 2$ for L(t), which we already tested.

FIG. 12.

View large Download slide

The log-log plot of the decaying energy as a function of decaying length scale L(t). Red dots are the DNS data, the green curve is an exact theoretical curve, and the dashed blue line is its asymptotic limit $E (t) \propto L {(t)}^{- 5 / 2}$ corresponding to $E \propto t^{- 5 / 4}$ ⁠.

Let us now turn to the energy spectrum. The data are not as good here as the energy decay data for two reasons. First, each wavevector 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 (L (L + 1) (L + 2) / 6)$ different values between 0 and $L \frac{\sqrt{3}}{2}$ ⁠. Unfortunately, the available data^42,48 aggregate the statistics at 512 equidistant bins in $| \vec{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 (k \sqrt{t + t_{0}})}{\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(t). 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. First, we consider the Mellin transform of the energy spectrum, normalized at p = 0 together with its derivative

M (p, t) = \frac{e^{- α p} \int d k E (k, t) {(k L (t))}^{p}}{\int d k E (k, t)},

(74)

α = {⟨ \frac{\int d k E (k, t) log [(k L (t))]}{\int d k E (k, t)} ⟩}_{t},

(75)

\partial_{p} {⟨ M (p, t) ⟩}_{t} = 0.

(76)

The coefficient α corresponds to the global scaling factor in the momentum k. Our definition of p differ from conventional

p' = - p - 1

⁠.

This function does not depend on time in the turbulent region described by our theory. So, all DNS data must collapse on a single universal curve. The K41 theory would correspond to a single pole at

p = \frac{2}{3}

⁠, which means for the normalized function

M_{K 41} (p, t) = \frac{exp (- \frac{3}{2} p)}{1 - \frac{3}{2} p} .

(77)

We computed this Mellin integral for the DNS data and compared it with our theory as well as K41 (see Fig. 13). We fitted the global scaling factor α in our theory to the DNS data by means of mean squares for the

log {⟨ M (p, t) ⟩}_{t}

(see Ref. 36). Our match with DNS is good but not perfect, perhaps due to finite size effects on 1 K lattice and finite Reynolds number

Rey \sim 80

⁠.

FIG. 13.

View large Download slide

The Mellin integral (74) for the DNS (red dots with error bars), our theory (green curve), and K41 pole term (77), shown as a blue curve.

The K41 is totally wrong here. There is no pole at $p = \frac{2}{3}$ as prescribed by K41. Large systematic errors invalidate the DNS Mellin transform with large positive p, dominated by lattice artifacts at large k.

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

⟨ Δ {\vec{v}}^{2} ⟩ (r) = \int \frac{d^{3} k}{4 π^{3}} \frac{(1 - e^{ı \vec{k} \cdot \vec{r}}) E (k, t)}{4 π {\vec{k}}^{2}},

(78)

⟨ Δ {\vec{v}}^{2} ⟩ (r) = ⟨ {(\vec{v} (\vec{r}, t) - \vec{v} (0, t))}^{2} ⟩ .

(79)

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

(80)

\begin{matrix} ξ_{2} (r, t) = r \partial_{r} log ⟨ Δ v^{2} ⟩ (r) \\ = \frac{\int_{0}^{\infty} d k (- cos (k r) + \frac{sin (k r)}{k r}) E (k, t)}{\int_{0}^{\infty} d k (1 - \frac{sin (k r)}{k r}) E (k, t)}, \end{matrix}

(80a)

{⟨ ξ_{2} (x L (t), t) ⟩}_{t} = f (x) .

(80b)

This universal function f(x) 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 14 shows what we found instead for the DNS data.^42,48

FIG. 14.

View large Download slide

The plot of f(x) in (80) as a function of $log x$ (red dots with error bars). The green curve is the theoretical index without adjustment of the coordinate scale, and the dashed blue line is the K41 constant value $\frac{2}{3}$ ⁠.

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 $log x$ 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.

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 Fig. 9).

The theoretical curve (green line) is shifted left by some amount $δ log x = 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 (t) \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.²⁴ Our theory assumes incompressibility, so it does not apply to these data in the air turbulence. Also, the magnitude of the errors in Ref. 24 is unclear.

Nevertheless, we still compared our prediction for the second moment with the air data from Ref. 24. The shape of the experimental curve for $⟨ Δ v^{2} ⟩ (r)$ (Fig. 1 in Ref. 24) 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 (Fig. 1 in Ref. 24; the slope linearly decreases with $log r$ ⁠, similar to our decreasing slope in Fig. 14.

However, the numerical values for the slope and curvature in Ref. 24 are quite different from those we derived above from the incompressible DNS,^42,48 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²⁴ without increasing errors, which are already large enough in the data in Refs. 42 and 48. Unlike numerical differentiation of the wind tunnel data, the Fourier integration of the DNS data suppress the noise, which makes the DNS data for effective index $ξ_{2} (r)$ 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.

IX. 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.

A. The energy spectrum

The spectral density

H (κ)

is in . The Mellin transform is the following integral, assuming our function decreases at infinity:

H (κ) = \int_{- ϵ - ı \infty}^{- ϵ + ı \infty} \frac{d p}{2 π ı} h (p) κ^{p},

(81)

h (p) = \int_{0}^{\infty} \frac{d κ}{κ} κ^{- p} H (κ) .

(82)

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(p) having a single pole in the left semi-plane. The position

p = - a

of this pole becomes an index of the power law that

H (κ) \propto κ^{- a} .

(83)

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 (n)

⁠. In application to the moments of velocity difference, this parameter n is the degree of the moment

M_{n} = ⟨ Δ v^{n} ⟩

⁠. This pole will produce multifractal scaling laws that

M_{n} \sim r^{a (n)} .

(84)

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

h (p) = \frac{f (p) ζ (p + \frac{15}{2}) Γ (- p)}{(2 p + 7) (2 p + 17) ζ (p + \frac{17}{2})},

(85)

f (p) = 20 \int_{Δ_{1}}^{Δ_{2}} d Δ (1 - Δ) C^{p - 1} (A C - B p),

(86)

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

This convergence makes this function an entire function without any finite singularities. The values of

C (Δ)

are bounded by two positive limits such that

0.054 198 4 < C < 0.063 075 5, \forall Δ_{1} < Δ < Δ_{2} .

(87)

Therefore, the entire function f(p) decreases as

e^{- 2.76342 p}

in the right semiplane and grows as

e^{- 2.9151 p}

in the left semiplane. It oscillates along an imaginary axis, which is our integration path.

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

The singularities of the Mellin transform h(p) for the energy spectrum

H (κ)

are given by the following table of simple poles:

| \begin{matrix} energy spectrum indexes \\ - 7 / 2 \\ - 13 / 2 \\ - 8 \pm ı t_{n} if n \in ℤ \\ - 17 / 2 - 2 n if n \in ℤ \land n \geq 0 \\ n if n \in ℤ \land n \geq 0 \end{matrix} | .

(88)

Here

\pm t_{n}

are imaginary parts of zeros of the ζ function on the critical line

z = 1 / 2 + ı t

⁠. About 10¹³ zeros are already known, though the Riemann hypothesis (no other complex zeros) still needs to be proven. These poles are shown in Fig. 15.

FIG. 15.

Complex poles of the Mellin transform h(p) of H ( κ ) in (85) defining the critical indexes of the energy spectrum as a function of κ = | k → | ν ̃ t .

View large Download slide

Complex poles of the Mellin transform h(p) of $H (κ)$ in (85) defining the critical indexes of the energy spectrum as a function of $κ = | \vec{k} | \sqrt{\tilde{ν} t}$ ⁠.

Only poles with negative real parts contribute to the power expansion at

κ \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.

| W_{n} | κ^{- 8} cos (t_{n} log κ + arg W_{n}) .

(89)

These slow oscillations are visible as regime change in the log-log plots of effective indexes in Figs. 5 and 6.

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

There is, however, an apparent contradiction between our fast decay $k^{- 7 / 2}$ of the energy spectrum and the bounds on this decay index established by Sulem and Fricsh.⁵⁷

At a closer look, this paradox is explained by an assumption made in Ref. 57 that the whole energy dissipation is caused by the Kolmogorov-like flux $\vec{v} \cdot (\vec{v} \cdot \vec{\nabla}) \vec{v}$ ⁠. The viscous contribution $ν {\vec{ω}}^{2}$ to the dissipation was neglected in that paper: they assumed that it “converges to zero as $ν \to 0$ ⁠.” Thus, their bound on decay index implies that our dissipation has zero Kolmogorov energy flux; therefore, it is dominated by the dissipative anomaly $E = lim_{ν \to 0} ν ⟨ {\vec{ω}}^{2} ⟩$ ⁠, missed in that old paper.

Our theory, based on the exact solution of the Navier–Stokes equation, does not need any assumptions about the “energy cascade” or “Kolmogorov energy flux.” In the physical picture described in Subsection IV A, there is no such flux: the energy pumped from the boundaries (such as an oscillating grid in the grid turbulence experiment¹⁰) is dissipated in the thin vortex filaments due to dissipative anomalies of such filaments.

B. The energy decay

The energy E(t) is related to the same function

\begin{matrix} E (t) = \int_{t}^{\infty} d t E (t) = 4 π \tilde{ν} \int_{t}^{\infty} \frac{d T}{T^{2}} \int_{k_{0} \sqrt{\tilde{ν} T}}^{\infty} H (x) x^{2} d x . \end{matrix}

(90)

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

(91)

E (t) = \int_{- 1 - ϵ - ı \infty}^{- 1 - ϵ + ı \infty} \frac{d q}{2 π ı} e (q) {(k_{0}^{2} \tilde{ν} t)}^{q},

(91a)

e (q) = 2 π k_{0}^{2} {\tilde{ν}}^{2} \frac{h (2 q - 1)}{q (q + 1)} .

(91b)

The table of complex poles of this function is

| \begin{matrix} energy decay indexes \\ - 5 / 4 \\ - 11 / 4 \\ - 7 / 2 \pm ı / 2 t_{n} if n \in ℤ \\ - 15 / 4 - n if n \in ℤ \land n \geq 0 \\ n / 2 if n \in ℤ \land n \geq 0 \end{matrix} | .

(92)

These poles are shown in Fig. 16. The leading pole is at

p = - 5 / 4

⁠, 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.

FIG. 16.

View large Download slide

Complex poles of the Mellin transform e(q) of E(t) in (91) defining the spectrum of critical indexes of remaining energy as a function of time.

C. The vorticity correlation function in coordinate space

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

\begin{matrix} ⟨ Δ {\vec{v}}^{2} ⟩ (r) = 2 \int \frac{d^{3} k}{{(2 π)}^{3}} (1 - exp (ı \vec{k} \cdot (\vec{r}))) \frac{E (k, t)}{4 π {\vec{k}}^{2}} \\ = \frac{{\tilde{ν}}^{2}}{π^{2} ν t} \int_{ϵ - ı \infty}^{ϵ + ı \infty} \frac{d p}{2 π ı} {(\frac{| \vec{r} |}{\sqrt{\tilde{ν} t}})}^{p} \\ cos (\frac{π p}{2}) Γ (- p - 1) h (- 1 - p), \end{matrix}

(93)

⟨ Δ {\vec{v}}^{2} ⟩ (r) = ⟨ {(\vec{v} (\vec{r}, t) - \vec{v} (0, t))}^{2} ⟩ .

(94)

The Mellin transform simplifies to

\begin{matrix} ⟨ Δ {\vec{v}}^{2} ⟩ (r) = \frac{{\tilde{ν}}^{2}}{ν t} \int_{ϵ - ı \infty}^{ϵ + ı \infty} \frac{d p}{2 π ı} V (p) {(\frac{| \vec{r} |}{\sqrt{\tilde{ν} t}})}^{p}, \end{matrix}

(95)

\begin{matrix} V (p) = - \frac{f (- 1 - p) ζ (\frac{13}{2} - p) csc (\frac{π p}{2})}{16 π^{2} (p + 1) (2 p - 15) (2 p - 5) ζ (\frac{15}{2} - p)} . \end{matrix}

(96)

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, \vec{r})

⁠. Should our theory be a CFT (which it is not), this spectrum would be related to the spectrum of anomalous dimensions Δ_n in the OPE:

CFT : \vec{v} (0) \cdot \vec{v} (r) \sim \sum_{n} O_{n} {| r |}^{p_{n}}; p_{n} = 2 Δ_{v} - Δ_{n} .

We have such an expansion with

\frac{| \vec{r} |}{\sqrt{\tilde{ν} t}}

in place of

| r |

and a factor

\frac{{\tilde{ν}}^{2}}{ν 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:

| \begin{matrix} indexes of velocity correlation \\ - 1 \\ 0 \\ 2 n if n \in ℤ \land n \geq 1 \\ 5 / 2 \\ 11 / 2 \\ 1 / 2 (15 + 4 n) if n \in ℤ \land n \geq 0 \\ 7 \pm ı t_{n} if n \in ℤ \end{matrix} | .

(97)

The poles in the complex p plane are shown in Fig. 17.

FIG. 17.

View large Download slide

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 $\vec{v} (0) \cdot \vec{v} (r) \sim \sum_{n} O_{n} {| \vec{r} |}^{p_{n}}$ ⁠.

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

| \vec{r} | \to 0 / \infty

⁠. The leading term at

r \to 0

is r², which is calculable in general form from its definition after expanding the exponential and averaging over directions of

\vec{k}

⁠,

⟨ Δ {\vec{v}}^{2} ⟩ (r) \to \frac{r^{2} E_{\infty} (t)}{24 π^{3} ν} = \frac{r^{2} \tilde{ν}}{20736 π ζ (3) ν t^{2}} .

(98)

The next term is

r^{5 / 2}

as it follows from the table in . As mentioned about the energy spectrum, there is no K41 scaling index

p = 2 / 3

⁠. 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 Fig. 18 for the oscillation of the index ξ₂ at the latest stage). This region is inaccessible with modern scale of the DNS.

FIG. 18.

View large Download slide

Oscillations of the effective index $ξ_{2} (r)$ 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.

The theoretical curve in Fig. 19 agrees with the Fourier-transformed data of Refs. 42 and 48 but deviates from Ref. 24. The probable reason is the compressibility of the air in real experiments in Ref. 24.

FIG. 19.

The universal function (95) as a function of ρ = | r → | / ν ̃ t . The turnover is caused by subleading terms in the power expansion starting with ρ 2 . The next terms involve quantum oscillations manifesting as a turnover from power growth to power decay. Asymptotic at large ρ is const as it follows from the table (97) of the decay indexes.

View large Download slide

The universal function (95) as a function of $ρ = | \vec{r} | / \sqrt{\tilde{ν} t}$ ⁠. The turnover is caused by subleading terms in the power expansion starting with $ρ^{2}$ ⁠. The next terms involve quantum oscillations manifesting as a turnover from power growth to power decay. Asymptotic at large ρ is $const$ as it follows from the table (97) of the decay indexes.

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.

X. DISCUSSION

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

A. Myth and reality of turbulent scaling laws

More than 80 years ago, Kolmogorov and Obukhov made a breakthrough in turbulence theory by establishing the relation 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

\begin{matrix} ⟨ v_{α} ({\vec{r}}_{0}) v_{β} ({\vec{r}}_{0}) ω_{ν} (\vec{r} + {\vec{r}}_{0}) ⟩ \\ = e_{μ γ ν} \partial_{μ} (\vec{r}) ⟨ v_{α} ({\vec{r}}_{0}) v_{β} ({\vec{r}}_{0}) v_{ν} (\vec{r} + {\vec{r}}_{0}) ⟩ \\ \propto e_{μ γ ν} (δ_{α γ} δ_{β ν} + δ_{α ν} δ_{β γ} - \frac{2}{d} δ_{α β} δ_{μ γ}) = 0. \end{matrix}

(99)

The last relation follows from the symmetry of the expression in the brackets with respect to the exchange of the indexes

γ \Leftrightarrow μ

⁠.

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. In Fourier space, this term becomes linear combinations of gradients of delta function $\vec{\nabla} δ^{3} (\vec{k})$ ⁠, rather than some constant flow through the whole spectrum. Such a constant flow would require an $| x | log | x |$ term rather than a linear term in coordinates. This linear term is an example of the harmonic terms added to the Biot–Savart integral for the velocity field.

Thus, we are challenging not only K41 scaling laws but also the concept of the “Kolmogorov energy flux” $v v \nabla v$ ⁠. Indeed, our fast decay of the energy spectrum, according to the old bounds established in Ref. 57, corresponds to a vanishing Kolmogorov energy flux in the turbulent limit. The dissipative anomaly, neglected in that old work, explains the paradox. As we argued in Subsection IV A, the energy is pumped through the boundaries and dissipated directly on viscous micro-structures (Burgers filaments). Our exact computations of the dissipation and energy spectrum in Sec. IX A are compatible with the dissipative anomaly but not with the energy cascade.

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 lognormal 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^55,59,60 have invalidated the K41 scaling laws (including the lognormal model) over the past 30 years. Regarding decaying turbulence, the experimental data⁴⁸ have diverged even further from Kolmogorov scaling laws despite all attempts to stretch these data or discard the non-fitting region as “erosion.” We highlight significant deviations—six orders of magnitude—from the $k^{- 5 / 3}$ scaling in Fig. 20. There are also recent measurements²⁴ with significant deviations of the log-log derivative of the second velocity moment $⟨ Δ v^{2} ⟩ (r)$ from the $\frac{2}{3}$ predicted by K41.

FIG. 20.

View large Download slide

The log-log plot of deviation from K41 energy spectrum in decaying turbulence. Reproduced from Panickacheril John, Donzis, and Sreenivasan, Philos. Trans. R. Soc. A. 380(2218), 20210089 (2022), licensed under a Creative Commons Attribution (CC BY) license. The observed curve approximately matches our theoretical curve in Fig. 7. The K41 spectrum would correspond to a horizontal line, a total mismatch.

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,⁴⁹ which, as a phenomenological model, successfully described observed deviations from the K41 laws in forced turbulence.^55,59

However, there are no theoretical grounds for conformal symmetry in turbulence; the “current conservation” conditions $\partial_{α} v_{α} = 0, \partial_{γ} ω_{γ} = 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 Fig. 2 (top) in Ref. 24. 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 $log k$ ⁠, spanning several decades.

Both the DNS and experimental papers^24,48 noted significant deviations from scaling laws (whether K41 or multifractal). The conclusion of Ref. 48 was cautiously negative: “it is somewhat disappointing that the results are not more closely aligned with theoretical arguments.” The most recent paper²⁴ 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³² 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^42,48 to compute the effective index of the velocity correlation by numerical Fourier transform of their energy spectrum (see Sec. VIII). Our effective index is plotted in Fig. 9. The K41 scaling law $Δ v^{2} \sim r^{2 / 3}$ 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 (97).

In this way, our theory extends the multifractal scaling laws by accommodating an infinite discrete spectrum of scaling dimensions. According to recent DNS results,⁵⁶ 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 Sec. VIII.

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¹⁸ and Chandrasekhar⁸ 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:⁵⁴

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 (88). 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 Fig. 5, 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 $⟨ Δ v^{2} ⟩ (r)$ with effective index (log derivative) shown in Fig. 14. 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.

B. Stochastic solution of the Navier–Stokes equation and ergodic hypothesis

Richard Feynman wrote about turbulence in his Lectures in Physics:¹³

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 60 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^{- 5 / 4}$ ⁠. Feynman did not distinguish between decaying turbulence and the steady state. However, his example of the steady flow of water through the pipe belongs to the grid turbulence category covered by our theory.

The steady state of the water flow past the grid is achieved by forcing at the position of the grid, but not down the stream, which in our theory would be the “initial condition” in the coordinate frame moving with the water. In this frame, the energy pumping occurred initially when the water passed the grid, after which the energy decayed with time at every point of the water flow.

This is a steady state in the original frame because the decaying pieces of the flow are constantly swept downstream. Thus, the energy does not depend on time at a fixed distance from the grid. However, it decays with this distance $z = \int v_{z} d t$ by the same $z^{- 5 / 4}$ law as the time decay of the energy of the decaying turbulence.

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 π. 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.^29,30

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 (ı S / ℏ)$ 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 $Γ_{C} [v] = \oint_{C} \vec{v} \cdot d \vec{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 some positive 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 of (12), (20), and (24) exactly satisfies the Navier–Stokes equation at a finite time for the loop functional in the turbulent limit $N \to \infty, ν \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 $\vec{P} (θ)$ ⁠. The continuous quantum mechanical integral over phase space with equal weight $1 / (2 π)$ 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.

C. 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,21,22,27,41} where the loop equations were first derived.^25,26

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 (ı \vec{p} \cdot \vec{x})$ ⁠. This superposition is the solution we have found, except the dot product $\vec{p} \cdot \vec{x}$ becomes a symplectic form $\oint \vec{P} (θ) \cdot d \vec{C} (θ)$ in the loop space.

Our momentum $\vec{P} (θ, t)$ is not an integral of motion, but simple scaling properties of the pure Navier–Stokes equation lead to the solution with $\vec{P} (θ, t) \propto \vec{F} (θ) / \sqrt{t}$ ⁠, with $\vec{F} (θ)$ 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 $\vec{F} (θ)$ ⁠.

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 Sec. VII.

D. Classical flow and quantum mechanics

Our computations rely significantly on number theory, particularly Jordan's multitotients, $φ_{l} (q)$ ⁠, which extend the Euler totient function.⁴³ 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.²⁹ 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 Sec. X C.

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,^29,30 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 cancelation 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 $\vec{P} (ξ + 1) = \vec{P} (ξ)$ ⁠.

From the conventional perspective, the fractal curve in complex momentum space, $\vec{P} (θ) \in ℂ_{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 and Liouville theory^23,51 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 $ℂ_{3}$ ⁠. Anomalous dissipation is achieved through numerous finite discontinuities of the fractal curve $\vec{P} (ξ) \in ℂ_{3}$ ⁠.

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

E. Renormalizability of the inviscid limit of the loop equation

Let us examine the relationship $ν N^{2} = const$ between diminishing viscosity and the increasing number N of discontinuities on the momentum loop $\vec{P} (θ)$ ⁠.

The Navier–Stokes (NS) 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 $ℝ_{3}$ but stochastic variables, with $Δ \vec{v} \sim {(Δ \vec{r})}^{α}$ ⁠.

Fractal Calculus¹⁴ 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 $⟨ {(Δ \vec{v})}^{n} ⟩ \sim {(Δ \vec{r})}^{ζ_{n}}$ cannot be explained in this manner.

In our theory, there is a universal function (95) for the velocity correlation function as a function of the separation $\vec{r}$ ⁠. However, the $\vec{r}$ dependence is not defined by a single power: the Mellin transform of this function reveals infinitely many poles (97) 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 $\vec{P} (θ)$ 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₀ and the lattice spacing a in QCD: $a \propto g_{0}^{α} exp (- β / g_{0}^{2})$ ⁠.

The naive local limit $a \to 0, g_{0} = 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{ν} = ν N^{2} \to const$ ⁠. Furthermore, the energy spectrum, expressed as a function of $k \sqrt{\tilde{ν} 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{ν}$ ⁠, defining the timescale.

The energy scale is related to the fluid's physical viscosity ν in the denominator. Formally setting ν = 0 would tend the physical energy scale to infinity. The fictitious limit $ν \to 0$ was only taken to compute the turbulent limit of the energy spectrum as a function of Reynolds number $Γ / ν$ ⁠.

In the real world, this number is large because of the large circulation Γ compared to finite viscosity ν. Still, we can compute these dimensionless functions by taking the inviscid limit, as described by Eq. (69) in Sec. VII B.

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

F. 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,²⁵ 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 $\vec{v}$ ⁠. As extensively discussed in Refs. 25, 26, and 30, any Stokes-type functional $Ψ (γ, C)$ that satisfies the boundary condition at a shrunk loop $Ψ [γ, 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 Refs. 29 and 30 (and also here, in Sec. IV B) how the velocity distribution for random uniform vorticity in the fluid was reproduced by a singular momentum loop $\vec{P} (θ)$ ⁠.

The solution for $\vec{P} (θ)$ at this specific fixed point of the loop equation was complex and had slowly decreasing Fourier coefficients, leading to a discontinuity $sign (θ - θ')$ 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 / ν$ 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.

XI. 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 $ν \to 0, N \to \infty$ ⁠. The nth moment of $Δ \vec{v} = \int d^{3} k \frac{ı \vec{k} \times {\vec{ω}}_{k}}{{\vec{k}}^{2}} (exp (ı \vec{k} \cdot \vec{r}) - exp (ı \vec{k} \cdot \vec{r}'))$ 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³¹ can be evaluated to compute vorticity correlation functions in the Navier–Stokes equation with an infinitesimal random force.

XII. 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 80 years, but it perfectly fits the available experimental data^9,10 and recent DNS data.^42,48,53 It offers an intriguing perspective, revealing unexpected connections between nonlinear classical physics and quantum mechanics.

DEDICATION

To dear friend and an outstanding physicist, Sreeni, my constant source of knowledge and support.

ACKNOWLEDGMENTS

Maxim Bulatov helped me with mathematical suggestions and discussions. He also participated in an unpublished paper⁷ where we simulated the Euler ensemble as a Markov chain. The strong cancellations in the vorticity correlation function prevented this simulation from producing reliable numerical results.

I filled the gaps in my knowledge of Number Theory in discussions of the Euler ensemble with Peter Sarnak, Alexandru Zaharesku, Konstantin Khanin, and Debmalya Basak.

Yang-Hui and I discussed Euler's totients at the Cambridge University workshop, where this theory was first reported in November 2023. His comments helped me derive the asymptotic distribution for scaling variables.

I am also grateful to the organizers and participants of the “Field Theory and Turbulence” workshop at ICTS in Bengaluru, India, where this work was advanced in December 2023. Discussions with Katepalli Sreenivasan, Rahul Pandit, and Gregory Falkovich were especially useful. They helped me understand the physics of decaying turbulence in a finite system and match my solution with the DNS data.

This work was discussed at the “Conformal Field Theory, Integrability, and Geometry” conference in Stony Brook on March 11–15, 2024. I am very grateful to Nikita Nekrasov, Sasha Polyakov, Sasha Zamolodchikov, Dennis Sullivan, and other participants for deep and inspiring discussions.

This work was also thoroughly discussed for a week at the Perimeter Institute for Theoretical Physics, where I spoke at a colloquium on April 3. I received many interesting questions and comments, reflected in Sec. X of this paper. I thank Davide Gaiotto, Luis Lehner, Sabrina Pastersky, Sergey Sibiryakov, Lee Smolin, and Andrey Shkerin for inspiring discussions. The help from Zechuan Zheng in Mathematica computations of the energy spectrum was also very useful.

I recently presented my theory at the 126th Statistical Mechanics Conference in Rutgers on May 19. There were very useful discussions, and I also got some unpublished data from John Panickacheril John and Diego A. Donzis.

After the Oxford Physics Department colloquium, Alex Schekochihin and I discussed the data and physics of decaying turbulence. I am grateful to Alex for these deep discussions and sharp comments.

The most important help and support came from Sreeni, who, for the last year, discussed the physics of decaying turbulence with me and guided me through the maze of inconsistent DNS and physical experiment data.

This research was supported by a Simons Foundation Award No. 686282 at NYU Abu Dhabi.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Alexander Migdal: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are openly available in the shared Google Drive, Ref. 42.

APPENDIX A: GLOBAL RANDOM ROTATION AND MOMENTUM LOOP SPACE

This path integral was computed in Refs. 30 and 31 for a special stochastic solution of the Navier–Stokes equation: the global rotation with Gaussian random rotation matrix. The initial velocity distribution is Gaussian, with a slowly varying correlation function. The corresponding loop field reads (we set γ = 1 for simplicity in this section)

\begin{matrix} Ψ_{0} [C] \equiv Ψ {(1, C)}_{t = 0} \\ = exp (- \frac{1}{2 ν^{2}} \int_{C} d \vec{C} (θ) \cdot d \vec{C} (θ') f (\vec{C} (θ) - \vec{C} (θ'))), \end{matrix}

(A1)

where

f (\vec{r})

is the velocity correlation function given by

⟨ v_{α} (r) v_{β} (r') ⟩ = (δ_{α β} - \partial_{α} \partial_{β} \partial_{μ}^{- 2}) f (r - r') .

(A2)

The potential part drops out in the closed loop integral. The correlation function varies at the macroscopic scale, which means that one could expand it in the Taylor series

f (r - r') \to f_{0} - f_{1} {(r - r')}^{2} + \dots .

(A3)

The first term f₀ is proportional to the initial energy density

\frac{1}{2} ⟨ v_{α}^{2} ⟩ = \frac{d - 1}{2} f_{0},

(A4)

and the second one is proportional to initial energy dissipation rate

E_{0}

f_{1} = \frac{E_{0}}{2 d (d - 1) ν},

(A5)

where d = 3 is the dimension of space. The constant term in as well as

r^{2} + r'^{2}

terms drop from the closed loop integral, so we are left with the cross term

r r'

⁠, which reduces to a full square

Ψ_{0} [C] \to exp (- \frac{f_{1}}{ν^{2}} {(\oint d C_{α} (θ) C_{β} (θ))}^{2}) .

(A6)

This distribution is almost Gaussian: it reduces to Gaussian one by extra integration

\begin{matrix} Ψ_{0} [C] \to const \int (d ϕ) exp (- ϕ_{α β}^{2}) exp (2 ı \frac{\sqrt{f_{1}}}{ν} ϕ_{μ ν} \oint d C_{μ} (θ) C_{ν} (θ)) . \end{matrix}

(A7)

The integration here involves all $\frac{d (d - 1)}{2} = 3$ independent $α < β$ components of the antisymmetric tensor $ϕ_{α β}$ ⁠. Note that this is ordinary integration, not the functional one.

This distribution can be translated into the momentum loop space. The resulting stochastic function

\vec{P} (θ)

is defined by the Fourier expansion on the circle as

P_{α} (θ) = \sum_{odd n = 1}^{\infty} P_{α, n} e^{ı n θ} + {\bar{P}}_{α, n} e^{- ı n θ},

(A8)

P_{α, n} = N (0, 1),

(A9)

{\bar{P}}_{α, n} = \frac{4 \sqrt{f_{1}}}{n ν} ϕ_{α β} P_{β, n},

(A10)

ϕ_{α β} = - ϕ_{β α},

(A11)

ϕ_{α β} = N (0, 1) \forall α < β .

(A12)

At fixed tensor

ϕ

⁠, the correlations are

{⟨ P_{α, n} {\bar{P}}_{β, m} ⟩}_{t = 0} = \frac{4 \sqrt{f_{1}}}{m ν} δ_{n m} ϕ_{α β},

(A13)

{⟨ P_{α} (θ) P_{β} (θ') ⟩}_{t = 0} = 2 ı \frac{\sqrt{f_{1}}}{ν} ϕ_{α β} sign (θ' - θ),

(A14)

Ψ_{0} [C] = {⟨ exp (\frac{ı}{ν} \oint d \vec{C} (θ) \vec{P} (θ)) ⟩}_{P, ϕ} .

(A15)

Though this special solution does not describe isotropic turbulence, it helps understand the mathematical properties of the loop technology. In particular, it shows the significance of the discontinuities of the momentum loop

\vec{P} (θ)

⁠.

APPENDIX B: THE MARKOV CHAIN AND ITS FERMIONIC REPRESENTATION

Here is a new representation of the Euler ensemble, leading us to the exact analytic solution below.

We start by replacing independent random variables σ with fixed sum by a Markov process, as suggested in Ref. 31. We start with n random values of $σ_{i} = 1$ and remaining N–n values of $σ_{i} = - 1$ ⁠. Instead of averaging over all of these values simultaneously, we follow a Markov process of picking $σ_{N}, \dots, σ_{1}$ one after another. At each step, there will be $M = N, \dots, 0$ remaining σ. We get a transition $n \Rightarrow n - 1$ with probability $\frac{n}{M}$ and $n \Rightarrow n$ with complementary probability.

Multiplying these probabilities and summing all histories of the Markov process is equivalent to the computation of the product of the Markov matrices

\prod_{M = 1}^{N} Q (M),

(B1)

Q (M) | n ⟩ = \frac{M - n}{M} | n ⟩ + \frac{n}{M} | n - 1 ⟩ .

(B2)

This Markov process will be random until n = 0. After that, all remaining σ_k will have negative signs and be taken with probability 1, keeping n = 0.

The expectation value of some function

\hat{X} ({σ})

reduces to the matrix product

\begin{matrix} P [\hat{X}] = \sum_{n = 0}^{N_{+}} \\ ⟨ n | (\prod_{M = 1}^{N} \hat{Q} (M)) \cdot \hat{X} \cdot | N_{+} ⟩, \end{matrix}

(B3)

\begin{matrix} \hat{Q} (M) \cdot \hat{X} | n ⟩ = \frac{n}{M} \hat{X} (σ_{M} \to 1) | n - 1 ⟩ + \frac{M - n}{M} \hat{X} (σ_{M} \to - 1) | n ⟩ . \end{matrix}

(B4)

Here $N_{+} = (N + \sum σ_{l}) / 2 = (N + q r) / 2$ is the number of positive sigmas. The operator $\hat{Q} (M)$ sets in $\hat{X} | n ⟩$ the variable σ_M to 1 with probability $\frac{n}{M}$ and to –1 with complementary probability. The generalization of the Markov matrix Q(M) to the operator $\hat{Q} (M)$ will be presented shortly.

Once the whole product is applied to $\hat{X}$ ⁠, all the sigma variables in all terms will be specified so that the result will be a number.

This Markov process is implemented as a computer code in Ref. 7, leading to a fast simulation with $O (N^{0})$ memory requirement.

Now, we observe that quantum Fermi statistics can represent the Markov chain of Ising variables. Let us construct the operator

\hat{Q} (M)

with Fermionic creation and annihilation operators, with occupation numbers

ν_{k} = (1 + σ_{k}) / 2 = (0, 1)

⁠. These operators obey (anti)commutation relations, and they create/annihilate σ = 1 state as follows (with Kronecker delta

δ [n] \equiv δ_{n, 0}

⁠):

{[a_{i}, a_{j}^{†}]}_{+} = δ_{i j},

(B5)

{[a_{i}, a_{j}]}_{+} = {[a_{i}^{†}, a_{j}^{†}]}_{+} = 0,

(B6)

\begin{matrix} a_{n}^{†} | σ_{1}, \dots, σ_{N} ⟩ = δ [σ_{n} + 1] | σ_{1}, \dots, σ_{n} \to 1, \dots, σ_{N} ⟩, \end{matrix}

(B7)

\begin{matrix} a_{n} | σ_{1}, \dots, σ_{N} ⟩ = δ [σ_{n} - 1] | σ_{1}, \dots, σ_{n} \to - 1, \dots, σ_{N} ⟩, \end{matrix}

(B8)

{\hat{ν}}_{n} = a_{n}^{†} a_{n},

(B9)

{\hat{ν}}_{n} | σ_{1}, \dots, σ_{N} ⟩ = δ [σ_{n} - 1] | σ_{1}, \dots, σ_{N} ⟩ .

(B10)

The number n(M) of positive sigmas

\sum_{l = 1}^{M} δ [σ_{l} - 1]

coincides with the occupation number of these Fermi particles such that

\hat{n} (M) = \sum_{l = 1}^{M} {\hat{ν}}_{l} .

(B11)

This relation leads to the representation

\hat{Q} (M) = {\hat{ν}}_{M} \frac{\hat{n} (M)}{M} + (1 - {\hat{ν}}_{M}) \frac{M - \hat{n} (M)}{M} .

(B12)

The variables σ_l can also be expressed in terms of this operator algebra by using

{\hat{σ}}_{l} = 2 {\hat{ν}}_{l} - 1.

(B13)

The Wilson loop in can now be represented as an average over the small Euler ensemble

E (N)

of a quantum trace expression

Ψ (γ, C) = \frac{{⟨ \hat{W} [C] ⟩}_{\hat{Ω}, E (N)}}{{⟨ \hat{W} [0] ⟩}_{\hat{Ω}, E (N)}},

(B14)

\hat{W} [C] = tr (\hat{Z} (q r) exp (\frac{i γ \hat{Γ}}{ν}) \prod_{M = 1}^{N} \hat{Q} (M)),

(B15)

\hat{Γ} = \sum_{l} Δ {\vec{C}}_{l} \cdot \hat{Ω} \cdot {\vec{P}}_{l} (t),

(B16)

\hat{Z} (s) = \oint \frac{d ω}{2 π} exp (ı ω (\sum_{l} {\hat{σ}}_{l} - s)),

(B17)

Δ {\vec{C}}_{l} = \vec{C} (\frac{l + 1}{N}) - \vec{C} (\frac{l}{N}),

(B18)

{\vec{P}}_{l} (t) = \sqrt{\frac{ν}{2 (t + t_{0})}} \frac{{\vec{F}}_{l}}{γ}, \hat{Ω} \in O (3),

(B19)

{\vec{F}}_{l} = \frac{{cos ({\hat{α}}_{l}), sin ({\hat{α}}_{l}), 0}}{2 sin (\frac{β}{2})},

(B20)

\begin{matrix} E (N) : p, q, r \in ℤ - N \leq q r \leq N, \end{matrix}

with

\begin{matrix} 0 < p < q < N, gcd (p, q) = 1, \end{matrix}

(B21)

{\hat{α}}_{l} = β \sum_{k = 1}^{l - 1} (2 {\hat{ν}}_{k} - 1) .

(B22)

The last component of the vector ${\vec{F}}_{l}$ is set to 0 as this component does not depend on l and yields zero in the sum over the loop $\sum_{l} Δ {\vec{C}}_{l} = 0$ ⁠.

The proof of equivalence to the combinatorial formula with an average over $σ_{l} = \pm 1$ can be given using the following Lemma (obvious for a physicist).

Lemma 1. The operators ${\hat{ν}}_{l}$ all commute with each other.

Proof. Using commutation relations, we can write

\begin{matrix} {\hat{ν}}_{l} {\hat{ν}}_{n} = a_{l}^{†} (δ_{l n} - a_{n}^{†} a_{l}) a_{n} = a_{l}^{†} a_{n} δ_{l n} - a_{l}^{†} a_{n}^{†} a_{l} a_{n} . \end{matrix}

(B23)

Interchanging indexes l, n in this relation, we see that the first term does not change due to Kronecker delta, and the second term does not change because

a_{l}^{†}, a_{n}^{†}

anti-commute, as well as a_l, a_n, so the second term is symmetric as well. Therefore,

{\hat{ν}}_{l} {\hat{ν}}_{n} = {\hat{ν}}_{n} {\hat{ν}}_{l}

⁠.

Quantum Trace Theorem. The trace formula (B14) equals the expectation value of the momentum loop ansatz (12), (20), and (24) in the big Euler ensemble.

Proof. As all the operators ${\hat{ν}}_{l}$ commute with each other, the operators $\hat{Q} (M)$ can be applied in arbitrary order to the states $Σ = | σ_{1}, \dots, σ_{N} ⟩$ involved in the trace. The same is true about individual terms in the circulation in the exponential of the Wilson loop. These terms ${\vec{F}}_{l}$ involve the operators ${\hat{α}}_{l}$ ⁠, which commute with each other and with each $\hat{Q} (M)$ ⁠. Thus, we can use the ordered product of the operators ${\hat{G}}_{l} = \hat{Q} (l) exp (ı ω {\hat{σ}}_{l} + \frac{i γ}{ν} Δ {\vec{C}}_{l} \cdot {\vec{P}}_{l} (t))$ ⁠. Each of the operators ${\hat{G}}_{l}$ acting in turn on arbitrary state Σ will create two terms with $δ [σ_{l} \pm 1]$ ⁠. The exponential in ${\hat{G}}_{l}$ will involve ${\hat{σ}}_{k}, k \leq l$ ⁠. As a result of the application of the operator ${\hat{Z}}_{l} = \prod_{k = 1}^{l} {\hat{G}}_{k}$ to the state vector Σ, we get $2^{l}$ terms with $Σ \prod_{k = 1}^{l} δ [σ_{k} - η_{k}], η_{k} = \pm 1$ ⁠. The factors ${\hat{Z}}_{l}$ will involve only ${\hat{σ}}_{k}, k \leq l$ ⁠, which are all reduced to $η_{k}, k \leq l$ in virtue of the product of the Kronecker deltas. Multiplying all operators ${\hat{G}}_{M}$ will lead to superposition ${\hat{Π}}_{N}$ of $2^{N}$ terms, each with product $\prod_{M = 1}^{N} δ [σ_{M} - η_{M}]$ with various choices of the signs $η_{i} = \pm 1$ for each i. Furthermore, the product of Kronecker deltas will project the total sum of $2^{N}$ combinations of the states Σ in the trace $tr \dots$ to a single term corresponding to a particular history $η_{1}, \dots, η_{N}$ of the Markov process. The product of Kronecker deltas in each history will be multiplied by the same state vector Σ, by the product of Markov transition probabilities, and by the exponential $exp (\frac{i γ}{ν} \sum_{l} Δ {\vec{C}}_{l} \cdot \hat{Ω} \cdot {\vec{P}}_{l} (t))$ with the operators $\hat{σ}$ in ${\vec{P}}_{k} (t)$ replaced by numbers η leading to the usual numeric ${\vec{P}}_{l} (t)$ ⁠. The transition probabilities of the Markov process are designed to reproduce combinatorial probabilities of random sigmas, adding up to one after summation over histories.⁴⁵ The integration over ω will produce $δ [\sum_{l} {\hat{η}}_{l} - s]$ ⁠. This delta function will reduce the trace to the required sum over all histories of the Markov process with a fixed $\sum_{l} η_{l}$ ⁠.

We have found a third vertex of the triangle of equivalent theories: the decaying turbulence in three-dimensional space, the fractal curve in complex space, and Fermi particles on a ring. By degrees of freedom, this is a one-dimensional Fermi-gas in the statistical limit $N \to \infty$ ⁠. However, there is no local Hamiltonian in this quantum partition function, just a trace of certain products of operators in Fock space. So, an algebraic (or quantum statistical) problem remains to find the continuum limit of this theory of the Fermion ring.

APPENDIX C: PATH INTEGRAL OVER MARKOV HISTORIES

Let us represent the product Π_N of the transitional probabilities of the particular history of the Markov processes as follows (with

n_{\pm} \equiv n_{\pm} (l), Δ n_{\pm} = - 1

⁠):

Π_{N} = exp (N Λ_{N}),

(C1)

Λ_{N} = \frac{1}{N} \sum_{l} G_{l},

(C2)

\begin{matrix} G_{l} = Δ n_{+} log (\frac{n_{+}}{n_{+} + n_{-}}) + Δ n_{-} log (\frac{n_{-}}{n_{+} + n_{-}}), \end{matrix}

(C3)

n_{+} = \sum_{k \leq l} ν_{k},

(C4)

n_{-} = \sum_{k \leq l} (1 - ν_{k}) .

(C5)

These

n_{\pm}

are net numbers of

η = \pm 1

in terms of Ising spins or occupation numbers

ν_{k} = (1, 0)

in the Fermi representation. There is an extra constraint on the Markov process

n_{+} + n_{-} = l; \forall l,

(C6)

which follows from the above definition in terms of the occupation numbers. We can redefine

n_{\pm}

as N times the piecewise constant functions as follows:

n_{\pm} = N f_{\pm} (ξ),

(C7)

f_{\pm} (ξ) = \sum_{k = 1}^{⌊ N ξ ⌋} \frac{ν_{k}}{N},

(C8)

f'_{\pm} (ξ) = \sum_{l = 1}^{N} δ (ξ - \frac{l}{N}) \sum_{k = 1}^{l} \frac{ν_{k}}{N},

(C9)

0 < ξ < 1.

(C10)

The sums can be rewritten as Lebesgue integrals

\begin{matrix} Λ_{N} = \int_{1}^{0} (d f_{+} (ξ) log (\frac{f_{+} (ξ)}{ξ}) + d f_{-} (ξ) log (\frac{f_{-} (ξ)}{ξ})) . \end{matrix}

(C11)

The sum over histories of the Markov process will become a path integral over the difference

ϕ = f_{+} (ξ) - f_{-} (ξ)

⁠, in that

\begin{matrix} \sum_{η_{.} = \pm 1} exp (N (Λ_{N} + ı Λ_{N}^{(1)})) \to \int D ϕ exp (N (Λ_{N} + ı Λ_{N}^{(1)})) . \end{matrix}

(C12)

This path integral will be dominated by the “classical history,” maximizing the product of transitional probabilities if such a classical trajectory exists.

The first term (without the circulation term) brings the variational problem

max_{ϕ} Λ_{N} [ϕ],

(C13)

\begin{matrix} Λ_{N} [ϕ] = \int_{1}^{0} d ξ (\frac{d f_{+}}{d ξ} log (\frac{f_{+}}{ξ}) + \frac{d f_{-}}{d ξ} log (\frac{f_{-}}{ξ})), \end{matrix}

(C14)

f_{\pm} (ξ) = \frac{1}{2} (ξ \pm ϕ (ξ)),

(C15)

f_{\pm} (ξ) \geq 0.

(C16)

This problem is, however, a degenerate one, as the functional reduces to the integral of the total derivative:

\frac{δ Λ_{N} [ϕ]}{δ ϕ (ξ)} = 0,

(C17)

\begin{matrix} Λ_{N} [ϕ] = \int_{1}^{0} d (f_{+} log f_{+} + f_{-} log f_{-}) + 1 + \int_{0}^{1} d ξ log ξ \\ = - \frac{1}{2} (1 - s) log (1 - s) - \frac{1}{2} (1 + s) log (1 + s) + log (2) . \end{matrix}

(C18)

It depends on the boundary condition $ϕ (0) = 0, ϕ (1) = s$ ⁠, but not on the shape of $ϕ (ξ)$ ⁠.

This expression matches the Stirling formula for the logarithm of the binomial coefficient in the combinatorial solution³¹ for the Euler ensemble

\begin{matrix} lim_{N \to \infty} \frac{log (\begin{matrix} N \\ N (1 + s) / 2 \end{matrix})}{N} = log (2) + - \frac{1}{2} (1 - s) log (1 - s) - \frac{1}{2} (1 + s) log (1 + s) . \end{matrix}

(C19)

This $Λ (s) = Λ_{\infty} (s)$ is a smooth even function of s taking positive values from $Λ (\pm 1) = 0$ to the maximal value $Λ (0) = log (2)$ (see Fig. 21).

FIG. 21.

The plot of the function Λ ( s ) . As required, it is positive, takes a maximal value log ( 2 ) at s = 0, and vanishes at both ends s = ± 1 of the physical region.

View large Download slide

The plot of the function $Λ (s)$ ⁠. As required, it is positive, takes a maximal value $log (2)$ at s = 0, and vanishes at both ends $s = \pm 1$ of the physical region.

Now, let us add the circulation term to the exponential of the partition function . This term can be directly expressed in terms of the difference between our two densities

N ϕ (ξ) = N f_{+} (ξ) - N f_{-} (ξ)

⁠:

ı N Λ_{N}^{(1)} [ϕ, C_{Ω}] = \frac{ı}{\sqrt{2 ν t}} \int_{0}^{1} d {\vec{C}}_{Ω} (ξ) \cdot \vec{F} (ξ),

(C20)

\vec{F} (ξ) = \frac{{sin (β N ϕ (ξ)), cos (β N ϕ (ξ)), 0}}{2 sin (β / 2)},

(C21)

{\vec{C}}_{Ω} (θ) = \hat{Ω} \cdot \vec{C} (θ) .

(C22)

We remember that the last component of the vector

\vec{F} (ξ)

does not contribute to the circulation integral in with the closed loop

{\vec{C}}_{Ω} (ξ)

⁠. This is why we replaced it with zero, not because it is small but because it drops. The key assumption is, of course, the existence of the smooth limit of the charge density

ϕ (ξ)

of these Fermions when they are densely covering this loop.

We are working with $α (ξ) = β N ϕ (ξ)$ in the following.

The measure for paths $[D α]$ is undetermined. The derivatives of these alphas were quantized in the original Fermi theory: each step $α' (ξ) \approx N Δ α = N β σ = \pm N β$ ⁠.

As we demonstrate below, in continuum theory, this discrete distribution can be replaced by a Gaussian distribution with the same mean square

\sum_{α^{'} = \pm N β} \leftrightarrow \int d α^{'} exp (- \frac{{(α^{'})}^{2}}{2 N^{2} β^{2}}) .

(C23)

To demonstrate that, we consider in the critical region

β^{2} \sim N^{- 1} \to 0

the most general term that arises in the moments of the circulation in (see Ref. 5 for some exact computations of these moments):

\begin{matrix} 2^{- N} \sum_{σ_{i} = \pm 1} exp (ı β \sum_{i} k_{i} σ_{i}) \\ = \prod_{i} cos β k_{i} \to exp (- β^{2} / 2 \sum_{i} k_{i}^{2}), \end{matrix}

(C24)

where k_i are some integers. With a large number N of these integers, the sum in the exponential becomes an integral, which is equivalent to a Gaussian integral,

\begin{matrix} exp (- β^{2} / 2 \sum_{i} k_{i}^{2}) \\ = \prod_{i} \int_{- \infty}^{\infty} \frac{d σ_{i}}{\sqrt{2 π}} exp (- σ_{i}^{2} / 2) exp (ı β k_{i} σ_{i}) \\ \to exp (- N β^{2} / 2 \int_{0}^{1} d ξ k {(ξ)}^{2}) . \end{matrix}

(C25)

We arrive at the standard path integral measure

\int [D α] = \int D α (ξ) exp (- \int_{0}^{1} d ξ \frac{{(α^{'})}^{2}}{2 N β^{2}}),

(C26)

\begin{matrix} \frac{\int [D α] exp (ı N \int_{0}^{1} d ξ α (ξ) K (ξ))}{\int [D α]} \\ = exp (- \frac{N^{2} \iint d ξ_{1} d ξ_{2} K (ξ_{1}) K (ξ_{2}) G_{1, 2}}{2}), \end{matrix}

(C27)

G_{1, 2} = ⟨ α (ξ_{1}) α (ξ_{2}) ⟩ .

(C28)

Appendix D will compute this Green's function

G_{1, 2} = G (ξ_{1}, ξ_{2})

⁠.

Thus, we arrive at the following path integral in the continuum limit:

(C29)

\begin{matrix} Ψ [C] = \frac{\sum_{p < q; (p, q)} \int_{Ω \in O (3)} d Ω \int [D α]}{\sum_{p < q; (p, q)} | O (3) | \int [D α]} \\ exp (ı \frac{\int_{0}^{1} d ξ Im (C_{Ω}^{'} (ξ) exp (ı α (ξ)))}{2 sin (π p / q) \sqrt{2 ν (t + t_{0})}}), \end{matrix}

(C29a)

C_{Ω} (θ) = \vec{C} (θ) \cdot \hat{Ω} \cdot {ı, 1, 0} .

(C29b)

We get the U(1) statistical model with the boundary condition $α (1) = α (0) + β N s$ ⁠. The period $β N s = 2 π p r$ is a multiple of $2 π$ ⁠, which is irrelevant at $N \to \infty$ ⁠. The effective potential for this theory is a linear function of the loop slope $\vec{C}' (ξ)$ ⁠.

This model is yet another representation of the Euler ensemble, suitable for the continuum limit.

APPENDIX D: MATCHING PATH INTEGRAL WITH COMBINATORIAL SUMS

The results of the path integration over α must match the combinatorial calculations with

σ_{l} = \pm 1

in the limit of large N. Without the interaction provided by the circulation term in , this path integral is dominated by a linear trajectory

α_{c l} (ξ) = β N ξ s .

(D1)

We already saw the match between the classical Action $Λ_{N} [ϕ (ξ) = ξ s]$ and the asymptotic value of the logarithm of the Binomial coefficient of the combinatorial solution for the sum over σ variables.

Let us verify some examples of the expectation values over σ. The simplest is (with

n \neq m

⁠)

{⟨ σ_{n} σ_{m} ⟩}_{\sum σ = N s} .

(D2)

The direct calculation using methods of Refs. 5 and 31 leads to

\begin{matrix} {⟨ σ_{n} σ_{m} ⟩}_{\sum σ = N s} = - \oint \frac{d ω}{2 π Z} e^{ı ω Ns} cos ((N - 2) ω) {sin}^{2} (ω) \\ = \frac{(1 - N s^{2})}{2^{N - 3} Z N (1 - s^{2})} (\begin{matrix} N - 2 \\ \frac{1}{2} (N + N s - 2) \end{matrix}), \end{matrix}

(D3)

Z = \oint \frac{d ω}{2 π} e^{ı ω Ns} cos (N ω) = 2^{- N} (\begin{matrix} N \\ \frac{1}{2} (s N + N) \end{matrix}) .

(D4)

Using Gamma function properties, this ratio simplifies to

\frac{N s^{2} - 1}{N - 1} .

(D5)

This result can be derived from symmetry without any integration.

\begin{matrix} {⟨ σ_{n} σ_{m} ⟩}_{\sum σ = N s} = A (N, s) (1 - δ_{n m}) + δ_{n m}, \end{matrix}

(D6)

\sum_{n, m} {⟨ σ_{n} σ_{m} ⟩}_{\sum σ = N s} = s^{2} N^{2} = A (N, s) N (N - 1) + N,

(D7)

A (N, s) = \frac{N s^{2} - 1}{N - 1} .

(D8)

The same limit

A (\infty, s) = s^{2}

follows from the classical trajectory

{⟨ σ_{n} σ_{m} ⟩}_{\sum σ = N s} \to \frac{α_{c l}^{'} (ξ)}{β N} \frac{α_{c l}^{'} (ξ')}{β N} = s^{2} .

(D9)

Let us consider less trivial example^5,31

U_{n, m} \to \sum_{k = n}^{m - 1} exp (ı α_{k, n}),

(D10)

α_{k, n} = β \sum_{\begin{matrix} l = 0 \\ l \neq n \end{matrix}}^{k} σ_{l} .

(D11)

We shall set s = 0, as this is the leading contribution to the partition function. The expectation value of

U_{n, m}

in our continuum limit becomes

\begin{matrix} ⟨ U_{n, m} ⟩ = N \int_{ξ_{1}}^{ξ_{2}} d ξ ⟨ exp (ı α (ξ)) ⟩ = N \int_{ξ_{1}}^{ξ_{2}} d ξ exp (- 1 / 2 G (ξ, ξ)) . \end{matrix}

(D12)

Here

G (ξ_{1}, ξ_{2})

is the Green's function corresponding to a 1D particle on a line interval

ξ \in (0, 1)

⁠, introduced in Appendix C. It satisfies the equation, which follows from our Gaussian Action,

\partial_{ξ}^{2} G (ξ, ξ') = - β^{2} N δ (ξ - ξ'),

(D13)

G (0, ξ') = G (ξ, 0) = 0.

(D14)

The solution is

G (ξ, ξ') = 1 / 2 β^{2} N (ξ + ξ' - | ξ - ξ' |) .

(D15)

Thus, we find

\begin{matrix} ⟨ U_{n, m} ⟩ = N \int_{ξ_{1}}^{ξ_{2}} d ξ exp (- 1 / 2 β^{2} N ξ) \\ = \frac{2}{β^{2}} (exp (- y / 2) - exp (- x / 2)), \end{matrix}

(D16)

x = β^{2} N ξ_{1},

(D17)

y = β^{2} N ξ_{2},

(D18)

in agreement with Refs. 5 and 31 in the critical region

N \to \infty, β^{2} \sim 1 / N

⁠. Finally, the expectation value

⟨ U_{n, m} {\bar{U}}_{n, m} ⟩ = \sum_{l = n}^{m - 1} \sum_{k = n}^{m - 1} ⟨ exp (ı α_{k n} - ı α_{l n}) ⟩ .

(D19)

Here, the Gaussian path integration yields

\begin{matrix} ⟨ U_{n, m} {\bar{U}}_{n, m} ⟩ \to N^{2} \int_{ξ_{1}}^{ξ_{2}} d ξ \int_{ξ_{1}}^{ξ_{2}} d ξ' \\ exp (- 1 / 2 (G (ξ, ξ) + G (ξ', ξ') - 2 G (ξ, ξ'))) \\ = N^{2} \int_{ξ_{1}}^{ξ_{2}} d ξ \int_{ξ_{1}}^{ξ_{2}} d ξ' exp (1 / 2 β^{2} N | ξ - ξ' |) \\ = \frac{4}{β^{4}} (2 exp ((y - x) / 2) + x - y - 2) . \end{matrix}

(D20)

This result also agrees with combinatorial computations in Refs. 5 and 31.

APPENDIX E: SMALL EULER ENSEMBLE IN STATISTICAL LIMIT

The remaining problem is averaging over the variables $N, p, q, r$ of the small Euler ensemble.

The variable

s = \frac{q r}{N}

is distributed between

- 1, 1

with the binomial weight^5,31

(\begin{matrix} N \\ N (1 + s) / 2 \end{matrix})

peaked at s = 0. There is a finite term coming from r = 0 plus a continuum spectrum coming from large r

W (r) = {\begin{matrix} 1, & if r = 0, \\ \frac{\sqrt{2 π N}}{q} exp (- \frac{{(q r)}^{2}}{2 N}), & otherwise . \end{matrix}

(E1)

As it was conjectured in Ref. 31 and supported by rigorous estimates in Ref. 5, the r = 0 term dominates the sums, after which the variables y, x can be treated as continuous variables.

The variable p at fixed q has a discrete distribution

f_{p} (p | q) = \frac{\sum_{\begin{matrix} p = 1 \\ (p, q) \end{matrix}}^{q - 1} δ (p - n)}{φ (q)} .

(E2)

As we shall see, rather than p, we would need an asymptotic distribution of a scaling variable

X (p, q) = \frac{1}{q^{2}} {cot}^{2} (\frac{π p}{q}) .

(E3)

This distribution for X(p, q) at fixed

q \to \infty

can be found analytically, using newly established relations for the cotangent sums (see the Appendix in Ref. 31 and Mathematica^® notebook³⁴). Asymptotically, at large q, these relations read

\begin{matrix} ⟨ X^{n} ⟩ \equiv lim_{q \to \infty} \frac{\sum_{\begin{matrix} p = 1 \\ (p, q) \end{matrix}}^{q - 1} X {(p, q)}^{n}}{q} = δ_{n, 0} + \frac{2 π^{- 2 n} ζ (2 n)}{(2 n + 1) ζ (2 n + 1)} . \end{matrix}

(E4)

This relation can be transformed further as

\begin{matrix} ⟨ X^{n} ⟩ = {\begin{matrix} 1, & if n = 0, \\ \frac{2 \sum_{k = 1}^{\infty} φ (k) k^{- (2 n + 1)}}{(2 n + 1) π^{2 n}}, & if n > 0. \end{matrix} \end{matrix}

(E5)

The Mellin transform of these moments leads to the following singular distribution:

⟨ X^{n} ⟩ = \int_{0}^{\infty} f_{X} (X) d X X^{n},

(E6)

f_{X} (X) = (1 - α) δ (X) + π X \sqrt{X} Φ (⌊ \frac{1}{π \sqrt{X}} ⌋),

(E7)

\begin{matrix} α = π \int_{0}^{\infty} X \sqrt{X} d X Φ (⌊ \frac{1}{π \sqrt{X}} ⌋) \\ = \frac{2 π}{5} \sum_{1}^{\infty} Φ (k) (\frac{1}{{(π k)}^{5}} - \frac{1}{{(π (k + 1))}^{5}}) \\ = \frac{2}{5 π^{4}} \sum_{1}^{\infty} \frac{φ (k)}{k^{5}} = \frac{1}{225 ζ (5)}, \end{matrix}

(E8)

where

Φ (n)

is the totient summatory function, defined as

Φ (q) = \sum_{n = 1}^{q} φ (n) .

(E9)

The distribution can also be rewritten as an infinite sum

\begin{matrix} \int d x f_{X} (x) F (x) = (1 - α) F (0) + π \sum_{n = 1}^{\infty} φ (n) \int_{0}^{\frac{1}{π^{2} n^{2}}} x^{\frac{3}{2}} dxF (x) . \end{matrix}

(E10)

The normalization of this distribution comes out 1 as it should, with factor $1 - α$ in front of the delta function.

The upper limit of X is

X_{max} = X (q - 1, q) \to \frac{1}{π^{2}} .

(E11)

Our distribution (E6) is consistent with this upper limit, as the argument $⌊ \frac{1}{π \sqrt{X}} ⌋$ becomes zero at $X π^{2} > 1$ ⁠. It is plotted in Fig. 22.

FIG. 22.

The log-log-plot of the distribution f X ( X ) . It is equals π Φ ( k ) X X at each internal 1 π 2 ( k + 1 ) 2 < X < 1 π 2 k 2 with positive integer k. Asymptotically f X ( X ) → 3 X π 3 at X → 0 .

View large Download slide

The log-log-plot of the distribution $f_{X} (X)$ ⁠. It is equals $π Φ (k) X \sqrt{X}$ at each internal $\frac{1}{π^{2} {(k + 1)}^{2}} < X < \frac{1}{π^{2} k^{2}}$ with positive integer k. Asymptotically $f_{X} (X) \to \frac{3 \sqrt{X}}{π^{3}}$ at $X \to 0$ ⁠.

Once we are zooming into the tails of the p, q distribution, we also must recall that

P (q < y N) = \frac{Φ (⌊ N y ⌋)}{Φ (N)} \to y^{2},

(E12)

f_{y} (y) = \frac{\sum_{q = 2}^{N} δ (y - \frac{q}{N}) φ (q)}{Φ (N)} .

(E13)

APPENDIX F: THE VORTICITY CORRELATION IN FOURIER SPACE

Let us outline an analytical solution. We shift the time variable by $t + t_{0} \Rightarrow t$ to simplify the formulas.

The correlation function reduces to the following average over the big Euler ensemble

E

of our random curves in complex space:³¹

\begin{matrix} ⟨ \vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{r}) ⟩ = \int_{O (3)} \frac{d Ω}{4 t^{2} | O (3) |} \sum_{0 \leq n < m < N} {⟨ {\vec{ω}}_{m} \cdot {\vec{ω}}_{n} e^{ı \vec{ρ} \cdot \hat{Ω} \cdot ({\vec{S}}_{n, m} - {\vec{S}}_{m, n})} ⟩}_{E}, \end{matrix}

(F1)

{\vec{S}}_{n, m} = \frac{\sum_{k = n}^{m - 1} {\vec{F}}_{k}}{m - n (mod N)},

(F2)

{\vec{ω}}_{k} = {0, 0, \frac{ı σ_{k}}{2} cot (\frac{β}{2})},

(F3)

\vec{ρ} = \frac{\vec{r}}{2 \sqrt{ν t}},

(F4)

{⟨ X [σ \dots] ⟩}_{p, q, r} \equiv \frac{\sum_{E} X [σ \dots] δ [q r - \sum σ]}{\sum_{E} δ [q r - \sum σ]} .

(F5)

Integrating the global rotation matrix O(3) is part of the ensemble averaging.

Let us apply our path integral to the expectation value over spins

σ = \pm 1

in the big Euler ensemble, with the distribution of q, X established in Appendix E. In the continuum limit, we replace summation with integration. We arrive at the following expression for the correlation function:

\begin{matrix} ⟨ \vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{r}) ⟩ \\ \propto \sum_{even q < N} \sum_{p; (p | q)} \frac{{cot}^{2} (π p / q)}{{(p / q)}^{2}} \int_{0 < ξ_{1} < ξ_{2} < 1} d ξ_{1} d ξ_{2} \int_{O (3)} d Ω \frac{\int [D α] α^{'} (ξ_{1}) α^{'} (ξ_{2}) e^{ı \frac{\vec{r} \cdot \hat{Ω} \cdot Im \vec{V} (ξ_{1}, ξ_{2})}{\sqrt{ν t}}}}{t^{2} Φ (N) | O (3) | \int [D α]}, \end{matrix}

(F6)

\begin{matrix} \vec{V} (ξ_{1}, ξ_{2}) = q \sqrt{X} {ı, 1, 0} (S (ξ_{1}, ξ_{2}) - S (ξ_{2}, 1 + ξ_{1})), \end{matrix}

(F7)

S (a, b) = \frac{\int_{a}^{b} d ξ e^{ı α (ξ)}}{b - a} .

(F8)

Here and in the following, we skip all positive constant factors, including powers of N. Ultimately, we restore the correct normalization of the vorticity correlation using its value at $\vec{r} = 0$ computed in previous work.³¹

The computations significantly simplify in Fourier space.

\begin{matrix} ⟨ \vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k}) ⟩ = \int d^{3} \vec{r} ⟨ \vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{r}) ⟩ exp (- ı \vec{k} \cdot \vec{r}) \\ \propto \sum_{even q < N} \sum_{p; (p | q)} {cot}^{2} (p / q) / {(p / q)}^{2} \int_{O (3)} d Ω \int_{0 < ξ_{1} < ξ_{2} < 1} d ξ_{1} d ξ_{2} \\ \frac{\int [D α] α^{'} (ξ_{1}) α^{'} (ξ_{2}) δ (\frac{\vec{r} \cdot \hat{Ω} \cdot Im \vec{V} (ξ_{1}, ξ_{2})}{\sqrt{ν t}} - \vec{k})}{t^{2} Φ (N) | O (3) | \int [D α]} . \end{matrix}

(F9)

The angular integration

\int d Ω

yields

\begin{matrix} \int_{O (3)} d Ω δ (\frac{\hat{Ω} \cdot Im \vec{V} (ξ_{1}, ξ_{2})}{\sqrt{ν t}} - \vec{k}) \\ \propto \frac{\sqrt{ν t}}{{\vec{k}}^{2}} δ (| Im \vec{V} (ξ_{1}, ξ_{2}) | - | \vec{k} | \sqrt{ν t}) . \end{matrix}

(F10)

Now, using the Lagrange multiplier λ for this condition, and shifting integration over λ to the real axis, we have to minimize effective action

\begin{matrix} A [α, λ] = \frac{π y^{2} X}{2} \int_{ξ_{1}}^{1 + ξ_{1}} {(α^{'})}^{2} + λ y \sqrt{X} | \frac{\int_{ξ_{1}}^{ξ_{2}} d ξ e^{ı α}}{ξ_{2} - ξ_{1}} - \frac{\int_{ξ_{2}}^{1 + ξ_{1}} d ξ e^{ı α}}{1 + ξ_{1} - ξ_{2}} |, \end{matrix}

(F11)

\partial_{λ} A [α, λ] = | \vec{k} | \sqrt{ν t} .

(F12)

This variational problem reduces to two pendulum equations:

α_{1}^{″} + \frac{r}{ξ_{2} - ξ_{1}} sin α_{1} = 0; \forall ξ_{1} < ξ < ξ_{2},

(F13)

α_{2}^{″} + \frac{r}{ξ_{2} - ξ_{1} - 1} sin α_{2} = 0; \forall ξ_{2} < ξ < 1 + ξ_{1},

(F14)

r = \frac{λ}{π y \sqrt{X} I (r)},

(F15)

I (r) = | \frac{\int_{ξ_{1}}^{ξ_{2}} d ξ e^{ı α_{1}}}{ξ_{2} - ξ_{1}} - \frac{\int_{ξ_{2}}^{1 + ξ_{1}} d ξ e^{ı α_{2}}}{1 + ξ_{1} - ξ_{2}} | .

(F16)

The well-known solution is Jacobi amplitude

am (x | u)

⁠,

\begin{matrix} α_{1} (ξ) = 2 am (\frac{ξ - α_{2}}{2} a_{1} | \frac{r}{a_{1}^{2} (ξ_{2} - ξ_{1})}), \end{matrix}

(F17a)

\begin{matrix} α_{2} (ξ) = 2 am (\frac{ξ - α_{2}}{2} a_{2} | - \frac{r}{a_{2}^{2} (1 - ξ_{2} + ξ_{1})}) . \end{matrix}

(F17b)

The free parameters

a_{1}, a_{2}, α_{1}, α_{2}

satisfy the following four equations:

α_{1} (ξ_{2}) = α_{2} (ξ_{2}),

(F18a)

α'_{1} (ξ_{2}) = α'_{2} (ξ_{2}),

(F18b)

α_{1} (ξ_{1}) = α_{2} (1 + ξ_{1}),

(F18c)

α'_{1} (ξ_{1}) = α'_{2} (1 + ξ_{1}),

(F18d)

together with the constraint following from the variation of the Lagrange multiplier λ:

I (r) = \frac{| \vec{k} | \sqrt{ν t}}{y \sqrt{X}} .

(F19)

APPENDIX G: TURBULENT VISCOSITY AND THE LOCAL LIMIT

These five equations, in general, are quite complex, but there is one simplifying property. In the local limit

N \to \infty

⁠, the remaining effective action at the extremum

\begin{matrix} N A [α_{c}, λ_{c}] = \frac{π N y^{2} X}{2} (\int_{ξ_{1}}^{ξ_{2}} d ξ {(α_{1}^{'} (ξ))}^{2} + \int_{ξ_{2}}^{1 + ξ_{1}} d ξ {(α_{2}^{'} (ξ))}^{2}) \end{matrix}

(G1)

grows as N, unless both

α_{1} (ξ) \sim α_{2} (ξ) \sim N^{- 1 / 2} \to 0

⁠. In this case, the above constraint can be expanded in

α_{1}, α_{2}

⁠. As we show in Ref. 39, the leading constant and linear terms both cancel so that the quadratic term remains

\begin{matrix} \frac{2 | \vec{k} | \sqrt{ν t}}{y \sqrt{X}} = | \int_{ξ_{1}}^{ξ_{2}} \frac{d ξ α_{1}^{2} (ξ)}{ξ_{2} - ξ_{1}} - \int_{ξ_{2}}^{1 + ξ_{1}} \frac{d ξ α_{2}^{2} (ξ)}{1 + ξ_{1} - ξ_{2}} | \sim \frac{1}{N} . \end{matrix}

(G2)

This estimate then requires vanishing viscosity in the local limit, at fixed turbulent viscosity

\tilde{ν} = ν N^{2} \to const .

(G3)

This phenomenon of renormalization of viscosity by a factor of N² was already observed in our first paper.³¹ Our Euler ensemble in the local limit $N \to \infty$ can only solve the inviscid limit of the Navier–Stokes decaying turbulence, with finite $\tilde{ν}$ acting as a turbulent viscosity.

The desired anomalous dissipation phenomenon takes place in this limit of our theory.

APPENDIX H: LINEARIZED CLASSICAL TRAJECTORY

Returning to the elliptic function solution, we rewrite it in the linearized case at

a_{1} \sim a_{2} \to 0

⁠. This linearization is equivalent to replacing

sin (α) \to α

in the differential equation and studying the resulting linear ODE as a boundary problem. We choose different parametrizations in this linear case

\begin{matrix} α_{1} (ξ) = a (cos (K_{1} (ξ - ξ_{2})) + \frac{b}{K_{1}} sin (K_{1} (ξ - ξ_{2}))), \end{matrix}

(H1)

\begin{matrix} α_{2} (ξ) = a (cos (K_{2} (ξ - ξ_{2})) + \frac{b}{K_{2}} sin (K_{2} (ξ - ξ_{2}))), \end{matrix}

(H2)

K_{1} = \sqrt{\frac{r}{Δ}},

(H3)

K_{2} = \sqrt{\frac{r}{Δ - 1}},

(H4)

Δ = ξ_{2} - ξ_{1} .

(H5)

In the physical region

0 < Δ < 1, r < 0

⁠, K₂ is real, and K₁ imaginary, but the solution stays real. The matching conditions at

α_{1} (ξ_{2}) = α_{2} (ξ_{2}), α'_{1} (ξ_{2}) = α'_{2} (ξ_{2})

are identically satisfied with this Anzatz. The derivative match

α'_{1} (ξ_{1}) = α'_{2} (1 + ξ_{1})

can be solved exactly for b as

b = \frac{P}{Q},

(H6)

\begin{matrix} P = \sqrt{\frac{r}{Δ - 1}} sin ((1 - Δ) \sqrt{\frac{r}{Δ - 1}}) + \sqrt{\frac{r}{Δ}} sin (Δ \sqrt{\frac{r}{Δ}}), \end{matrix}

(H7)

Q = cos ((Δ - 1) \sqrt{\frac{r}{Δ - 1}}) - cos (Δ \sqrt{\frac{r}{Δ}}) .

(H8)

The remaining matching condition

α_{1} (ξ_{1}) = α_{2} (1 + ξ_{1})

reduces to the root of the function

\begin{matrix} g (r, Δ) = \frac{(2 Δ - 1) sin (\sqrt{(Δ - 1) r}) sin (\sqrt{Δ r})}{\sqrt{(Δ - 1) Δ}} \\ + 2 cos (\sqrt{(Δ - 1) r}) cos (\sqrt{Δ r}) - 2 . \end{matrix}

(H9)

This function has multiple roots, but we are looking for the real root

r_{0} (Δ)

with minimal value of the action at given Δ

A_{c} (r, Δ) = \int_{ξ_{1}}^{ξ_{2}} d ξ α'_{1} {(ξ)}^{2} + \int_{ξ_{2}}^{1 + ξ_{1}} d ξ α'_{2} {(ξ)}^{2} .

(H10)

This integral is elementary, but the expression is too long to be presented here. It can be found in the Mathematica^® notebook,³⁹ where it is used to select the roots

r_{0} (Δ)

of

g (r, Δ)

⁠, minimizing

A_{c} (r_{0} (Δ), Δ)

for a given value of

Δ \in (0, 1)

⁠. All the solutions for Action for different roots of

r_{0} (Δ)

are shown in Fig. 23.

FIG. 23.

View large Download slide

Plot of $A_{c} (r_{0} (Δ), Δ)$ for the three real solutions of the equation $g (r, Δ) = 0$ ⁠. At $Δ = Δ_{1}$ and $Δ = Δ_{2}$ ⁠, the action curves intersect; at $Δ = 1 / 2$ ⁠, there is a gap between the lowest action (⁠ $= 0$ ⁠) and the lowest of the other two. So, there are second-order phase transitions at $Δ_{1}, Δ_{2}$ and the first-order phase transition at $Δ = 1 / 2$ ⁠.

This lowest action root is plotted in Fig. 24. The corresponding value of minimal action $L (Δ) = A_{c} (r_{0} (Δ), Δ)$ is plotted in Fig. 25.

FIG. 24.

View large Download slide

Log plot of $r_{0} (Δ)$ ⁠.

FIG. 25.

View large Download slide

Log plot of $L (Δ) = A_{c} (r_{0} (Δ), Δ)$ ⁠.

There are phase transitions at

Δ = {Δ_{1}, Δ_{2}, Δ_{3}},

(H11)

Δ_{1} = 0.157 143,

(H12)

Δ_{2} = 0.430 15,

(H13)

Δ_{3} = 1 / 2.

(H14)

These branch points in Δ correspond to the switch of the lowest action solution. At small positive

Δ - 1 / 2

⁠,

r_{0} (Δ) \to - 48 (Δ - 1 / 2) - \frac{1536 {(Δ - 1 / 2)}^{3}}{7} + \dots,

(H15)

A_{0} (r_{0} (Δ), Δ) \to 256 {(Δ - 1 / 2)}^{2} + \dots .

(H16)

At

Δ \to 1

⁠, all solutions go to

- \infty

as

r_{n} (Δ) \to - \frac{6}{1 - Δ} + O (1) .

(H17)

This behavior matches numerical computations in Mathematica^®.³⁹

The constraint (G2) is also reduced to elementary functions, too lengthy to quote here (see Ref. 39).

This constraint yields the quadratic relation for the last unknown parameter a in our solution

a^{2} = \frac{| \vec{k} | \sqrt{ν t}}{S (Δ) y \sqrt{X}}

(H18)

with universal function

S (Δ)

presented in Ref. 39 and shown in Fig. 26.

FIG. 26.

View large Download slide

Plot of $S (Δ)$ ⁠.

The resulting integral (up to the pre-exponential factor Q) is equal to

\begin{matrix} \int [D α] α^{'} (ξ_{1}) α^{'} (ξ_{2}) δ (| Im \vec{V} (ξ_{1}, ξ_{2}) | - | \vec{k} | \sqrt{ν t}) \\ \propto Q ⟨ α' (ξ_{1}) α' (ξ_{2}) ⟩ exp (- y \sqrt{X} | \vec{k} | \sqrt{\tilde{ν} t} \frac{π L (Δ)}{2 | S (Δ) |}) . \end{matrix}

(H19)

The factor

⟨ α' (ξ_{1}) α' (ξ_{2}) ⟩

contains two terms:

⟨ α' (ξ_{1}) α' (ξ_{2}) ⟩ = α_{1}^{'} (ξ_{1}) α_{1}^{'} (ξ_{2}) + ⟨ δ α' (ξ_{1}) δ α' (ξ_{2}) ⟩ .

(H20)

The first term is the contribution of the classical solution we have just found, and the second term comes from Gaussian fluctuations

δ α (ξ)

around this solution.

The classical term is calculated as (see Ref. 39)

α_{1}^{'} (ξ_{1}) α_{1}^{'} (ξ_{2}) = \frac{| \vec{k} | \sqrt{\tilde{ν} t}}{N y \sqrt{X}} \frac{J (Δ)}{| S (Δ) |},

(H21)

J (Δ) = {\frac{r A (r, Δ) B (r, Δ)}{Δ (Δ - 1) C {(r, Δ)}^{2}} |}_{r = r_{0} (Δ)},

(H22)

\begin{matrix} A (r, Δ) = Δ sin (\sqrt{(Δ - 1) r}) - \sqrt{(Δ - 1) Δ} sin (\sqrt{Δ r}), \end{matrix}

(H23)

\begin{matrix} B (r, Δ) = Δ sin (\sqrt{(Δ - 1) r}) cos (\sqrt{Δ r}) \\ - \sqrt{(Δ - 1) Δ} sin (\sqrt{Δ r}) cos (\sqrt{(Δ - 1) r}), \end{matrix}

(H24)

C (r, Δ) = cos (\sqrt{(Δ - 1) r}) - cos (\sqrt{Δ r})

(H25)

(see Fig. 27).

FIG. 27.

View large Download slide

Plot of universal function $J (Δ)$ ⁠. The four corves correspond to four phases [solutions for $r_{0} (Δ)$ ].

The fluctuation term $⟨ δ α' (ξ_{1}) δ α' (ξ_{2}) ⟩$ is also proportional to $1 / N$ ⁠; therefore, we must keep this term as well.

As for the pre-exponential factor Q in the saddle point integral, it is given by the functional determinant of the operator

\hat{L}

corresponding to linearized effective action in the vicinity of the saddle point

λ_{c}, α_{c} (ξ)

⁠.

\begin{matrix} A [α_{c} + δ α, λ_{c} + δ λ] \to A [α_{c}, λ_{c}] + 1 / 2 tr (V^{†} | \hat{L} | V), \end{matrix}

(H26)

V = {δ λ, δ α},

(H27)

Q (Δ, τ) = exp {(- 1 / 2 tr (\frac{log \hat{L}}{{\hat{L}}^{α}}))}_{α \to 0},

(H28)

τ = y \sqrt{X} .

(H29)

The fluctuation correction reduces to the inverse operator

\hat{L}

⁠, which we compute in Appendix I.

Now, we can reduce multiple sum/integral in to the following:

⟨ \vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k}) ⟩ = \frac{{\tilde{ν}}^{3 / 2} H (k \sqrt{\tilde{ν} t})}{\sqrt{t}},

(H30)

\begin{matrix} H (κ) = \frac{1}{Z} \sum_{n = 1}^{\infty} φ (n) \int_{0}^{1 / n} d τ (τ^{5} / n^{5} - τ^{10}) \int_{0}^{1} d Δ (1 - Δ) G (Δ, τ, κ), \end{matrix}

(H31)

\begin{matrix} G (Δ, τ, κ) = exp (- \frac{τ κ L (Δ)}{2 π | S (Δ) |}) \\ Q (Δ, τ) (\frac{κ}{N τ} \frac{J (Δ)}{| S (Δ) |} + ⟨ δ α' (ξ_{1}) δ α' (ξ_{2}) ⟩), \end{matrix}

(H32)

where

Z

is the normalization constant to be determined later.

APPENDIX I: FUNCTIONAL DETERMINANT IN THE PATH INTEGRAL

As we have discussed in Appendix H, in the limit $a \to 0$ the classical solution $α_{1, 2} (ξ) \propto a \to 0$ ⁠.

This observation simplifies the linearized theory corresponding to this quadratic form $⟨ V | \hat{L} | V ⟩$ ⁠. First, integrate the fluctuations $δ λ$ of λ around the saddle point solution.

The Lagrange multiplier at the saddle point vanishes, as we show in Ref. 39,

λ_{0} = τ r_{0} (Δ) I (r_{0} (Δ)) = 0 .

(I1)

The quadratic term comes from the first derivatives

I_{λ} = \partial_{λ} I, I_{r} = \partial_{r} I, λ_{r} = \partial_{r} λ

⁠, which can be simplified by switching to

λ (r) = τ t I (r)

as

A_{λ λ} = τ I_{λ} = \frac{τ I_{r}}{λ_{r}} = \frac{τ I_{r}}{τ r I_{r}} = \frac{1}{r} .

(I2)

The bilinear term

λ δ α

also simplifies

A_{α λ} (δ α) = ı τ δ λ F [δ α],

(I3)

F [δ α] = \int_{ξ_{1}}^{ξ_{2}} \frac{d ξ δ α (ξ)}{ξ_{2} - ξ_{1}} - \int_{ξ_{2}}^{1 + ξ_{1}} \frac{d ξ δ α (ξ)}{1 + ξ_{1} - ξ_{2}} .

(I4)

We can integrate out λ, producing the extra pre-exponential factor $Q_{λ} = \sqrt{| r_{0} (Δ) |} / \sqrt{N}$ ⁠.

The bilinear term in the exponential after λ integration leads to the following effective quadratic Action for

δ α

⁠:

A_{eff} [δ α] = \frac{τ^{2}}{2} \int_{ξ_{1}}^{1 + ξ_{1}} δ α_{2}^{'} + \frac{r_{0} (Δ) τ^{2}}{2} F {[δ α]}^{2} .

(I5)

There is a zero-mode

δ α (ξ) = const

⁠, related to translational invariance of

A_{eff} [δ α]

⁠. Naturally, this zero-mode must be eliminated from the spectrum when we compute the functional determinant and the resolvent below.

After discarding the zero-mode, this effective action becomes a positive definite functional of $δ α$ only in the region of Δ where $r_{0} (Δ) > 0$ ⁠, i.e., for $Δ_{1} < Δ < Δ_{2}$ ⁠.

As we shall see below, the spectrum of fluctuations is positive only in this region. Therefore, we restrict our integration to this region.

The

{(δ α)}^{2}

term corresponds to the linear eigenvalue equation with

f_{1, 2} = δ α_{1, 2}

⁠,

f_{1, 2}^{″} (ξ) - μ_{1, 2} F [f] = - ω f_{1, 2} (ξ),

(I6)

μ_{1, 2} = {\frac{r}{Δ}, \frac{r}{Δ - 1}},

(I7)

ϵ = ω τ^{2},

(I8)

r = r_{0} (Δ) .

(I9)

The solution matching with first derivative at

ξ = ξ_{2}, ξ = (ξ_{1}, 1 + ξ_{1})

is built the same way as in . Equations for

f_{1, 2}

being linear homogeneous, we can fix the normalization as

F [f] = 1

⁠,

\begin{matrix} f_{1} (x) = a sin (\sqrt{ω} (ξ_{-} ξ_{2})) + B_{1} cos (\sqrt{ω} (ξ_{-} ξ_{2})) + \frac{r}{ω Δ}, \end{matrix}

(I10)

\begin{matrix} f_{2} (ξ) = a sin (\sqrt{ω} (ξ_{-} ξ_{2})) + B_{2} cos (\sqrt{ω} (ξ_{-} ξ_{2})) + \frac{r}{ω (Δ - 1)} . \end{matrix}

(I11)

The spectrum

ω = ω_{n}

is defined by the transcendental equation (the discriminant of this linear system of equations), which we found in Ref. 39,

f (ω_{n}, Δ) = 0,

(I12)

\begin{matrix} f (ω, Δ) = (Δ - 1) Δ \sqrt{ω} sin (\frac{\sqrt{ω}}{2}) ((Δ - 1) Δ ω + r) \\ + r cos (\frac{1}{2} (1 - 2 Δ) \sqrt{ω}) - r cos (\frac{\sqrt{ω}}{2}), \end{matrix}

(I13)

r = r_{0} (Δ) .

(I14)

The spectrum is positive in the interval

Δ_{1} < Δ < Δ_{2}

where

r_{0} (Δ) > 0

so that the solution for

ω_{n} (Δ)

is stable. In the following, we only select the stable region with positive

r_{0} (Δ)

⁠.

The first levels of the spectrum satisfying equation $f (ω, Δ) = 0$ are shown in Fig. 28.

FIG. 28.

View large Download slide

The first levels of the spectrum satisfying equation $f (ω, Δ) = 0$ ⁠. The colored lines correspond to four phases: red: $0 < Δ < Δ_{1}$ ⁠, green: $Δ_{1} < Δ < Δ_{2}$ ⁠, blue: $Δ_{2} < Δ < 1 / 2$ ⁠, brown: $1 / 2 < Δ < 1$ ⁠. The green zone is left as stable, and others are eliminated because $r_{0} (Δ) < 0$ in these zones. Naturally, we eliminate the zero-mode $ϵ_{0} = 0$ corresponding to translational invariance of the effective Action.

The functional determinant, resulting from the WKB approximation to the α path integral, would be related to the infinite product of positive eigenvalues

ϵ_{n} = τ^{2} ω_{n}

⁠, which can be written using a contour integral

\begin{matrix} Q_{α} (Δ, τ) = \prod_{ω_{n} > 0} {(τ^{2} ω_{n})}^{- 1 / 2} \\ = exp (1 / 2 \partial_{x} I m \oint_{Γ} \frac{f' (ω)}{f (ω)} \frac{d ω}{2 π {(ω τ^{2})}^{x}}) | x \to 0, \end{matrix}

(I15)

and the integration contour Γ encircles anticlockwise the positive real poles of the meromorphic function

f' (ω) / f (ω)

⁠. The integral converges at

x > 1 / 2

and should be analytically continued to x = 0.

For this purpose, let us introduce another function

Φ (ω) = \frac{f (ω)}{cos (\sqrt{ω} / 2)} ω^{- 3 / 2} .

(I16)

We show in Ref. 39 that at large

ω = ı y

⁠, this function reaches finite limits

Φ (ı y) \to ı sign y {(Δ - 1)}^{2} Δ^{2} + \frac{(Δ - 1) Δ r}{| y |} .

(I17)

The logarithmic derivative of the original function differs from

\frac{Φ' (ω)}{Φ (ω)}

by the following meromorphic function:

\frac{f' (ω)}{f (ω)} - \frac{Φ' (ω)}{Φ (ω)} = - \frac{tan (\frac{\sqrt{ω}}{2})}{4 \sqrt{ω}} + \frac{3}{2 ω} .

(I18)

This difference produces a calculable contribution to our integral. By summing residues of the poles of the tangent, we get

\begin{matrix} \oint_{Γ} (- \frac{tan (\frac{\sqrt{ω}}{2})}{4 \sqrt{ω}} + \frac{3}{2 ω}) \frac{d ω}{2 π {(ω τ^{2})}^{x}} = ı (1 - 2^{2 x}) {(2 π τ)}^{- 2 x} ζ (2 x) . \end{matrix}

(I19)

The derivative at x = 0 yields a constant

1 / 2 \partial_{x} I m ı (1 - 2^{2 x}) {(2 π τ)}^{- 2 x} ζ (2 x) \to \frac{log (2)}{2},

(I20)

leading to an irrelevant renormalization of

Q_{α} (Δ, τ)

by a factor

\sqrt{2}

⁠.

The remaining integral with

f (ω) \Rightarrow Φ (ω)

already converges at

Re x > - 1

⁠, so that we can set x = 0 and rotate the integration contour Γ parallel to the imaginary axis at

Re Γ = ϵ > 0

⁠:

\begin{matrix} Q_{α} (Δ, τ) = exp (1 / 2 I m \int_{ϵ - ı \infty}^{ϵ + ı \infty} \frac{Φ' (ω)}{Φ (ω)} \frac{log (ω τ^{2}) d ω}{2 π}) . \end{matrix}

(I21)

The remarkable property of this functional determinant is the factorization of the τ dependence

Q_{α} (Δ, τ) = τ^{μ (Δ)} Q_{α} (Δ, 1),

(I22)

μ (Δ) = Im \int_{ϵ - ı \infty}^{ϵ + ı \infty} \frac{Φ' (ω)}{Φ (ω)} \frac{d ω}{2 π} .

(I23)

The index

μ (Δ)

has a topological origin and can be computed analytically as

μ (Δ) = \frac{arg Φ (ı \infty) - arg Φ (- ı \infty)}{2 π} = 1 / 2 .

(I24)

Our result for the correlation function is given by with

Q (Δ, τ) = Q_{α} (Δ, 1) τ^{1 / 2} \sqrt{r_{0} (Δ)}

(I25)

and

Q_{α} (Δ, 1)

given by . All the constant factors we have omitted here are absorbed by the normalization factor

Z

⁠, which we determine at the end of Appendix J.

APPENDIX J: THE FLUCTUATION TERM IN $α' (ξ_{1}) α' (ξ_{2})$

The last missing term is the fluctuation contribution to

α' (ξ_{1}) α' (ξ_{2})

⁠. In the Gaussian approximation, valid at

N \to \infty

⁠, this term equals

⟨ δ α' (ξ_{1}) δ α' (ξ_{2}) ⟩ = \frac{- 1}{N τ^{2}} {[\partial_{ξ} \partial_{ξ'} G (ξ, ξ')]}_{ξ = ξ_{1}, ξ' = ξ_{2}},

(J1)

where

G (ξ, ξ')

is a resolvent for the effective quadratic Action . This resolvent satisfies the equation

\partial_{ξ}^{2} G (ξ, ξ') - μ (ξ) F [G] = δ (ξ - ξ'),

(J2)

G (ξ_{1}, ξ') = G (ξ_{1} + 1, ξ') = 0,

(J3)

\begin{matrix} F [G] = \int_{ξ_{1}}^{ξ_{2}} \frac{d ξ G (ξ, ξ')}{ξ_{2} - ξ_{1}} - \int_{ξ_{2}}^{1 + ξ_{1}} \frac{d ξ G (ξ, ξ')}{1 + ξ_{1} - ξ_{2}}, \end{matrix}

(J4)

μ (ξ) = {\begin{matrix} \frac{r}{Δ}, & if ξ_{1} \leq ξ < ξ_{2}, \\ \frac{r}{Δ - 1}, & if ξ_{2} \leq ξ < 1 + ξ_{1} . \end{matrix}

(J5)

The solution of this equation, matching with the first derivative at

ξ = ξ_{2}

⁠, is

\begin{matrix} G {(ξ, ξ')}_{ξ_{1} \leq ξ < ξ_{2}} = A + \frac{| ξ - ξ' |}{2} + B (ξ - ξ_{2}) + \frac{F [G] r {(ξ - ξ_{2})}^{2}}{2 Δ}, \end{matrix}

(J6)

\begin{matrix} G {(ξ, ξ')}_{ξ_{2} \leq ξ < 1 + ξ_{1}} = A + \frac{| ξ - ξ' |}{2} + B (ξ - ξ_{2}) + \frac{F [G] r {(ξ - ξ_{2})}^{2}}{2 (Δ - 1)} . \end{matrix}

(J7)

The linear functional

F [G]

on this solution becomes a linear function of these unknown parameters A, B. Two boundary conditions

G (ξ_{1}, ξ') = G (ξ_{1} + 1, ξ') = 0

fix these parameters as functions of

ξ_{1}, ξ_{2}, ξ'

⁠.

The result derived in Ref. 39 is too lengthy to present here. The desired quantity is quite simple such that

\frac{1}{N τ^{2}} {[\partial_{ξ} \partial_{ξ'} G (ξ, ξ')]}_{ξ = ξ_{1}, ξ' = ξ_{2}} = \frac{2 (r - 6)}{(r + 12) N τ^{2}} .

(J8)

Finally, we get the following correlation (H30) (absorbing the constant factors in $Z$ ⁠):

(J9)

ν ⟨ \vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k}) ⟩ = \frac{{\tilde{ν}}^{3 / 2} H (k \sqrt{\tilde{ν} t})}{\sqrt{t}},

(J9a)

\begin{matrix} H (κ) = \int_{Δ_{1}}^{Δ_{2}} d Δ (1 - Δ) \sum_{n = 1}^{\infty} φ (n) \int_{0}^{1 / n} \frac{d τ}{τ^{\frac{5}{2}}} (τ^{5} / n^{5} - τ^{10}) G (Δ, τ κ), \end{matrix}

(J9b)

\begin{matrix} G (Δ, x) = \frac{Q_{α} (Δ, 1) \sqrt{r_{0} (Δ)}}{Z} (x \frac{J (Δ)}{S (Δ)} + \frac{2 (r_{0} (Δ) - 6)}{(12 + r_{0} (Δ))}) exp (- \frac{x L (Δ)}{2 π S (Δ)}) . \end{matrix}

(J9c)

APPENDIX K: MELLIN INTEGRAL FOR THE ENERGY SPECTRUM AND ENERGY DECAY

The spectral function

H (κ)

can be computed as follows. Let us represent the theta function as inverse Mellin transform (at

n > 0, τ > 0

⁠)

θ (1 - n τ) = \int_{+ 0 - ı \infty}^{+ 0 + ı \infty} \frac{d p}{2 π ı p} n^{- p} τ^{- p} .

(K1)

Substituting this representation into our integral in , and interchanging summations/integrations, we find

\begin{matrix} H (κ) = \int_{Δ_{1}}^{Δ_{2}} d Δ (1 - Δ) \int_{+ 0 - ı \infty}^{+ 0 + ı \infty} \frac{d p}{2 π ı p} \sum_{n = 1}^{\infty} \frac{φ (n)}{n^{p}} \\ \int_{0}^{\infty} d t (\frac{t^{5 / 2}}{n^{5}} - t^{15 / 2}) t^{- p} (A + B t κ) e^{- C t κ}, \end{matrix}

(K2)

A = \frac{Q_{α} (Δ, 1) \sqrt{r_{0} (Δ)}}{Z} \frac{2 (r_{0} (Δ) - 6)}{(12 + r_{0} (Δ))},

(K3)

B = \frac{Q_{α} (Δ, 1) \sqrt{r_{0} (Δ)}}{Z} \frac{J (Δ)}{S (Δ)},

(K4)

C = \frac{L (Δ)}{2 π S (Δ)} .

(K5)

The last integral reduces to Gamma functions and powers of n, after which the sum over n reduces to the ratio of two zeta functions

\begin{matrix} H (κ) = \int_{Δ_{1}}^{Δ_{2}} d Δ (1 - Δ) \int_{- 2 - ı \infty}^{- 2 + ı \infty} \frac{d z κ^{z}}{2 π ı} \\ \frac{20 C^{z - 1} (A C - B z) ζ (z + \frac{15}{2}) Γ (- z)}{(2 z + 7) (2 z + 17) ζ (z + \frac{17}{2})} . \end{matrix}

(K6)

According to the Riemann hypothesis, the function $ζ (w), w = z + \frac{17}{2}$ in the denominator has trivial zeros at negative even points $z = - 2, - 4, \dots$ and nontrivial zeros at the critical line $Re w = 1 / 2$ ⁠.

There is also a contribution from the poles at $z + \frac{15}{2} = 1, (2 z + 7) = 0, (2 z + 17) = 0$ ⁠.

We are not going to sum the residues in these poles, but rather integrate along the line $z = - 2 + ı y$ ⁠, where the integrand exponentially decreases as $exp (- \frac{π | y |}{2})$ in both directions.

There are oscillations on top of this exponential decrease. The integrals can be computed with arbitrary accuracy in Mathematica^® using the “DoubleExponentialOscillatory” integration method. Each integration takes a fraction of a second. The remaining integral over Δ was performed by computing the table of (positive) values of integrand $I (Δ)$ on a grid with step $δ Δ = \frac{Δ_{2} - Δ_{1}}{100}$ ⁠, and interpolating $log (I (Δ))$ between these values by fifth order polynomial $P_{5} (Δ)$ ⁠. Integral of the exponential of this polynomial $\int d Δ (1 - Δ) exp (P_{5} (Δ))$ yields good enough accuracy (floating precision $10^{- 8}$ ⁠).

This process was repeated for all M = 1000 values of $κ = 0.1 {(1000 / 0.1)}^{n / M}, n = 0, \dots, M$ ⁠, corresponding to equidistant log grid. This table was computed in parallel in Mathematica^® in Ref. 33, which took just a few minutes. Here is the resulting function $H (κ) κ^{5 / 3}$ (see Fig. 7). It shows strong deviations from K41 spectrum, in agreement with the DNS at Fig. 20.

The asymptotic behavior of this function is

H (κ) \propto κ^{- 7 / 2} .

(K7)

The effective index

μ (κ) = \frac{d log H (κ)}{d log κ}

slowly decreases reaching

μ (\infty) = - 7 / 2

⁠, as shown in Fig. 5.

The function $F (κ) = \int_{κ}^{\infty} x^{2} H (x) d x$ was obtained by numerical integration of exponential of interpolated function as before, with asymptotic tail $H (κ) \to const κ^{- 7 / 2}$ integrated analytically.

The log-log plot of the energy dissipation as a function of time is shown in Fig. 29, and it is curved in the log-log scale due to the quantum effects (complex poles of Mellin transform at Riemann ζ function zeros).

FIG. 29.

View large Download slide

The log-log plot energy dissipation $t^{2} E (t)$ as a function of t has a regime change due to quantum effects. It starts as constant at small t and asymptotically decays as $t^{2} E \propto t^{- 1 / 4}$ ⁠.

REFERENCES

1.

M. E.

Agishtein

and

A. A.

Migdal

, “

Simulations of four-dimensional simplicial quantum gravity as dynamical triangulation

,”

Mod. Phys. Lett. A

07

(

12

),

1039

–

1061

(

1992

).

https://doi.org/10.1142/S0217732392000938

Google Scholar

Crossref

2.

P. D.

Anderson

and

M.

Kruczenski

, “

Loop equations and bootstrap methods in the lattice

,”

Nucl. Phys. B

921

,

702

–

726

(

2017

).

https://doi.org/10.1016/j.nuclphysb.2017.06.009

Google Scholar

Crossref

3.

G.

Apolinario

,

L.

Moriconi

,

R.

Pereira

, and

V.

valadão

, “

Vortex gas modeling of turbulent circulation statistics

,”

Phys. Rev. E

102

,

041102

(

2020

).

https://doi.org/10.1103/PhysRevE.102.041102

Google Scholar

Crossref

PubMed

4.

A.

Ashtekar

, “

New variables for classical and quantum gravity

,”

Phys. Rev. Lett.

57

,

2244

–

2247

(

1986

).

https://doi.org/10.1103/PhysRevLett.57.2244

Google Scholar

Crossref

PubMed

5.

D.

Basak

and

A.

Zaharescu

, “

Connections between number theory and the theory of turbulence

,” (unpublished).

6.

E.

Brézin

and

V. A.

Kazakov

, “

Exactly solvable field theories of closed strings

,”

Phys. Lett. B

236

(

2

),

144

–

150

(

1990

).

https://doi.org/10.1016/0370-2693(90)90818-Q

Google Scholar

Crossref

7.

M.

Bulatov

and

A.

Migdal

, “

Numerical simulation of Euler ensemble

,” (unpublished).

8.

S.

Chandrasekhar

, “

On Heisenberg's elementary theory of turbulence

,”

Proc. R. Soc. London Ser. A. Math. Phys. Sci.

200

(

1060

),

20

–

33

(

1949

).

Google Scholar

9.

G.

Comte-Bellot

and

S.

Corrsin

, “

Simple Eulerian time correlation of full-and narrow-band velocity signals in grid-generated, ‘isotropic’ turbulence

,”

J. Fluid Mech.

48

(

2

),

273

–

337

(

1971

).

https://doi.org/10.1017/S0022112071001599

Google Scholar

Crossref

10.

G.

Comte-Bellot

and

S.

Corrsin

, “

The use of a contraction to improve the isotropy of grid-generated turbulence

,”

J. Fluid Mech.

25

(

4

),

657

–

682

(

1966

).

https://doi.org/10.1017/S0022112066000338

Google Scholar

Crossref

11.

M. R.

Douglas

and

S. H.

Shenker

, “

Strings in less than one dimension

,”

Nucl. Phys. B

335

(

3

),

635

–

654

(

1990

).

https://doi.org/10.1016/0550-3213(90)90522-F

Google Scholar

Crossref

12.

T. D.

Drivas

and

G. L.

Eyink

, “

An Onsager singularity theorem for Leray solutions of incompressible Navier–Stokes

,”

Nonlinearity

32

(

11

),

4465

(

2019

).

https://doi.org/10.1088/1361-6544/ab2f42

Google Scholar

Crossref

13.

R. P.

Feynman

,

R. B.

Leighton

, and

M.

Sands

, “

3-7 How did it get that way?

,” in

The Feynman Lectures on Physics

, Mainly Mechanics, Radiation, and Heat Vol.

1

, edited by

M. A.

Gottlieb

and

R.

Pfeiffer

(

Basic Books

,

2011

).

Google Scholar

14.

A. K.

Golmankhaneh

,

Fractal Calculus and Its Applications

(

World Scientific

,

2022

).

Google Scholar

Crossref

15.

J.-B.

Gorce

and

E.

Falcon

, “

Freely decaying Saffman turbulence experimentally generated by magnetic stirrers

,”

Phys. Rev. Lett.

132

,

264001

(

2024

).

https://doi.org/10.1103/PhysRevLett.132.264001

Google Scholar

Crossref

PubMed

16.

D. J.

Gross

and

A. A.

Migdal

, “

Nonperturbative two-dimensional quantum gravity

,”

Phys. Rev. Lett.

64

(

2

),

127

–

130

(

1990

).

https://doi.org/10.1103/PhysRevLett.64.127

Google Scholar

Crossref

PubMed

17.

G. H.

Hardy

and

E. M.

Wright

,

An Introduction to the Theory of Numbers

, 6th ed. (

Oxford University Press

,

Oxford

,

2008

). Revised by D. R. Heath-Brown and J. H. Silverman, with a foreword by Andrew Wiles.

Google Scholar

Crossref

18.

W.

Heisenberg

, “

Die bedeutung des schönen in der exakten naturwissenschaft

,”

Phys. Bl.

27

(

3

),

97

–

197

(

1971

).

https://doi.org/10.1002/phbl.19710270301

Google Scholar

Crossref

19.

K. P.

Iyer

,

S. S.

Bharadwaj

, and

K. R.

Sreenivasan

, “

The area rule for circulation in three-dimensional turbulence

,”

Proc. Natl. Acad. Sci. U. S. A.

118

(

43

),

e2114679118

(

2021

).

https://doi.org/10.1073/pnas.2114679118

Google Scholar

Crossref

PubMed

20.

K. P.

Iyer

,

K. R.

Sreenivasan

, and

P. K.

Yeung

, “

Circulation in high Reynolds number isotropic turbulence is a bifractal

,”

Phys. Rev. X

9

,

041006

(

2019

).

https://doi.org/10.1103/PhysRevX.9.041006

Google Scholar

Crossref

21.

V.

Kazakov

and

Z.

Zheng

, “

Bootstrap for lattice Yang-Mills theory

,”

Phys. Rev. D

107

(

5

),

L051501

(

2023

).

https://doi.org/10.1103/PhysRevD.107.L051501

Google Scholar

Crossref

22.

V.

Kazakov

and

Z.

Zheng

, “

Bootstrap for finite n lattice Yang-Mills theory

,” arXiv:2404.16925 (

2024

).

Google Scholar

23.

V. G.

Knizhnik

,

A. M.

Polyakov

, and

A. B.

Zamolodchikov

, “

Fractal structure of 2D–quantum gravity

,”

Mod. Phys. Lett. A

03

(

08

),

819

–

826

(

1988

).

https://doi.org/10.1142/S0217732388000982

Google Scholar

Crossref

24.

C.

Küchler

,

G. P.

Bewley

, and

E.

Bodenschatz

, “

Universal velocity statistics in decaying turbulence

,”

Phys. Rev. Lett.

131

,

024001

(

2023

).

https://doi.org/10.1103/PhysRevLett.131.024001

Google Scholar

Crossref

PubMed

25.

Y.

Makeenko

and

A. A.

Migdal

, “

Exact equation for the loop average in multicolor QCD

,”

Phys. Lett. B

88

(

1

),

135

–

137

(

1979

).

https://doi.org/10.1016/0370-2693(79)90131-X

Google Scholar

Crossref

26.

A.

Migdal

, “

Loop equations and

\frac{1}{N}

expansion

,”

Phys. Rep.

201

,

199

–

290

(

1983

).

Google Scholar

27.

A. A.

Migdal

, “

Momentum loop dynamics and random surfaces in QCD

,”

Nucl. Phys. B

265

(

4

),

594

–

614

(

1986

).

https://doi.org/10.1016/0550-3213(86)90331-7

Google Scholar

Crossref

28.

A. A.

Migdal

, “

Second quantization of the Wilson loop

,”

Nucl. Phys. B Proc. Suppl.

41

(

1

),

151

–

183

(

1995

).

https://doi.org/10.1016/0920-5632(95)00433-A

Google Scholar

Crossref

29.

A.

Migdal

, “

Loop equation and area law in turbulence

,” in

Quantum Field Theory and String Theory

, edited by

L.

Baulieu

,

V.

Dotsenko

,

V.

Kazakov

, and

P.

Windey

(

Springer US

,

1995

), pp.

193

–

231

.

Google Scholar

Crossref

30.

A.

Migdal

, “

Statistical equilibrium of circulating fluids

,”

Phys. Rep.

1011

,

1

–

117

(

2023

).

https://doi.org/10.1016/j.physrep.2023.02.001

Google Scholar

Crossref

31.

A.

Migdal

, “

To the theory of decaying turbulence

,”

Fractal Fract.

7

(

10

),

754

(

2023

).

https://doi.org/10.3390/fractalfract7100754

Google Scholar

Crossref

32.

A.

Migdal

, “

Topological vortexes, asymptotic freedom, and multifractals

,”

Fractals Fract.

7

,

351

(

2023

).

https://doi.org/10.3390/fractalfract7050351

Google Scholar

Crossref

33.

A.

Migdal

, “

abc

,” (

2024

), https://www.wolframcloud.com/obj/sasha.migdal/Published/ABCInterpolator.m.nb.

Google Scholar

34.

A.

Migdal

, “

bernsum

,” (

2024

), https://www.wolframcloud.com/obj/sasha.migdal/Published/BernSum.nb.

Google Scholar

35.

A.

Migdal

, “

corroscillations

,” (

2024

), https://www.wolframcloud.com/obj/sasha.migdal/Published/CorrelationOscillation.nb.

Google Scholar

36.

A.

Migdal

, “

dnsfit

,” (

2024

), https://www.wolframcloud.com/obj/sasha.migdal/Published/DataDNSFitFromTable.nb.

Google Scholar

37.

A.

Migdal

, “

Hierarchical structure of quantum solution

,” (

2024

), https://sashamigdal.github.io/QuantumSolution/.

Google Scholar

38.

A.

Migdal

, “

instantoncomputations

,” (

2024

), https://www.wolframcloud.com/obj/sasha.migdal/Published/DualTheory3.nb.

Google Scholar

39.

A.

Migdal

, “

instantoncomputations

,” (

2024

), https://www.wolframcloud.com/obj/sasha.migdal/Published/DualTheory1.nb.

Google Scholar

40.

A.

Migdal

, “

instantoncomputations

,” (

2024

), https://www.wolframcloud.com/obj/sasha.migdal/Published/DualTheory2.nb.

Google Scholar

41.

A. A.

Migdal

, “

Hidden symmetries of large N QCD

,”

Prog. Theor. Phys. Suppl.

131

,

269

–

307

(

1998

).

https://doi.org/10.1143/PTPS.131.269

Google Scholar

Crossref

42.

S.

Migdal

, “

Decaying turbulence experiment

,” (

2024

), https://drive.google.com/drive/folders/1e9piiWYvupbAv4gEss2 Mm5Kh-51MXJYm?q=sharedwith:public\%20parent:1e9piiWYvupbAv4gEss2 Mm5Kh-51MXJYm, accessed 16 May 2024.

Google Scholar

43.

P.

Moree

,

S.

Saad Eddin

,

A.

Sedunova

, and

Y.

Suzuki

, “

Jordan totient quotients

,”

J. Number Theory

209

,

147

–

166

(

2020

).

https://doi.org/10.1016/j.jnt.2019.08.014

Google Scholar

Crossref

44.

N. P.

Müller

,

J. I.

Polanco

, and

G.

Krstulovic

, “

Intermittency of velocity circulation in quantum turbulence

,”

Phys. Rev. X

11

,

011053

(

2021

).

https://doi.org/10.1103/PhysRevX.11.011053

Google Scholar

Crossref

45.

J. R.

Norris

,

Markov Chains

(

Cambridge University Press

,

2007

).

Google Scholar

46.

N.

Nekrasov

private communication (March 2024) suggested to me an algorithm of generating this solution as a set of adjacent triangles in complex three-space and pointed out an invariant measure in phase space, made of lengths of shared sides and angles between them. Unfortunately, this beautiful construction does not guarantee real circulation, requiring further study.

47.

K.

Ohkitani

, “

Study of the Hopf functional equation for turbulence: Duhamel principle and dynamical scaling

,”

Phys. Rev. E

101

(

1

),

013104

(

2020

).

https://doi.org/10.1103/PhysRevE.101.013104

Google Scholar

Crossref

PubMed

48.

J.

Panickacheril John

,

D. A.

Donzis

, and

K. R.

Sreenivasan

, “

Laws of turbulence decay from direct numerical simulations

,”

Philos. Trans. A Math. Phys. Eng. Sci.

380

(

2218

),

20210089

(

2022

).

https://doi.org/10.1098/rsta.2021.0089

Google Scholar

Crossref

PubMed

49.

G.

Parisi

and

U.

Frisch

, “

On the singularity structure of fully developed turbulence

,” in

Turbulence and predictability in Geophysical Fluid Dynamics

, edited by

R.

Benzi

,

M.

Ghil

, and

G.

Parisi

(

Scientific Research Publishing

,

North-Holland

,

1985

), pp.

84

–

88

.

Google Scholar

50.

J. I.

Polanco

,

N. P.

Müller

, and

G.

Krstulovic

, “

Vortex clustering, polarisation and circulation intermittency in classical and quantum turbulence

,”

Nat. Commun.

12

(

1

),

7090

(

2021

).

https://doi.org/10.1038/s41467-021-27382-6

Google Scholar

Crossref

PubMed

51.

A. M.

Polyakov

, “

Quantum geometry of bosonic strings

,”

Phys. Lett. B

103

,

207

–

210

(

1981

).

https://doi.org/10.1016/0370-2693(81)90743-7

Google Scholar

Crossref

52.

C.

Rovelli

and

L.

Smolin

, “

Knot theory and quantum gravity

,”

Phys. Rev. Lett.

61

,

1155

–

1158

(

1988

).

https://doi.org/10.1103/PhysRevLett.61.1155

Google Scholar

Crossref

PubMed

53.

M.

Sinhuber

,

E.

Bodenschatz

, and

G. P.

Bewley

, “

Decay of turbulence at high Reynolds numbers

,”

Phys. Rev. Lett.

114

,

034501

(

2015

).

https://doi.org/10.1103/PhysRevLett.114.034501

Google Scholar

Crossref

PubMed

54.

K. R.

Sreenivasan

, “

Chandrasekhar's fluid dynamics

,”

Annu. Rev. Fluid Mech.

51

,

1

–

24

(

2019

).

https://doi.org/10.1146/annurev-fluid-010518-040537

Google Scholar

Crossref

55.

K. R.

Sreenivasan

and

V.

Yakhot

, “

Dynamics of three-dimensional turbulence from Navier-Stokes equations

,”

Phys. Rev. Fluids

6

,

104604

(

2021

).

https://doi.org/10.1103/PhysRevFluids.6.104604

Google Scholar

Crossref

56.

K. R.

Sreenivasan

, “

Decaying turbulence

,” (

2023

), https://youtu.be/Ghr3-mZJ9I8?si=DrRdu4rDQkaUHomH, accessed 4 February 2024.

Google Scholar

57.

P. L.

Sulem

and

U.

Frisch

, “

Bounds on energy flux for finite energy turbulence

,”

J. Fluid Mech.

72

(

3

),

417

–

423

(

1975

).

https://doi.org/10.1017/S0022112075003059

Google Scholar

Crossref

58.

Wikipedia

, “

Burgers vortex

,” (

2022

), https://en.wikipedia.org/wiki/Burgers_vortex, accessed 27 April 2022.

59.

V.

Yakhot

and

K. R.

Sreenivasan

, “

Towards a dynamical theory of multifractals in turbulence

,”

Phys. A: Stat. Mech. Appl.

343

,

147

–

155

(

2004

).

https://doi.org/10.1016/j.physa.2004.07.037

Google Scholar

Crossref

60.

V.

Yakhot

and

V.

Zakharov

, “

Hidden conservation laws in hydrodynamics; energy and dissipation rate fluctuation spectra in strong turbulence

,”

Phys. D: Nonlinear Phenom.

64

(

4

),

379

–

394

(

1993

).

https://doi.org/10.1016/0167-2789(93)90050-B

Google Scholar

Crossref

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/).

Quantum solution of classical turbulence: Decaying energy spectrum

I. PROLOGUE

II. DEFINITIONS AND NOTATIONS

III. SUMMARY

IV. INTRODUCTION

A. Physical introduction: The energy flow and random vorticity structures

B. Mathematical introduction: The loop equation and its solution

C. The big and small Euler ensembles

V. THE FERMI RING AND ITS CONTINUUM LIMIT

VI. INSTANTON IN THE PATH INTEGRAL

VII. THE DECAYING ENERGY IN FINITE SYSTEM

A. Computation of the energy dissipation

B. The energy normalization problem

VIII. COMPARING OUR THEORY WITH THE DNS

IX. THE SPECTRUM OF SCALING DIMENSIONS

A. The energy spectrum

B. The energy decay

C. The vorticity correlation function in coordinate space

X. DISCUSSION

A. Myth and reality of turbulent scaling laws

B. Stochastic solution of the Navier–Stokes equation and ergodic hypothesis

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

D. Classical flow and quantum mechanics

E. Renormalizability of the inviscid limit of the loop equation

F. Relation of our solution to weak turbulence

XI. REMAINING PROBLEMS

XII. CONCLUSION

DEDICATION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: GLOBAL RANDOM ROTATION AND MOMENTUM LOOP SPACE

APPENDIX B: THE MARKOV CHAIN AND ITS FERMIONIC REPRESENTATION

APPENDIX C: PATH INTEGRAL OVER MARKOV HISTORIES

APPENDIX D: MATCHING PATH INTEGRAL WITH COMBINATORIAL SUMS

APPENDIX E: SMALL EULER ENSEMBLE IN STATISTICAL LIMIT

APPENDIX F: THE VORTICITY CORRELATION IN FOURIER SPACE

APPENDIX G: TURBULENT VISCOSITY AND THE LOCAL LIMIT

APPENDIX H: LINEARIZED CLASSICAL TRAJECTORY

APPENDIX I: FUNCTIONAL DETERMINANT IN THE PATH INTEGRAL

APPENDIX J: THE FLUCTUATION TERM IN α ′ ( ξ 1 ) α ′ ( ξ 2 )

APPENDIX K: MELLIN INTEGRAL FOR THE ENERGY SPECTRUM AND ENERGY DECAY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

APPENDIX J: THE FLUCTUATION TERM IN $α' (ξ_{1}) α' (ξ_{2})$