Exact Solution of Decaying Turbulence

1. Loop equation

1.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}$$

(1)

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

(2)

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}$$

(3)

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}$$

(4a)

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

(4b)

$$\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}$$

(4c)

$$\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}$$

(4d)

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

(4e)

$$\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}$$

(4f)

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 $〈\dots 〉$ 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 & & {\mathsf{\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}$$

(5)

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

(6)

In the review, [2,17], 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 )}$$

(7)

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}$$

(8)

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 );$$

(9)

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}$$

(10)

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

1.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)〉$$

(11)

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 [20] and any other Navier-Stokes stationary solutions with a singularity at fixed loop C in space.

The functional derivative, acting on the exponential in (11), 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)$$

(12)

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 }$$

(13)

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 \mathrm{sgn}\left(\mu \right)P_{\alpha }\left(\theta +ϵ\mu \right)\hfill \end{array}$$

(14)

The vorticity (8) and velocity (9) 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 (14).

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}$$

(15)

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}$$

(16)

$$\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}$$

(17)

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

(18)

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}$$

(19)

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}$$

(20)

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}$$

(21)

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

1.3. Random global vorticity

Possible initial data for the reduced dynamics were suggested in the original papers [2,17]. 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}$$

(22)

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 })$$

(23)

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+\dots$$

(24)

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

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

(25)

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

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

(26)

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

The constant term in (24) 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}$$

(27)

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}$$

(28)

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 [17]. 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 (27).

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,17] 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}$$

(29)

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

(30)

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

(31)

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

(32)

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

(33)

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}$$

(34)

$$\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}$$

(35)

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$$

(36)

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}$$

(37)

$$\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 }{\mathsf{\Sigma }}_{\alpha \beta }\left(C\right);\hfill \end{array}$$

(38)

The tensor area $\mathsf{\Sigma }$ was defined in (6). 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.

1.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}$$

(39)

$$\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}$$

(40)

The Decay scenario in the nonlinear ODE (21) 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}{c}\overrightarrow{P}=\sqrt{\frac{\nu }{2(t+t_0)}}\frac{\overrightarrow{F}}{\gamma };\hfill \end{array}$$

(41)

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

(42)

The Fixed Point would correspond to the vanishing right side of the momentum loop equation (21). 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}$$

(43)

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 (43) 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}$$

(44)

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

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 (16), leading to a trivial solution of the loop equation $\mathsf{\Psi }[\gamma ,C]=1$.

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

2. Fractal curve in complex space

2.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}$$

(45)

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

(46)

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}$$

(47a)

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

(47b)

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}$$

(48)

Note that this relation is parametrically 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 (47) 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 discontinuities must be present at every $\theta$. 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 $\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 (47).

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

(49a)

$$\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}$$

(49b)

2.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}$$

(50a)

$$\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}$$

(50b)

after which we can find the next value

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

(51)

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}$$

(52)

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}$$

(53)

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}$$

(54)

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

(55)

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

2.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,\dots {\overrightarrow{q}}_{N-1}$

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

(56)

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.

2.4. Mirror pairs of solutions

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

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}$$

(57)

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

(58)

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}$$

(59)

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.

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

3. Exact analytic solution

3.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 [21] of our recurrent equation (49) 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}$$

(60)

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}$$

(61)

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

(62)

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

(63)

This sequence with arbitrary signs ${\sigma }_k=\pm 1$ solves recurrent equation (49) 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}$$

(64)

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,\dots ,1,-1,\dots ,-1\}}_{perm};\hfill \end{array}$$

(65)

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

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

(66)

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

(67)

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

(68)

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}$$

(69)

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}$$

(70)

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}$$

(71)

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

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

(72)

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

(73)

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}$$

(74)

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}$$

(75)

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

(76)

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}$$

(77)

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}$$

(78)

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}$$

(79)

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 that 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 [22,23] is dual to the continuum theory of quantum gravity.

3.2. The Euler ensemble

The discrete set of fractions $\frac{p}{q}$ with denominator $q<N$ is well known in the number theory [24], 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\dots {\sigma }_N\right\}\hfill \end{array}$$

(80)

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

(81)

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

(82)

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

(83)

$$\begin{array}{ccc}& & W_N\left(p,q,r,\left\{{\sigma }_1\dots {\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}$$

(84)

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}$$

(85)

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

(86)

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

(87)

$$\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}$$

(88)

$$\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}$$

(89)

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,\dots 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\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)$ [24] counts such numbers.

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

(90)

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

(91)

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

(92)

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

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

(93)

$$\begin{array}{ccc}& & {〈F\left(p,q,r,\left\{{\sigma }_1\dots {\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\dots {\sigma }_N\right\}\right);\hfill \end{array}$$

(94)

$$\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}$$

(95)

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}$$

(96)

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}$$

(97)

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}$$

(98)

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}$$

(99)

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}$$

(100)

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 can compute this sum numerically and extrapolate 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}$$

(101)

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

Recently, the number theory answered this challenge [25]. The result of their computation confirms the $\sqrt{N}$ factor but rejects the ${log}^{\alpha }\left(N\right)$ 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. Still, it does not contribute to the leading term of the grand canonical ensemble we use in the rest of the paper.

The results of [25] correspond to $Z{(N,N)}_{\mathrm{odd}N}\to \sqrt{N}/\sqrt{2\pi }$ without any factors of $logN$. This reminds us that mathematical laws must be derived theoretically rather than fitted to numerical data.

3.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 that 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}$$

(102)

$$\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}$$

(103)

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}$$

(104)

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.

With the grand canonical ensemble, the ergodic hypothesis becomes a weaker restriction: any smooth weight function $W\left(N\right)$ in the distribution will cancel in expectation values in the limit $\mu \to 0$. For example, the saddle point calculation yields

$$\begin{array}{c}\hfill \frac{{\int }_2^{\infty }dNW\left(N\right)N^{\alpha }F\left(N\right)exp\left(-\mu N\right)}{{\int }_2^{\infty }dNW\left(N\right)N^{\alpha }exp\left(-\mu N\right)}\to F(\alpha /\mu )\end{array}$$

(105)

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.

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

4.1. General formulas

The simplest observable quantities we can extract from the loop functional are the vorticity correlation functions [17], 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(\mu \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}$$

(106)

$$\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}$$

(107)

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

(108)

The averaging $〈\dots 〉$ 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}$$

(109)

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 (21) 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}$$

(110)

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

(111)

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}$$

(112)

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}$$

(113)

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 (60) 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 (106).

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}$$

(114)

The terms $S_{m,n},S_{n,m}$ in (106) 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}$$

(115)

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

(116)

$$\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}$$

(117)

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

(118)

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

(119)

4.2. Critical phenomena in statistical limit

Now we are prepared to average over spin variables ${\sigma }_l,l=0\dots 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}$$

(120)

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}$$

(121)

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{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}$$

(122)

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}{c}{〈{\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 \\ \mathsf{\Phi }\left(A_{n,m},\omega ,\beta \right)\Phi \left(-A_{m,n+N},\omega ,\beta \right);\hfill \end{array}$$

(123)

$$\begin{array}{ccc}& & \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}$$

(124)

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}$$

(125)

$$\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}$$

(126)

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

(127)

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}$$

(128)

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 computation represents another challenge to the number theorists.

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.

4.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 (106) (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}$$

(129)

where the big Euler ensemble average ${〈\dots 〉}_{\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 3.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 (114) 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}$$

(130)

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}$$

(131)

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}$$

(132)

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}$$

(133)

This integral is calculable (see [26]):

$$\begin{array}{ccc}& & {〈{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n〉}_{\sigma }=\frac{2(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}$$

(134)

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-4}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}$$

(135)

$$\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}$$

(136)

We investigated this new function $S\left(q\right)$ in [19], and represented it in terms of so-called multitotients [27]. For the reader’s convenience, we present the computations leading to this representation in Appendix A.

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

(137)

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

(138)

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

(139)

Here $p|q$ are prime factors of q. This remarkable identity can be directly verified using Mathematica^®[28]. 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}$$

(140)

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}$$

(141)

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 a 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-1)}$ 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$.

4.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}$$

(142)

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}$$

(143)

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}$$

(144)

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 3.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)}}{16\sqrt{2\pi }Z\left(\mu \right)}\hfill \end{array}$$

(145)

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

$$\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}$$

(146)

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}\nu }{{\mu }^2{(t+t_0)}^2};\hfill \end{array}$$

(147)

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

(148)

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

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

(149)

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 contribution from odd N to the enstrophy

$$\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}$$

(150)

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

(151)

We have the following estimate for the odd N contribution to the enstrophy in 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}$$

(152)

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,20](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.

4.5. The higher moments of the enstrophy

The higher moments of the distribution of enstrophy are also calculable [29].

We take the $2n$ point correlation function in the limit of the vanishing loop

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

(153)

$$\begin{array}{ccc}& & {\omega }_k=\frac{ı{\sigma }_k}{2}cot\left(\frac{\beta }{2}\right)\hfill \end{array}$$

(154)

With all different $m_1<\dots <m_{2n}$, the averaging over the Ising variables ${\sigma }_l$ is straightforward. It leads to the integral over $\omega$

$$\begin{array}{c}\hfill J(n,qr,N)=\oint \frac{d\omega }{2\pi }exp\left(-ıqr\right){cos}^{N-2n}\left(\omega \right){sin}^{2n}\left(\omega \right)\end{array}$$

(155)

This integral is reduced to the hypergeometric function in [26]:

$$\begin{array}{ccc}& & J(n,qr,N)=2{(-1)}^n\left({\displaystyle \genfrac{}{}{0pt}{}{N-2n}{\frac{1}{2}(N-qr)}}\right)\hfill \\ & & {}_2{F}_1\left(-2n+N+1,\frac{1}{2}(N+qr+2);\frac{1}{2}(-4n+N+qr+2);-1\right)\hfill \end{array}$$

(156)

In particular,

$$\begin{array}{ccc}& & J(n,0,N)=\frac{{(-1)}^n\Gamma \left(\frac{1}{2}(-2n+N+1)\right)}{\Gamma \left(\frac{1}{2}-n\right)\Gamma \left(\frac{N}{2}+1\right)}\hfill \end{array}$$

(157)

This term does not depend on q; it dominates in the local limit $N\to \infty$ at fixed n.

The general formula for the numerator of the moments at fixed N reads

$$\begin{array}{ccc}& & {〈\omega {\left(\overrightarrow{0}\right)}^{2n}〉}_{\mathcal{E}\left(N\right)}\propto \left({\displaystyle \genfrac{}{}{0pt}{}{N}{2n}}\right){\left(4t\right)}^{-2n}\sum _{2<q<N}S(n,q)\sum _{\begin{array}{c}r=-N/q\\ 2\left|\right(N-qr)\end{array}}^{N/q}J(n,qr,N);\hfill \end{array}$$

(158)

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

(159)

This cotangent sum is reduced to the superposition of multi-totients in Appendix.A

$$\begin{array}{ccc}& & S(n,q)={(-1)}^n\phi \left(q\right)-{(-4)}^n\sum _{j=1}^n\frac{B_{2j}{\phi }_{2j}\left(q\right)\mathrm{BernSum}(n,n-j)}{\left(2j\right)!};\hfill \end{array}$$

(160)

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

(161)

The rational coefficients $\mathrm{BernSum}(n,n-j)$ are given by recurrent relations

$$\begin{array}{ccc}& & \mathrm{BernSum}(0,0)=1;\hfill \end{array}$$

(162)

$$\begin{array}{ccc}& & \mathrm{BernSum}(n,m)=0\mathrm{if}n<m;\hfill \end{array}$$

(163)

$$\begin{array}{ccc}& & \mathrm{BernSum}(n,m)=\sum _{j_1=0}^m\sum _{j_2=0}^{m-j_1}\frac{B_{2j_1}{B}_{2j_2}\mathrm{BernSum}(n-1,m-j_1-j_2)}{\left(2j_1\right)!\left(2j_2\right)!};\hfill \end{array}$$

(164)

The specific cases are considered in Appendix A. We are interested in the local limit, even $N\to \infty$.

In this limit, the highest totient ${\phi }_{2n}\left(q\right)$ dominates the sum, and we find for the sum

$$\begin{array}{ccc}& & {〈\omega {\left(\overrightarrow{0}\right)}^{2n}〉}_{\mathcal{E}\left(\mu \right)}=\frac{{\displaystyle \sum _{evenN}exp\left(-\mu N\right)\mathsf{\Omega }(N,n)}}{{(t+t_0)}^{2n}Z\left(\mu \right)};\hfill \end{array}$$

(165)

$$\begin{array}{ccc}& & \mathsf{\Omega }(N,n)=-\frac{{\left(-\frac{1}{4}\right)}^n{B}_{2n}{N}^{2n+1}\Gamma \left(n+\frac{1}{2}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{N}{2n}}\right)\Gamma \left(-n+\frac{N}{2}+\frac{1}{2}\right)}{\pi \zeta (2n+1)\Gamma (2n+2)\Gamma \left(\frac{N}{2}+1\right)}\hfill \end{array}$$

(166)

Finally, in the limit $\mu \to +0$, replacing the sum over even N by $1/2$ of the integral of the asymptotic of this product of gamma functions, we find

$$\begin{array}{ccc}& & 〈\omega {\left(\overrightarrow{0}\right)}^{2n}〉\to \frac{{\Xi }_n}{{(t+t_0)}^{2n}{\mu }^{3n-1}}\mathrm{if}n>0;\hfill \end{array}$$

(167)

$$\begin{array}{ccc}& & {\Xi }_n=\frac{\pi \sqrt{\pi }{(-1)}^{n-1}{2}^{2-3n}{B}_{2n}\Gamma \left(3n+\frac{3}{2}\right)}{9\zeta (2n+1)\Gamma (n+1)\Gamma (2n+2)}\hfill \end{array}$$

(168)

These coefficients ${\Xi }_n$ are all positive. Here are the first five values

(169)

4.6. 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[29].

Let us consider the Fourier transform of the PDF corresponding to the vorticity moments,

$$\begin{array}{c}\hfill H(\overrightarrow{J},\mu ,t)={〈exp\left(ı\overrightarrow{J}·\widehat{\mathsf{\Omega }}·\overrightarrow{\omega }\left(0\right)\right)〉}_{\widehat{\mathsf{\Omega }}\in O\left(3\right),\mathbb{E}\left(\mu \right)}\end{array}$$

(170)

The averaging over the global rotation matrix $\widehat{\mathsf{\Omega }}\in O\left(3\right)$ yields

$$\begin{array}{c}\hfill H(\overrightarrow{J},\mu ,t)={〈\frac{sin\left(\sqrt{{\overrightarrow{J}}^2\overrightarrow{\omega }{\left(0\right)}^2}\right)}{\sqrt{{\overrightarrow{J}}^2\overrightarrow{\omega }{\left(0\right)}^2}}〉}_{\mathbb{E}\left(\mu \right)}=\sum _{n=0}^{\infty }\frac{{(-{\overrightarrow{J}}^2)}^n}{(2n+1)!}{〈\overrightarrow{\omega }{\left(0\right)}^{2n}〉}_{\mathbb{E}\left(\mu \right)}\end{array}$$

(171)

Using our solution for the moments, we get a universal scaling function

$$\begin{array}{ccc}& & H(\overrightarrow{J},\mu ,t)=1+\mu \mathcal{F}\left(\frac{{\overrightarrow{J}}^2}{{(t+t_0)}^2{\mu }^3}\right);\hfill \end{array}$$

(172)

$$\begin{array}{ccc}& & \mathcal{F}\left(z\right)=\sum _{n=1}^{\infty }\frac{z^n{(-1)}^n{\Xi }_n}{(2n+1)!}\hfill \end{array}$$

(173)

This expansion for the scaling function has a finite radius of convergence, as it follows from the asymptotic

$$\begin{array}{ccc}& & \frac{{\Xi }_n}{(2n+1)!}\to \frac{\pi }{3n^2{R}^n};\hfill \end{array}$$

(174)

$$\begin{array}{ccc}& & R=\frac{128{\pi }^2}{27};\hfill \end{array}$$

(175)

The singularity is located at negative $z=-R$. The expansion can be truncated inside the convergence region $\left|z\right|<R$. Here is the corresponding plot for expansion truncated at $n=200$ (see Figure 7).

At large positive z, it behaves as

$$\begin{array}{ccc}& & \mathcal{F}\left(z\right)\sim \pi /3\sum _{n=1}^{\infty }\frac{{(-z/R)}^n}{n^2}=\pi /3{\mathrm{Li}}_2\left(-\frac{1}{27}\left(128{\pi }^{2}z\right)\right)\sim -\frac{\pi }{6}{log}^2\left(z\right)\hfill \end{array}$$

(176)

This singularity of the Fourier transform at the imaginary axis corresponds to the exponential decrease of the PDF

$$\begin{array}{ccc}& & W\left(\right|\overrightarrow{\omega }\left(0\right)|>x)\propto exp\left(-ax(t+t_0){\mu }^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right);\hfill \end{array}$$

(177)

$$\begin{array}{ccc}& & a=\sqrt{R}=\frac{8}{3}\sqrt{\frac{2}{3}}\pi ;\hfill \end{array}$$

(178)

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

$$\begin{array}{c}\hfill {\overrightarrow{j}}_R={\mu }^{-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\overrightarrow{j}\end{array}$$

(179)

eliminates the singularity in the same way as the redefinition of viscosity

$$\begin{array}{c}\hfill {\nu }_R={\mu }^{-2}\nu \end{array}$$

(180)

eliminated the singularity in the energy dissipation.

5. The decay index spectrum

The vicinity of the fixed point in nonlinear dynamic systems provides the most interesting physical phenomena, such as anomalous dimensions of various local operators in the theory of renormalization group.

Our system is no exception. Let us perturb the momentum loop equation(21) and study the general properties of this perturbation $\delta \overrightarrow{P}(t,\theta )=\overrightarrow{Q}(t,\theta )$.

5.1. Linearized loop equation

We get a linearized equation for $\overrightarrow{Q}(t,\theta )$ of the form

$$\begin{array}{c}\hfill \nu {\partial }_t\overrightarrow{Q}={\widehat{H}}_1\left[\overrightarrow{P}\right]·\overrightarrow{Q}+{\widehat{H}}_2\left[\overrightarrow{P}\right]·\Delta \overrightarrow{Q}\end{array}$$

(181)

These matrices ${\widehat{H}}_{1,2}\left[\overrightarrow{P}\right]$ for the decaying solution (41) have an explicit time dependence

$$\begin{array}{c}\hfill {\widehat{H}}_{1,2}\left[\overrightarrow{P}\right]=\frac{\nu }{2(t+t_0){\gamma }^2}{\widehat{H}}_{1,2}\left[\overrightarrow{F}\right]\end{array}$$

(182)

which follows from the fact that the RHS of the original loop equation (21) represents the third-degree homogeneous functional of $\overrightarrow{P}$.

We find the linear equation of the form

$$\begin{array}{c}\hfill (t+t_0){\partial }_t\overrightarrow{Q}={\widehat{L}}_1·\overrightarrow{Q}+{\widehat{L}}_2·\Delta \overrightarrow{Q}\end{array}$$

(183)

with $L_{1,2}=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\widehat{H}}_{1,2}\left[\overrightarrow{F}\right]$ being some tensors in ${\mathbb{R}}_d$. There is no explicit $\theta$ dependence other than through the fixed point solution $\overrightarrow{F}\left(\theta \right)$. We shall write explicit equations in a minute.

This equation has power-like solutions

$$\begin{array}{c}\hfill \overrightarrow{Q}(t,\theta )={(t+t_0)}^{-\lambda }\overrightarrow{G}\left(\theta \right)\end{array}$$

(184)

with the index $\lambda$ determined from the eigenvalue problem

$$\begin{array}{c}\hfill -\lambda \overrightarrow{G}={\widehat{L}}_1·\overrightarrow{G}+{\widehat{L}}_2·\Delta \overrightarrow{G}\end{array}$$

(185)

We can go one step forward before specifying the parameters of this equation. Let us use our polygonal approximation for $\overrightarrow{F},\overrightarrow{G}$. This equation becomes a recurrent equation

$$\begin{array}{ccc}& & (2\lambda \widehat{I}-{\widehat{A}}_k)·{\overrightarrow{G}}_{k+1}=(2\lambda \widehat{I}-{\widehat{B}}_k)·{\overrightarrow{G}}_k;\hfill \end{array}$$

(186)

$$\begin{array}{ccc}& & {\widehat{A}}_k={\widehat{L}}_1+2{\widehat{L}}_2;\hfill \end{array}$$

(187)

$$\begin{array}{ccc}& & {\widehat{B}}_k=-{\widehat{L}}_1+2{\widehat{L}}_2;\hfill \end{array}$$

(188)

We can solve it as a matrix product (in reverse order)

$$\begin{array}{ccc}\hfill & & {\overrightarrow{G}}_{k+1}={\widehat{M}}_k\left(\lambda \right)·{\overrightarrow{G}}_0;\hfill \end{array}$$

(189)

$$\begin{array}{ccc}& & {\widehat{M}}_k\left(\lambda \right)=\prod _{i=k}^{i=0}{(\widehat{I}-{\widehat{A}}_k/\left(2\lambda \right))}^{-1}(\widehat{I}-{\widehat{B}}_k/\left(2\lambda \right))\hfill \end{array}$$

(190)

The periodicity requires $G_N=G_0$ which leads to the eigenvalue equation (this is already $d\otimes d$ matrix equation)

$$\begin{array}{ccc}& & {\widehat{M}}_N\left(\lambda \right)·{\overrightarrow{G}}_0={\overrightarrow{G}}_0\hfill \end{array}$$

(191)

As a result, we arrive at the spectral equation for the fractal dimensions of decaying turbulence

$$\begin{array}{c}\hfill det\left[{\widehat{M}}_N\left(\lambda \right)-\widehat{I}\right]=0\end{array}$$

(192)

Here is the explicit form of these two matrices ( with $\Delta {\overrightarrow{F}}_k={\overrightarrow{F}}_{k+1}-{\overrightarrow{F}}_k$, see [30])

$$\begin{array}{ccc}& & {\widehat{A}}_k=\gamma \widehat{I}+(2\gamma -4i{\overrightarrow{F}}_k·\Delta {\overrightarrow{F}}_k+i)\Delta {\overrightarrow{F}}_k\otimes {\overrightarrow{F}}_k\hfill \end{array}$$

$$\begin{array}{c}\left(-\gamma +2i({\overrightarrow{F}}_k·\Delta {\overrightarrow{F}}_k)(2{\overrightarrow{F}}_k·\Delta {\overrightarrow{F}}_k-1)\right)\Delta {\overrightarrow{F}}_k\otimes \Delta {\overrightarrow{F}}_k+i{\overrightarrow{F}}_k\otimes \Delta {\overrightarrow{F}}_k;\hfill \end{array}$$

(193)

$$\begin{array}{c}{\widehat{B}}_k=-\gamma \widehat{I}+\left(2\gamma -4i{\overrightarrow{F}}_k·\Delta {\overrightarrow{F}}_k-i\right)\Delta {\overrightarrow{F}}_k\otimes {\overrightarrow{F}}_k+\hfill \\ \left(\gamma +2i({\overrightarrow{F}}_k·\Delta {\overrightarrow{F}}_k)(2{\overrightarrow{F}}_k·\Delta {\overrightarrow{F}}_k+1)\right)\Delta {\overrightarrow{F}}_k\otimes \Delta {\overrightarrow{F}}_k-i{\overrightarrow{F}}_k\otimes \Delta {\overrightarrow{F}}_k\hfill \end{array}$$

(194)

Note that these two matrices are functions of $\gamma$, unlike the fixed point solution (60),(61), where the $\gamma$ dependence dropped.

This fact will be important for the distribution of the velocity circulation.

Another important fact. The vector ${\overrightarrow{G}}_0$ is determined up to total normalization

$$\begin{array}{c}\hfill {\overrightarrow{G}}_0=g{\overrightarrow{\psi }}_1\left(M\right)\end{array}$$

(195)

where the eigenvector ${\overrightarrow{\psi }}_1\left(M\right)$ is a unit $3D$ vector, and g is an arbitrary parameter, having dimension

$$\begin{array}{c}\hfill \left[g\right]=\left[L\right]{\left[T\right]}^{\lambda -1}\end{array}$$

(196)

This parameter g is not universal and depends on the initial data specific to the flow, like the energy dissipation in the K41 model.

This issue deserves further numerical and analytical study, which is outside our scope.

5.2. The circulation distribution

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

$$\begin{array}{c}\hfill W[C,\Gamma ]=〈\delta \left(\Gamma -{\oint }_{C}d\overrightarrow{r}·\overrightarrow{v}\left(\overrightarrow{r}\right)\right)〉={\int }_{-\infty }^{\infty }\frac{d\gamma }{2\pi \nu }exp\left(-ı\frac{\gamma }{\nu }\Gamma \right)\mathsf{\Psi }[\gamma ,C];\end{array}$$

(197)

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 $\gamma$. Formally, we get the delta function $\delta (\Gamma )$, 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:

$$\begin{array}{ccc}\hfill & & W[C,\Gamma ]=\frac{ı}{\nu }{\int }_{-\infty }^{\infty }\frac{d\gamma }{2\pi }exp\left(-ı\frac{\gamma }{\nu }\Gamma \right)\hfill \\ & & {〈exp\left(+\frac{ı}{\sqrt{2\nu (t+t_0)}}\sum _k\Delta {\overrightarrow{C}}_k·{\overrightarrow{F}}_k\right){(t+t_0)}^{-\lambda }\sum _k\Delta {\overrightarrow{C}}_k·{\overrightarrow{G}}_k〉}_{\mathbb{E}\left(N\right)};\hfill \end{array}$$

(198)

The solution (189) for $G_k$ depends upon $\gamma$, and so does the eigenvalue $\lambda$, which satisfies the spectral equation (192). This distribution becomes highly nontrivial for any large finite N. The hope is that there is a limit $N\to \infty$, where some simple fractal (or multifractal) laws will emerge.

5.3. The spectral identity and Wilson loop asymptotics

The distribution of the eigenvalues $\lambda$ can be studied using the resolvent

$$\begin{array}{c}\hfill R_N(\gamma ,\lambda )=\frac{1}{d}\mathrm{tr}{\widehat{M}}_N^{\prime }\left(\lambda \right){\left({\widehat{M}}_N\left(\lambda \right)-\widehat{I}\right)}^{-1}\end{array}$$

(199)

By construction, at finite N this is a rational function of $\lambda$ with simple poles at the spectrum.

In the continuum limit, it becomes a holomorphic function in a complex plane cut along the spectrum.

Using exact representation (189) we can prove that

$$\begin{array}{ccc}& & {\widehat{M}}_N(\lambda \to \infty )\to \widehat{I}+\frac{\widehat{\mathsf{\Sigma }}}{2\lambda };\hfill \end{array}$$

(200)

$$\begin{array}{ccc}& & \widehat{\mathsf{\Sigma }}=\sum _k({\widehat{A}}_k-{\widehat{B}}_k);\hfill \end{array}$$

(201)

$$\begin{array}{ccc}& & R_N(\gamma ,\lambda \to \infty )\to \frac{1}{d}\mathrm{tr}\frac{-\widehat{\mathsf{\Sigma }}}{2{\lambda }^2}\frac{2\lambda }{\widehat{\mathsf{\Sigma }}}=\frac{-1}{\lambda }\hfill \end{array}$$

(202)

Consider the clockwise contour $\omega$ surrounding the spectrum of $\lambda$ in a complex plane (not necessarily a real axis).

Then, we can use the residue at infinity to compute the contour integral

$$\begin{array}{c}\hfill {\int }_{\omega }\frac{d\lambda }{2\pi ı}{R}_N\left(\lambda \right)=1\end{array}$$

(203)

The left side of this identity can be inserted inside the statistical average, such as the circulation PDF integral (198):

$$\begin{array}{ccc}\hfill & & W[C,\Gamma ]={〈exp\left(\frac{ı}{\sqrt{2\nu (t+t_0)}}\sum _k\Delta {\overrightarrow{C}}_k·{\overrightarrow{F}}_k\right)J\left(t\right)〉}_{\mathbb{E}\left(N\right)};\hfill \end{array}$$

(204)

$$\begin{array}{ccc}& & J\left(t\right)={\int }_{-\infty }^{\infty }\frac{d\gamma }{2\pi }exp\left(-ı\frac{\gamma }{\nu }\Gamma \right){\int }_{\omega }\frac{d\lambda }{2\pi }{R}_N(\gamma ,\lambda )\frac{{\sum }_k\Delta {\overrightarrow{C}}_k·{\overrightarrow{G}}_k}{\nu {(t+t_0)}^{\lambda }}\hfill \end{array}$$

(205)

By taking residues in the poles of the resolvent, we can relate this integral to the sum over the spectrum of poles. However, it is better 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 asymptotic behavior of $W[C,\Gamma ]$ at a large time t is determined by the saddle point equations for two variables $\gamma ,\lambda$

$$\begin{array}{ccc}& & logt={〈\frac{\partial log{R}_N(\gamma ,\lambda )}{\partial \lambda }〉}_{\mathbb{E}\left(N\right)}+{〈\frac{\partial log{\sum }_k\Delta {\overrightarrow{C}}_k·{\overrightarrow{G}}_k}{\partial \lambda }〉}_{\mathbb{E}\left(N\right)};\hfill \end{array}$$

(206)

$$\begin{array}{ccc}& & \frac{ı\Gamma }{\nu }={〈\frac{\partial log{R}_N(\gamma ,\lambda )}{\partial \gamma }〉}_{\mathbb{E}\left(N\right)}+{〈\frac{\partial log{\sum }_k\Delta {\overrightarrow{C}}_k·{\overrightarrow{G}}_k}{\partial \gamma }〉}_{\mathbb{E}\left(N\right)};\hfill \end{array}$$

(207)

The stability of our fixed-point solution requires the condition of the saddle point ${\gamma }_c,{\lambda }_c$

$$\begin{array}{c}\hfill \mathbf{Re}{\lambda }_c\ge 0\end{array}$$

(208)

We plan to study the spectrum of fractal dimensions and asymptotic distribution of the circulation in the forthcoming paper [19], using the NYU AD supercomputer cluster.

6. Discussion

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

In our view, the most advanced research in this field so far has been performed by Tao (see his recent review [31]). Based on a simplified model of the averaged Navier-Stokes, Tao conjectured [18] 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 [17,20], 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 $\nu \to 0$ limit of the Navier-Stokes equation.

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: the singularity we have found exists at finite viscosity, in a spirit of Tao’s conjecture, though we are not studying a particular solution with a finite time blow-up.

Again, let us stress it: we are not claiming any theorems about finite time singularities of the Navier-Stokes equation. We are studying a different phenomenon: turbulent flow, a stochastic solution of the Navier-Stokes equation, covering a certain manifold.

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

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 the averaging over 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 $\mu \to 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 enter the game, leading 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.

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

The long-term evolution of Newton’s dynamical system with many particles eventually covers the energy surface (microcanonical ensemble). The ergodic hypothesis (accepted in Physics but still not proven mathematically) states that this energy surface is covered uniformly.

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\left(ı\overrightarrow{p}·\overrightarrow{x}\right)$.

This is the solution we have found, except the dot product $\overrightarrow{p}·\overrightarrow{x}$ becomes a symplectic form $\oint \overrightarrow{P}\left(\theta \right)·d\overrightarrow{C}\left(\theta \right)$ in the loop space.

Our momentum $\overrightarrow{P}(\theta ,t)$ is not an integral of motion, but simple scaling properties of the pure Navier-Stokes equation lead to the solution with $\overrightarrow{P}(\theta ,t)\propto \overrightarrow{F}\left(\theta \right)/\sqrt{t}$, with $\overrightarrow{F}\left(\theta \right)$ being the integral of motion.

The rest is a purely technical task: substituting this scaling solution into the Navier-Stokes equation and solving the resulting universal equation for a fixed point $\overrightarrow{F}\left(\theta \right)$.

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\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 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\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 (104) and the expectation value of the energy dissipation (147). As discussed in the previous section, the singularities at $\mu \to 0$ indicate inconsistencies in the Navier-Stokes equation as an idealization of molecular dynamics.

6.4. Classical flow and quantum geometry

Our computations heavily rely on the number theory, particularly Lehmer’s multitotients ${\phi }_l\left(q\right)$, (137), generalizing [27] 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[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.

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\ll 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 [18]; however, it has a different meaning. Our singularity is a divergence of 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 divergencies in these correlation functions, like QED, though these divergencies emerge in the exact solution beyond perturbation theory.

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

There is no gauge invariance regarding the velocity field in fluid dynamics (though there is such invariance in the Clebsch variables [17]). 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 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.

6.6. 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,17], 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,17] (and also here, in Section 1.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 (35). The corresponding Wilson loop was equal to the Stokes-type functional (27).

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

6.7. Continuum spectrum of anomalous dimensions and multifractality

We studied the linear perturbations of our fixed manifold, which decay with time as $t^{-\lambda }$, with a calculable spectrum of decay indexes $\lambda$. At fixed parameters $N,p,q,r,{\sigma }_i$ of the Euler ensemble, this spectrum is discrete, but with random values of these parameters, in the local limit $\mu \to 0,〈N〉\to \infty$, the spectrum becomes continuous.

The explicit equation (192) for the spectrum of anomalous dimensions in decaying turbulence is a nice surprise. Future investigation will likely relate it to some multifractal properties.

We derived the representation for the averages over this continuum spectrum for the Wilson loop PDF (204).

This kind of saddle point integral is used to justify the multifractal models. In our case, the weight in that integral is calculable in terms of the Euler ensemble. The analytic result for the large-time asymptotic would require more work, however.

The most urgent task is to study the stability of our spectrum of decay indexes $\mathbf{Re}\lambda >0$. We do not see any alternative to numerical study of this spectrum [19].

6.8. 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].

The most striking nonperturbative effect 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 discrete statistical system similar 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 (148) is related to the prime factorization of large integers used in so-called multi-totients [27].

Whether this quantum solution to the classical problem is realized in Nature remains to be seen.