Decaying Turbulence as a Fractal Curve

1. Loop Equation

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

Here is a summary for the new generation.

We write the Navier-Stokes equation as follows

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

(4)

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

(5)

The Wilson loop average for the Turbulence

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

(6)

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

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

(7a)

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

(7b)

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

(7c)

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

(7d)

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

(7e)

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

(7f)

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

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

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

(8)

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

(9)

In the review, [1,2], we present the explicit limiting procedure needed to define these functional derivatives in terms of finite variations of the loop while keeping it closed.

All the operators ${\partial }_{\mu },{\widehat{\omega }}_{\alpha \beta },{\widehat{v}}_{\alpha }$ are expressed in terms of the spike operator

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

(10)

The area derivative operator can be regularized as

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

(11)

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

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

(12)

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

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

(13)

In three dimensions, it follows from identity $\overrightarrow{\nabla }·\overrightarrow{\omega }=0$; in general dimension $d>3$, the dual vorticity $\tilde{\omega }$ is an antisymmetric tensor with $d-2$ components. The divergence of this tensor equals zero identically.

However, for the loop functional, this restriction is not an identity; it reflects that this functional is a function of a circulation of some vector field, averaged by some set of parameters.

This constraint was analyzed in [2] 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 customary Gaussian random forces, choosing instead to randomize the initial data for the velocity field.

These random forces would lead to the potential term [2] in the loop Hamiltonian ${\mathcal{H}}_C$, breaking certain symmetries needed for the dimensional reduction we study below.

With random initial data instead of time-dependent delta-correlated random forcing, we no longer describe the steady state (i.e., statistical equilibrium) but decaying Turbulence, which is also an interesting process, manifesting the same critical phenomena.

The energy is pumped in at the initial moment $t=0$ and slowly dissipates over time, provided the viscosity is small enough, corresponding to the large Reynolds number we are studying.

2. Dimensional Reduction

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

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

(14)

The difference with the original definition of $\Psi \left[C\right]$ is that our new function $P_{\alpha }\left(\theta \right)$ depends directly on $\theta$ rather then through the function $v_{\alpha }\left(r\right)$ taken at $r_{\alpha }=C_{\alpha }\left(\theta \right)$.

This transformation is the dimensional reduction $d\Rightarrow 1$ we mentioned above. From the point of view of the loop functional, there is no need to deal with field $v\left(r\right)$; one could take a shortcut.

The reduced dynamics must be fitted to the Navier-Stokes dynamics of the original field. With the loop calculus developed above, we have all the necessary tools to build these reduced dynamics.

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 in the Introduction to the Kelvinon theory and any other Navier-Stokes stationary solutions with a singularity at fixed loop C in space.

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

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

(15)

The 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){\Omega }_{\beta \alpha }$$

(16)

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

The spike derivative D in the above equation

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

(17)

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

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

As we shall see, the viscosity enters initial data so that at any finite time t, the solution for P still depends on viscosity.

Another observation is that the spike derivative $D(\theta ,ϵ)$ turns to the discontinuity $\Delta P\left(\theta \right)=P\left({\theta }^{+}\right)-P\left({\theta }^{-}\right)$ in the limit $ϵ\to 0^{+}$

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

(18)

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

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}{c}{\Omega }_{\alpha \beta }(\theta ,0^{+})=\frac{-ı}{\nu }{P}_{\alpha \beta }\left(\theta \right);\hfill \end{array}$$

(19)

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

(20)

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

(21)

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

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

(22)

The Bianchi constraint is identically satisfied as it should

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

(23)

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

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

(24)

This equation is complex due to the irreversible dissipation effects in the Navier-Stokes equation.

The viscosity dropped from the right side of this equation; it can be absorbed in units of time. Viscosity also enters the initial data, as we shall see in the next Section on the example of the random rotation.

However, the large-time asymptotic behavior of the solution would be universal, as it should be in the Turbulent flow.

We are looking for a degenerate fixed point [2], 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.

3. Random Global Rotation

Possible initial data for the reduced dynamics were suggested in the original papers [1,2]. The initial velocity field’s simplest meaningful distribution is the Gaussian one, with energy concentrated in the macroscopic motions. The corresponding loop field reads

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

(25)

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

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

(26)

The potential part drops out in the closed loop integral.

The correlation function varies at the macroscopic scale, which means that we could expand it in the Taylor series

$$f(r-r^{\prime })\to f_0-f_1{(r-r^{\prime })}^2+\cdots$$

(27)

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

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

(28)

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

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

(29)

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

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

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

(30)

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

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

(31)

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

The physical meaning of this $\varphi$ is the random uniform vorticity 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 [2]. 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 time derivative at this special initial data vanishes so that the exact solution of the loop equation with this initial data equals its initial value (30).

Naturally, the time derivative of the momentum loop with the corresponding initial data will vanish as well.

It is instructive to look at the momentum trajectory $P_{\alpha }\left(\theta \right)$ for this fixed point.

The functional Fourier transform [1,2] 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}{c}{P}_{\alpha }\left(\theta \right)=\sum _{n=1}^{\infty }{P}_{\alpha ,n}exp\left(ın\theta \right)+{\overline{P}}_{\alpha ,n}exp\left(-ın\theta \right);\hfill \end{array}$$

(32)

$$\begin{array}{c}{P}_{\alpha ,n}=\mathcal{N}(0,1)\forall \alpha ,n>0;\hfill \end{array}$$

(33)

$$\begin{array}{c}{\overline{P}}_{\alpha ,n}=\frac{4\sqrt{f_1}}{n}{\varphi }_{\alpha \beta }{P}_{\beta ,n};\forall \beta ,n>0;\hfill \end{array}$$

(34)

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

(35)

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

(36)

As for the constant part $P_{\alpha ,0}$ of $P_{\alpha }\left(\theta \right)$ , it 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}{c}{〈P_{\alpha ,n}{P}_{\beta ,m}〉}_{t=0}=\frac{4\sqrt{f_1}}{m}{\delta }_{-nm}{\varphi }_{\alpha \beta };\hfill \end{array}$$

(37)

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

(38)

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

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

(39)

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

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

Here are the steps involved

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

(40)

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

(41)

Here the tensor area $\Sigma$ was defined in (9). It is an antisymmetric tensor; therefore its trace with a symmetric tensor ${\widehat{\Omega }}_{\alpha \beta }{\widehat{\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 some distribution for every point on the loop.

4. Decay or Fixed Point

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

$$\begin{array}{cc}\hfill \mathbf{Decay}:& \overrightarrow{P}\to 0;\Psi \left[C\right]\to 1;\hfill \end{array}$$

(42)

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

(43)

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

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

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

(44)

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

(45)

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

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

(46)

The fixed point could mean self-sustained Turbulence, which would be too good to be true, violating the second law of Thermodynamics. Indeed, it is easy to see that this fixed point cannot exist.

The fixed point equation (46) is a linear relation between two vectors $\overrightarrow{P},\Delta \overrightarrow{P}$ with coefficients depending on various scalar products. The generic solution is simply

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

(47)

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

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

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

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

5. Markov Process in Complex Space

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}{c}\Delta \overrightarrow{F}\stackrel{?}{=}\lambda \overrightarrow{F};\hfill \end{array}$$

(48)

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

(49)

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

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

(50a)

$$\begin{array}{c}{(\overrightarrow{F}·\Delta \overrightarrow{F})}^2={\overrightarrow{F}}^2+ı\overrightarrow{F}·\Delta \overrightarrow{F}\hfill \end{array}$$

(50b)

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

Also, note that this complex vector $\overrightarrow{F}\left(\theta \right)$ is dimensionless, and the fixed point equation (50) is completely universal!

One can build this solution as a 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 following discrete version of the discontinuity equations (50).

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

(51a)

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

(51b)

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

$$\begin{array}{c}\hfill \left\{{\overrightarrow{q}}_k^2=1,{\overrightarrow{F}}_k·{\overrightarrow{q}}_k=\frac{1}{2}\left(ı\pm \sqrt{4{\overrightarrow{F}}_k^2-1+2ı}\right)\right\}\end{array}$$

(52)

after which we can find the next value

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

(53)

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$. On top of this Markov process, there is an extra periodicity requirement ${\overrightarrow{F}}_N={\overrightarrow{F}}_0$.

This requirement represents a nonlinear restriction on initial position ${\overrightarrow{F}}_0$.

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

(54)

Note that the complex unit vector is not defined with the Euclidean metric in six 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}$$

(55)

which leads to two conditions between real and imaginary parts

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

(56)

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

(57)

In d dimensions, there are $d-1$ complex parameters of the unit vector; with an extra linear constraint in (52), 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.

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

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

(58)

We use this complex vector constraint to fix the arbitrary initial complex vector ${\overrightarrow{F}}_0$ as a function of all remaining parameters.

Looking ahead into the rest of our investigation, it turns out that the closure conditions fix only half of the $2d$ real parameters in the initial point ${\overrightarrow{F}}_0$. The remaining parameters are free zero modes of our fixed manifold.

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; therefore, the translational zero modes of the momentum loop do not lead to ambiguities.

However, the missing d out of $2d$ parameters in ${\overrightarrow{F}}_0$ mean that some other d parameters should be adjusted to provide the momentum loop closure.

We discuss this issue in the next Section, where we derive the SDE for the closed momentum loop in three dimensions. This SDE has explicit terms, which we computed in Appendix A and coded in Mathematica^® in [32].

The adjustment of parameters we mentioned earlier yields three constraints on the Wiener process we derived.

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

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

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

(59)

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

(60)

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

(61)

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.

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.

In the best-known example, the microcanonical Gibbs distribution covers the energy surface with a uniform measure (ergodic hypothesis, widely accepted in Physics).

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’s $N\to \infty$ limit, is parametrized by N arbitrary local rotations, 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.

6. The Probability Measure

The simplest case where these equations have nontrivial solutions is the three-dimensional space. For smaller dimensions of space, there is only a degenerate solution with zero vorticity (a vanishing cross product $\widehat{\Omega }\propto \overrightarrow{P}\times \Delta \overrightarrow{P}$ ). Thus, we only consider $d>2$ in the rest of the paper.

The complex unit vector in d dimensions can be parametrized by rotation matrix and a unit real vector in $d-2$ dimensions

$$\begin{array}{c}\overrightarrow{q}=\widehat{O}·\overrightarrow{u}({\alpha }_1,{\alpha }_2,\overrightarrow{w},\beta );\hfill \end{array}$$

(62a)

$$\begin{array}{c}\overrightarrow{u}({\alpha }_1,{\alpha }_2,\overrightarrow{w},\beta )=\{{\alpha }_1,{\alpha }_2\overrightarrow{w},ı\beta \};\hfill \end{array}$$

(62b)

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

(62c)

$$\begin{array}{c}{\alpha }_1^2+{\alpha }_2^2=1+{\beta }^2;\hfill \end{array}$$

(62d)

The following steps lead to this canonical form. Take a general complex d-vector $\overrightarrow{q}$ and choose the rotation $\widehat{O}\in O\left(d\right)$ to direct its imaginary part at the last axis d.

The imaginary part of the condition ${\overrightarrow{q}}^2=1$ implies the real part of this vector has zero component d. This real vector in $d-1$ dimensions can be parametrized as $\{{\alpha }_1,{\alpha }_2\overrightarrow{w}\}$ with the unit vector $\overrightarrow{w}\in {\mathbb{S}}_{d-3}$ and arbitrary real parameters ${\alpha }_1,{\alpha }_2$.

There is a multiple counting of the same unit vector with this parametrization: the rotation matrix space $O\left(d\right)$ must be factored by rotations $O(d-3)$ of the unit vector $\overrightarrow{w}$.

$$\begin{array}{c}\hfill \widehat{O}\in {(}^{O\left(d\right)}{/}_{O(d-3)})\end{array}$$

(63)

Also, the sign change of ${\alpha }_2$ is equivalent to the reflection of the vector $\overrightarrow{w}$, so we have to factor out such reflections and keep the sign of ${\alpha }_2$ arbitrary

$$\begin{array}{c}\overrightarrow{w}\in {(}^{{\mathbb{S}}^{d-3}}{/}_{{\mathbb{Z}}^2});\hfill \end{array}$$

(64)

$$\begin{array}{c}\{{\alpha }_1,{\alpha }_2\}\in {\mathbb{R}}_2\hfill \end{array}$$

(65)

The complex constraint for ${\overrightarrow{F}}_k·{\overrightarrow{q}}_k$ can be used to fix these ${\alpha }_1,{\alpha }_2$ as a linear function of $\beta$ given a complex vector

$$\begin{array}{c}{\overrightarrow{f}}_k={\widehat{O}}_k^T·{\overrightarrow{F}}_k;\hfill \end{array}$$

(66)

as follows:

$$\begin{array}{c}\{{\alpha }_1,{\alpha }_2\}={\widehat{M}}^{-1}.\{\mathbf{Re}\left(R\right)-\beta \mathbf{Im}\left(c\right),\beta \mathbf{Re}\left(c\right)+\mathbf{Im}\left(R\right)\};\hfill \end{array}$$

(67)

$$\begin{array}{c}R=\frac{1}{2}\left(ı\pm \sqrt{4{\overrightarrow{f}}_k^2-1+2ı}\right)\hfill \end{array}$$

(68)

where ${\overrightarrow{f}}_k=\{a,\overrightarrow{b},c\}$ and

$$\widehat{M}=\left(\begin{array}{cc}\mathbf{Re}\left(a\right)& \mathbf{Re}(\overrightarrow{b}·\overrightarrow{w})\\ \mathbf{Im}\left(a\right)& \mathbf{Im}(\overrightarrow{b}·\overrightarrow{w})\end{array}\right)$$

(69)

After that, ${\alpha }_1^2+{\alpha }_2^2=1+{\beta }^2$ yields a quadratic equation for $\beta$.

Note in passing that $\overrightarrow{u}$ belongs to De Sitter space $d{S}_{d-1}$. However, this is where an analogy with the ADS/CFT duality ends.

There are, in general, four solutions for $\beta$: two signs for R in (68) and two more signs in a solution of the quadratic equation for $\beta$ (62a).

We have to choose a particular real solution for $\beta$. A universal option is to choose the step with the smallest Euclidean distance ${\left(\mathbf{Re}\overrightarrow{q}\right)}^2+{\left(\mathbf{Im}\overrightarrow{q}\right)}^2$. We used this choice in our initial simulations [32], but later we switched to another method, using the SDE we describe later in this work. The SDE guarantees the closure condition, unlike the naive random walk approach.

We arrive at the invariant distribution for our fractal curve. At a fixed N

$$\begin{array}{c}d{\Omega }_d\left(N\right)=\prod _{k=0}^{N-1}\frac{\left(d{\widehat{O}}_k\right)\left(d{\overrightarrow{w}}_k\right)}{2\left|O(d-3)\right|}\int d^{2d}{\overrightarrow{F}}_0{\delta }^{2d}\left(\overrightarrow{Q}\right);\hfill \end{array}$$

(70a)

$$\begin{array}{c}\overrightarrow{Q}={\overrightarrow{F}}_N-{\overrightarrow{F}}_0=\sum _{k=0}^{N-1}{\overrightarrow{q}}_k;\hfill \end{array}$$

(70b)

$$\begin{array}{c}{\overrightarrow{F}}_{k+1}={\overrightarrow{F}}_k+{\overrightarrow{q}}_k;\forall k=0,\cdots N-1;\hfill \end{array}$$

(70c)

where ${\overrightarrow{q}}_k$ are complex vectors, parametrized by ${\widehat{O}}_0,\cdots {\widehat{O}}_{N-1},{\overrightarrow{w}}_0,\cdots {\overrightarrow{w}}_{N-1}$ via recurrent equations (62),(67).

The complex vector’s integration and delta function is understood as a product of its real and imaginary parts.

We conclude that the fixed manifold ${\mathbb{T}}_d\left(N\right)$ of the decaying Turbulence is a subset of the tensor product of rotational and spherical spaces.

$$\begin{array}{c}\hfill {\mathbb{T}}_d\left(N\right)\in {\left({(}^{O\left(d\right)}{/}_{O(d-3)})\otimes {(}^{{\mathbb{S}}^{d-3}}{/}_{{\mathbb{Z}}^2})\right)}^{\otimes N}\end{array}$$

(71)

This subset is selected by imposing the closure condition ${\sum }_k{\overrightarrow{q}}_k=0$, which in general provides $2d$ nonlinear relations between all parameters: $\overrightarrow{\lambda },{\widehat{O}}_k,{\overrightarrow{w}}_k$ for each choice of the real solutions for $\beta$ on each step.

This closure condition is sufficient by parameter count to eliminate a complex vector ${\overrightarrow{F}}_0$; however, we discovered in 3D that half of the parameters in complex vector ${\overrightarrow{F}}_0$ are left undetermined, leading instead to three real constraints on the parameters of the rotation matrices.

We cannot resolve the global closure conditions, but we found a way to achieve the same goal by an SDE describing the evolution of our curve from one of the two exact symmetric solutions we have found; this method preserves the closure condition at each infinitesimal step of the stochastic process.

7. Symmetric Fixed Point

The above formal definition of the probability measure does not offer a practical simulation method for covering this manifold.

We attempted to simulate a random walk ${\overrightarrow{F}}_k\Rightarrow {\overrightarrow{F}}_{k+1}$ step by step, taking random rotation matrices. Unfortunately, there was a rapidly diminishing probability of the return to the vicinity of the initial point ${\overrightarrow{F}}_N={\overrightarrow{F}}_0$ after N steps.

We could not numerically solve the resulting transcendental equation for the initial position ${\overrightarrow{F}}_0$ at large N, neither by analytical nor by Monte Carlo methods.

Instead, we have found an alternative algorithm for covering this manifold, preserving the closed curve.

First, we have found a symmetric solution [32] of our recurrent equation (51) for arbitrary N

$$\begin{array}{c}{\overrightarrow{\Phi }}_k=\frac{1}{2sin\left(\frac{\pi }{N}\right)}\left\{cos\left(\frac{2\pi k}{N}\right),sin\left(\frac{2\pi k}{N}\right)\overrightarrow{w},ıcos\left(\frac{\pi }{N}\right)\right\}\hfill \end{array}$$

(72)

Here $\overrightarrow{w}\in {\mathbb{S}}^{d-3}$ is a unit vector. As we pointed out in Section 5, a reflected sequence ${\overrightarrow{\Phi }}_{N-k}^{*}$ also represents a solution to the recurrent equations (51).

The random walk step ${\overrightarrow{\varphi }}_k={\overrightarrow{\Phi }}_{k+1}-{\overrightarrow{\Phi }}_k$ is a real unit vector in this case, corresponding to $\widehat{O}=1,\beta =0$ in (62) and

$$\begin{array}{c}{\alpha }_1+ı{\alpha }_2=ıexp\left(\frac{\pi ı(2k+1)}{N}\right);\hfill \end{array}$$

(73)

Our covering algorithm will use a symmetric fixed point or its reflection as a starting element.

Due to the global $O\left(d\right)$ symmetry, the rotated curve $\{\widehat{O}·{\overrightarrow{\Phi }}_0,\cdots \widehat{O}·{\overrightarrow{\Phi }}_{N-1}\}$ with arbitrary orthogonal matrix $\widehat{O}$ is also a solution. We use random global rotation in numerical simulations.

The imaginary parts of the steps ${\overrightarrow{q}}_k$ are zero vectors at the start. Still, the evolution below will involve complex infinitesimal rotations $\delta {\overrightarrow{q}}_k={\overrightarrow{\mu }}_k\times {\overrightarrow{q}}_k$ so that the imaginary parts appear later in the evolution.

Some constraints are left in the solution for the vectors ${\overrightarrow{q}}_k$ even after the complex rotations. In three dimensions, three scalar constraints remain on the imaginary parts of the complex rotation vectors ${\overrightarrow{\mu }}_k$.

These constraints are needed to provide the closure condition. There are only three scalar constraints among N real vectors, which leads to a nontrivial $3N-3$ dimensional quotient space.

8. Infinitesimal Complex Rotations in 3D

Let us assume we already know a particular solution ${\overrightarrow{F}}_0,{\overrightarrow{q}}_0,\cdots {\overrightarrow{q}}_{N-1}$ of the recurrent equations (51) in $d=3$ and perturb it by an infinitesimal transformation of the complex vectors ${\overrightarrow{q}}_k$, preserving their square.

We also shift the initial point ${\overrightarrow{F}}_0$ to keep the loop closed after infinitesimal transformations of all the steps ${\overrightarrow{q}}_k$.

$$\begin{array}{c}\delta {\overrightarrow{q}}_k={\overrightarrow{\mu }}_k\times {\overrightarrow{q}}_k;\hfill \end{array}$$

(74)

$$\begin{array}{c}\delta {\overrightarrow{F}}_0=\overrightarrow{\lambda };\hfill \end{array}$$

(75)

Here ${\overrightarrow{\mu }}_k,\overrightarrow{\lambda }$ are infinitesimal complex 3D vectors.

The real part $\mathbf{Re}{\overrightarrow{\mu }}_k\in {\mathbb{R}}^3$ comes from the infinitesimal group transformation ${\delta }_L$ of rotation matrices in our canonical form (62)

$$\begin{array}{c}{\delta }_L{\widehat{O}}_k={\widehat{\Omega }}_k·{\widehat{O}}_k;\hfill \end{array}$$

(76)

$$\begin{array}{c}{\widehat{\Omega }}_k^{\alpha \beta }=e^{\alpha \beta \gamma }{\Omega }_k^{\gamma };\hfill \end{array}$$

(77)

$$\begin{array}{c}{\overrightarrow{\Omega }}_k=-\mathbf{Re}{\overrightarrow{\mu }}_k;\hfill \end{array}$$

(78)

$$\begin{array}{c}{\delta }_L{\overrightarrow{q}}_k=\mathbf{Re}{\overrightarrow{\mu }}_k\times {\overrightarrow{q}}_k;\hfill \end{array}$$

(79)

The imaginary part $\mathbf{Im}{\overrightarrow{\mu }}_k$ leads to the infinitesimal transformation of parameters ${\alpha }_1,{\alpha }_2,\beta$ in two-dimensional de Sitter space $d{S}_2$; therefore, there are only two independent components of $\mathbf{Im}{\overrightarrow{\mu }}_k$.

We do not need an explicit split of the parameters of ${\overrightarrow{\mu }}_k$ into these two transformations; it is sufficient to know that cross product ${\overrightarrow{\mu }}_k\times {\overrightarrow{q}}_k$ with any complex vector ${\overrightarrow{\mu }}_k$ is orthogonal to ${\overrightarrow{q}}_k$, as we need it in our random walk with ${\overrightarrow{q}}_k^2=1$.

Below, we will parameterize $\mathbf{Im}{\overrightarrow{\mu }}_k$ by two scalar parameters.

There are two contributions to the variation of each position ${\overrightarrow{F}}_k$. One variation comes from the rotation of the step from the previous position, and another comes from the variation of the previous position.

$$\begin{array}{c}\delta {\overrightarrow{F}}_k=\delta {\overrightarrow{F}}_{k-1}+{\overrightarrow{\mu }}_{k-1}\times {\overrightarrow{q}}_{k-1}=\lambda +\sum _0^{k-1}{\overrightarrow{\mu }}_l\times {\overrightarrow{q}}_l;\hfill \end{array}$$

(80)

By variation of the second of the constraints in (51), which we rewrite as

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

(81)

we find the following set of relations between infinitesimal ${\overrightarrow{\mu }}_k,\overrightarrow{\lambda }$

$$\begin{array}{c}{G}_k{\overrightarrow{F}}_k·{\overrightarrow{\mu }}_k\times {\overrightarrow{q}}_k=\left(2{\overrightarrow{F}}_k-G_k{\overrightarrow{q}}_k\right)·\left(\overrightarrow{\lambda }+\sum _0^{k-1}{\overrightarrow{\mu }}_l\times {\overrightarrow{q}}_l\right);\hfill \end{array}$$

(82a)

$$\begin{array}{c}{G}_k=(2{\overrightarrow{F}}_k·{\overrightarrow{q}}_k-ı);\hfill \end{array}$$

(82b)

We treat it as a recurrent system of equations for $\mathbf{Im}{\overrightarrow{\mu }}_k$, assuming known values of $\mathbf{Re}{\mu }_l,\overrightarrow{\lambda }$.

After that, the complex vector $\overrightarrow{\lambda }$ is supposed to be determined from the closure equation

$$\begin{array}{c}\hfill \sum _{l=0}^{N-1}{\overrightarrow{\mu }}_l\times {\overrightarrow{q}}_l=0\end{array}$$

(83)

assuming all ${\overrightarrow{\mu }}_l$ expressed as linear combinations of $\mathbf{Re}\overrightarrow{\lambda },\mathbf{Im}\overrightarrow{\lambda },\mathbf{Re}{\overrightarrow{\mu }}_i$.

As we found in [32] in three dimensions, this system of equations for $\overrightarrow{\lambda }$ is degenerate: three parameters in $\overrightarrow{\lambda }$ are left undetermined.

The solution for $\overrightarrow{\lambda }$ exists only if N vectors $\{\mathbf{Re}{\overrightarrow{\mu }}_0,\cdots \mathbf{Re}{\overrightarrow{\mu }}_{N-1}\}$ obey three scalar constraints.

In other words, the complex vector equation (83) reduces to three constraints for $\overrightarrow{\lambda }$ and another three constraints for $\mathbf{Re}{\overrightarrow{\mu }}_k$. The complex vector $\overrightarrow{\lambda }$ is left with three free components, and the vectors $\{\mathbf{Re}{\overrightarrow{\mu }}_0,\cdots \mathbf{Re}{\mu }_{N-1}\}$ are left with $3N-3$ free components out of $3N$.

The solution of these equations, which we find in Appendix A, has the form

$$\begin{array}{c}\overrightarrow{\lambda }=\sum _{l=0}^{N-1}{\widehat{\Lambda }}_l·\mathbf{Re}{\overrightarrow{\mu }}_l;\hfill \end{array}$$

(84)

$$\begin{array}{c}\mathbf{Im}{\overrightarrow{\mu }}_k=\sum _{l=0}^{N-1}{\widehat{S}}_{kl}·\mathbf{Re}{\overrightarrow{\mu }}_l\hfill \end{array}$$

(85)

with real $3\times 3$ matrices ${\widehat{S}}_{kl}$, and complex $3\times 3$ matrices ${\widehat{\Lambda }}_l$; these matrices depend on the current values of all the vectors ${\overrightarrow{F}}_k$.

In addition, we have found three linear constraints on $\mathbf{Re}{\overrightarrow{\mu }}_l$, related to three complex null vectors of a block matrix $\widehat{H}$ involved in the equation for $\overrightarrow{\lambda }$.

The vector $\overrightarrow{\lambda }$ is defined modulo this null space

$$\begin{array}{c}\hfill \overrightarrow{\lambda }\Rightarrow \overrightarrow{\lambda }+\sum _{i=1}^3{c}_i{\overrightarrow{\psi }}_i\end{array}$$

(86)

Due to the closure of the original loop $C\in b{R}^d$, the translation of $\overrightarrow{\lambda }$ by arbitrary complex vector does not change the circulation in (54). This translation of $\overrightarrow{\lambda }$ leads to the global translation of our momentum curve $\overrightarrow{P}\left(\theta \right)$, preserving the circulation over the closed loop in space.

We resolved this ambiguity of $\overrightarrow{\lambda }$ by choosing the pseudo-inverse of the degenerate matrix $\widehat{H}$ when computing the coefficients ${\widehat{\Lambda }}_l$.

The three constraints on infinitesimal rotations have a form (A16). These constraints define a subspace $\mathcal{S}$ of the whole space ${\mathbb{R}}_3^{\otimes N}$ of our rotation vectors $\mathbf{Re}{\overrightarrow{\mu }}_k$ (dual to elements of Lie algebra on each $SO\left(3\right)$)

$$\begin{array}{c}\mathcal{S}:\sum _k{\overrightarrow{\Theta }}_{ik}·\mathbf{Re}{\overrightarrow{\mu }}_k=0;i=1,2,3;\hfill \end{array}$$

(87)

$$\begin{array}{c}{\overrightarrow{\Theta }}_{ik}=\mathbf{Re}\left({\overrightarrow{\psi }}_i^{*}{\widehat{W}}_k\right);\hfill \end{array}$$

(88)

The rotation vectors $\mathbf{Re}{\overrightarrow{\mu }}_k$ vary in the quotient space

$$\begin{array}{c}\hfill \mathcal{F}={(}^{{\mathbb{R}}_3^{\otimes N}}{/}_{\mathcal{S}})\end{array}$$

(89)

The null-vectors ${\overrightarrow{\psi }}_i$ and coefficients ${\widehat{W}}_k$, depending on the current positions ${\overrightarrow{F}}_k$ are computed using recurrent equations in [32].

We get numerical results on a laptop for arbitrary $N<100$. Larger values of N would require a supercomputer.

9. Brownian Motion on Fixed Manifold

Now, we are ready to write down the SDE for the evolution of our complex curve using the stochastic process $d{\overrightarrow{\xi }}_l=\mathbf{Re}{\overrightarrow{\mu }}_l$:

$$\begin{array}{c}d{\overrightarrow{q}}_k=\sum _{l=0}^{N-1}{\widehat{T}}_{kl}·d{\overrightarrow{\xi }}_l;\hfill \end{array}$$

(90a)

$$\begin{array}{c}{\widehat{T}}_{kl}={\delta }_{kl}{\widehat{q}}_k+ı{\widehat{q}}_k·{\widehat{S}}_{kl};\hfill \end{array}$$

(90b)

$$\begin{array}{c}d{\overrightarrow{F}}_0=\sum _{l=0}^{N-1}{\widehat{\Lambda }}_l·d{\overrightarrow{\xi }}_l;\hfill \end{array}$$

(90c)

$$\begin{array}{c}{\overrightarrow{F}}_k={\overrightarrow{F}}_0+\sum _{l=0}^{k-1}{\overrightarrow{q}}_l;\hfill \end{array}$$

(90d)

$$\begin{array}{c}\sum _{k=0}^{N-1}{\overrightarrow{\Theta }}_{ik}·d{\overrightarrow{\xi }}_k=0;i=1,2,3;\hfill \end{array}$$

(90e)

These constrained stochastic differential equations describe the evolution of the point on our fixed manifold ${\mathbb{T}}_3\left(N\right)$ of closed complex curves subject to the loop equations (51), starting with one of the symmetric fixed points (72)

$$\begin{array}{c}\hfill {\left.{\overrightarrow{F}}_k\right|}_{\tau =0}={\overrightarrow{\Phi }}_k\mathit{or}{\overrightarrow{\Phi }}_{N-k}^{*}\end{array}$$

(91)

The constrained SDE was studied in the mathematical literature [33].

We use a standard method of the projection of the Brownian motion to a quotient space. Let us introduce new stochastic real vector variables $d\widehat{\eta }=\{d{\overrightarrow{\eta }}_0,\cdots d{\overrightarrow{\eta }}_{N-1}\}\in {\mathbb{R}}_3^{\otimes N}$ and project out the constraints as follows (in matrix notations)

$$\begin{array}{c}d\widehat{\xi }=d\widehat{\eta }-\mathcal{P}·d\widehat{\eta };\hfill \end{array}$$

(92)

$$\begin{array}{c}\mathcal{P}={\widehat{\Theta }}^{†}·{\left(\widehat{\Theta }·{\widehat{\Theta }}^{†}\right)}^{-1}·\widehat{\Theta };\hfill \end{array}$$

(93)

The variables $d\widehat{\eta }$ are assumed to be delta correlated (in proper units of stochastic time)

$$\begin{array}{c}\hfill 〈\frac{d\widehat{\eta }}{d\tau }\otimes \frac{d\widehat{\eta }}{d{\tau }^{\prime }}〉=\widehat{I}\delta (\tau -{\tau }^{\prime })\end{array}$$

(94)

Here $\widehat{I}$ is a unit matrix in $3N$ dimensions.

It is straightforward to check that $d\xi$ satisfies the constraints for arbitrary $d\widehat{\eta }$

$$\begin{array}{c}\hfill \widehat{\Theta }·d\widehat{\xi }=0\end{array}$$

(95)

The variables $d\widehat{\xi }$ do not change when variables $d\widehat{\eta }$ are shifted by superposition of transposed constraints

$$\begin{array}{c}\delta d\widehat{\eta }={\widehat{\Theta }}^{†}·dw;\hfill \end{array}$$

(96)

$$\begin{array}{c}\delta d\widehat{\xi }=0\hfill \end{array}$$

(97)

So, our stochastic process $d\widehat{\eta }$ has some redundant (gauge) degrees of freedom $dw$.

The variables $d\widehat{\xi }$ evolve in the quotient space $\mathcal{F}$, covering it with an $O\left(3\right)$ invariant measure. This invariance is easy to check by noticing that all the matrices $\widehat{\Theta },\widehat{\Lambda },\widehat{S},\widehat{T}$ in our equations are made of rotation-covariant parameters in the linearized recurrent equations. These parameters are direct products of vectors times some dot products of other vectors.

In mathematical terms, $d\widehat{\eta }$ is a Wiener process in ${\mathbb{R}}_3^{\otimes N}$ with a unit variance matrix, and $d\widehat{\xi }$ is a Brownian motion in the quotient space $\mathcal{F}$. This quotient space evolves with stochastic time, as the constraint matrix $\widehat{\Theta }$ depends on current values of all vectors ${\overrightarrow{F}}_{.}$.

The projection can be used to redefine the matrices

$$\begin{array}{c}\mathcal{T}=\widehat{T}-\mathcal{P}·\widehat{T};\hfill \end{array}$$

(98)

$$\begin{array}{c}\mathcal{L}=\widehat{\Lambda }-\mathcal{P}·\widehat{\Lambda };\hfill \end{array}$$

(99)

after which our SDE takes a usual form

$$\begin{array}{c}d{\overrightarrow{q}}_k=\sum _{l=0}^{N-1}{\mathcal{T}}_{kl}·d{\overrightarrow{\eta }}_l;\hfill \end{array}$$

(100)

$$\begin{array}{c}d{\overrightarrow{F}}_0=\sum _{l=0}^{N-1}{\mathcal{L}}_l·d{\overrightarrow{\eta }}_l;\hfill \end{array}$$

(101)

$$\begin{array}{c}{\overrightarrow{F}}_k={\overrightarrow{F}}_0+\sum _{l=0}^{k-1}{\overrightarrow{q}}_l;\hfill \end{array}$$

(102)

We propose this stochastic process as a definition of the fixed manifold of decaying Turbulence.

The proof of this conjecture and extension to higher dimensions is left for a detailed mathematical study, which is beyond the scope of this work.

We coded these SDE in [32] using Mathematica^®. This code is convenient for theoretical development, but the optimized computations should be translated into Python and C++ and run on a supercomputer or a Tensorflow cluster.

Once we fix the initial value at one of the two mirror fixed points ${\overrightarrow{\Phi }}_k,{\overrightarrow{\Phi }}_{N-k}^{*}$, the evolution is unambiguous, unlike the global description of the manifold in Section 6, where we had to choose between four solutions of two quadratic equations for the point $\{{\alpha }_1,{\alpha }_2,\beta \}$ in de Sitter space $d{S}_2$.

We are still left with a choice of one of the two mirror solutions or, in the general case, the coefficients of their linear superposition in the Wilson loop.

Such linear superposition will still solve the loop equation (7a), as this equation is linear in loop space.

Section Mirror solution and inequality for the Wilson loop

There is an obvious problem with the solution we have found. The loop equation for $\overrightarrow{P}\left(\theta \right)$ is complex, and so is the solution, particularly the vorticity in (19). Since the equation for $\overrightarrow{P}$ is nonlinear, we cannot take a real part of $\overrightarrow{P}$.

The negative imaginary part of the circulation in momentum space may lead to violation of inequality $\left|\Psi \left[C\right]\right|\le 1$. The following restriction must be imposed to preserve this inequality

$$\begin{array}{c}\hfill \oint d\overrightarrow{C}\left(\theta \right)·\mathbf{Im}\overrightarrow{P}\left(\theta \right)\ge 0;\end{array}$$

(103)

Such restriction on a fractal curve $\overrightarrow{P}\left(\theta \right)$ would contradict its assumed universality. If this curve is related to $C\left(\theta \right)$, it would be involved in variations when the loop C is varied in the loop equation.

Here is a resolution of this paradox that we have found.

In the previous Sections, we described two mirror solutions, originating in (72) and evolving by an SDE (90a).

For any particular loop, we have to choose the solution with the positive imaginary part of the circulation

$$\begin{array}{c}\Psi \left[C\right]=〈exp\left(\frac{ı\Gamma }{\nu }\right)\Theta (\mathbf{Im}\Gamma )〉+\{\Gamma \iff \tilde{\Gamma }\};\hfill \end{array}$$

(104a)

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

(104b)

$$\begin{array}{c}\tilde{\Gamma }=\oint d\overrightarrow{C}\left(\theta \right)·{\overrightarrow{P}}^{*}(2\pi -\theta );\hfill \end{array}$$

(104c)

The averaging $〈\cdots 〉$ corresponds to averaging over the stochastic process or, equivalently, over the stochastic time $\tau$.

At any moment of stochastic time, the inequality restricts the loop C, but not the momentum loop $\overrightarrow{P}$: for some loops C, the circulation $\Gamma$ has a positive imaginary part; for other loops, the reflected circulation $\tilde{\Gamma }$ does.

This choice is like selecting a decaying wave function for the bound state in the Schrödinger equation for a quantum potential problem.

The theta functions in this solution represent certain boundary conditions for the loop functional in the areas (if they exist) where $\mathbf{Im}\Gamma =0$ or $\mathbf{Im}\tilde{\Gamma }=0$.

One could imagine a situation where both circulations have negative imaginary parts; the Wilson loop will be zero in these regions of parameter space.

We leave investigating these boundary conditions for a future complete mathematical theory. This exploratory paper has different goals.

Section Vorticity Distribution and Energy Dissipation

The simplest quantity to compute in our theory is the local vorticity distribution.

As we shall see, it determines the energy dissipation rate.

The local vorticity for our decaying solution of the loop equation

$$\begin{array}{c}\overrightarrow{\omega }=\frac{-ı\overrightarrow{F}\left(\theta \right)\times \Delta \overrightarrow{F}\left(\overrightarrow{\theta }\right)}{2(t+t_0)};\hfill \end{array}$$

(105)

Here $\theta$ is an arbitrary point at the loop, which makes this expression a random variable.

Note that viscosity is canceled here, as it should be by dimensional counting (vorticity has the dimension of $1/t$).

In our random walk representation, the complex vorticity operator

$$\begin{array}{c}{\overrightarrow{\omega }}_k=\frac{-ı{\overrightarrow{G}}_k}{2(t+t_0)};\hfill \end{array}$$

(106)

$$\begin{array}{c}{\overrightarrow{G}}_k={\overrightarrow{F}}_k\times {\overrightarrow{F}}_{k+1};\hfill \end{array}$$

(107)

The time derivative of energy density in our theory is

$$\begin{array}{c}{E}^{\prime }\left(t\right)=-\frac{\nu \kappa }{4{\left(t_0+t\right)}^2};\hfill \end{array}$$

(108)

$$\begin{array}{c}\kappa =\frac{1}{N}\sum _0^{N-1}\mathbf{Re}\left(-{\overrightarrow{G}}_k^2\right);\hfill \end{array}$$

(109)

Solving this equation with boundary value $E(t=\infty )=0$ we relate $t_0$ to mean initial energy

$$\begin{array}{c}〈E\left(t\right)〉=\frac{\nu 〈\kappa 〉}{4\left(t_0+t\right)};\hfill \end{array}$$

(110)

$$\begin{array}{c}{t}_0=\frac{\nu 〈\kappa 〉}{4〈E_0〉}\hfill \end{array}$$

(111)

The probability distribution of $\kappa$ and its mean value $〈\kappa 〉$ can be computed using our random walk. For the anomalous dissipation, we need the mean enstrophy to diverge [2,18] so that viscosity is compensated in the extreme turbulent limit.

$$\begin{array}{c}\hfill 〈\kappa 〉=\infty \end{array}$$

(112)

As we shall see later, in Section 11, this happens in our numerical simulations.

The microscopic picture of this infinite enstrophy differs from the singular vortex line.

In the Euler theory, divergence came from the singularity of the classical field. However, in our dual theory, it comes from the large fluctuation of the fractal curve in momentum space.

We can now write down our result for the Wilson loop in decaying Turbulence as a functional of the contour C. Let us start with a three-dimensional case

$$\begin{array}{c}\Psi {\left[C\right]}_{t\to \infty }\to {〈exp\left(\frac{-ı\oint d\theta \overrightarrow{C}\left(\theta \right)·{\overrightarrow{F}}^{\prime }\left(\theta \right)}{\sqrt{2\nu (t+t_0)}}\right)\Theta (\Gamma )〉}_F+\mathbf{reflected};\hfill \end{array}$$

(113)

The finite steps approximation we considered above

$$\begin{array}{c}〈exp\left(-\frac{ı\oint d\theta \overrightarrow{C}\left(\theta \right)·{\overrightarrow{F}}^{\prime }\left(\theta \right)}{\sqrt{2\nu (t+t_0)}}\right)〉=\hfill \end{array}$$

$$\begin{array}{c}\underset{N\to \infty }{lim}〈exp\left(-\frac{ı}{2}\frac{{\sum }_{k=0}^{N-1}\left({\overrightarrow{C}}_{k+1}+{\overrightarrow{C}}_k\right)·{\overrightarrow{q}}_k}{2\sqrt{2\nu (t+t_0)}}\right)〉;\hfill \end{array}$$

(114)

$$\begin{array}{c}{\overrightarrow{C}}_k=\overrightarrow{C}\left(\frac{2\pi k\pi }{N}\right);\hfill \end{array}$$

(115)

For the simplest circular loop in an $xy$ plane, we have

$$\begin{array}{c}\overrightarrow{C}\left(\overrightarrow{\theta }\right)=R\{cos\theta ,sin\theta ,0\};\hfill \end{array}$$

(116)

$$\begin{array}{c}\frac{{\overrightarrow{C}}_k+{\overrightarrow{C}}_{k+1}}{2}=Rcos\left(\frac{\pi }{N}\right)\left\{cos\left(\frac{\pi (2k+1)}{N}\right),sin\left(\frac{\pi (2k+1)}{N}\right),0\right\}\hfill \end{array}$$

(117)

We observe that even at the large time $t\gg t_0$ when the asymptotic fractal curve is already in place, there is a region of parameters

$$\begin{array}{c}\hfill R\sim \sqrt{\nu t}\end{array}$$

(118)

where the Wilson loop is a nontrivial universal function of a single variable

$$\begin{array}{c}\hfill \Psi \left[C\right]\to \psi \left(\frac{R}{\sqrt{2\nu t}}\right)\end{array}$$

(119)

We can compute our prediction for this function by numerically simulating the SDE for our vectors ${\overrightarrow{F}}_k$ and wait for the results with physical or numerical experiments in conventional three-dimensional decaying Turbulence.

10. Correlation Functions

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

The correlation function reduces to a random walk with a complex weight

$$\begin{array}{c}〈\overrightarrow{\omega }\left(\overrightarrow{0}\right)\otimes \overrightarrow{\omega }\left(\overrightarrow{r}\right)〉=\frac{1}{4{(t+t_0)}^2}\hfill \end{array}$$

$$\begin{array}{c}〈\sum _{n,m}\frac{{\overrightarrow{F}}_m\times {\overrightarrow{F}}_{m+1}\otimes {\overrightarrow{F}}_n\times {\overrightarrow{F}}_{n+1}}{N^2}exp\left(\frac{ı\overrightarrow{r}·\left({\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right)}{2\sqrt{\nu (t+t_0)}}\right)〉+\mathbf{reflected};\hfill \end{array}$$

(120)

$$\begin{array}{c}{\overrightarrow{S}}_{m,n}=\frac{{\sum }_m^n{\overrightarrow{F}}_k}{(n-mmodN)};\hfill \end{array}$$

(121)

In Fourier space, this reduces to the restricted random walk:

$$\begin{array}{c}〈\overrightarrow{\omega }\otimes \overrightarrow{\omega }\left(\overrightarrow{k}\right)〉=\frac{1}{4{(t+t_0)}^2}\hfill \\ 〈\sum _{n,m}\frac{{\overrightarrow{F}}_m\times {\overrightarrow{F}}_{m+1}\otimes {\overrightarrow{F}}_n\times {\overrightarrow{F}}_{n+1}}{N^2}{\delta }^d\left(\overrightarrow{k}+\frac{{\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}}{2\sqrt{\nu (t+t_0)}}\right)〉+\mathbf{reflected};\hfill \end{array}$$

(122)

Presumably, the vorticity vectors ${\overrightarrow{F}}_m\times {\overrightarrow{F}}_{m+1}$ as well as the vectors ${\overrightarrow{S}}_{m,n}$ are distributed by some power laws in our random walk on a fixed manifold; this would lead to scaling laws with some fractal dimensions.

The numerical simulation of this correlation function would require significant computer resources.

Still, these resources are much more modest than those for full d dimensional simulations of the Navier-Stokes equation.

In our theory, the dimension of space enters as the number of components of the one-dimensional fluctuating field $\overrightarrow{F}\left(\theta \right)$ rather than the number of variables $\overrightarrow{r}\in {\mathbb{R}}_d$ in the fluctuating velocity field $\overrightarrow{v}\left(\overrightarrow{r}\right)$.

Also, note that our quantum problem of the complex random walk naturally fits quantum computer architecture. Thus, in the future, when large quantum computers would become available for researchers, we can expect a real breakthrough in numerical simulations of the loop equation.

11. Open Loop Computations

We wrote a Mathematica^® program [34] generating our random walk, starting with a random complex vector ${\overrightarrow{F}}_0$ and using random orthogonal $SO\left(3\right)$ matrix $O_k$ at every step. If there is more than one real solution for $\beta$, we have chosen the shortest step in Euclidean metric ${\left|\overrightarrow{q}\right|}^2=1+2{\beta }^2$, i.e., the one with minimal $\left|\beta \right|$.

We have chosen the simplest circular coordinate loop C in (116) and imposed the inequality (103) on the last step.

For 1000 steps, it takes a few seconds on a laptop to compute the whole path. We generate a parallel table of 100000 paths, each with 1000 steps, with various random initial vectors ${\overrightarrow{F}}_0$ with a random set of rotation matrices for each step.

The path’s closure requires a numerical solution of the SDE (90a), which we plan to implement later on a supercomputer. This open path simulation only covers the big space of the direct product of rotation matrices at every step; the true turbulent fixed point corresponds to the projection of this space onto the closure conditions. We cannot do it at the global level, only at the level of the SDE described in the previous Section.

Thus, this open-path simulation cannot be used for predictions of fractal dimensions in the scaling laws; this has to wait until the SDE simulation is performed at the supercomputer.

With these comments in mind, let us analyze the open curves’ fractal properties, discarding the closure conditions.

The simplest quantity to compute is a fractal dimension $d_f$ of this random walk, defined as

$$\begin{array}{c}\hfill \frac{1}{d_f}=\underset{N\to \infty }{lim}\frac{dlog\left|{\overrightarrow{F}}_N-{\overrightarrow{F}}_0\right|}{dlogN}\end{array}$$

(123)

The ordinary Brownian motion (linear random walk) has $d_f=2$, but our random walk is very different, mainly because the Euclidean distance of an elementary step $\left|{\overrightarrow{F}}_{k+1}-{\overrightarrow{F}}_k\right|$ in De Sitter space is unlimited from above (though it is limited by 1 from below).

Here is the plot of $log\left|{\overrightarrow{F}}_N-{\overrightarrow{F}}_0\right|$ vs $logN$ (Figure 2).

The statistical data for parameters

$$\begin{array}{ccccc}& \mathrm{Estimate}\hfill & \mathrm{Standard}\mathrm{Error}\hfill & \mathrm{t}-\mathrm{Statistic}\hfill & \mathrm{P}-\mathrm{Value}\hfill \\ 1\hfill & -2.33876\hfill & 0.0241623\hfill & -96.7941\hfill & 4.3388931215906555{*}^{\wedge }-42\hfill \\ \xi \hfill & 0.976443\hfill & 0.00457433\hfill & 213.461\hfill & 2.1004008453460444{*}^{\wedge }-53\hfill \end{array}$$

(124)

This data is compatible with $d_f=1.02412\pm 0.005$.

The distribution of the Euclidean length of each step $\left|{\overrightarrow{F}}_{k+1}-{\overrightarrow{F}}_k\right|$ (Figure 3).

The statistical table for the parameters of this fit

$$\begin{array}{ccccc}& \mathrm{Estimate}\hfill & \mathrm{Standard}\mathrm{Error}\hfill & \mathrm{t}-\mathrm{Statistic}\hfill & \mathrm{P}-\mathrm{Value}\hfill \\ 1\hfill & 13.0661\hfill & 0.00203512\hfill & 6420.3\hfill & 0.\hfill \\ log\left(\mathrm{step}\right)\hfill & -2.00076\hfill & 0.000737843\hfill & -2711.64\hfill & 0.\hfill \end{array}$$

(125)

The mean is finite, $〈step〉=1.07056$, but the variance of the step is divergent.

Such a low decay of the step distribution undermines the concept of a finite fractal dimension as defined in (123). The linear fit is inadequate for such large statistics.

With large statistics, one can reach a perfect fit by adding the next correction to the linear log-log law

$$\begin{array}{c}\hfill log\left|{\overrightarrow{F}}_N-{\overrightarrow{F}}_0\right|\approx a+blogN+cloglogN\end{array}$$

(126)

The fit at larger interval of N becomes perfect, with a very different coefficient b in front of $logN$:

$$\begin{array}{c}\hfill \begin{array}{ccccc}& \mathrm{Estimate}\hfill & \mathrm{Standard}\mathrm{Error}\hfill & \mathrm{t}-\mathrm{Statistic}\hfill & \mathrm{P}-\mathrm{Value}\hfill \\ 1\hfill & 11.4051\hfill & 0.067611\hfill & 168.688\hfill & 1.201777144099837{*}^{\wedge }-128\hfill \\ \xi \hfill & 4.85103\hfill & 0.0218164\hfill & 222.357\hfill & 4.3433787030416533{*}^{\wedge }-141\hfill \\ log\left(\xi \right)\hfill & -20.5553\hfill & 0.109809\hfill & -187.192\hfill & 2.4775697694919745{*}^{\wedge }-133\hfill \end{array}\end{array}$$

(127)

This data fit is shown at Figure 4.

Our random walk with unbounded step size differs from an ordinary fractal curve. The fractal dimension does not properly describe this random object as the distance grows by a more complex law than a pure power of the number of steps.

Another interesting distribution is the enstrophy density $\kappa$ in (108).

The CDF is shown in Figure 5. The tail is compatible with ${\kappa }^{-0.936266}$ decay, corresponding to the ${\kappa }^{-1.936266}$ decay of the PDF. The mean value and all higher moments diverge, leading to anomalous dissipation.

The statistical table for the parameters of this fit

$$\begin{array}{ccccc}& \mathrm{Estimate}\hfill & \mathrm{Standard}\mathrm{Error}\hfill & \mathrm{t}-\mathrm{Statistic}\hfill & \mathrm{P}-\mathrm{Value}\hfill \\ 1\hfill & 17.1851\hfill & 0.00330287\hfill & 5203.09\hfill & 0.\hfill \\ log\left(\kappa \right)\hfill & -0.936266\hfill & 0.000320989\hfill & -2916.81\hfill & 0.\hfill \end{array}$$

(128)

The computation of the Wilson loop and related correlation functions of vorticity needs an ensemble of closed fractal loops with various sets of random matrices.

The closure condition for the loop would require some computational effort because the probability of the random curve with fractal dimension $d_f\sim 1.$ returning to an initial point goes to zero with the increased number of steps.

An alternative approach of starting with a large closed loop ${\overrightarrow{F}}_k;{\overrightarrow{F}}_N={\overrightarrow{F}}_0$ and randomizing it point by point while preserving its closure.

This approach would replace an SDE (90a) with a Monte-Carlo process in a closed polygon space. Each step would correspond to a small rotation of some randomly chosen link ${\overrightarrow{q}}_k$ of the polygon and compensating shift of the initial position ${\overrightarrow{F}}_0$

The zeroth approximation to this shift would be our solution in Appendix A of the linearized equations. Then, a couple of Newton iterations will finalize this shift to provide the closure condition.

These extra layers of computational complexity would require a supercomputer, which we plan to do later.

12. Discussion

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

The solution expresses the probability distribution for the Wilson loop at any given moment of the time t in terms of a nonlinear SDE in auxiliary time $\tau$. The loop is approximated as a polygon with $N\to \infty$ sides.

Compared to the original Navier-Stokes equation, this is the reduction of $d\Rightarrow 1$ of the dimension of space. This SDE is straightforward to simulate by a Monte Carlo method.

The equivalence of a strong coupling phase of the fluctuating vector field to quantum geometry is a well-known phenomenon in gauge theory (the ADS/CFT duality).

In our case, this is a much simpler quantum geometry: a fractal curve in complex space.

An expert in the traditional approach to Turbulence may wonder why the solutions of the Loop equation have anything in common with the statistical distribution of the velocity field in a decaying turbulent flow.

Such questions were asked and answered in the gauge theories, including QCD[6,8,9,10] where the loop equations were derived first [4,5].

The solution of the loop equation with finite area derivative, satisfying Bianchi constraint, belongs to the so-called Stokes-type functionals [4], the same as the Wilson loop for Gauge theory and fluid dynamics.

Extra complications in the gauge theory, which are fortunately absent in fluid dynamics, are the short-distance singularities related to the infinite number of fluctuating degrees of freedom in quantum field theory.

As we discussed in detail in [2,4,5], any Stokes-type functional $\Psi \left[C\right]$ satisfying boundary condition at shrunk loop $\Psi \left[0\right]=1$, and solving the loop equation can be iterated in the nonlinear term in the Navier-Stokes equations (which applies at large viscosity).

The resulting expansion in inverse powers of viscosity (weak Turbulence) exactly coincides with the ordinary perturbation expansion of the Navier-Stokes equations for the velocity field, averaged over the distribution of initial data or boundary conditions at infinity.

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

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

Our general Anzatz (14) satisfies the loop equation and boundary condition at $\Psi [C=0]=1$. It has a finite area derivative, which obeys the Bianchi constraint, making it a Stokes-type functional.

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

Therefore, it represents a statistical distribution in a turbulent Navier-Stokes flow, a degenerate fixed point of the Hopf equation for 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. Whether this exact solution is realized in Nature remains to be seen.

It is infinitely more complex than the randomly rotated fluid, as our curve $\overrightarrow{P}\left(\theta \right)$ has a discontinuity at every $\theta$, corresponding to a distributed random vorticity.

This solution is described by a fractal curve in complex d dimensional space, a limit of a random walk with nonlinear algebraic relations between the previous position ${\overrightarrow{F}}_k$ and the next one ${\overrightarrow{F}}_{k+1}$. These relations are degenerate: each step ${\overrightarrow{q}}_k={\overrightarrow{F}}_{k+1}-{\overrightarrow{F}}_k$ is characterized by an arbitrary element ${\widehat{O}}_k\in {(}^{O\left(d\right)}{/}_{O(d-3)})$ and an arbitrary element ${\overrightarrow{w}}_k\in {(}^{{\mathbb{S}}^{d-3}}{/}_{{\mathbb{Z}}^2})$. This step also depends upon the previous position ${\overrightarrow{F}}_k$, making this process a Markov chain.

The periodicity condition ${\overrightarrow{F}}_N={\overrightarrow{F}}_0$ provides a nonlinear equation for an initial position ${\overrightarrow{F}}_0$ as a function of the above free parameters ${\widehat{O}}_k,{\overrightarrow{w}}_k$.

This periodicity condition presents a hard problem, particularly in the limit $N\to \infty$, when the probability of our random walk returning to the initial point afterN steps rapidly diminishes with N.

We found a way around this problem by finding an exact periodic solution to the momentum loop equations (51). This analytical solution for arbitrary N can serve as initial data for the SDE, which preserves the periodicity. In Appendix A, we constructed this SDE for $d=3$, leaving the generalizations to mathematicians.

This SDE describes the Brownian motion of the rotation matrices ${\widehat{O}}_k\in SO\left(3\right)$ in our canonical representation (62) of the solution to the discrete loop equations (51). Each matrix moves independently, while the remaining parameters $\{{\alpha }_1,{\alpha }_2,\beta \}$ move around de Sitter space $d{S}_2$ to satisfy the loop equation (51).

The closure condition further restricts the set of N infinitesimal rotations $\delta {\overrightarrow{q}}_k=\delta {\overrightarrow{\theta }}_k\times {\overrightarrow{q}}_k$: there are three linear relations between N vector parameters $\delta {\overrightarrow{\theta }}_k$ of these rotations. We found the projection matrix required to project the whole array of vector rotations $\delta {\overrightarrow{\theta }}_k$ onto the quotient space, satisfying the closure condition.

After this projection, we obtain the motion in the quotient space, where the closure condition is satisfied at every step.

Presumably, this SDE uniformly covers our fixed manifold ${\mathbb{T}}_3\left(N\right)$ for arbitrary N. The limit $N\to \infty$ presents a computational challenge, and we are planning to address this challenge in the next publication using a supercomputer.

We simulated the open random walk (without the closure condition) in three dimensions and studied its statistical properties. The distribution of lengths of steps in Euclidean six-dimensional space $\mathbf{Re}\overrightarrow{F},\mathbf{Im}\overrightarrow{F}$ has a long tail $\mathit{PDF}\propto x^{-2}$.

The fractal dimension is not well defined for a random walk with such an intermittent step size, unbounded from above. The linear log-log fit as in (123) yields $d_f\approx 1.20$, but this fit is imperfect with our large statistics.

As for the distribution of an enstrophy density, it has a power tail $x^{-1.9}$ corresponding to an infinite mean value and all higher moments. This infinity is how anomalous dissipation manifests in our solution.

These numerical simulations must be repeated on a supercomputer with better statistics and more steps. There are many things to do next with this conjectured solution to the decaying turbulence problem; the first is to look for unnoticed inconsistencies.

One important step is yet to be made: the MC simulation of the SDE (90a). Let us assume that the qualitative properties and fractal dimensions we have found for the open fractal curves will stay the same or at least close.

As a first test of this hypothesis, let us compare it with various experimental data and those from DNS [35].

There is no agreement between these data, they vary in Reynolds number, and they have other differences related to the experimental setup. No value n for the decay power $t^{-n}$ would fit all that data. However, a consensus seems to be around $n\approx 1.2-1.4$, which means faster decay than we have.

We are skeptical about these data. As we recently learned [36], there is a regime change at large Reynolds numbers; the numbers achievable in modern DNS may belong to such a transitional regime.

Besides, fitting powers is not a reliable method of deriving physical laws.

For example, we took a formula $1/(t-0.5)$, added random noise between $(-0.1,0.1)$ and fitted this data to $b{t}^{-n}$. The best fit produced some fake power $n\approx 1.43$ and some fake coefficient $b\approx 1.88$ in front (see Figure 6)

Instead, one should compare a hypothetical theory with a null hypothesis by estimating the log-likelihood of both fits. In case the new theory is more likely as an explanation of the data, you may temporarily accept it until a better theory or better data will appear.

A good history lesson is fitting the power n in Newton’s gravity law to explain the astronomic data for the Mercury perihelion before the General Relativity theory. A small correction to $n=1$ "explained" the data, but this was useless without a theory.

Presumably, our fixed point corresponds to a true infinite Reynolds limit, as it is completely universal and does not depend on the Reynolds scales.

If you assume no hidden scales are left, our $E\propto \nu /t$ law follows from dimensional analysis. Observed or simulated data with $n>1$ all have the powers of some other dimensional parameters related to the Reynolds number. They rely on (multifractal versions of) K41 spectra and other intermediate turbulent phenomena.

We have an anomalous dissipation rate: the mean value of the vorticity square diverges, compensating for the viscosity factor in the energy decay in extreme turbulent limit.

This mechanism of anomalous dissipation differs from the one we studied in the Kelvinon [2,18]. In those fixed points, the viscosity canceled in the dissipation rate due to the singular vorticity configurations with the thin vortex line resolved as a core of a Burgers vortex.

Here, in the dual theory of fractal momentum loop, the large fluctuations of this momentum loop would lead to the divergent expectation value of the enstrophy.

Our solution is universal, rotational, and translational invariant. It has the expected properties of extreme isotropic Turbulence. Is it THE solution? Time will tell.