Exact Solution of Decaying Turbulence

1. Introduction

A while ago, we derived [1,2] a functional equation for the so-called loop average [3,4] or Wilson loop in turbulence. The path to an exact solution by a dimensional reduction in this equation was proposed in the ’93 paper [2] 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 [2], could not be verified in DNS at the time with existing computer power.

The situation has changed over the last decades. No alternative microscopic theory based on the Navier-Stokes equation emerged, but our understanding of the strong turbulence phenomena grew significantly.

On the other hand, the loop equations technology in the gauge theory also advanced over the last decades. The correspondence between the loop space functionals and the original vector fields was better understood, and various solutions to the gauge loop equations were found.

In particular, the momentum loop equation was developed, similar to our momentum loop used below [5,6,7]. Recently, some numerical methods were found to solve loop equations beyond perturbation theory [8,9].

The loop dynamics was extended to quantum gravity, where it was used to study nonperturbative phenomena [10,11].

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 [12,13,14,15] 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 [16].

However, as we argued in [17,18], at those Reynolds numbers, the DNS cannot yet distinguish between pure scaling laws with anomalous dimension $\zeta \left(n\right)$ and some algebraic function of the logarithm of scale $\zeta (n,logr)$ modifying the K41 scaling.

Theoretically, we studied the loop equation in the confinement region (large circulation over large loop C), and we have justified the Area law, suggested back in ’93 on heuristic arguments [2].

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 [12] which triggered the further development of the geometric theory of turbulence[13,14,15,17,18,19,20,21,22,23,24,25,26,27,28,29,30].

In particular, the Area law was justified for flat and quadratic minimal surfaces [22], and an exact scaling law in confinement region $\Gamma \propto \sqrt{Area}$ was derived [21]. The area law was verified with better precision in [13].

It was later conjectured in [17] that the dominant field configurations in extreme turbulence are so-called Kelvinons, which were shown to solve stationary Navier-Stokes equations assuming the sparse distribution of vorticity structures.

These topological solitons of the Euler theory are built around a vortex sheet bounded by a singular vortex line. This vortex line is locally equivalent to the cylindrical Burgers vortex [31], with infinitesimal thickness in the limit of a large Reynolds number.

As we argued in [17,18], the Kelvinon has an anomalous dissipation, surviving the strong turbulent limit. This dissipation is proportional to the square of constant circulation of the Burgers vortex times a line integral of the tangent component of the strain along the loop.

The Kelvinon minimizes the energy functional, with anomalous terms coming from the Burgers core of the vortex line. There is also a constant scale factor Z in the representation of the Kelvinon vorticity in terms of spherical Clebsch variables:

$$\begin{array}{ccc}\hfill & & \overrightarrow{\omega }=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.Z{e}_{abc}{S}_a\overrightarrow{\nabla }{S}_b\times \overrightarrow{\nabla }{S}_c=\overrightarrow{\nabla }{\varphi }_1\times \overrightarrow{\nabla }{\varphi }_2;\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & S_1^2+S_2^2+S_3^2=1;\hfill \end{array}$$

(2)

$$\begin{array}{ccc}& & {\varphi }_2=arg(S_1+ıS_2);{\varphi }_1=Z{S}_3;\hfill \end{array}$$

(3)

In that paper, the constant Z was related to the Kolmogorov energy dissipation density and the boundary value of the $S_3$ variable at the loop C.

The anomalous Hamiltonian [17,18] explicitly violated the K41 scaling by the logarithmic terms $logZ/\nu$ in the region of small loops C. This region resembles the asymptotically free QCD. The logarithmic terms were summed up by RG equation with running coupling constant logarithmically small in this region.

These exciting developments explain and quantitatively describe many interesting phenomena [18] but do not provide a complete microscopic theory covering the full inertial range of turbulence without simplifying assumptions of the sparsity of vortex structures.

Moreover, while the Kelvinon (presumably) solves the stationary Navier-Stokes equations, it does not solve the loop equations for the following reason.

The loop equation assumes that the velocity field is independent of the loop C. In this case, the circulation ${\oint }_C{v}_{\alpha }d{r}_{\alpha }$ variations in the loop functional by the shape C of the loop can be reduced to the Navier-Stokes equation.

Otherwise, the variation would also involve the variation of the velocity field ${\oint }_C\delta v_{\alpha }d{r}_{\alpha }$.

This problem does not invalidate the Kelvinon theory as an ideal gas of random vortex rings sparsely distributed in a turbulent flow.

The loop functional is not needed for that statistical theory, and the stationary solution of the Navier-Stokes equation is sufficient. The shape of the loop and the vortex sheet inside would become random variables influenced by a background strain like in the pure vortex sheet solutions [18].

These objections, however, prevent the Kelvinon gas model from being a complete theory of strong isotropic turbulence. This model is merely an approximation of the full theory.

In the present work, we develop the theory free of these assumptions and approximations by analytically solving the loop equations for decaying turbulence.

2. Loop equation

2.1. Loop operators

We introduced the loop equation in the Lecture Series at Cargese and Chernogolovka Summer Schools [2].

Here is a summary for the new generation.

We write the Navier-Stokes equation as follows

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

(4)

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

(5)

The Wilson loop average for the turbulence

$$\begin{array}{c}\hfill \mathsf{\Psi }[\gamma ,C]=〈exp\left(\frac{ı\gamma }{\nu }{\oint }_C{v}_{\alpha }d{r}_{\alpha }\right)〉\end{array}$$

(6)

treated as a function of time and a functional of the periodic function $C:r_{\alpha }=C_{\alpha }\left(\theta \right);\theta \in (0,2\pi )$ (not necessarily a single closed loop), satisfies the following functional equation

$$\begin{array}{ccc}\hfill & & ı\nu {\partial }_t\mathsf{\Psi }={\mathcal{H}}_C\mathsf{\Psi };\hfill \end{array}$$

(7a)

$$\begin{array}{ccc}& & {\mathcal{H}}_C={\mathcal{H}}_C^{\left(1\right)}+{\mathcal{H}}_C^{\left(2\right)}\hfill \end{array}$$

(7b)

$$\begin{array}{ccc}& & {\mathcal{H}}_C^{\left(1\right)}=\nu \gamma {\oint }_{C}d{r}_{\alpha }{\partial }_{\beta }{\widehat{\omega }}_{\alpha \beta }\left(r\right);\hfill \end{array}$$

(7c)

$$\begin{array}{ccc}& & {\mathcal{H}}_C^{\left(2\right)}=\gamma {\oint }_{C}d{r}_{\alpha }{\widehat{\omega }}_{\alpha \beta }\left(r\right){\widehat{v}}_{\beta }\left(r\right);\hfill \end{array}$$

(7d)

$$\begin{array}{ccc}& & {\widehat{\omega }}_{\alpha \beta }\equiv -ı\frac{\nu }{\gamma }\frac{\delta }{\delta {\sigma }_{\alpha \beta }}\hfill \end{array}$$

(7e)

$$\begin{array}{ccc}& & {\widehat{v}}_{\beta }\left(r\right)=\frac{1}{{\partial }_{\mu }^2}{\partial }_{\alpha }{\widehat{\omega }}_{\beta \alpha }\left(r\right)\hfill \end{array}$$

(7f)

We added a dimensionless factor $\gamma$ in the exponential compared to some previous definitions as an extra parameter of the Wilson loop. Without loss of generality, we shall assume that $\gamma >0$. The negative $\gamma$ corresponds to a complex conjugation of the Wilson loop.

In Abelian gauge theory, this parameter would be the continuous electric charge.

The statistical averaging $〈\cdots 〉$ corresponds to initial randomized data, to be specified later.

The area derivative $\frac{\delta }{\delta {\sigma }_{\alpha \beta }}$ is related to the variation of the functional when the little closed loop $\delta C$ is added

$$\begin{array}{ccc}\hfill & & {\Sigma }_{\alpha \beta }\left(\delta C\right)\frac{\delta F\left[C\right]}{\delta {\sigma }_{\alpha \beta }\left(r\right)}=F[C+\delta C]-F\left[C\right];\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill & & {\Sigma }_{\alpha \beta }\left(\delta C\right)=\frac{1}{2}{\oint }_{\delta C}{r}_{\alpha }d{r}_{\beta }\hfill \end{array}$$

(9)

In the review, [2,18], 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 ${\partial }_{\mu },{\widehat{\omega }}_{\alpha \beta },{\widehat{v}}_{\alpha }$ are expressed in terms of the spike operator

$$D_{\alpha }(\theta ,ϵ)={\int }_{-ϵ}^{+ϵ}d\xi \left(1-\frac{\left|\xi \right|}{ϵ}\right)\frac{\delta }{\delta C_{\alpha }(\theta +\xi )}$$

(10)

The area derivative operator can be regularized as

$$\begin{array}{ccc}\hfill & & {\mathsf{\Omega }}_{\alpha \beta }(\theta ,ϵ)=-ı\frac{\nu }{\gamma }\frac{\delta }{\delta C_{\alpha }^{\prime }\left(\theta \right)}{\int }_{-ϵ}^{ϵ}d\xi \frac{\delta }{\delta C_{\beta }(\theta +\xi )}-\{\alpha \leftrightarrow \beta \};\hfill \end{array}$$

(11)

and velocity operator (with $\delta ,ϵ\to 0^{+}$)

$$V_{\alpha }(\theta ,ϵ,\delta )=\frac{1}{D_{\mu }^2(\theta ,ϵ)}{D}_{\beta }(\theta ,ϵ){\mathsf{\Omega }}_{\beta \alpha }(\theta ,\delta );$$

(12)

In addition to the loop equation, every valid loop functional $F\left[C\right]$ must satisfy the Bianchi constraint [3,4]

$$\begin{array}{c}\hfill {\partial }_{\alpha }\frac{\delta F\left[C\right]}{\delta {\sigma }_{\beta \gamma }\left(r\right)}+\mathbf{cyclic}=0\end{array}$$

(13)

In three dimensions, it follows from identity $\overrightarrow{\nabla }·\overrightarrow{\omega }=0$; in general dimension $d>3$, the dual vorticity $\tilde{\omega }$ 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 [18] 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 [18] in the loop Hamiltonian ${\mathcal{H}}_C$, breaking certain symmetries needed for the dimensional reduction we study below.

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 at the initial moment $t=0$ and slowly dissipates over time, provided the viscosity is small enough, corresponding to the large Reynolds number we are studying.

2.2. Dimensional reduction

The crucial observation in [2] 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

$$\mathsf{\Psi }[\gamma ,C]=〈exp\left(\frac{ı\gamma }{\nu }\oint d{C}_{\alpha }\left(\theta \right)P_{\alpha }\left(\theta \right)\right)〉$$

(14)

The difference with the original definition of $\mathsf{\Psi }[\gamma ,C]$ is that our new vector function $\overrightarrow{P}\left(\theta \right)$ depends directly on $\theta$ rather then through the vector field $\overrightarrow{v}\left(\overrightarrow{r}\right)$ in ${\mathbb{R}}_d$ projected at $\overrightarrow{r}=\overrightarrow{C}\left(\theta \right)$.

This transformation is the dimensional reduction $d\Rightarrow 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, depending on our position C in the loop space).

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 $\overrightarrow{P}(\theta ,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 our objection (see the Introduction) to the Kelvinon theory and any other Navier-Stokes stationary solutions with a singularity at fixed loop C in space.

The functional derivative, acting on the exponential in (14), could be replaced by the derivative $P^{\prime }$ as follows

$$\frac{\delta }{\delta C_{\alpha }\left(\theta \right)}\leftrightarrow -\frac{ı\gamma }{\nu }{P}_{\alpha }^{\prime }\left(\theta \right)$$

(15)

The loop equation for $P\left(\theta \right)$ as a function of $\theta$ and also a function of time, reads:

$${\partial }_t{P}_{\alpha }=\left(\nu D_{\beta }-V_{\beta }\right){\mathsf{\Omega }}_{\beta \alpha }$$

(16)

where the operators $V,D,\mathsf{\Omega }$ should be regarded as ordinary numbers, with the following definitions.

The spike derivative D in the above equation

$$\begin{array}{ccc}\hfill & & D_{\alpha }(\theta ,ϵ)=-\frac{ı\gamma }{\nu }{\int }_{-1}^{1}d\mu sgn\left(\mu \right)P_{\alpha }\left(\theta +ϵ\mu \right)\hfill \end{array}$$

(17)

The vorticity (11) and velocity (12) also become singular functionals of the trajectory $P\left(\theta \right)$.

The first observation about this equation is that the viscosity factor cancels after the substitution (17).

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(\theta ,ϵ)$ turns to the discontinuity $\Delta P\left(\theta \right)=P\left({\theta }^{+}\right)-P\left({\theta }^{-}\right)$ in the limit $ϵ\to 0^{+}$

$$\begin{array}{c}\hfill D(\theta ,0^{+})=-\frac{ı\gamma }{\nu }\Delta P\left(\theta \right)\end{array}$$

(18)

The relation of the operators in the QCD loop equation to the discontinuities of the momentum loop was noticed, justified, and investigated in [6,7].

The momentum loop could be piecewise constant with an arbitrary number of such discontinuities.

In the Navier-Stokes theory, this relation provides the key to the exact solution.

In the same way, we find the limit for vorticity

$$\begin{array}{ccc}\hfill & & {\mathsf{\Omega }}_{\alpha \beta }(\theta ,0^{+})=\frac{-ı\gamma }{\nu }{P}_{\alpha \beta }\left(\theta \right);\hfill \end{array}$$

(19)

$$\begin{array}{ccc}& & P_{\alpha \beta }\left(\theta \right)=\Delta P_{\alpha }\left(\theta \right)P_{\beta }\left(\theta \right)-\{\alpha \leftrightarrow \beta \};\hfill \end{array}$$

(20)

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

(21)

and velocity (skipping the common argument $\theta$ )

$$\begin{array}{ccc}& & V_{\alpha }=\frac{\Delta P_{\beta }}{\Delta P_{\mu }^2}{P}_{\beta \alpha }=P_{\alpha }-\frac{\Delta P_{\alpha }\Delta P_{\beta }{P}_{\beta }}{\Delta P^2}\hfill \end{array}$$

(22)

The Bianchi constraint is identically satisfied as it should

$$\begin{array}{c}\hfill \Delta P_{\alpha }\left(\Delta P_{\beta }{P}_{\gamma }-\{\beta \leftrightarrow \gamma \}\right)+\mathbf{cyclic}=0\end{array}$$

(23)

We arrive at a singular loop equation for $P_{\alpha }\left(\theta \right)$

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

(24)

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 [18], 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 $\overrightarrow{P}(\theta ,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.

2.3. Random global vorticity

Possible initial data for the reduced dynamics were suggested in the original papers [2,18]. The initial velocity field’s simplest meaningful distribution is the Gaussian one, with energy concentrated in the macroscopic motions. The corresponding loop field reads (we set $\gamma =1$ for simplicity in this section)

$$\begin{array}{ccc}& & {\mathsf{\Psi }}_0\left[C\right]=exp\left(-\frac{1}{2{\nu }^2}{\int }_{C}d\overrightarrow{C}\left(\theta \right)·d\overrightarrow{C}\left({\theta }^{\prime }\right)f\left(\overrightarrow{C}\left(\theta \right)-\overrightarrow{C}\left({\theta }^{\prime }\right)\right)\right)\hfill \end{array}$$

(25)

where $f\left(\overrightarrow{r}\right)$ is the velocity correlation function

$$〈v_{\alpha }\left(r\right)v_{\beta }\left(r^{\prime }\right)〉=\left({\delta }_{\alpha \beta }-{\partial }_{\alpha }{\partial }_{\beta }{\partial }_{\mu }^{-2}\right)f(r-r^{\prime })$$

(26)

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^{\prime })\to f_0-f_1{(r-r^{\prime })}^2+\cdots$$

(27)

The first term $f_0$ is proportional to initial energy density,

$$\frac{1}{2}〈v_{\alpha }^2〉=\frac{d-1}{2}{f}_0$$

(28)

and the second one is proportional to initial energy dissipation rate ${\mathcal{E}}_0$

$$f_1=\frac{{\mathcal{E}}_0}{2d(d-1)\nu }$$

(29)

where $d=3$ is the dimension of space.

The constant term in (27) as well as $r^2+r^{\prime 2}$ terms drop from the closed loop integral, so we are left with the cross-term $r{r}^{\prime }$, which reduces to a full square

$$\begin{array}{ccc}\hfill & & {\mathsf{\Psi }}_0\left[C\right]\to exp\left(-\frac{f_1}{{\nu }^2}{\left(\oint d{C}_{\alpha }\left(\theta \right)C_{\beta }\left(\theta \right)\right)}^2\right)\hfill \end{array}$$

(30)

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

$$\begin{array}{ccc}& & {\mathsf{\Psi }}_0\left[C\right]\to \mathit{const}\int \left(d\varphi \right)exp\left(-{\varphi }_{\alpha \beta }^2\right)\hfill \\ & & exp\left(2ı\frac{\sqrt{f_1}}{\nu }{\varphi }_{\mu \nu }\oint d{C}_{\mu }\left(\theta \right)C_{\nu }\left(\theta \right)\right)\hfill \end{array}$$

(31)

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

The physical meaning of this $\varphi$ is the random uniform vorticity $\widehat{\omega }=\sqrt{f_1}\widehat{\varphi }$ 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 [18]. The uniform global rotation represents a fixed point of the Navier-Stokes equation for arbitrary uniform vorticity tensor.

Gaussian integration by $\varphi$ 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 (30).

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_{\alpha }\left(\theta \right)$ for this fixed point.

The functional Fourier transform [2,18] leads to the following simple result for the initial values of $P_{\alpha }\left(\theta \right)$.

In terms of Fourier harmonics, this initial data read

$$\begin{array}{ccc}& & P_{\alpha }\left(\theta \right)=\sum _{\mathrm{odd}n=1}^{\infty }{P}_{\alpha ,n}exp\left(ın\theta \right)+{\overline{P}}_{\alpha ,n}exp\left(-ın\theta \right);\hfill \end{array}$$

(32)

$$\begin{array}{ccc}& & P_{\alpha ,n}=\mathcal{N}(0,1);\hfill \end{array}$$

(33)

$$\begin{array}{ccc}& & {\overline{P}}_{\alpha ,n}=\frac{4\sqrt{f_1}}{n\nu }{\varphi }_{\alpha \beta }{P}_{\beta ,n};\hfill \end{array}$$

(34)

$$\begin{array}{ccc}& & {\varphi }_{\alpha \beta }=-{\varphi }_{\beta \alpha };\hfill \end{array}$$

(35)

$$\begin{array}{ccc}& & {\varphi }_{\alpha \beta }=\mathcal{N}(0,1)\forall \alpha <\beta ;\hfill \end{array}$$

(36)

The constant part $P_{\alpha ,0}$ of $P_{\alpha }\left(\theta \right)$ is not defined, but it drops from equations by translational invariance.

Note that this initial data is not real, as ${\overline{P}}_{\alpha ,n}\ne P_{\alpha ,n}^{★}$ . Positive and negative harmonics are real but unequal, leading to a complex Fourier transform. At fixed tensor $\varphi$ the correlations are

$$\begin{array}{ccc}& & {〈P_{\alpha ,n}{P}_{\beta ,m}〉}_{t=0}=\frac{4\sqrt{f_1}}{m\nu }{\delta }_{nm}{\varphi }_{\alpha \beta };\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill & & {〈P_{\alpha }\left(\theta \right)P_{\beta }\left({\theta }^{\prime }\right)〉}_{t=0}=2ı\frac{\sqrt{f_1}}{\nu }{\varphi }_{\alpha \beta }sign({\theta }^{\prime }-\theta );\hfill \end{array}$$

(38)

This correlation function immediately leads to the uniform expectation value of the vorticity

$$〈P_{\alpha }\left(\theta \right)\Delta P_{\beta }\left(\theta \right)〉=4ı\frac{\sqrt{f_1}}{\nu }{\varphi }_{\alpha \beta };\forall \theta$$

(39)

The uniform constant vorticity kills the linear term of the Navier-Stokes equation in the original loop space, involving ${\partial }_{\alpha }{\widehat{\mathsf{\Omega }}}_{\alpha \beta }=0$.

The nonlinear term ${\widehat{V}}_{\alpha }{\widehat{\mathsf{\Omega }}}_{\alpha \beta }$ vanishes in the coordinate loop space only after integration around the loop.

Here are the steps involved

$$\begin{array}{ccc}& & {\widehat{V}}_{\beta }=\frac{1}{2}{\widehat{\mathsf{\Omega }}}_{\alpha \beta }{C}_{\beta };\hfill \end{array}$$

(40)

$$\begin{array}{ccc}& & \oint {\widehat{\mathsf{\Omega }}}_{\alpha \beta }{C}_{\beta }{\widehat{\mathsf{\Omega }}}_{\beta \gamma }d{C}_{\alpha }\propto {\widehat{\mathsf{\Omega }}}_{\alpha \beta }{\widehat{\mathsf{\Omega }}}_{\beta \gamma }{\Sigma }_{\alpha \beta }\left(C\right);\hfill \end{array}$$

(41)

The tensor area $\Sigma$ was defined in (). It is an antisymmetric tensor; therefore its trace with a symmetric tensor ${\widehat{\mathsf{\Omega }}}_{\alpha \beta }{\widehat{\mathsf{\Omega }}}_{\beta \gamma }$ 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.

2.4. Decay or fixed point

The absolute value of loop average $\mathsf{\Psi }[\gamma ,C]$ stays below 1 at any time, which leaves two possible scenarios for its behavior at a large time.

$$\begin{array}{ccc}\hfill \mathbf{Decay}:& & \overrightarrow{P}\to 0;\mathsf{\Psi }[\gamma ,C]\to 1;\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill \mathbf{Fixed}\mathbf{Point}:& & \overrightarrow{P}\to {\overrightarrow{P}}_{\infty };\mathsf{\Psi }[\gamma ,C]\to {\mathsf{\Psi }}_{\infty }\left[C\right];\hfill \end{array}$$

(43)

The Decay scenario in the nonlinear ODE (24) corresponds to the $1/\sqrt{t}$ decrease of $\overrightarrow{P}$.

Omitting the common argument $\theta$, we get the following exact time-dependent solution (not just asymptotically, at $t\to +\infty$).

$$\begin{array}{ccc}\hfill & & \overrightarrow{P}=\sqrt{\frac{\nu }{2(t+t_0)}}\frac{\overrightarrow{F}}{\gamma };\hfill \\ \hfill & & \left({(\Delta \overrightarrow{F})}^2-1\right)\overrightarrow{F}=\hfill \end{array}$$

(44)

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

(45)

The Fixed Point would correspond to the vanishing right side of the momentum loop equation (24). Multiplying by ${(\Delta \overrightarrow{P})}^2$ and reducing the terms, we find a singular algebraic equation

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

(46)

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 (46) is a linear relation between two vectors $\overrightarrow{P},\Delta \overrightarrow{P}$ with coefficients depending on various scalar products. The generic solution is simply

$$\begin{array}{c}\hfill \Delta \overrightarrow{P}=\lambda \overrightarrow{P};\end{array}$$

(47)

with the complex parameter $\lambda$ to be determined from the equation (46).

This solution is degenerate: the fixed point equation is satisfied for arbitrary complex $\lambda$.

The discontinuity vector $\Delta \overrightarrow{P}$ aligned with the principal value $\overrightarrow{P}$ corresponds to vanishing vorticity in (19), leading to a trivial solution of the loop equation $\mathsf{\Psi }[\gamma ,C]=1$.

We are left with the decaying turbulence scenario () as the only remaining physical solution.

3. Fractal curve in complex space

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

$$\begin{array}{ccc}& & \Delta \overrightarrow{F}\stackrel{?}{=}\lambda \overrightarrow{F};\hfill \end{array}$$

(48)

$$\begin{array}{ccc}& & {\lambda }^2{\overrightarrow{F}}^2-1\stackrel{?}{=}{\lambda }^2{\overrightarrow{F}}^2;\hfill \end{array}$$

(49)

There is, however, another solution where the vectors $\Delta \overrightarrow{F},\overrightarrow{F}$ are not aligned. This solution requires the following relations

$$\begin{array}{ccc}& & {(\Delta \overrightarrow{F})}^2=1;\hfill \end{array}$$

(50a)

$$\begin{array}{ccc}& & {(2\overrightarrow{F}·\Delta \overrightarrow{F}-ı\gamma )}^2+{\gamma }^2=4{\overrightarrow{F}}^2\hfill \end{array}$$

(50b)

These relations are very interesting. The complex numbers reflect irreversibility, and lack of alignment leads to vorticity distributed along the loop.

$$\begin{array}{ccc}& & {\mathsf{\Omega }}_{\alpha \beta }=\frac{ı}{2(t+t_0)}\left(F_{\alpha }\Delta F_{\beta }-F_{\beta }\Delta F_{\alpha }\right);\hfill \end{array}$$

(51)

$$$$

(52)

Note that this relation is parametricaly invariant (being local and independent of $\theta$).

Also, note that this complex vector $\overrightarrow{F}\left(\theta \right)$ is dimensionless, and the fixed point equation (Section 3.1) is completely universal, up to a single dimensionless parameter $\gamma$. 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 $\Delta \overrightarrow{F}=0$, so these discontinuities must be present at every $\theta$. We approximate this fractal curve by a polygon, with piecewise constant $\overrightarrow{F}\left(\theta \right)$ and N gaps $\Delta \overrightarrow{F}\left({\theta }_k\right)$ at equidistant angles ${\theta }_k=2\pi k/N$. If the limit $N\to \infty$ exists, we get the desired fractal curve $\overrightarrow{F}\left(\theta \right)$ with a discontinuity at every $\theta$.

One can build such a solution as a limit of the Markov process by the following method. Start with a complex vector $\overrightarrow{F}(\theta =0)={\overrightarrow{F}}_0$.

We compute the next values ${\overrightarrow{F}}_k=\overrightarrow{F}\left(\frac{2\pi k}{N}\right)$ from the discontinuity equations (Section 3.1).

$$\begin{array}{ccc}& & {\left({\overrightarrow{F}}_{k+1}-{\overrightarrow{F}}_k\right)}^2=1;\hfill \end{array}$$

(53a)

$$\begin{array}{ccc}& & {\left({\overrightarrow{F}}_{k+1}^2-{\overrightarrow{F}}_k^2-ı\gamma \right)}^2+{\gamma }^2={\left({\overrightarrow{F}}_{k+1}+{\overrightarrow{F}}_k\right)}^2\hfill \end{array}$$

(53b)

Naturally, this polygon at finite N does not correspond to the decaying turbulence; one must go to the local limit $N\to \infty$. As we shall see, this is where we find interesting critical phenomena.

3.2. Constraints imposed on a random step

A solution to these equations can be represented using a complex vector ${\overrightarrow{q}}_k$ subject to two complex constraints

$$\begin{array}{ccc}& & {\overrightarrow{q}}_k^2=1;\hfill \end{array}$$

(54a)

$$\begin{array}{ccc}\hfill & & {\left(2{\overrightarrow{F}}_k·{\overrightarrow{q}}_k-ı\gamma \right)}^2=4{\overrightarrow{F}}_k^2+\gamma (2ı-\gamma )\hfill \end{array}$$

(54b)

after which we can find the next value

$$\begin{array}{ccc}& & {\overrightarrow{F}}_{k+1}={\overrightarrow{F}}_k+{\overrightarrow{q}}_k;\hfill \end{array}$$

(55)

We assume N steps, each with the angle shift $\Delta \theta =\frac{2\pi }{N}$.

This recurrent sequence is a Markov process because each step only depends on the current position ${\overrightarrow{F}}_k$. The closure requirement ${\overrightarrow{F}}_N={\overrightarrow{F}}_0$ makes it a periodic Markov process, with the period N.

This requirement represents a nonlinear restriction on all the variables ${\overrightarrow{F}}_k$, which we discuss below.

With this discretization, the circulation can be expressed in terms of these steps

$$\begin{array}{c}\hfill \oint \overrightarrow{F}\left(\theta \right)·d\overrightarrow{C}\left(\theta \right)=-\oint \overrightarrow{C}\left(\theta \right)·d\overrightarrow{F}\left(\theta \right)\Rightarrow -\sum _{k=0}^{N-1}\frac{{\overrightarrow{C}}_{k+1}+{\overrightarrow{C}}_k}{2}·{\overrightarrow{q}}_k\end{array}$$

(56)

Note that the complex unit vector is not defined with the Euclidean metric in $2d$ dimensions $〈\overrightarrow{A},\overrightarrow{B}〉=\mathbf{Re}\overrightarrow{A}·\mathbf{Re}\overrightarrow{B}+\mathbf{Im}\overrightarrow{A}·\mathbf{Im}\overrightarrow{B}$. Instead, we have a complex condition

$$\begin{array}{c}\hfill {\overrightarrow{q}}^2=1\end{array}$$

(57)

which leads to two conditions between real and imaginary parts

$$\begin{array}{ccc}& & {\left(\mathbf{Re}\overrightarrow{q}\right)}^2=1+{\left(\mathbf{Im}\overrightarrow{q}\right)}^2;\hfill \end{array}$$

(58)

$$\begin{array}{ccc}& & \mathbf{Re}\overrightarrow{q}·\mathbf{Im}\overrightarrow{q}=0;\hfill \end{array}$$

(59)

In d dimensions, there are $d-1$ complex parameters of the unit vector; with an extra constraint in (53a), 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 (53a) for the projection ${\overrightarrow{F}}_k·{\overrightarrow{q}}_k$.

3.3. Closure condition

At the last step, $k=N-1$, we need to get a closed loop ${\overrightarrow{F}}_N={\overrightarrow{F}}_0$. This is one more constraint on the complex vectors ${\overrightarrow{q}}_0,\cdots {\overrightarrow{q}}_{N-1}$

$$\begin{array}{ccc}& & \sum _0^{N-1}{\overrightarrow{q}}_k=0;\hfill \end{array}$$

(60)

We use this complex vector constraint to fix some of the remaining parameters.

Due to the closure of the space loop $\overrightarrow{C}\left(\theta \right)$, the global translation of the momentum loop $\overrightarrow{P}\left(\theta \right)$ 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.

3.4. Mirror pairs of solutions

Return to the general study of the discrete loop equations (Section 3.2).

There is a trivial solution to these equations at any even N

$$\begin{array}{ccc}& & {\overrightarrow{f}}_k=\frac{{(-1)}^k\overrightarrow{q}}{2};\hfill \end{array}$$

(61)

$$\begin{array}{ccc}& & {\overrightarrow{q}}^2=1;\hfill \end{array}$$

(62)

We reject this solution as unphysical: the corresponding vorticity equals zero, as all the vectors ${\overrightarrow{f}}_k$ are aligned.

Our set of equations has certain mirror reflection symmetry

$$\begin{array}{c}\hfill {\overrightarrow{F}}_k\leftrightarrow {\overrightarrow{F}}_{N-k}^{*}\end{array}$$

(63)

Thus, the complex solutions come in mirror pairs ${\overrightarrow{F}}_k,{\overrightarrow{F}}_{N-k}^{*}$. The real solutions are only a particular case of the above trivial solution with real $\overrightarrow{q}$.

Each nontrivial solution represents a periodic random walk in complex vector space ${\mathbb{C}}^d$. The complex unit step ${\overrightarrow{q}}_k\in {\mathbb{C}}^d$ depends on the current position ${\overrightarrow{F}}_k\in {\mathbb{C}}^d$, or, equivalently, on the initial position ${\overrightarrow{F}}_0$ plus the sum of the preceding steps.

We are interested in the limit of infinitely many steps $N\to \infty$, corresponding to a closed fractal curve with a discontinuity at every point.

3.5. 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\to \infty$ sign variables, like an Ising model, plus an arbitrary global rotation matrix $\widehat{O}\in SO\left(d\right)$ and a global parameter $\beta$, 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.

4. Exact analytic solution

4.1. Random walk on a circle

How could a complex curve describe real circulation? This magic trick is possible if the imaginary part of $\overrightarrow{P}\left(\theta \right)$ does not depend on $\theta$.

Such an imaginary term will drop after integration over closed loop $\overrightarrow{C}\left(\theta \right)$.

We have found a family of such solutions [32] of our recurrent equation (Section 3.1) for arbitrary N

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

(64)

Here $\overrightarrow{w}\in {\mathbb{S}}^{d-3}$ is a unit vector.

The angles ${\alpha }_k$ must satisfy recurrent relation

$$\begin{array}{ccc}\hfill & & {\alpha }_{k+1}={\alpha }_k+{\sigma }_k\beta ;\hfill \end{array}$$

(65)

$$\begin{array}{ccc}& & {\alpha }_N-{\alpha }_0=0\left(\mathrm{mod}2\pi \right);\hfill \end{array}$$

(66)

$$\begin{array}{ccc}& & {\sigma }_k^2=1\hfill \end{array}$$

(67)

This sequence with arbitrary signs ${\sigma }_k=\pm 1$ solves recurrent equation (Section 3.1) independently of $\gamma$.

The closure condition requires certain relations between these numbers.

The main condition is that $\beta$ must be a rational fraction of $2\pi$:

$$\begin{array}{c}\hfill \beta =\frac{2\pi p}{q};0<p<q<N;q\ne 2;\end{array}$$

(68)

The $q=2$ case is eliminated. It corresponds to $p=1,\beta =\pi ,{\overrightarrow{F}}_k=\infty$.

In that case, the periodic solution for ${\alpha }_k$ will have correspond to the following set of ${\sigma }_k$

$$\begin{array}{ccc}& & \sigma ={\{1,\cdots ,1,-1,\cdots ,-1\}}_{perm};\hfill \end{array}$$

(69)

This array has $N_{+}$ positive values and $N_{-}$ negative values where

$$\begin{array}{ccc}& & N_{\pm }=\lfloor \frac{N\pm rq}{2}\rfloor ;\hfill \end{array}$$

(70)

$$\begin{array}{ccc}& & N_{+}+N_{-}=N;\hfill \end{array}$$

(71)

$$\begin{array}{ccc}& & N_{+}-N_{-}=rq-(N-rq)\left(\mathrm{mod}2\right);\hfill \end{array}$$

(72)

The symbol $perm$ stands for a random permutation of the array, which preserves its sum

$$\begin{array}{c}\hfill \sum {\sigma }_k=(N_{+}-N_{-})\end{array}$$

(73)

This sum must be a multiple of q for periodicity, which leads to another restriction

$$\begin{array}{c}\hfill (N-rq)\left(\mathrm{mod}2\right)=0\end{array}$$

(74)

In other words, $rq$ must have the same parity as N.

These properties lead to periodicity

$$\begin{array}{c}\hfill {\alpha }_N-{\alpha }_0=\beta \sum {\sigma }_k=2\pi pr\end{array}$$

(75)

At fixed denominator q, the winding number r can take the values

$$\begin{array}{ccc}& & r=\left\{-r_{max},\cdots ,0,\cdots ,r_{max}\right\};\hfill \end{array}$$

(76)

$$\begin{array}{ccc}& & r_{max}=N/q\hfill \end{array}$$

(77)

The sequence with all spins flipped: ${\sigma }_k\Rightarrow -{\sigma }_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\Rightarrow -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=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$

$$\begin{array}{c}\hfill W(N,p,q,r)=2^{-N}\left({\displaystyle \genfrac{}{}{0pt}{}{N}{N_{+}}}\right)=\frac{N!}{2^N{N}_{+}!N_{-}!}\end{array}$$

(78)

The next section discusses this ensemble of rational numbers and its statistics at $N\to \infty$.

Given the rational number $\frac{p}{q}$, we can generate the sequence of angles $\sum \pm \beta$, adding to a $2\pi$ multiple.

The random walk step ${\overrightarrow{q}}_k={\overrightarrow{F}}_{k+1}-{\overrightarrow{F}}_k$ is a real unit vector in this solution

$$\begin{array}{ccc}& & {\overrightarrow{q}}_k={\sigma }_k\{-sin{\delta }_k,\overrightarrow{w}cos{\delta }_k,0\};\hfill \end{array}$$

(79)

$$\begin{array}{ccc}& & {\delta }_k={\alpha }_k+\frac{\beta {\sigma }_k}{2}\hfill \end{array}$$

(80)

The direction of this vector is not random, though; in addition to the random sign ${\sigma }_k$ and random unit vector $\overrightarrow{w}$ in $d>3$ dimensions, its direction depends on the previous position ${\alpha }_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 $\beta =\frac{2\pi 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 ${\overrightarrow{F}}_k$ by a constant vector:

$$\begin{array}{c}\hfill {\overrightarrow{F}}_k\Rightarrow {\overrightarrow{F}}_k+\overrightarrow{V}\end{array}$$

(81)

This gauge invariance follows from the closure of the loop C: any constant term in $\overrightarrow{F}\left(\theta \right)$ yields zero after integration $\oint d\overrightarrow{C}=0$, or summation $\sum \Delta \overrightarrow{C}=0$.

Using this invariance, we can drop the last component of ${\overrightarrow{F}}_k$ so that they become real vectors

$$\begin{array}{ccc}& & {\overrightarrow{F}}_k\Rightarrow \frac{1}{2}csc\left(\frac{\beta }{2}\right)\left\{cos\left({\alpha }_k\right),sin\left({\alpha }_k\right)\overrightarrow{w},0\right\};\hfill \end{array}$$

(82)

The vorticity operator in this gauge will become a purely imaginary vector in the z direction:

$$\begin{array}{c}\hfill {\overrightarrow{\omega }}_k=\left\{0,0,\frac{ı{\sigma }_k}{2}cot\left(\frac{\beta }{2}\right)\right\}\end{array}$$

(83)

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 $\overrightarrow{r}$ in an "inertial interval," is finite and real after integration over the global rotation matrix. Its limit at $\overrightarrow{r}\to 0$ may be singular, so the anomalous dissipation may emerge.

This discrete set $N,p,q,r,{\sigma }_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,{\sigma }_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\to \infty$. 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 the properties of such an ensemble can be studied by a number theory methods.

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〉\to \infty$.

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 [33,34] is dual to the continuum theory of quantum Gravity.

4.2. The Euler ensemble

The discrete set of fractions $\frac{p}{q}$ with denominator $q<N$ is well known in the number theory [35], starting with Euler.

However, our extra restriction ${\sum }_l{\sigma }_l\left(\mathrm{mod}q\right)=0$ ties the set of fractions to the set $Z_2^{\otimes N}$ of N Ising variables. We will study the statistics of the corresponding ensemble, which we call the Euler ensemble, honoring my idol, Leonard Euler.

He never thought that his $\phi$ functions and his equation for ideal fluid would meet in theoretical physics, but a genius creates subconsciously.

We distinguish between big Euler ensembles and small Euler ensembles.

The big ensemble $\mathbb{E}\left(N\right)$ assigns equal weight for each element of the large set

$$\begin{array}{ccc}& & \mathrm{variables}\text{:}p,q,r,\left\{{\sigma }_1\cdots {\sigma }_N\right\}\hfill \end{array}$$

(84)

$$\begin{array}{ccc}& & 0<p<q<N;\hfill \end{array}$$

(85)

$$\begin{array}{ccc}& & -N\le qr\le N;\hfill \end{array}$$

(86)

$$\begin{array}{ccc}& & {\sigma }_1,\cdots {\sigma }_N=\pm 1;\hfill \end{array}$$

(87)

$$\begin{array}{ccc}& & W_N\left(p,q,r,\left\{{\sigma }_1\cdots {\sigma }_N\right\}\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}(p,q)=1,\sum _i{\sigma }_i=qr\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(88)

The small ensemble $\mathcal{E}\left(N\right)$ results from averaging the big ensemble over the $\sigma$ variables

$$\begin{array}{ccc}& & \mathrm{variables}\text{:}p,q,r,\hfill \end{array}$$

(89)

$$\begin{array}{ccc}& & 0<p<q<N;\hfill \end{array}$$

(90)

$$\begin{array}{ccc}& & -N\le qr\le N;\hfill \end{array}$$

(91)

$$\begin{array}{ccc}& & w_N\left(p,q,r\right)=\left\{\begin{array}{cc}{w}_N(q,r)\hfill & \mathrm{if}(p,q)=1\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(92)

$$\begin{array}{ccc}& & w_N(q,r)=\left\{\begin{array}{cc}{2}^{-N}\left({\displaystyle \genfrac{}{}{0pt}{}{N}{(N+qr)/2}}\right)\hfill & \mathrm{if}2\left|\right(N-qr)\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(93)

The binomial coefficients count the number of states with ${\sum }_i{\sigma }_i=qr$ among all $2^N$ states of the set of ${\sigma }_i=\pm 1,i=1,\cdots N$.

Note that we divided the statistical weights (the number of allowed configurations) by the total number $2^N$ of spin configurations. This normalization makes $w_N\left(q,r\right)$ the probability of finding the values $q,r$ in the big Euler ensemble with random ${\sigma }_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 ${\sum }_r{w}_N(q,r)$.

At given $q<N$ the allowed numbers of p are all integers $0<p<q$ such that $\mathbf{gcd}(p,q)=1$. The famous Euler’s totient $\phi \left(q\right)$ [35] counts such numbers.

$$\begin{array}{ccc}\hfill & & Prob[q<Q]=\frac{Z(Q,N)}{Z(N,N)}\hfill \end{array}$$

(94)

$$\begin{array}{ccc}& & Z(Q,N)=\sum _{2<q<Q}\phi \left(q\right)\sum _r{w}_N(q,r);\hfill \end{array}$$

(95)

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

(96)

Thus, we can relate the big Euler ensemble average of some function to the small ensemble average as follows

$$\begin{array}{ccc}\hfill & & {〈F(\cdots )〉}_{\mathbb{E}\left(N\right)}={〈{〈F(\cdots )〉}_{\sigma }〉}_{\mathcal{E}\left(N\right)};\hfill \end{array}$$

(97)

$$\begin{array}{ccc}& & {〈F\left(p,q,r,\left\{{\sigma }_1\cdots {\sigma }_N\right\}\right)〉}_{\sigma }\equiv 2^{-N}\sum _{\begin{array}{c}{\sigma }_i=\pm 1\\ {\sum }_i{\sigma }_i=qr\end{array}}F\left(p,q,r,\left\{{\sigma }_1\cdots {\sigma }_N\right\}\right);\hfill \end{array}$$

(98)

$$\begin{array}{ccc}& & {〈F(p,q,r,N)〉}_{\mathcal{E}\left(N\right)}=\frac{{\displaystyle \sum _{2<q<N}\sum _{\begin{array}{c}p=1\\ (p,q)=1\end{array}}^{q-1}\sum _{\begin{array}{c}r\\ 2\left|\right(N-qr)\end{array}}F(p,q,r,N)}}{Z(N,N)};\hfill \end{array}$$

(99)

The advantage of this representation is that we can first average by the $\sigma$ 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\to \infty$ 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

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\mathrm{odd}r,\mathrm{odd}q\hfill & \mathrm{odd}N\hfill \\ \mathrm{integer}r,\mathrm{even}q\hfill & \mathrm{even}N\hfill \\ \mathrm{even}r,\mathrm{odd}q\hfill & \mathrm{even}N\hfill \end{array}\right.\end{array}$$

(100)

The ratios

$$\begin{array}{c}\hfill R_N\left(q\right)=\frac{{\sum }_{r\ne 0}{w}_N(q,r)}{w_N(q,0)}\end{array}$$

(101)

were computed for $N=1000,1001$ (see Figure 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\ne 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 q- independent limit

$$\begin{array}{c}\hfill w_N(q,0)=2^{-N}\left({\displaystyle \genfrac{}{}{0pt}{}{N}{\frac{N}{2}}}\right)\to \sqrt{\frac{2}{\pi N}}\end{array}$$

(102)

which does not restrict the sum of Euler totients ${\sum }_q\phi \left(q\right)$.

Therefore, for even N

$$\begin{array}{c}\hfill Z{(N,N)}_{\mathrm{even}N}\to \sqrt{\frac{2}{\pi N}}\sum _{2<q<N}\phi \left(q\right)\to \frac{3\sqrt{2}{N}^{3/2}}{{\pi }^{5/2}}\end{array}$$

(103)

For odd N, the leading term is missing, so we have to sum the terms with $r\ne 0$.

$$\begin{array}{c}\hfill Z{(N,N)}_{\mathrm{odd}N}=\sum _{2<q<N}\phi \left(q\right)\sum _{r>0}2w_N(q,r)\end{array}$$

(104)

This term converges at $q\ll N$ where ${\sum }_{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 are left to direct numerical computation of this sum and extrapolation to $N\to \infty$.

We computed it numerically in Mathematica^®and fitted it to $\sqrt{N}$ times a power of $logN$ (see Figure 2).

$$\begin{array}{c}\hfill Z{(N,N)}_{\mathrm{odd}N}\approx 1.2072\sqrt{N}{log}^{0.078776}\left(N\right)\end{array}$$

(105)

It is a challenge to the number theory to produce exact asymptotic behavior, replacing my fitted "law."

4.3. 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^{-\mu N}$, then N is treated as the rest of the variables. The chemical potential $\mu$ 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 $\mu$ compared to the historical definition: our $\mu \to +0$, so the opposite sign would be inconvenient.

The partition function and the ensemble averages in the grand canonical ensemble are

$$\begin{array}{ccc}& & Z\left(\mu \right)=\sum _{N}Z(N,N)e^{-\mu N};\hfill \end{array}$$

(106)

$$\begin{array}{ccc}& & {〈F(p,q,r,N)〉}_{\mathcal{E}\left(\mu \right)}=\frac{{\displaystyle \sum _N{e}^{-\mu N}\sum _{2<q<N}\sum _{\begin{array}{c}p=1\\ (p,q)=1\end{array}}^{q-1}\sum _{\begin{array}{c}r\\ 2\left|\right(N-qr)\end{array}}F(p,q,r,N)}}{Z\left(\mu \right)};\hfill \end{array}$$

(107)

With the grand canonical ensemble, the ambiguity of the local limit disappears. At the critical point $\mu \to 0$, the even N dominate, with the following result

$$\begin{array}{c}\hfill Z\left(\mu \right)\to \frac{1}{2}\frac{9}{2\sqrt{2}{\pi }^2{\mu }^{5/2}}\end{array}$$

(108)

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 $\mathbb{E}\left(\mu \right),\mathcal{E}\left(\mu \right)$ in the thermodynamic limit $\mu \to 0$, where we find the critical phenomena.

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 ${\mu }_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.

5. Correlation functions

In this section and the rest of the paper, we only consider the three-dimensional space we live in. 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.

5.1. General formulas

The simplest observable quantities we can extract from the loop functional are the vorticity correlation functions [18], corresponding to the loop C backtracking between two points in space ${\overrightarrow{r}}_1=0,{\overrightarrow{r}}_2=\overrightarrow{r}$, see Figure 3. The vorticity operators are inserted at these two points.

The correlation function reduces to the following average over the ensemble $\mathbb{E}\left(N\right)$ of our random curves in complex space.

$$\begin{array}{ccc}\hfill & & 〈\overrightarrow{\omega }\left(\overrightarrow{0}\right)·\overrightarrow{\omega }\left(\overrightarrow{r}\right)〉=\frac{{〈{\displaystyle \sum _{0\le n<m<N}{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_{n}exp\left(ı\overrightarrow{\rho }·\left({\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right)\right)}〉}_{\mathbb{E}\left(\mu \right)}}{4{(t+t_0)}^2};\hfill \end{array}$$

(109)

$$\begin{array}{ccc}& & {\overrightarrow{S}}_{n,m}=\frac{{\sum }_{k=n}^{m-1}{\overrightarrow{F}}_k}{(m-n)\left(\mathrm{mod}N\right)};\hfill \end{array}$$

(110)

$$\begin{array}{ccc}& & \overrightarrow{\rho }=\frac{\overrightarrow{r}}{2\sqrt{\nu (t+t_0)}};\hfill \end{array}$$

(111)

The averaging $〈\cdots 〉$ in these formulas involves group integration ${\int }_{O\left(3\right)}d\widehat{O}$ with $\overrightarrow{\rho }\Rightarrow \widehat{O}·\overrightarrow{\rho }$.

$$\begin{array}{c}\hfill {〈H(\overrightarrow{\rho }·\widehat{\mathsf{\Omega }}·\overrightarrow{F})〉}_{O\left(3\right)}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\int }_{-1}^{1}dzH\left(|\overrightarrow{\rho }\left|\right|\overrightarrow{F}|z\right)\end{array}$$

(112)

Let us explain the origin of summation over two positions $n,m$ of the points ${\overrightarrow{r}}_1=0,{\overrightarrow{r}}_2=\overrightarrow{r}$ on the discreet loop $\overrightarrow{C}\left(\theta \right)$.

There is a degree of freedom we did not specify until now, namely, the reparametrization of the momentum loop $\overrightarrow{P}\left(\theta \right)$.

The loop equations (24) are invariant under the one-dimensional diffeomorphisms ( or reparametrizations)

$$\begin{array}{ccc}& & \overrightarrow{P}\left(\theta \right)\Rightarrow \overrightarrow{P}\left(f\left(\theta \right)\right);\hfill \end{array}$$

(113)

$$\begin{array}{ccc}& & f^{\prime }\left(\theta \right)>0,f\left(2\pi \right)=f\left(0\right)+2\pi \hfill \end{array}$$

(114)

Thus, the general solution involves an arbitrary monotonous function $f\left(\theta \right)$, 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 ${\theta }_1,{\theta }_2$ of the vorticity insertion points on a curve $\overrightarrow{C}\left(\theta \right)$. In the continuum theory, this would be an ordered Lebesgue integration (diffeomorphisms preserve the ordering of points on a curve)

$$\begin{array}{c}\hfill \underset{-\pi }{\overset{\pi }{\int \int }}\mathsf{\Theta }({\theta }_2-{\theta }_1)\left(\overrightarrow{P}\left({\theta }_1\right)\times d\overrightarrow{P}\left({\theta }_1\right)\right)\otimes \left(\overrightarrow{P}\left({\theta }_2\right)\times d\overrightarrow{P}\left({\theta }_2\right)\right)\end{array}$$

(115)

and in our case of piecewise constant curve with discontinuities $\Delta {\overrightarrow{P}}_k$, it becomes an ordered sum

$$\begin{array}{c}\hfill \sum _{0\le n<m<N}\left({\overrightarrow{P}}_m\times \Delta {\overrightarrow{P}}_m\right)\otimes \left({\overrightarrow{P}}_n\times \Delta {\overrightarrow{P}}_n\right)\end{array}$$

(116)

The discontinuity $\Delta \overrightarrow{P}\left(\theta \right)$ stays finite in the continuum limit $N\to \infty$. 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 (63) does not depend on the point on a circle. Therefore, it contributes a constant term into ${\overrightarrow{S}}_{m,n}$ which cancels in the difference ${\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}$ in the exponential, as it should.

Let us look at the correlation function (108).

First, we expand and simplify the dot product involved

$$\begin{array}{c}\hfill {\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n=\frac{-{\sigma }_m{\sigma }_n}{4}{cot}^2\left(\frac{\beta }{2}\right)\end{array}$$

(117)

The terms $S_{m,n},S_{n,m}$ in (108) have the following form

$$\begin{array}{ccc}& & \overrightarrow{\rho }·{\overrightarrow{S}}_{n,m}=exp\left(ı\beta {\sigma }_n\right)A_{n,m};\hfill \end{array}$$

(118)

$$\begin{array}{ccc}& & \overrightarrow{\rho }·{\overrightarrow{S}}_{m,n}=exp\left(ı\beta {\sigma }_m\right)A_{m,n+N};\hfill \end{array}$$

(119)

$$\begin{array}{ccc}& & A_{n,m}=\mathbf{Re}{R}^{★}\frac{{\displaystyle \sum _{k=n}^{m-1}exp\left(ı{\alpha }_{k,n}\right)}}{2sin(\beta /2)(m-n)};\hfill \end{array}$$

(120)

$$\begin{array}{ccc}& & {\alpha }_{k,n}=\beta \sum _{\begin{array}{c}l=0\\ l\ne n\end{array}}^k{\sigma }_l;\hfill \end{array}$$

(121)

$$\begin{array}{ccc}& & R={\rho }_x+ı{\rho }_y\hfill \end{array}$$

(122)

5.2. Critical phenomena in statistical limit

Now we are prepared to average over spin variables ${\sigma }_l,l=0\cdots N-1$.

This expression singles out the variables ${\sigma }_n,{\sigma }_m$ so we can sum over these two variables, leaving the rest of ${\sigma }_l$ free, except for a constraint ${\sum }_{\sigma }=rq$. This constraint can be implemented as a discrete Fourier integral:

$$\begin{array}{c}\hfill \delta [qr-\sum \sigma ]=\oint \frac{d\omega }{2\pi }{e}^{ı\omega \left(qr-\sum \sigma \right)}\end{array}$$

(123)

We start with

$$\begin{array}{ccc}& & \frac{-1}{4}{cot}^2\left(\frac{\beta }{2}\right)\hfill \\ & & {〈{\sigma }_m{\sigma }_n\oint \frac{d\omega }{2\pi }{e}^{ı\omega \left(qr-\sum \sigma \right)}exp\left(ıexp\left(ı\beta {\sigma }_n\right)A_{n,m}-ıexp\left(ı\beta {\sigma }_m\right)A_{m,n+N}\right)〉}_{{\sigma }_l=\pm 1};\hfill \end{array}$$

(124)

Then we take the $\omega$ integral out of the sum:

$$\begin{array}{ccc}& & \oint \frac{d\omega }{2\pi }{e}^{ı\omega qr}\hfill \\ & & {〈{\sigma }_m{\sigma }_n\oint \frac{d\omega }{2\pi }{e}^{-ı\omega \sum \sigma }exp\left(ıexp\left(ı\beta {\sigma }_n\right)A_{n,m}-ıexp\left(ı\beta {\sigma }_m\right)A_{m,n+N}\right)〉}_{{\sigma }_l=\pm 1};\hfill \end{array}$$

(125)

This expression can be readily averaged over two variables ${\sigma }_m,{\sigma }_n$, which reduces to the sum of four terms with ${\sigma }_n,{\sigma }_m=\pm 1$. Keeping only the factors depending on ${\sigma }_n,{\sigma }_m$

$$\begin{array}{ccc}& & {〈{\sigma }_m{\sigma }_n{e}^{-ı\omega ({\sigma }_n+{\sigma }_m)}exp\left(ıexp\left(ı\beta {\sigma }_n\right)A_{n,m}-ıexp\left(ı\beta {\sigma }_m\right)A_{m,n+N}\right)〉}_{{\sigma }_n,{\sigma }_m=\pm 1}=\hfill \end{array}$$

$$\begin{array}{ccc}& & \mathsf{\Phi }\left(A_{n,m},\omega ,\beta \right)\mathsf{\Phi }\left(-A_{m,n+N},\omega ,\beta \right);\hfill \end{array}$$

(126)

$$\begin{array}{ccc}& & \mathsf{\Phi }(A,\omega ,\beta )=sin(\omega -isin(A\beta \left)\right)(sin(cos\left(A\beta \right))-icos(cos\left(A\beta \right)\left)\right)\hfill \end{array}$$

(127)

The next step would be averaging over the remaining variables ${\sigma }_l,l\ne n,l\ne 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 $\sim m\sim n\sim N\to \infty$. We also have $\beta \sim 1/\sqrt{N}\to 0$ in that limit.

We are now considering the unconstrained distribution over ${\sigma }_l$, as the constraint is implemented via the discrete Fourier integral.

Let us compute the mean and variance of $U_{n,m}={\sum }_{k=n}^{m-1}exp\left(ı{\alpha }_{k,n}\right)$

$$\begin{array}{ccc}& & 〈U_{nm}〉=\sum _{k=n}^{m-1}{cos}^{k-1}\beta =\frac{cos{\left(\beta \right)}^n-cos{\left(\beta \right)}^m}{cos\beta (1-cos\beta )};\hfill \end{array}$$

(128)

$$\begin{array}{ccc}& & 〈{\left|U_{nm}-〈U_{nm}〉\right|}^2〉=\sum _{l=n}^{m-1}\sum _{k=n}^{m-1}〈{\xi }_k{\xi }_l^{★}〉=\sum _{l=n}^{m-1}\sum _{k=n}^{m-1}\left(cos{\beta }^{|l-k|}-cos{\beta }^{l}cos{\beta }^k\right);\hfill \end{array}$$

(129)

$$\begin{array}{ccc}& & {\xi }_r=exp\left(ı{\alpha }_{r,n}\right)-〈exp\left(ı{\alpha }_{r,n}\right)〉\hfill \end{array}$$

(130)

In the relevant critical region $\beta \to 0,m=x/{\beta }^2,n=y/{\beta }^2,x>y$, the ratio of the variance of $U_{n,m}$ to the square of its mean tends to a finite limit

$$\begin{array}{c}\hfill \frac{〈{\left|U_{nm}-〈U_{nm}〉\right|}^2〉}{{〈U_{nm}〉}^2}\to \frac{2e^{\frac{y-x}{2}}+x-y-2}{{\left(e^{-\frac{x}{2}}-e^{-\frac{y}{2}}\right)}^2}-1\end{array}$$

(131)

This function looks singular, but it is positive and finite in the allowed region $0<y<x$ ( see Figure 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.

This complication is bad news for attempted direct computation but 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 $\mu$, relating it to the functions of the number theory.

5.3. 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 $\overrightarrow{r}=0$ in (108) (which is possible at finite N) we find the following relation

$$\begin{array}{ccc}\hfill & & 〈\overrightarrow{\omega }{\left(\overrightarrow{0}\right)}^2〉=\frac{1}{4{(t+t_0)}^2}\sum _{0\le n<m<N}{〈{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n〉}_{\mathbb{E}\left(N\right)}\hfill \end{array}$$

(132)

where the big Euler ensemble average ${〈\cdots 〉}_{\mathbb{E}\left(N\right)}$ denotes averaging over $p,q,r$ and the Ising variables $\sigma$ subject to the global constraint $\sum {\sigma }_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 ${\sigma }_k$, which is rather simple.

Let us use the explicit expression (116) for the dot product ${\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n$.

Now, using the one-dimensional Fourier integral for the global constraint on $\sigma$, we get the unconstrained average:

$$\begin{array}{ccc}& & {〈{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n〉}_{\sigma }=-\frac{1}{4}{cot}^2\left(\frac{\beta }{2}\right){\int }_{-\pi }^{\pi }\frac{d\omega }{2\pi }{e}^{ı\omega qr}{〈{\sigma }_m{\sigma }_{n}exp\left(-ı\omega \sum _0^{N-1}{\sigma }_l\right)〉}_{\sigma }\hfill \end{array}$$

(133)

Averaging it over ${\sigma }_m,{\sigma }_n$ we find

$$\begin{array}{ccc}& & {〈exp\left(-ı\omega ({\sigma }_n+{\sigma }_m)\right){\sigma }_m{\sigma }_n〉}_{{\sigma }_m,{\sigma }_n}=-{sin}^2\left(\omega \right)\hfill \end{array}$$

(134)

Now, averaging over the remaining ${\sigma }_r,r\ne m,r\ne n$ is straightforward.

$$\begin{array}{c}\hfill {〈exp\left(-ı\omega \sum _{\begin{array}{c}l=0\\ l\ne n,l\ne m\end{array}}^{N-1}{\sigma }_l\right)〉}_{\sigma }={cos}^{N-2}\omega \end{array}$$

(135)

We arrive at the integral

$$\begin{array}{ccc}& & {〈{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n〉}_{\sigma }=\frac{{cot}^2\left(\frac{\beta }{2}\right)}{4}{\int }_{-\pi }^{\pi }\frac{d\omega }{2\pi }{e}^{ı\omega qr}{cos}^{N-2}\omega {sin}^2\omega \hfill \end{array}$$

(136)

This integral is calculable (see [36]):

$$\begin{array}{ccc}& & {〈{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n〉}_{\sigma }=\frac{(N-q^2{r}^2){cot}^2\left(\frac{\beta }{2}\right)}{2^N(N^2-q^2{r}^2)}\left({\displaystyle \genfrac{}{}{0pt}{}{N-2}{\frac{1}{2}(N+qr-2)}}\right)\hfill \end{array}$$

(137)

Summing over $0\le n<m<N$ yields $N(N-1)/2$, leading to our final answer

$$\begin{array}{ccc}\hfill & & 〈\overrightarrow{\omega }{\left(\overrightarrow{0}\right)}^2〉=\frac{1}{{(t+t_0)}^2}{〈2^{-N-5}S\left(q\right)\left(N-q^2{r}^2\right)\left({\displaystyle \genfrac{}{}{0pt}{}{N}{(N+qr)/2}}\right)〉}_{\mathcal{E}\left(\mu \right)};\hfill \end{array}$$

(138)

$$\begin{array}{ccc}& & S\left(q\right)=\sum _{\begin{array}{c}p=1\\ (p,q)=1\end{array}}^{q-1}{cot}^2\left(\frac{\pi p}{q}\right);\hfill \end{array}$$

(139)

We investigated this new function $S\left(q\right)$ in [37] and represented it in terms of so-called multitotients [38]

$$\begin{array}{ccc}\hfill & & S\left(q\right)=\frac{{\phi }_2\left(q\right)}{3}-{\phi }_1\left(q\right);\hfill \end{array}$$

(40)

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

(141)

$$\begin{array}{ccc}& & {\phi }_1\left(q\right)=\phi \left(q\right)\hfill \end{array}$$

(142)

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\left(100\right)=2360$. The same result using the multitotient formula takes 140 microseconds.

Here is the table of $S\left(q\right)$ for $2\le q\le 10$

$$\begin{array}{c}\hfill \left(\begin{array}{cc}2& 0\\ 3& \frac{2}{3}\\ 4& 2\\ 5& 4\\ 6& 6\\ 7& 10\\ 8& 12\\ 9& 18\\ 10& 20\end{array}\right)\end{array}$$

(143)

This $S\left(q\right)$ is an even integer for $q>3$, and here is a simple proof.

Proof. 

The second term $-\phi \left(q\right)$ in $S\left(q\right)$ is an even integer, and the first term can be rewritten as

$$\begin{array}{c}\hfill 1/3q^2\prod _{p|q}\left(1-\frac{1}{p^2}\right)=1/3\prod _{p|q}{p}^{2({\alpha }_p-1)}(p^2-1);\end{array}$$

(144)

where ${\alpha }_p\ge 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 one factor divisible by three, and it cannot be p as it is a prime number, with the only exception being $q=p=3$. Starting with $q=4$ there is always a factor of 3 either in $(2^2-1)$ or in $3^{2({\alpha }_3-2)}$ or in any other prime factors $(p^2-1)$ analyzed above. The divisibility by 2 is similar: both factors $p\pm 1$ for any prime $p>2$ are even, and any power $2^{2({\alpha }_2-1)}$ with ${\alpha }_2>1$ is divisible by 4. Therefore $p^2-1$ is divisible by 12, so that $\frac{{\phi }_2\left(q\right)}{3}$ is divisible by 4. Thus $S\left(q\right)$ is an even integer for $q>3$. □

The plot of the first 100000 values of $S\left(q\right)$ looks as follows (see Figure 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.

This formula 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 $\mu \to 0$.

5.4. The local limit of the energy dissipation

In the limit $N\to \infty$, the term with zero winding $r=0$ dominates the sum ${\sum }_r$, which is only possible for even N. Its asymptotic limit is

$$\begin{array}{c}\hfill 2^{-N-5}N\left({\displaystyle \genfrac{}{}{0pt}{}{N}{N/2}}\right)\to \frac{\sqrt{N}}{16\sqrt{2\pi }}\end{array}$$

(145)

The remaining terms, with $r\ne 0$ can be replaced by an integral of the asymptotic expansion of the exact expression

$$\begin{array}{c}\hfill \sqrt{N}{\int }_{-\infty }^{\infty }dz\frac{e^{-\frac{z^2}{2}}\left(1-z^2\right)}{16\sqrt{2\pi }}=0\end{array}$$

(146)

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

$$\begin{array}{c}\hfill {\int }_{-\infty }^{\infty }dz\frac{e^{-\frac{z^2}{2}}\left(z^2-1\right)\left(z^4-6z^2+3\right)}{192\sqrt{2\pi }\sqrt{N}}=0\end{array}$$

(147)

We do not need this term $O(1/\sqrt{N})$, as the dominant $r=0$ term grows as $\sqrt{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.

$$\begin{array}{ccc}& & {(t+t_0)}^2〈\overrightarrow{\omega }{\left(\overrightarrow{0}\right)}^2〉\to \frac{1}{2}\frac{{\displaystyle {\int }_0^{\infty }dN\sqrt{N}{e}^{-\mu N}\sum _{2<q<N}S\left(q\right)}}{64\sqrt{2\pi }Z\left(\mu \right)}\hfill \end{array}$$

(148)

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 [38]

$$\begin{array}{c}\hfill \sum _{m=2}^N{\phi }_l\left(m\right)\to \frac{N^{l+1}}{(l+1)\zeta (l+1)}\end{array}$$

(149)

where $\zeta \left(n\right)$ is the Riemann’s zeta function.

We obtain the following local limit of the energy decay in the grand canonical Euler ensemble

$$\begin{array}{ccc}\hfill & & {\partial }_{t}E=-\nu 〈\overrightarrow{\omega }{\left(\overrightarrow{0}\right)}^2〉=-\frac{\mathcal{B}}{{(t+t_0)}^2};\hfill \end{array}$$

(150)

$$\begin{array}{ccc}\hfill & & \mathcal{B}=\frac{\nu }{{\mu }^2}\frac{35{\pi }^2}{13824\zeta \left(3\right)};\hfill \end{array}$$

(151)

This energy dissipation stays finite in a turbulent limit $\nu \to 0$ provided

$$\mu \propto \sqrt{\nu }$$

(152)

Let us now estimate the odd N contribution to the enstrophy. Here, the sum of the spin average is dominated by $q\ll 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 average

$$\begin{array}{ccc}& & {\displaystyle \sum _{\mathrm{odd}N}{e}^{-\mu N}\sum _{\begin{array}{c}2<q<N\\ \mathrm{odd}q\end{array}}\left(\frac{{\phi }_2\left(q\right)}{3}-\phi \left(q\right)\right)\sum _{\begin{array}{c}r>0\\ \mathrm{odd}r\end{array}}{2}^{-N-4}\left(N-q^2{r}^2\right)\left({\displaystyle \genfrac{}{}{0pt}{}{N}{(N+qr)/2}}\right)}\hfill \end{array}$$

(153)

This contribution $\Delta E$ comes out negative. We fitted these results as $log(-\Delta E)\approx a+blogN+cloglogN$ (see Figure 6).

Here are the fit statistics

(154)

We have the following estimate for the odd N contribution to the enstrophy of the grand canonical Euler ensemble in the local limit:

$$\begin{array}{c}\hfill \frac{\Delta E}{Z\left(\mu \right)}\approx -\frac{0.068618{log}^{0.44420}\left(\frac{1}{\mu }\right)}{{\mu }^{0.45581}}\end{array}$$

(155)

The leading term grows $1/{\mu }^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 [17,18](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.

In the next publication with Maxim Bulatov [37], we present the simulations of the vorticity correlation in this theory using the supercomputer cluster in NYU Abu Dhabi.

6. Discussion

6.1. The Duality of Turbulence

We have presented an analytical solution of the Navier-Stokes loop equations for the Wilson loop in decaying turbulence as a functional of the shape and size of the loop in arbitrary dimension $d>2$.

The 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\to \infty$ sides.

This statistical system is similar to a one-dimensional Ising ring in an imaginary magnetic field $ı\beta =\frac{2ı\pi p}{q}$ and zero coupling constant. Some global observables, such as the 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\to \infty$ 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\left(-\mu N\right)$. 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 $\mu \to 0$. In this limit, we compute the partition function (107) and the expectation value of the energy dissipation (149).

These computations heavily rely on the number theory, particularly Lehmer’s multitotients ${\phi }_l\left(q\right)$, (139), generalizing [38] the Euler totient function.

6.2. Classical flow and quantum geometry

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 ground state 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[5,7,8,9] where the loop equations were derived first [3,4]. 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.

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.

Fortunately, 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 short-distance singularities. The anomalous dissipation is achieved in our solution via a completely different mechanism: large fluctuations of the fractal curve at $p\ll q$.

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[3,4].

6.3. Stokes-type functionals and vorticity correlations

There is no gauge invariance regarding the velocity field in fluid dynamics (though there is such invariance in the Clebsch variables [18]). 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 ${\overrightarrow{v}}_k=\frac{ı\overrightarrow{k}\times {\overrightarrow{\omega }}_k}{{\overrightarrow{k}}^2}$. Thus, the general solution for the Wilson loop functional $\mathsf{\Psi }[\gamma ,C]$ allows computing both vorticity and velocity correlation functions.

We demonstrated that in the last two sections by computing the correlation function of vorticity in coordinate space and the anomalous dissipation related to the singular limit of this correlation function at small distances.

6.4. 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 [3], 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 $\overrightarrow{v}$.

As we discussed in detail in [3,4,18], any Stokes-type functional $\mathsf{\Psi }[\gamma ,C]$ satisfying boundary condition at shrunk loop $\mathsf{\Psi }\left[0\right]=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 [2,18] (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 $\overrightarrow{P}\left(\theta \right)$.

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

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/\nu$ expansion in the Navier-Stokes equation.

In our solution, 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.

6.5. Conclusion

The exact solution for $\overrightarrow{P}\left(\theta \right)$ in decaying turbulence, which we have found in this paper, leads to the Stokes functional $\mathsf{\Psi }[\gamma ,C]$ satisfying the boundary value $\mathsf{\Psi }[\gamma ,0]=1$ at the shrunk loop.

Therefore, it describes a stochastic Navier-Stokes flow, corresponding to the degenerate fixed point of the Hopf equation for the generating functional of velocity circulation.

It sums up all the Wylde diagrams in the limit of vanishing random forces plus nonperturbative effects, which are missed in the Wylde functional integral, particularly the quantization of parameters of the solution.

These quantum effects follow from the exact equivalence of the Navier-Stokes statistics to the quantum mechanics in loop space [2].

Compared to the other critical phenomena, this theory is amazingly simple: It is not a field theory but a discrete statistical system isomorphic to the one-dimensional Ising ring in the presence of the imaginary quantized magnetic field.

Still, the solution is far from trivial and reveals unexpected relations between turbulence and number theory. In particular, the anomalous energy dissipation constant () is related to the prime factorization of large integers used in so-called multi-totients [38].

Whether this exact solution is realized in Nature remains to be seen.