To the theory of decaying turbulence

We have found an infinite dimensional manifold of exact solutions of the Navier-Stokes loop equation for the Wilson loop in decaying Turbulence in arbitrary dimension $d>2$. This solution family is equivalent to a fractal curve in complex space $\mathbb C^d$ with random steps parametrized by $N$ Ising variables $\sigma_i=\pm 1$, in addition to a rational number $\frac{p}{q}$ and an integer winding number $r$, related by $\sum \sigma_i = q r$. This equivalence provides a dual theory describing a strong turbulent phase of the Navier-Stokes flow in $\mathbb R_d$ space as a random geometry in a different space, like ADS/CFT correspondence in gauge theory. From a mathematical point of view, this theory implements a stochastic solution of the unforced Navier-Stokes equations. For a theoretical physicist, this is a quantum statistical system with integer-valued parameters, satisfying some number theory constraints. Its long-range interaction leads to critical phenomena when its size $N \rightarrow \infty$ or its chemical potential $\mu \rightarrow 0$. The system with fixed $N$ has different asymptotics at odd and even $N\rightarrow \infty$, but the limit $\mu \rightarrow 0$ is well defined. The energy dissipation rate is analytically calculated as a function of $\mu$ using methods of number theory. It grows as $\nu/\mu^2$ in the continuum limit $\mu \rightarrow 0$, leading to anomalous dissipation at $\mu \propto \sqrt{\nu} \to 0$. The same method is used to compute all the local vorticity distribution, which has no continuum limit but is renormalizable in the sense that infinities can be absorbed into the redefinition of the parameters. The small perturbation of the fixed manifold satisfies the linear equation we solved in a general form. This perturbation decays as $t^{-\lambda}$, with a continuous spectrum of indexes $\lambda$ in the local limit $\mu \to 0$.


Introduction
We derived a functional equation for the so-called loop average or Wilson loop in turbulence a while ago.All the references to our previous works can be found in a recent review paper [1].
The path to an exact solution by a dimensional reduction in this equation was proposed in the 1993 paper but has just been explored.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.

Fluctuating geometry in gauge theory and turbulent flow
The situation has changed over the last decades.No alternative microscopic theory based on the Navier-Stokes equation emerged, but our understanding of the strong turbulence phenomena grew significantly.On the other hand, the loop equations technology in the gauge theory also advanced over the last decades.The correspondence between the loop space functionals and the original vector fields was better understood, and various solutions to the gauge loop equations were found.In particular, the momentum loop equation was developed, similar to our momentum loop used below [2,3].Recently, some numerical methods were found to solve loop equations beyond perturbation theory [4,5].The loop dynamics was extended to quantum gravity, where it was used to study nonperturbative phenomena [6,7].
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 [8][9][10][11] with Reynolds numbers of tens of thousands revealed and quantified violations of the K41 scaling laws.These numerical experiments are in agreement with so-called multifractal scaling laws [12].
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 four years ago [8], which triggered the further development of the geometric theory of turbulence [1,10,11].In particular, the Area law was justified for flat and quadratic minimal surfaces, and an exact scaling law in confinement region Γ ∝ √ Area was derived [1].The area law was verified with better precision in [9].

Sparse vorticity picture
We know (or at least assume) that the vortex structures in this extreme turbulent flow collapse into thin clusters in physical space.Snapshots of vorticity in numerical simulations [13,14] show a collection of tube-like structures relatively sparsely distributed in space.The large vorticity domains and the large strain domain lead to anomalous dissipation.These domains dominate the enstrophy integral so that the viscosity factor in front of this integral is compensated, leading to a finite dissipation rate.It was observed years ago and studied in the DNS [15].An excellent recent DNS [8,9] studied statistics of vorticity structures in isotropic turbulence with a high Reynolds number.Their distribution of velocity circulation appeared compatible with circular vortex structures (which we now call Kelvinons).
The authors of [8,9] confirmed the area law and compared it with the tensor area law, which would correspond to a constant uniform vortex, irrelevant to the turbulence.These two laws are indistinguishable for simple loops like a square, and they both are solutions to the loop equation but are very different for twisted or non-planar loops.Their data supported the area law, inspiring our search for the relevant vorticity structures behind it.Some recent works also modeled sparse vortex structures in classical [10] and quantum [11] turbulence.
We know such singular structures in Euler dynamics: vortex surfaces and lines.Vorticity collapses into a thin boundary layer around the surface, or the core surrounding the vortex line, moving in a self-generated velocity field.The vortex surface motion is known to be unstable against the Kelvin-Helmholtz instabilities, which undermines the whole idea of random vortex surfaces.However, the exact solutions of the Navier-Stokes equations discovered in the previous century by Burgers [16,17] and Townsend [18] show stable planar sheets with Gaussian profile of the vorticity in the normal direction, peaked at the plane.Thus, the viscosity effects in certain cases suppress the Kelvin-Helmholtz instabilities, leading to stable, steady vortex surfaces.The same applies to the Burgers vortex, corresponding to a cylindrical core, regularizing a singular vortex line.The Burgers cylindrical vortex has a constant strain, which is not the most general; there is only one independent eigenvalue instead of two.The general asymmetric solution approaches the symmetric Burgers vortex with the potential background flow as it was subsequently found in [19] in the turbulent limit |Γ| ≫ ν with asymmetric strain.Some singular (weak ) solutions of Euler equations with nontrivial topology are described by the so-called Clebsch field, defined modulo gauge transformations.These transformations manifestly preserve vorticity and, therefore, velocity.These variables and their ambiguity were known for centuries [20].They were first utilized within hydrodynamics in the pioneering work of Khalatnikov [21] and subsequent publications of Kuznetzov and Mikhailov [22], and Levich [23] in the early 80-ties.The modern mathematical formu-lation using symplectomorphisms was initiated in [24].Yakhot and Zakharov [25] derived the K41 spectrum in weak turbulence using some approximation to the kinetic equations in Clebsch variables.The topological singular Euler solutions (vortex sheets and lines) are similar to those studied theoretically and experimentally in superfluids [22,[26][27][28][29][30][31].
In a recent paper, [35], we studied the topological circular vortex in a classical turbulent flow.In particular, we found an effective Euler Hamiltonian, with extra anomalous terms coming from the line singularity resolved as a Burgers vortex in the local tangent direction.These anomalous terms have logarithmic dependence on the Reynolds number (circulation divided by viscosity).We summed up the leading log terms, which resulted in extra powers of the logarithm of the scale, modifying the K41 scaling laws.These exciting developments explain and quantitatively describe many interesting phenomena but do not provide a complete microscopic theory covering the full inertial range of turbulence without simplifying assumptions of the sparsity of vortex structures.In the present work, we develop the theory free of these assumptions and approximations by solving the loop equations for decaying turbulence.Our solution is irregular (local vorticity and local enstrophy limit do not exist), thus resembling the Tao conjecture [32].
There is an important difference, however.We are not studying singular solutions of the Navier-Stokes equations; rather, we are solving the Hopf equations for a generating functional for the distribution of velocity circulation.We are looking for a statistical ensemble of Navier-Stokes solutions, and we arrive at the probability distribution.This distribution is singular in that expectation values of powers of local vorticity diverge in a local limit.At the same time, the correlation functions have singularities at coinciding points.The correlation functions of vorticity at separated points are studied in the forthcoming paper by numerical simulation of our theory.These correlations have power singularities at coinciding points, showing new fractal dimensions.

Loop equation 2.1. Wilson loop as reduced Hopf functional
We introduced the loop equation in the Lecture Series at Cargese and Chernogolovka Summer Schools (see references in [1]).Here is a summary for the new generation.We write the Navier-Stokes equation as follows (with Einstein's notation of summation over repeated Greek indexes) (1) The Wilson loop average for the turbulence is defined as We added a dimensionless factor γ in the exponential compared to some previous definitions as an extra parameter of the Wilson loop.Without loss of generality, we shall assume that γ > 0. The negative γ corresponds to a complex conjugation of the Wilson loop.In Abelian gauge theory, this parameter would be the continuous electric charge.The statistical averaging ⟨. . .⟩ corresponds to initial randomized data, to be specified later.
Comparing the Wilson loop with the Hopf functional (see [33] for recent review) we observe that the Wilson loop is a particular case of the Hopf functional with the source ⃗ J(⃗ r) concentrated on a fixed loop in space Three comments are necessary here.

•
The Wilson loop is a dimensional reduction of the Hopf functional.Instead of a source ⃗ J(⃗ r) mapping R d → R d we have a source ⃗ C(θ) mapping a circle to physical space S 1 → R d .This mapping makes the loop equation a one-dimensional problem with physical space becoming a target space, like in string theory.The functional differential equations in more than one dimension are never solved except for the linear (or linearizable) ones.One-dimensional problems, even nonlinear ones, are sometimes solvable, and the loop equation is no exception, as we shall see.Rather than being an abstract notion, the Hopf equation becomes an analytical tool after this projection to one dimension.

•
We have lost most of the information about vector field ⃗ v(⃗ r) by reducing the Hopf source ⃗ J(⃗ r) to a one-dimensional subset ⃗ C of the physical space.This information would be irrecoverable if we only consider the set of smooth closed loops in space.However, with the loop functional defined for arbitrary non-smooth loop C, the information about velocity correlation functions and other statistical observables can be recovered [1] as follows.The small variation of the loop corresponding to attaching a closed loop at some point (which necessarily makes the loop non-smooth) brings down vorticity at this point.It allows to computation of the vorticity correlation functions.After that, velocity correlations can be recovered from incompressibility (the Biot-Savart integral in coordinate space or orthogonal projection in Fourier space).

•
The imaginary factor ıγ in the source is necessary for the mathematical existence of this functional.As the DNS shows ( [8,9,34] ), and the Kelvinon theory supports [35], the PDF of circulation decays only exponentially.The Fourier transform of such PDF exists but in the case of real instead of imaginary factor ıγ = γ R , the integral over the circulation distribution would diverge either on the positive or the negative circulation, depending upon the sign of γ R ; it should then be defined as an analytic continuation (Laplace transform).
With imaginary factor ıγ, this Wilson loop becomes the wave function of the quantum mechanical system in loop space.When treated as a function of time and a functional of the periodic function C : r α = C α (θ); θ ∈ (0, 2π), the Wilson loop satisfies the following functional Schrödinger equation (with operators ω, v defined below for reader's convenience)

H
(2) The area derivative δ δσ αβ is related to the variation of the functional when the little closed loop δC is added In the review, [1], we present the explicit limiting procedure needed to define these functional derivatives in terms of finite variations of the loop while keeping it closed.All the operators ∂ µ , ωαβ , vα are expressed in terms of the spike operator The area derivative operator can be regularized as − {α ↔ β}; (10) and velocity operator (with δ, ϵ → 0 + ) In addition to the loop equation, every valid loop functional F[C] must satisfy the Bianchi constraint [36,37] In three dimensions, it follows from identity ⃗ ∇ • ⃗ ω = 0; in general dimension d > 3, the dual vorticity ω is an antisymmetric tensor with d − 2 components.The divergence of this tensor equals zero identically.However, for the loop functional, this restriction is not an identity; it reflects that this functional is a function of a circulation of some vector field, averaged by some set of parameters.This constraint was analyzed in [1] in the confinement region of large loops, where it was used to predict the Area law.The area derivative of the area of some smooth surface inside a large loop reduces to a local normal vector.The Bianchi constraint is equivalent to the Plateau equation for a minimal surface (mean external curvature equals zero).
In the Navier-Stokes equation, we did NOT add artificial random forces, choosing instead to randomize the initial data for the velocity field.These ad hoc random forces would lead to the potential term [1] in the loop Hamiltonian H C , breaking the translational symmetry in the loop space needed for the dimensional reduction we study below.What is worse, the random forces pollute the genuine turbulent dynamics with long-range artifacts, violating the universality of turbulence as a critical phenomenon.With random initial data instead of time-dependent delta-correlated random forcing, we no longer describe the steady state (i.e., statistical equilibrium) but decaying turbulence, which is also an interesting process, manifesting the same critical phenomena.The energy is pumped in before the initial moment, at the interval −t 0 < t < 0 and slowly dissipates over time, provided the viscosity is small enough, corresponding to the large Reynolds number we are studying.

Plane wave in loop space
Our crucial observation in '93 was that the right side of the Loop equation, without random forcing, dramatically simplifies in functional Fourier space.The dynamics of the loop field can be reproduced in an Ansatz The difference with the original definition of Ψ[γ, C] is that our new vector function ⃗ P(θ) depends directly on θ rather then through the vector field ⃗ v(⃗ r) in R d projected at ⃗ r = ⃗ C(θ).This transformation manifests the dimensional reduction d ⇒ 1 we mentioned above.
From the point of view of the quantum analogy, this Anzatz is a plane wave in the Loop space, solving the Schrödinger equation in the absence of forces (potential terms in the Hamiltonian, breaking the translation invariance in the loop space).At the technical level, with the Anzatz (13) would pass through the loop equation provided the momentum loop ⃗ P(θ) does not depend on the space loop ⃗ C(θ) In this case, the variations of ( 13) bring down the product of momentum or, in a more general case, an arbitrary functional of the loop derivative reduces to that of the momentum loop This reduction will make the loop equation an algebraic rather than a functional differential equation.Instead of functional derivatives, we shall have ordinary functions.This equation appears solvable by an extra piece of luck.

Random global vorticity
Let us study the momentum loop dynamics on a simple example from original papers [1].The simplest meaningful distribution of the velocity field is the Gaussian one, with energy concentrated in the macroscopic motions.The corresponding loop field reads (we set γ = 1 for simplicity in this section) where f (⃗ r) is the velocity correlation function 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 The first term f 0 is proportional to initial energy density, and the second one is proportional to initial energy dissipation rate E 0 where d = 3 is the dimension of space.The constant term in (19) as well as r 2 + r ′2 terms drop from the closed loop integral, so we are left with the cross-term rr ′ , which reduces to a full square This distribution is almost Gaussian: it reduces to Gaussian one by extra integration The integration here involves all = 3 independent α < β components of the antisymmetric tensor ϕ αβ .Note that this is ordinary integration, not the functional one.
The physical meaning of this ϕ is the random uniform vorticity ω = f 1 φ at the initial moment.However, as we see it now, this initial data represents a spurious fixed point unrelated to the turbulence problem.It was discussed in our review paper [1].The uniform global rotation represents a fixed point of the Navier-Stokes equation for arbitrary uniform vorticity tensor.Gaussian integration by ϕ keeps it as a fixed point of the Loop equation.The right side of the Navier-Stokes equation vanishes at this special initial data so that the exact solution of the loop equation with this initial data equals its initial value (22).Naturally, the time derivative of the momentum loop with the corresponding initial data will vanish as well.
It is instructive to look at the momentum trajectory P α (θ) for this fixed point.The functional Fourier transform [1] leads to the following simple result for the initial values of P α (θ).In terms of Fourier harmonics, this initial data read The constant part P α,0 of P α (θ) is not defined, but it drops from equations by translational invariance.Note that this initial data is not real, as Pα,n ̸ = P ⋆ α,n .Positive and negative harmonics are real but unequal, leading to a complex Fourier transform.At fixed tensor ϕ the correlations are This correlation function immediately leads to the uniform expectation value of the vorticity The uniform constant vorticity kills the linear term of the Navier-Stokes equation in the original loop space, involving ∂ α Ωαβ = 0.The nonlinear term Vα Ωαβ vanishes in the coordinate loop space only after integration around the loop.Here are the steps involved The tensor area Σ was defined in (8).It is an antisymmetric tensor; therefore its trace with a symmetric tensor Ωαβ Ωβγ vanishes.
This calculation demonstrates how an arbitrary uniform vorticity tensor satisfies the loop equation in coordinate loop space.We expect the turbulent solution of the loop equation to be more general, with the local vorticity tensor at the loop becoming a random variable with a non-Gaussian distribution for every point on the loop.

Reduced Loop Equation
Let us use the Navier-Stokes loop equation to derive a dynamical equation for the momentum loop ⃗ P(t, θ).The reduced dynamics must be equivalent to the Navier-Stokes dynamics of the original field.With the loop calculus developed above, we have all the necessary tools to ensure this equivalence.Let us stress an important point: the function ⃗ P(θ, t) is independent of the loop C. As we shall see later, it is a random variable with a universal distribution in functional space.
This independence removes objection to the Kelvinon theory [35] and any other Navier-Stokes stationary solutions with a singularity at fixed loop C in space.Such solutions do not satisfy the loop equations because the functional derivatives would also act on the velocity field by varying its boundary conditions at the loop C.These variations would lead to extra terms to the Navier-Stokes equation for the circulation of the velocity field.
The loop equation for P(θ) as a function of θ and also a function of time, reads: where the operators V, D, Ω should be regarded as ordinary numbers, with the following definitions.The spike derivative D in the above equation The vorticity (10) and velocity (11) also become singular functionals of the trajectory P(θ).
The first observation about this equation is that the viscosity factor cancels after the substitution (35).As we shall see, the viscosity enters initial data so that at any finite time t, the solution for P still depends on viscosity.Another observation is that the spike derivative D(θ, ϵ) turns to the discontinuity ∆P(θ) = P(θ + ) − P(θ − ) in the limit ϵ → 0 The relation of the operators in the QCD loop equation to the discontinuities of the momentum loop was noticed, justified, and investigated in [3,38].The momentum loop in QCD could be piecewise constant with an arbitrary number of such discontinuities.In the Navier-Stokes theory, the discontinuities must be present at every point on a parametric circle S 1 .
We find the local limit for vorticity and velocity (skipping the common argument θ ) The Bianchi constraint is identically satisfied as it should We arrive at a singular loop equation for P α (θ) This equation is complex due to the irreversible dissipation effects in the Navier-Stokes equation.The viscosity dropped from the right side of this equation; it can be absorbed in units of time.Viscosity also enters the initial data, as we shall see in the next section on the example of the random rotation.However, the large-time asymptotic behavior of the solution would be universal, as it should be in the Turbulent flow.
We are looking for a degenerate fixed point [1], a fixed manifold with some internal degrees of freedom.The spontaneous stochastization corresponds to random values of these hidden internal parameters.Starting with different initial data, the trajectory ⃗ P(θ, t) would approach this fixed manifold at some arbitrary point and then keep moving around it, covering it with some probability measure.The Turbulence problem is to find this manifold and determine this probability measure.

Decay or fixed point
The absolute value of loop average Ψ[γ, C] stays below 1 at any time, which leaves two possible scenarios for its behavior at a large time.
FixedPoint : The Decay scenario in the nonlinear ODE (42) corresponds to the 1/ √ t decrease of ⃗ P. Omitting the common argument θ, we get the following exact time-dependent solution (not just asymptotically, at t → +∞).
The fixed point would correspond to the vanishing right side of the momentum loop equation (42).Multiplying by (∆ ⃗ P) 2 and reducing the terms, we find a singular algebraic equation The fixed point could mean self-sustained turbulence, which would be too good to be true, violating the second law of Thermodynamics.Indeed, it is easy to see that this fixed point cannot exist.The fixed point equation ( 47) is a linear relation between two vectors ⃗ P, ∆ ⃗ P with coefficients depending on various scalar products.The generic solution is simply with the complex parameter λ to be determined from the equation ( 47).This solution is degenerate: the fixed point equation is satisfied for arbitrary complex λ.The discontinuity vector ∆ ⃗ P aligned with the principal value ⃗ P corresponds to vanishing vorticity in (37), leading to a trivial solution of the loop equation Ψ[γ, C] = 1.We are left with the decaying turbulence scenario (46) as the only remaining physical solution.

Random walk
One may try the solution where the discontinuity vector is proportional to the principal value.However, in this case, such a solution does not exist.Assuming such a solution, we are led to the following contradictory equations There is, however, another solution where the vectors ∆ ⃗ F, ⃗ F are not aligned.This solution requires the following relations A simple algebra using these two relations shows that both sides of the decay equation ( 46) vanish: 0∆ ⃗ F = 0 ⃗ F. At the same time, vectors ∆ ⃗ F and ⃗ F will not be aligned, so the vorticity remains finite.Conversely, without these relations, the vectors ∆ ⃗ F and ⃗ F will be aligned, leading to vanishing vorticity in (37).
These relations are very interesting.The complex numbers reflect irreversibility, and lack of alignment leads to vorticity distributed along the loop.
Note that the system of equations ( 51) is parametrically invariant (being local and independent of θ).Also, note that this complex vector ⃗ F(θ) is dimensionless, and this system is completely universal, up to a single dimensionless parameter γ.The viscosity dropped from this equation, and the dimension of vorticity is inverse time, which explains the factor 1/(t + t 0 ).These equations do not allow ∆ ⃗ F = 0, so discontinuities must be present at every θ.In other words, our solution cannot be smooth: it is a fractal curve rather than a piecewise smooth curve with discontinuities at several points.
We approximate this fractal curve by a polygon, with piecewise constant ⃗ F(θ) and N gaps ∆ ⃗ F(θ k ) at equidistant angles θ k = 2πk/N.If the limit N → ∞ exists, we get the desired fractal curve ⃗ F(θ) with a discontinuity at every θ.One can build such a solution as a limit of the Markov process by the following method.Start with a complex vector

Constraints imposed on a random step
A solution to these equations can be represented using a complex vector ⃗ q k subject to two complex constraints after which we can find the next value We assume N steps, each with the angle shift ∆θ = 2π N .This recurrent sequence is a Markov process because each step only depends on the current position ⃗ F k .The closure requirement ⃗ F N = ⃗ F 0 makes it a periodic Markov process, with the period N.This requirement represents a nonlinear restriction on all the variables ⃗ F k , which we discuss below.With this discretization, the circulation can be expressed in terms of these steps Note that the complex unit vector is not defined with the Euclidean metric in 2d which leads to two conditions between real and imaginary parts In d dimensions, there are d − 1 complex parameters of the unit vector; with an extra constraint in (54b), there are now d − 2 free complex parameters at every step of our iteration, plus the discrete choice of the sign of the root in the solution of the quadratic equation (54b) for the projection ⃗ F k • ⃗ q k .

Closure condition
At the last step, k = N − 1, we need to get a closed loop ⃗ F N = ⃗ F 0 .This is one more constraint on the complex vectors ⃗ q 0 , . ..⃗ q N−1 (60) We use this complex vector constraint to fix some of the remaining parameters.Due to the closure of the space loop ⃗ C(θ), the global translation of the momentum loop ⃗ P(θ) leaves invariant the Wilson loop; this leads to certain gauge invariance (see below).The circulation must correspond to a real number, though the Wilson loop is not real, as there is an asymmetry in the distribution of signs of circulation.We discuss this issue in the following sections.

Mirror pairs of solutions
Return to the general study of the discrete loop equations (54).There is a trivial solution to these equations at any even We reject this solution as unphysical: the corresponding vorticity equals zero, as all the vectors ⃗ f k are aligned.Our set of equations has certain mirror reflection symmetry Thus, the complex solutions come in mirror pairs ⃗ F k , ⃗ F * N−k .The real solutions are only a particular case of the above trivial solution with real ⃗ q.Each nontrivial solution represents a periodic random walk in complex vector space C d .The complex unit step ⃗ q k ∈ C d depends on the current position ⃗ F k ∈ C d , or, equivalently, on the initial position ⃗ F 0 plus the sum of the preceding steps.We are interested in the limit of infinitely many steps N → ∞, corresponding to a closed fractal curve with a discontinuity at every point.

The degenerate fixed point and its statistical meaning
This solution's degeneracy (fewer restrictions than the number of free parameters) is a welcome feature.One would expect this from a fixed point of the Hopf equation for the probability distribution.This degeneracy leads to stochastization of the Navier-Stokes flow at large Reynolds numbers.The solution comes as a manifold, and the flow covers this manifold with some invariant measure.In the best-known example, the microcanonical Gibbs distribution for Newton's mechanics covers the energy surface with a uniform measure (ergodic hypothesis, widely accepted in Physics, though still unproven in Mathematics).The parameters describing a point on this energy surface are not specified-in the case of an ideal Maxwell gas, these are arbitrary velocities of particles.
Likewise, the fixed manifold, corresponding to our fractal curve, is parametrized by N → ∞ sign variables, like an Ising model, plus an arbitrary global rotation matrix Ô ∈ SO(d) and a global parameter β, as discussed in the next section.This rich internal random structure of our fixed manifold, combined with its rotation and translation invariance in loop space C, makes it an acceptable candidate for extreme isotropic turbulence.

Exact analytic solution 4.1. Random walk on a circle
How could a complex curve describe real circulation?This remarkable cancellation of the imaginary part of the circulation is possible if the imaginary part of ⃗ P(θ) does not depend on θ.Such an imaginary term will drop after integration over closed loop ⃗ C(θ).We have found a family of such solutions [39] of our recurrent equation (53) for arbitrary N Here ⃗ w ∈ S d−3 is a unit vector.The angles α k must satisfy recurrent relation This sequence with arbitrary signs σ k = ±1 solves recurrent equation (53) independently of γ.The closure condition requires certain relations between these numbers.The main condition is that β must be a rational fraction of 2π: The q = 2 case is eliminated.
Otherwise, the periodic solution for α k will correspond to the following set of σ k σ = {1, . . ., 1, −1, . . ., −1} perm ; (69) This array has N + positive values and N − negative values where (70) The symbol perm stands for a random permutation of the array, which preserves its sum This sum must be a multiple of q for periodicity, which leads to another restriction In other words, rq must have the same parity as N.
These properties lead to periodicity At fixed denominator q, the winding number r can take the values r = {−r max , . . . , 0,. . . ,r max }; (76) The sequence with all spins flipped: σ k ⇒ −σ k also solves the loop equation.This sequence is a reflected solution we mentioned, so we include it in the statistical samples with equal probability.It corresponds to the winding number's reflection r ⇒ −r.The number of states with given N, p, q, r is a partition of N + positive and N − negative spins into N boxes.The probabilities are given by binomial distribution with w = 1 /2 The next section discusses this ensemble of rational numbers and its statistics at N → ∞.Given the rational number p q , we can generate the sequence of angles ∑ ±β, adding to a 2π multiple.The random walk step ⃗ q k = ⃗ F k+1 − ⃗ F k is a real unit vector in this solution The direction of this vector is not random, though; in addition to the random sign σ k and random unit vector ⃗ w in d > 3 dimensions, its direction depends on the previous position α k on a circle.So, this is a perfect example of a periodic Markov chain, with the periodicity condition analytically solved by quantizing the angular step to a rational number β = 2π p q .This solution corresponds to the real value of velocity circulation on each of these two solutions; however, the reflection changes this value.Thus, the arithmetic average of two Wilson loops with two reflected solutions is reflection-symmetric, but it is still a complex number.
Our solution has a peculiar gauge invariance.The circulation and, therefore, all observables are invariant under the shift of all ⃗ F k by a constant vector: This gauge invariance follows from the closure of the loop C: any constant term in ⃗ F(θ) yields zero after integration d ⃗ C = 0, or summation ∑ ∆ ⃗ C = 0. Using this invariance, we can drop the last component of ⃗ F k so that they become real vectors The vorticity operator in this gauge will become a purely imaginary vector in the z direction: As we shall see, this does not lead to complex numbers for the correlation functions of vorticity in physical space.The correlation function of two vorticities, separated by a finite distance⃗ r in an "inertial interval," is finite and real after integration over the global rotation matrix.Its limit at⃗ r → 0 may be singular so that the anomalous dissipation may emerge.
This discrete set N, p, q, r, σ k describes a particular solution of the loop equation for the Wilson loop in decaying turbulence.Here is an important point to keep in mind.Unlike the Navier-Stokes equation, the loop equation is linear.Therefore, any superposition of its solutions parametrized by N, p, q, r, σ k also solves the loop equation.We need to find a particular superposition that has the correct physical properties.In particular, this superposition must describe a continuum limit of our fractal curve when N → ∞.In the next two sections, we study an ensemble of such solutions, the Euler ensemble.This ensemble corresponds to adding every solution with equal weight.Naturally, the question arises: Why equal weight?Why not give the even N more weight or exclude prime numbers?This is our ergodic hypothesis.Equal weight is the most symmetric option from the mathematical point of view, plus methods of numbers theory can study the properties of such an ensemble.One argument favoring the equal weight hypothesis is that weight distribution may become irrelevant in the local limit.This limit, as we shall see in the rest of the paper, is determined by the statistical weight of the configurations of the ensemble.Therefore, the singular behavior of the ensemble when the mean number ⟨N⟩ of vertices in the fractal curve goes to infinity may be universal.In other words, the continuum fractal curve corresponding to the limit of this solution may be a unique mathematical object independent of the method used to approximate it by a polygon.
In that case, there may be an alternative way to describe this fractal curve without taking a limit of a large polygon.This description would require more sophisticated mathematical methods than those we use in this paper.The advantage of our polygonal approach to the fractal curve is that it is well-defined and is solvable (the ensemble averages are calculable) at every finite N and can be analytically extrapolated to the ensemble with ⟨N⟩ → ∞.In terms of modern QFT, the continuum limit of this statistical ensemble is dual to the statistical theory of the velocity field.In the same way, the discrete model of well-known dynamical triangulations is dual to the continuum theory of quantum gravity.

The Euler ensemble
The discrete set of fractions p q with denominator q < N is well known in the number theory [40], starting with Euler.However, our extra restriction ∑ l σ l (mod q) = 0 ties the set of fractions to the set Z ⊗N 2 of N Ising variables.We will study the statistics of the corresponding ensemble, which we call the Euler ensemble, honoring great Leonard Euler.He never thought that his φ functions and his equation for ideal fluid would meet in theoretical physics, but good theories have a life of their own.
We distinguish between big Euler ensembles and small Euler ensembles.The big ensemble E(N) assigns equal weight for each element of the large set variables: p, q, r, {σ 1 . . .σ N } (84) −N ≤ qr ≤ N; (86) The small ensemble E (N) results from averaging the big ensemble over the σ variables variables: p, q, r, (89) −N ≤ qr ≤ N; (91) The binomial coefficients count the number of states with ∑ i σ i = qr among all 2 N states of the set of σ i = ±1, i = 1, . . .N. We divided the statistical weights (the number of allowed configurations) by the total number 2 N of spin configurations.This normalization makes w N (q, r) the probability of finding the values q, r in the big Euler ensemble with random σ i .
Let us count all fractions with denominator q < N and proper parity, same as N.All the integers between 2 and N are allowed for q, and each such number would enter with the weight ∑ r w N (q, r).At given q < N the allowed numbers of p are all integers 0 < p < q such that gcd(p, q) = 1.The famous Euler's totient φ(q) [40] counts such numbers.
Thus, we can relate the big Euler ensemble average of some function to the small ensemble average as follows ⟨F(. . .)⟩ E(N) = ⟨⟨F(. . .)⟩ σ ⟩ E (N) ; (97) The advantage of this representation is that we can first average by the σ variables, which is a rather simple arithmetic mean.After that, we have to average over the small Euler ensemble, which involves only three variables.Resulting triple sums are calculable numerically with Mathematica ® and also can be studied in the local limit N → ∞ by the methods of number theory.
The odd and even ensembles have different asymptotic behavior with N.Here are the allowed parity of q, r for odd/even N      odd r, odd q odd N integer r, even q even N even r, odd q even N (100) The ratios were computed for N = 1000, 1001 (see Fig. 1).These ratios are finite, up to some value of q, after which they quickly go to zero (faster than exponentially).This fast decrease suppresses these r ̸ = 0 terms in the sum over q, which otherwise diverges at the upper limit and grows as N 2 .At large N the r = 0 weight tends to the qindependent limit Log-log plots of four ratios R N (q) for even/odd 2 < q < N, even/odd N = 1000, 1001.The larger N led to astronomically small ratios, so we did not use them.which does not restrict the sum of Euler totients ∑ q φ(q).Therefore, for even For odd N, the leading term is missing, so we have to sum the terms with r ̸ = 0.
This term converges at q ≪ N where ∑ r>0 2w N (q, r) is not exponentially small.The asymptotic formulas for summators of the Euler totients do not apply here, so we can compute this sum numerically and extrapolate to N → ∞.We computed it numerically in Mathematica ® and fitted it to √ N times a power of log N (see Fig. 2).
It is a challenge to the number theory to produce exact asymptotic behavior, replacing my fitted "law."Recently, the number theory answered this challenge [41].The result of their computation confirms the √ N factor but rejects the log α (N) correction factor as an overfit due to insufficiently large N.The actual pre-exponential coefficient is required for the full theory of the Euler ensemble.However, it does not contribute to the leading term of the grand canonical ensemble we use in the rest of the paper.The results of [41] correspond to Z(N, N) odd N → √ N/ √ 2π without any factors of log N.This mistake reminds us that mathematical laws must be derived theoretically rather than fitted to numerical data.

Grand canonical ensemble
The original Euler ensemble with fixed N can be regarded as a microcanonical ensemble of statistical mechanics.The grand canonical ensemble is more appropriate if the number of states can fluctuate, as in our case, where probabilities drastically change when N is shifted by 1.The weights for varying numbers N of degrees of freedom are multiplied by e −µN , then N is treated as the rest of the variables.The chemical potential µ will have to tend to zero in the thermodynamic limit, and the resulting singularity becomes a thermodynamic singularity corresponding to the critical phenomena.Note that we changed the sign of µ compared to the historical definition: our µ → +0 so that the opposite sign would be inconvenient.
The partition function and the ensemble averages in the grand canonical ensemble are With the grand canonical ensemble, the ambiguity of the local limit disappears.At the critical point µ → 0, the even N dominate, with the following result The extra factor of 1/2 came from skipping all the odd values of N. The remaining sum over even N tends to be half the asymptotic expression's integral for these even N.In the rest of the paper, we shall use the grand canonical ensembles E(µ), E (µ) in the thermodynamic limit µ → 0, where we find the critical phenomena.
With the grand canonical ensemble, the ergodic hypothesis becomes a weaker restriction: any smooth weight function W(N) in the distribution will cancel in expectation values in the limit µ → 0. For example, the saddle point calculation yields This phenomenon is "self-organized criticality," unlike conventional second-order phase transitions in statistical physics, which require tuning of thermodynamical potentials (temperature, pressure, chemical potential, etc.).The critical point µ c = 0 does not require fine-tuning.The self-organized criticality, or spontaneous stochastization of turbulence, is an expected property that was never proven theoretically but observed in numerical and real experiments.In our theory, this spontaneous stochastization naturally emerges in the solution of the loop equations in the local limit.

Correlation functions
In this section and the rest of the paper, we only consider the three-dimensional space we inhabit.There are interesting mathematical problems related to decaying turbulence in higher dimensions, which we leave to future researchers.In less than three dimensions, our solutions do not exist.

General formulas
The simplest observable quantities we can extract from the loop functional are the vorticity correlation functions [1], corresponding to the loop C backtracking between two points in space⃗ r 1 = 0,⃗ r 2 = ⃗ r, see Fig. 3.The vorticity operators are inserted at these two points.
The correlation function reduces to the following average over the ensemble E(µ) of our random curves in complex space.
The averaging ⟨. . .⟩ in these formulas involves group integration O(3 Let us explain the origin of summation over two positions n, m of the points ⃗ r 1 = 0,⃗ r 2 = ⃗ r on the discreet loop ⃗ C(θ).There is a degree of freedom we did not specify until now, namely, the reparametrization of the momentum loop ⃗ P(θ).The loop equations (42) are invariant under the one-dimensional diffeomorphisms ( or reparametrizations) ⃗ P(θ) ⇒ ⃗ P( f (θ)); (114) Thus, the general solution involves an arbitrary monotonous function f (θ), and averaging over the fixed manifold of the solutions of the Navier-Stokes equations involves functional integration over all such functions.This integration includes summation over the positions θ 1 , θ 2 of the vorticity insertion points on a curve ⃗ C(θ).In the continuum theory, this would be an ordered Lebesgue integration (diffeomorphisms preserve the ordering of points on a curve) and in our case of piecewise constant curve with discontinuities ∆ ⃗ P k , it becomes an ordered sum The discontinuity ∆ ⃗ P(θ) stays finite in the continuum limit N → ∞.The continuum limit can be taken only after integrating(summing) the internal degrees of freedom of the fixed manifold of the Loop equations.The imaginary part of our solution (64) does not depend on the point on a circle.Therefore, it contributes a constant term into ⃗ S m,n which cancels in the difference ⃗ S n,m − ⃗ S m,n in the exponential, as it should.
Let us look at the correlation function (110).First, we expand and simplify the dot product involved The terms S m,n , S n,m in (110) have the following form

Critical phenomena in statistical limit
Now we are prepared to average over spin variables σ l , l = 0 . . .N − 1.This expression singles out the variables σ n , σ m so we can sum over these two variables, leaving the rest of σ l free, except for a constraint ∑ σ = rq.This constraint can be implemented as a discrete Fourier integral: We start with ; (125) Then we take the ω integral out of the sum: This expression can be readily averaged over two variables σ m , σ n , which reduces to the sum of four terms with σ n , σ m = ±1.Keeping only the factors depending on The next step would be averaging over the remaining variables σ l , l ̸ = n, l ̸ = m.These variables are split into two sets: one is used in the A n,m , and the other is used in A m,n+N .The variables A n,m have a certain distribution in the statistical limit when both ∼ m ∼ n ∼ N → ∞.We also have β ∼ 1/ √ N → 0 in that limit.We are now considering the unconstrained distribution over σ l , as the constraint is implemented via the discrete Fourier integral.Let us compute the mean and variance of In the relevant critical region β → 0, m = x/β 2 , n = y/β 2 , x > y, the ratio of the variance of U n,m to the square of its mean tends to a finite limit This function looks singular, but it is positive and finite in the allowed region 0 < y < x ( see Fig. 4).Therefore, the CLT does not apply to the distribution of A n,m , and all the moments of this distribution must be computed to obtain the probability distribution as a Mellin transform.
Calculating these moments represents another challenge to the number theorists.This complication is good news for our theory: we have critical phenomena with a non-Gaussian distribution in the statistical limit.Still, in the next section, we derive an exact formula for the mean value of the enstrophy as a function of the chemical potential µ, relating it to the functions of the number theory.

Analytic solution for the enstrophy
Let us find analytic formulas for observables of this remarkable statistical distribution, which is isomorphic to the Ising model tied to the ensemble of fractions.The one-dimensional Ising models are usually solvable, and this one is no exception.The simplest quantity is the enstrophy related to our random variables in the previous section.Setting⃗ r = 0 in (110) (which is possible at finite N) we find the following relation where the big Euler ensemble average ⟨. . .⟩ E(N) denotes averaging over p, q, r and the Ising variables σ subject to the global constraint ∑ σ k = qr .The general theory in section 4.2 expressed the big Euler ensemble average in terms of the small Euler ensemble average of the average over spins.So, we start the computation of the enstrophy by averaging over spins σ k , which is rather simple.Let us use the explicit expression (118) for the dot product ⃗ ω m • ⃗ ω n .Now, using the one-dimensional Fourier integral for the global constraint on σ, we get the unconstrained average: Averaging it over σ m , σ n we find Now, averaging over the remaining σ r , r ̸ = m, r ̸ = n is straightforward. exp We arrive at the integral This integral is calculable (see [39]): Summing over 0 ≤ n < m < N yields N(N − 1)/2, leading to our final answer We investigated this new function S(q) in [42] and represented it in terms of so-called multitotients [43].For the reader's convenience, we present the computations leading to this representation in AppendixA.
Here p|q are prime factors of q.This remarkable identity can be directly verified using Mathematica ® [39].It takes over a minute of CPU to compute and simplify S(100) = 2360.The same result using the multitotient formula takes 140 microseconds.Here is the table of S(q) for 2 ≤ q ≤ 10 This S(q) is an even integer for q > 3, and here is a simple proof.
Proof.The second term −φ(q) in S(q) is an even integer, and the first term can be rewritten as where α p ≥ 1 is a multiplicity of the prime factor p. The remaining factors (p 2 − 1) = (p − 1)(p + 1).The sequence of three integers (p − 1), p, (p + 1) has a factor divisible by three, and it cannot be p as it is a prime number, with the only exception being q = p = 3.
The plot of the first 100, 000 values of S(q) looks as follows (see Fig. 5).The appearance of prime numbers in the fluid dynamics problem is exciting; it reveals hidden relations between these different branches of modern mathematics.It reminds me of similar unexpected relations between matrix models, orthogonal polynomials, and 2D quantum gravity.φ 2 (q) 3 -φ(q)

Out[ ]=
Figure 5. S(q) for 10 < q < 100000.It does not reach any smooth limit; several bands persist up to infinity, similar to the Euler totient φ 1 (q).
The formula (139) is our exact solution for enstrophy, expressing it as a calculable average over the small Euler ensemble.In the next section, we compute it in the local limit µ → 0.

The local limit of the energy dissipation
In the limit N → ∞, the term with zero winding r = 0 dominates the sum ∑ r , which is only possible for even N.Its asymptotic limit is The remaining terms, with r ̸ = 0 can be replaced by an integral of the asymptotic expansion of the exact expression The leading term cancels after integration, so the r = 0 contribution dominates the sum.The next term of expansion in 1/N cancels as well We do not need this term O(1/ √ N), as the dominant r = 0 term grows as √ N. As seen in the section 4.2, the remaining terms decrease faster than exponential with q, making it a nontrivial exercise in number theory to find an analytic formula for the partition function and the enstrophy in the odd Euler ensemble.Let us concentrate on even N, as this term dominates the grand canonical ensemble we study.
The extra factor of 1/2 came from skipping all the odd values of N. The remaining sum over even N tends to be half the asymptotic expression's integral for these even N.The asymptotic behavior of multitotient summators has been known for a century [43] where ζ(n) is the Riemann's zeta function.We obtain the following local limit of the energy decay in the grand canonical Euler ensemble This energy dissipation stays finite in a turbulent limit ν → 0 provided Let us now estimate the odd N contribution to the enstrophy.Here, the sum of the spin average is dominated by q ≪ N, where we cannot use the asymptotic formula of the number theory for totients.Thus, we used numerical results of direct computations of the small Euler ensemble contribution from odd N to the enstrophy This contribution ∆E comes out negative.We fitted these results as log(−∆E) ≈ a + b log N + c log log N (see Fig. 6).Here are the fit statistics We have the following estimate for the odd N contribution to the enstrophy in the grand canonical Euler ensemble in the local limit: The leading term grows 1/µ 2 ; therefore, this correction is negligible.
The enstrophy divergence corresponds to anomalous dissipation in our theory.This divergence is the dual version of the original anomalous dissipation in the Navier-Stokes velocity field, coming from singular vorticity regions [1,35](Burgers vortexes).Here, it comes from large fluctuations of the fractal curve in the grand canonical Euler ensemble.These are quantum effects related to the prime factorization of large integers.

The higher moments of the enstrophy
The higher moments of the distribution of enstrophy are also calculable [39].We take the 2n point correlation function in the limit of the vanishing loop (157) With all different m 1 < • • • < m 2n , the averaging over the Ising variables σ l is straightforward.It leads to the integral over ω This integral is reduced to the hypergeometric function in [39]: In particular, This term does not depend on q; it dominates in the local limit N → ∞ at fixed n.
The general formula for the numerator of the moments at fixed N reads This cotangent sum is reduced to the superposition of multi-totients in Appendix.A The rational coefficients BernSum(n, n − j) are given by recurrent relations The specific cases are considered in AppendixA.We are interested in the local limit, even N → ∞.In this limit, the highest totient φ 2n (q) dominates the sum, and we find for the sum Finally, in the limit µ → +0, replacing the sum over even N by 1/2 of the integral of the asymptotic of this product of gamma functions, we find (171)

The vorticity distribution
Using the analytic formulas for the moments of the enstrophy, we can study nonperturbative effects in our solution, such as the PDf of the vorticity [39].Let us consider the Fourier transform of the PDF corresponding to the vorticity moments, (t+t 0 ) 2 µ 3 , truncated at n = 200.
The averaging over the global rotation matrix Ω ∈ O(3) yields Using our solution for the moments, we get a universal scaling function This expansion for the scaling function has a finite radius of convergence, as it follows from the asymptotic The singularity is located at negative z = −R.The expansion can be truncated inside the convergence region |z| < R.Here is the corresponding plot for expansion truncated at n = 200 (see Fig. 7).At large positive z, it behaves as This singularity of the Fourier transform at the imaginary axis corresponds to the exponential decrease of the PDF There is no continuum limit for the distribution of vorticity.However, this statistical system is renormalizable in the sense that the redefinition of the source eliminates the singularity in the same way as the redefinition of viscosity eliminated the singularity in the energy dissipation.

The decay index spectrum
The vicinity of the fixed point in nonlinear dynamic systems provides the most interesting physical parameters, such as anomalous dimensions of various local operators in the theory of renormalization group.Our system is no exception.

Linearized loop equation
Let us perturb the momentum loop equation( 42) and study the general properties of this perturbation δ ⃗ P(t, θ) = ⃗ Q(t, θ).We get a linearized equation for ⃗ Q(t, θ) of the form These matrices Ĥ1,2 [ ⃗ P] for the decaying solution (45) have an explicit time dependence which follows from the fact that the RHS of the original loop equation (42) represents the third-degree homogeneous functional of ⃗ P. We find the linear equation of the form There is no explicit θ dependence other than through the fixed point solution ⃗ F(θ).We shall write explicit equations in a minute.This equation has power-like solutions with the index λ determined from the eigenvalue problem We found scaling laws with anomalous dimensions without any renormalization group.We studied linearized equations near the fixed point.In the vicinity of a conventional fixed point, the perturbations would decay or grow exponentially with time by the Lyapunov indexes of the linearized equation.In our decay fixed point, the base solution decays as 1/ √ t, making the linearized operator decay ∝ 1/t.This decay of the operator converts the exponential decay of perturbations into a power decay without having scale invariant field theory.
We can go one step forward before specifying the parameters of this equation.Let us use our polygonal approximation for ⃗ F, ⃗ G.This equation becomes a recurrent equation We can solve it as a matrix product (in reverse order) The periodicity requires G N = G 0 which leads to the eigenvalue equation (this is As a result, we arrive at the spectral equation for the fractal dimensions of decaying turbulence det MN (λ) Here is the explicit form of these two matrices ( with Note that these two matrices are functions of γ, unlike the fixed point solution (64),(65), where the γ dependence dropped.This fact will be important for the distribution of the velocity circulation.
The spectral equations are simpler than they look.The key to simplification is a simple algebra satisfied by the two vectors Using this algebra, we may expand the inverse matrix ( Îλ + Âk ) −1 (− Îλ + Bk ) on five various tensor products and solve linear equations for expansion coefficients.We get (see [39]) Thus, we have factored out the poles from the matrix Mk Mn = Ĥn Mn has a n-th degree poles at λ = ± γ 2 and no other poles.The spectrum is determined by the polynomial equation of degree 6N spectrum : det ĤN − Î((4λ 2 − 1)γ 2 ) N = 0 (208) In three dimensions, this equation can be written as a cubic characteristic polynomial equation x = ((4λ 2 − 1)γ 2 ) N (210)

The spectral identity and Wilson loop asymptotics
The Wilson loop (3) doubles as a Fourier transform of the PDF for the velocity circulation.This PDF can be obtained by inverse Fourier transform However, substituting the leading term of the solution of the loop equation into this general formula leads to a paradox: this leading term does not depend on γ.Formally, we get the delta function δ(Γ), which means that this leading term is not sufficient to get the PDF; we need the next correction.
We have found the linear correction to the fixed point, and now we can fix these formulas.Expanding the correction and keeping the leading linear term, we find: The distribution of the eigenvalues λ can be studied using the resolvent By construction, at finite N, this is a rational function of λ with simple poles at the spectrum, as it follows from the representation The asymptotic behavior of this rational function at infinity (for even N) follows from (193): Combining these two properties, we get the sum over the spectrum plus an irrelevant constant The number of zeros z i can be counted from the rational function At infinity, this becomes (with 3N poles at λ = − 1 /2γ and 3N poles at On the other hand, we can write this rational function as Consider the anticlockwise contour ω surrounding the whole spectrum of poles z i in a complex plane.Then, we can use the unit residues at z i to compute the contour integral, with arbitrary holomorphic function F(λ) The left side of this identity can be inserted inside the statistical average, such as the circulation PDF integral: This integral, by design, equals the sum over the decay spectrum and will be dominated at large times by the zero z 0 with the smallest real part.After averaging over the random parameters of the Euler ensemble, the spectrum will likely become continuous.However, this relation can still be used to define the large-time behavior of the Wilson loop.
In this case, it is appropriate to deform the integration contour to the steepest descent from the saddle point.This contour will allow us an analytic investigation of the large-time asymptotic.The resolvent poles will become discontinuous across the cut in the complex λ plane.The presumed finite density of these condensed poles will result in logarithmic factors modifying the power laws of the time decay.The asymptotic behavior of W[C, Γ] at a large time t is determined by the saddle point equations for two variables γ, λ The stability of our fixed-point solution requires the condition of the saddle point We plan to study the spectrum of fractal dimensions and asymptotic distribution of the circulation in the forthcoming paper [42], using the NYU AD supercomputer cluster.

Singular solutions of the Navier-Stokes equation?
The issue of "finite time singularity" of the Navier-Stokes equation, particularly that without random forcing, has attracted much interest from mathematicians.The Millennium Prize Problem of proving or disproving the smoothness in the solution of the Navier-Stokes equation remains unsolved.The most advanced research in this field so far has been performed by Tao (see his recent review [44]).Based on a simplified model of the averaged Navier-Stokes, Tao conjectured [32] that irregular behavior occurs in finite time.One of his arguments is based on anomalous dissipation, coming from divergent enstrophy due to singular velocity gradients.The anomalous dissipation occurs only in the vanishing viscosity limit in the Navier-Stokes equation.
We investigated anomalous dissipation in the previous papers [1,35], where we have found the anomalous terms in the Euler Hamiltonian related to the Burgers vortex.This vortex corresponds to a singular Euler solution, with vorticity becoming the delta function at the infinitely thin vortex line in the ν → 0 limit of the Navier-Stokes equation.The analysis in [1,35] says nothing about the finite time singularity; this solution was suggested there as a stationary solution of the Navier-Stokes equation in the limit of vanishing viscosity without specifying the time evolution preceding this stationary solution.In the dual theory developed in this paper, the anomalous dissipation comes from large fluctuations of our fractal curve, leading to divergent enstrophy.Surprisingly, we have more than anomalous dissipation in the present theory: the singularity we have found exists at finite viscosity, in the spirit of Tao's conjecture.However, we are not studying a particular solution with a finite time blow-up.On the contrary, we have a time evolution in our solution (time decay) such that vorticity distribution is singular at every moment.Again, let us stress it: we are not claiming any theorems about finite time singularities of solutions of Navier-Stokes equation with some initial data.

Stochastic solution of the Navier-Stokes equation and ergodic hypothesis
Statistical "analysis of circulating or turbulent fluids" was defined by Feynman [45] as the last unsolved problem of classical physics.We are pursuing this problem by finding a stochastic solution of the Navier-Stokes equation covering a certain manifold.Our singularities arise in correlation functions after averaging over this manifold of solutions.We find this manifold (Euler ensemble) by solving the loop equation (a subset of the Hopf functional equation for the generating functional of velocity field probability).None of the particular solutions in this manifold has a finite time blow-up.The singularity emerges from averaging the distribution of these solutions over the Euler ensemble.
We take the most natural invariant measure from the point of the number theory: each element of the Euler ensemble enters with equal weight.We call it our ergodic hypothesis.This hypothesis is not necessary to solve the loop (i.e., Hopf) equation, as any linear superposition of the found solutions would satisfy the loop equation.The singularities of our Euler grand canonical ensemble at µ → 0 remain in local variables such as enstrophy and its PDF, indicating an intrinsic problem of the Navier-Stokes equation.It has to be regularized at small distances, and viscosity does not provide enough regularization in our solution.
The situation reminds that of the QED in the mid-twentieth century.The continuum limit of the theory showed unexpected divergent integrals due to an infinite number of local degrees of freedom of a continuous electromagnetic field.In both cases, the continuum theory was an idealization of the real physical system: in the case of QED, there were other forces at small distances, eventually merging QED into the Standard Model, which is still inconsistent at Planck's ultra-small distances where the quantum gravity enters the game.In the case of the Navier-Stokes equation, it ignores the molecular structure of fluid.The incompressibility is also an idealization: at large gradients of velocity, the sound waves related to compressibility lead to some physical cutoff of infinite vorticity.
In other words, the Navier-Stokes equation has limited applicability in the physical world and needs to be modified at large gradients.After all, this is a phenomenological equation resulting from the truncated expansion of friction forces in velocity gradients.The common presumption that "viscosity regularizes the velocity field" must be tested beyond the perturbation expansion, which we did in this work.We present a singular decay solution of the Navier-Stokes loop equations in arbitrary dimension d > 2. However, we did not prove that the solution of the Navier-Stokes equation as a PDE with smooth initial data would eventually approach our stochastic solution as an asymptotic regime.

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.The Turbulence theory aims to find a replacement of the microcanonical ensemble for the Navier-Stokes equation.This surface would also take part 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 specification of the mechanism of the stochastization.Indeed, the Gibbs, as well as the 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 as a fixed point manifold.
In the case of the decaying turbulence, the loop equations represent a closed subset of the Hopf equations, which are still sufficient to generate the dynamics of vorticity.In this case, the exact solution we have found for the Hopf functional also follows from the integrals of motion, this time, the conservation laws in the loop space.The loop space Hamiltonian we have 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 the plane wave exp(ı⃗ p • ⃗ x).This plane wave is the solution we have found, except the dot product ⃗ p • ⃗ x becomes a symplectic form ⃗ P(θ) • d ⃗ C(θ) in the loop space.Our momentum ⃗ P(θ, t) is not an integral of motion, but simple scaling properties of the pure Navier-Stokes equation lead to the solution with ⃗ P(θ, t) ∝ ⃗ F(θ)/ √ t, with ⃗ F(θ) being the integral of motion.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 ⃗ F(θ).
Our solution expresses the probability distribution and expected value for the Wilson loop at any given moment t in terms of an ensemble of fractal loops in complex momentum space.The loop is represented by a polygon with N → ∞ sides.This statistical system is similar to a one-dimensional Ising ring in an imaginary magnetic field ıβ = 2ıπ p q and zero coupling constant.Some global observables, such as the moments of enstrophy, are calculable for arbitrary N as an analytic function of N, p, q, relating it to the Euler totients and similar functions of the number theory.The continuum limit N → ∞ differs for odd and even N, which means this limit does not exist.This ambiguity disappears in the grand canonical ensemble.In this ensemble, the number N of degrees of freedom is not fixed but can also fluctuate, with the weight exp(−µN).These fluctuations smooth out the difference between odd and even ensembles so that the grand canonical ensemble is unambiguous in the continuum limit.The continuum limit in the grand canonical ensemble corresponds to µ → 0. In this limit, we compute the partition function (108) and the expectation value of the energy dissipation (151).As discussed in the previous section, the singularities at µ → 0 indicate inconsistencies in the Navier-Stokes equation as an idealization of molecular dynamics.

Classical flow and quantum geometry
Our computations heavily rely on the number theory, particularly Lehmer's multitotients φ l (q), (141), generalizing [43] the Euler totient function.What could the number theory have in common with the turbulent flow?The quantization of parameters of the fixed manifold of decaying turbulence stems from the deep quantum correspondence we have discovered.The statistical distribution of a nonlinear classical Navier-Stokes PDE is exactly related to the wave functional of a linear Schrödinger equation in the loop space.The quantization mechanism is the same as in ordinary quantum mechanics: this is a requirement of the periodicity of the solution.The equivalence of a strong coupling phase of the fluctuating vector field to quantum geometry is a well-known duality phenomenon in gauge theory (the ADS/CFT duality), ringing a bell to the modern theoretical physicist.In our case, this is a simpler quantum geometry: a fractal curve in complex space.
An expert in the traditional approach to turbulence may wonder why the loop equation's solutions have any relation to the velocity field's statistics in a decaying turbulent flow.Such questions were raised and answered in the last few decades in the gauge theories, including QCD [2,4,5,38] where the loop equations were derived first [36,37].The short answer is that duality only applies to the correlation functions of two theories with different dynamical variables; there is no correspondence between these variables, but the correlation functions are identical.Mathematical physics sometimes has alternative languages for the same phenomena; examples are the duality between Schrödinger's wave equation and Heisenberg's matrix mechanics, between dynamical triangulation and Liouville theory in 2D quantum gravity.
Extra complications in the gauge theory are the short-distance singularities related to the infinite number of fluctuating degrees of freedom in quantum field theory.The Wilson loop functionals in coordinate space are singular in the gauge field theory and cannot be multiplicatively renormalized.Perturbatively, there is no short-distance divergence in the Navier-Stokes equations nor the Navier-Stokes loop dynamics.The Euler equations represent the singular limit, which, as we argued, should be resolved using singular topological solitons regularized by the Burgers vortex.In the present theory, we keep viscosity constant and do not encounter any singularities in coordinate space.The anomalous dissipation is achieved in our solution via a completely different mechanism: large fluctuations of the fractal curve at p ≪ q.However, we have found the singularity from large vorticity fluctuations at any finite viscosity value.It cannot be attributed to the Euler singularities such as line vortexes.Those vortexes are regularized by finite viscosity, unlike our singularities.
This singularity resembles the conjecture by Tao [32]; however, it has a different meaning.Our singularity displays itself in correlation functions of vorticity, not the local vorticity as a function of coordinates.Our vorticity is not a smooth field, developing finite-time singularity at some region of physical space, like a vortex line or a vortex sheet; it is a stochastic field with the singularities (discontinuities) distributed all over the physical space with some multi-dimensional distribution.We compute the correlation functions for this distribution from the dual theory of momentum loop dynamics.We have singularities in these correlation functions, like QED, though these singularities emerge in the exact solution beyond perturbation theory.

Stokes-type functionals and vorticity correlations
The loop equation describes the gauge invariant sector of the gauge field theory.Therefore, the gauge degrees of freedom are lost in the loop functional.However, the gauge-invariant correlations of the field strength are recoverable from the solutions of the loop equation [36,37].There is no gauge invariance regarding the velocity field in fluid dynamics (though there is such invariance in the Clebsch variables [1]).The longitudinal, i.e., a potential part of the velocity, has a physical meaning -it is responsible for pressure and energy pumping.This part is lost in the loop functional but is recoverable from the rotational part (the vorticity) using the Biot-Savart integral.In the Fourier space, the correlation functions of the velocity field are algebraically related to those of vorticity Thus, the general solution for the Wilson loop functional Ψ[γ, C] allows computing both vorticity and velocity correlation functions.We demonstrated that in the last two sections by computing the moments of the enstrophy and resulting anomalous dissipation.This computation is nonperturbative: it corresponds to the extreme turbulent limit and cannot be expanded in inverse powers of viscosity.

Relation of our solution to the weak turbulence
The solution of the loop equation with finite area derivative, satisfying Bianchi constraint, belongs to the so-called Stokes-type functionals [36], the same as the Wilson loop for Gauge theory and fluid dynamics.The Navier-Stokes Wilson loop is a case of the Abelian loop functional, with commuting components of the vector field ⃗ v.As we discussed in detail in [1,36,37], any Stokes-type functional Ψ[γ, C] satisfying boundary condition at shrunk loop Ψ[0] = 1, and solving the loop equation can be iterated in the nonlinear term in the Navier-Stokes equations (which iterations would apply at large viscosity).
The resulting expansion in inverse powers of viscosity (weak turbulence) exactly coincides with the ordinary perturbation expansion of the Navier-Stokes equations for the velocity field, averaged over the distribution of initial data or boundary conditions at infinity.We have demonstrated in [1,46] (and also here, in Section 2.3) how the velocity distribution for the random uniform vorticity in the fluid was reproduced by a singular momentum loop ⃗ P(θ).The solution for ⃗ P(θ) in this special fixed point of the loop equation was random complex and had slowly decreasing Fourier coefficients, leading to a discontinuity sign(θ − θ ′ ) in a pair correlation function (30).The corresponding Wilson loop was equal to the Stokes-type functional (22).Using this example, we demonstrated how a discontinuous momentum loop describes the vorticity distribution in the stochastic Navier-Stokes flow.In this example, the vorticity is a global random variable corresponding to a random uniform fluid rotation: a well-known exact solution of the Navier-Stokes equation.This example corresponds to a special fixed point for the loop equation, not general enough to describe the turbulent flow but mathematically ideal as a toy model for the loop technology.It demonstrates how the momentum loop solution sums up all the terms of the 1/ν expansion in the Navier-Stokes 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 independent of viscosity.At a given initial condition, after a finite time, the solution will still depend on viscosity and initial condition.At large time, though, it will approach our universal fixed manifold and (supposedly, for random initial data) cover it uniformly, according to the Euler ensemble measure.The vorticity will become a random variable with a singular distribution in the local limit, suggesting intrinsic inconsistency of the Navier-Stokes equation.

Continuum spectrum of anomalous dimensions and multifractality
We studied the linear perturbations of our fixed manifold, which decay with time as t −λ with the spectrum of λ determined by the large N polynomial equation (208).The calculable time decay spectrum at finite N is a nice surprise, but the limit N → ∞ remains an unsolved problem.Future investigation will relate it to some "multifractal" properties in the probability distribution emerging after the saddle point computation.The effective decay spectrum is continuous, meaning the modification of naive scaling laws for decaying turbulence by some powers of the logarithm of time.The analytic result for the largetime asymptotic would require more work, however.Meanwhile, we are involved in the numerical study of this asymptotic behavior [42].

•
The solution for ⃗ P(θ) in the loop equation for decaying turbulence, which we have found in this paper, exhibits nonperturbative features, particularly the quantization of parameters.These quantum effects follow from the exact equivalence of the Navier-Stokes statistics to the quantum mechanics in loop space [1].The most striking nonperturbative feature is the singularity of vorticity distribution in the local limit for a finite viscosity.

•
Compared to the other critical phenomena, this theory is amazingly simple: It is not a field theory but a quantum statistical system similar to the one-dimensional Ising ring in the presence of the imaginary quantized magnetic field.The moments of the distribution of vorticity are analytically calculable.

•
Still, the solution is far from trivial and reveals unexpected relations between turbulence and number theory.The analytic computation of the decay spectrum from the general equation (214) for the resolvent remains a challenge to the number theory.• This work does not claim a full solution to the decaying turbulence problem.At best, this is the new path to the solution; following this path will require joint efforts of mathematical physicists, number theorists, and computer scientists.
Funding: This research was supported by a Simons Foundation award ID 686282 at NYU Abu Dhabi.

Data Availability Statement:
The Mathematica ® notebooks used to verify the equations and compute some functions are available for download in [39].

Figure 2 .
Figure 2. Partition function Z(N, N) for odd N, fitted as a √ N log(N) b

Figure 3 .
Figure 3. Backtracking wires corresponding to vorticity correlation function in (110).With these backtracking wires, the correlation function reduces to the closed loop functional, which is represented by our solution with fractal momentum loop ⃗ P(θ).

Figure 4 .
Figure 4.The relative variance computed in the text in the statistical limit.

Figure 6 .
Figure 6.The direct computation of the odd N contribution ∆E to the enstrophy with 20 digits fitted as ∆E ≈ −e a N b log c N.
These coefficients Ξ n are all positive.Here are the first five values