Decaying Turbulence as a Fractal Curve

1. Introduction

A while ago, we derived [1,3] a functional equation for the so-called loop average [4,5] or Wilson loop in Turbulence. The path to an exact solution by a dimensional reduction in this equation was proposed in the ’93 paper [1], but we have only started exploring this path.

At the time, we could not compare a theory with anything but crude measurements in physical and numerical experiments at modest Reynolds numbers. All these experiments agreed with the K41 scaling, so the exotic equation based on unjustified methods of quantum field theory was premature.

The specific prediction of the Loop equation, namely the Area law [1], could not be verified in DNS at the time with existing computer power.

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

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

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

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

All these old and new developments made loop equations a major nonperturbative approach to gauge field theory.

So, it is time to revive the hibernating theory of the loop equations in Turbulence, where these equations are much simpler.

The latest DNS [13,14,15,16] with Reynolds numbers of tens of thousands revealed and quantified violations of the K41 scaling laws. These numerical experiments are in agreement with so-called multifractal scaling laws [17].

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

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

This law says that the PDF tails of velocity circulation PDF in the confinement region are functions of the minimal area inside this loop.

It was verified in DNS four years ago [13] which triggered the further development of the geometric theory of turbulence[2,14,15,16,18,19,20,21,22,23,24,25,26,27,28,29,30].

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

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

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

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

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

$$\begin{array}{ccc}\hfill & & \overrightarrow{\omega }=\frac{1}{2}Z{e}_{abc}{S}_a\overrightarrow{\nabla }{S}_b\times \overrightarrow{\nabla }{S}_c=\overrightarrow{\nabla }{\varphi }_1\times \overrightarrow{\nabla }{\varphi }_2;\hfill \end{array}$$

(1)

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

(2)

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

(3)

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

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

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

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

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

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

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

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

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

In the present work, we develop the theory free of these assumptions by exactly solving the loop equations for decaying turbulence. Our representation of velocity circulation does not run into the problems of the Kelvinon gas model; nor do we make any approximations in the loop equations we are solving.

2. 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}{ccc}& & {\partial }_t{v}_{\alpha }=\nu {\partial }_{\beta }{\omega }_{\beta \alpha }-v_{\beta }{\omega }_{\beta \alpha }-{\partial }_{\alpha }\left(p+\frac{v_{\beta }^2}{2}\right);\hfill \end{array}$$

(4)

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

(5)

The Wilson loop average for the Turbulence

$$\begin{array}{c}\hfill \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}{ccc}\hfill & & ı\nu {\partial }_t\Psi ={\mathcal{H}}_C\Psi ;\hfill \end{array}$$

(7a)

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

(7b)

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

(7c)

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

(7d)

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

(7e)

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

(7f)

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

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

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

(8)

$$\begin{array}{ccc}& & {\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}{ccc}\hfill & & {\Omega }_{\alpha \beta }(\theta ,ϵ)=\hfill \\ & & -ı\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 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.

3. 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 a simple Ansatz

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

(13)

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 removes an objection we made in the Introduction to the Kelvinon theory as well as any other Navier-Stokes stationary solutions, with a singularity at fixed loop C in space.

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

(14)

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

(15)

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

This transformation yields the following formula for the spike derivative D in the above equation

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

(16)

The vorticity operator (11) and velocity (12) become ordinary numbers, which are some singular functionals of the trajectory $P\left(\theta \right)$.

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

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

(17)

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

(18)

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

(19)

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

(20)

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

(21)

We arrive at a singular equation

$$\begin{array}{ccc}\hfill & & \nu {\partial }_t\overrightarrow{P}=-{\left(\Delta \overrightarrow{P}\right)}^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}$$

(22)

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.

4. Random Global Rotation

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

$$\begin{array}{ccc}& & {\Psi }_0\left[C\right]=\hfill \\ & & 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}$$

(23)

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

(24)

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

(25)

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

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

(26)

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

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

(27)

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

The constant term in (25) 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 & & {\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}$$

(28)

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

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

(29)

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

The physical meaning of this $\omega$ is the random uniform vorticity at the initial moment.

Unfortunately, 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 $\omega$ 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 (28).

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

(30)

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

(31)

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

(32)

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

(33)

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

(34)

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 $\omega$ the correlations are

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

(35)

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

(36)

5. 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}{ccc}\hfill \mathbf{Decay}:& & \overrightarrow{P}\to 0;\Psi \left[C\right]\to 1;\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill \mathbf{Fixed}\mathbf{Point}:& & \overrightarrow{P}\to {\overrightarrow{P}}^{*};\Psi \left[C\right]\to {\Psi }^{*}\left[C\right];\hfill \end{array}$$

(38)

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

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

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

(39)

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

(40)

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

$$\begin{array}{ccc}\hfill & & {\left(\Delta \overrightarrow{P}\right)}^2\left({\left(\Delta \overrightarrow{P}\right)}^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{\left(\Delta \overrightarrow{P}\right)}^2\right);\hfill \end{array}$$

(41)

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

(42)

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

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 (18), leading to a trivial solution of the loop equation $\Psi \left[C\right]=1$.

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

6. Random Walk 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}{ccc}& & \Delta \overrightarrow{F}\stackrel{?}{=}\lambda \overrightarrow{F};\hfill \end{array}$$

(43)

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

(44)

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}& & {\left(\Delta \overrightarrow{F}\right)}^2=1;\hfill \end{array}$$

(45a)

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

(45b)

These relations are 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 (Section 6) is completely universal!

One can iteratively build this solution as a fractal curve 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 (Section 6).

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

(46a)

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

(46b)

The solution of 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}$$

(47)

after which we can find the next value

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

(48)

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

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

(49)

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

(50)

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

(51)

In d dimensions, there are $d-1$ complex parameters of the unit vector; with an extra linear constraint in (47), 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}{ccc}& & \sum _0^{N-1}{\overrightarrow{q}}_k=0;\hfill \end{array}$$

(52)

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

We have here a nonlinear 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$.

The initial position ${\overrightarrow{F}}_0\in {\mathbb{C}}^d$ is adjusted to close the trajectory after N steps at given rotations.

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.

7. The Probability measure

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

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

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

(53a)

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

(53b)

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

(53c)

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

(53d)

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}{ccc}& & {\overrightarrow{f}}_k={\widehat{O}}_k^T·{\overrightarrow{F}}_k;\hfill \end{array}$$

(54)

as follows:

$$\begin{array}{ccc}\hfill & & \{{\alpha }_1,{\alpha }_2\}={\widehat{M}}^{-1}.\hfill \end{array}$$

$$\begin{array}{ccc}& & \left\{\mathbf{Re}\right(R)-\beta \mathbf{Im}(c),\beta \mathbf{Re}(c)+\mathbf{Im}(R\left)\right\};\hfill \end{array}$$

(55)

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

(56)

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

(57)

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 () and two more signs in a solution of the quadratic equation for $\beta$ (53a). We have to choose one of these four solutions at each step, also selecting the real solution for $\beta$.

The simple universal choice, which works in numerical simulations, is the solution with the smallest Euclidean length ${\left|\overrightarrow{q}\right|}^2=\mathbf{Re}{\overrightarrow{q}}^2+\mathbf{Im}{\overrightarrow{q}}^2=1+2{\beta }^2$.

The motivation for this minimization is that there is no upper limit on the solution for $\beta$. The unminimized random walk may rapidly go to infinity due to instabilities.

With this minimization restriction, the step size s distribution we find below in numerical simulation decays as $s^{-3.27}$, sufficient for finite mean and variance.

We also found that the mean Euclidean distance $\left|{\overrightarrow{F}}_N-{\overrightarrow{F}}_0\right|$ of unconstrained random walk grows approximately linearly with the number N of steps, not showing any sign of instabilities.

We conclude that the fixed manifold of the decaying Turbulence is a tensor product of $SO\left(d\right)\otimes {\mathbb{S}}^{d-3}$ spaces

$$\begin{array}{c}\hfill \mathbf{Fixed}\mathbf{Manifold}\left(d\right)={\left(SO\left(d\right)\otimes {\mathbb{S}}^{d-3}\right)}^{\otimes N}\end{array}$$

(58)

We arrive at the invariant distribution for our fractal curve. At a fixed number of steps, N

$$\begin{array}{ccc}& & d{\Omega }_d\left[F\right]=\prod _{k=0}^{N-1}\left(d{\widehat{O}}_k\right)\left(d{\overrightarrow{w}}_k\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & Z\left({\widehat{O}}_1,\cdots {\widehat{O}}_{N-1},{\overrightarrow{w}}_1,\cdots {\overrightarrow{w}}_{N-1}\right);\hfill \end{array}$$

(59a)

$$\begin{array}{ccc}& & {\overrightarrow{F}}_{k+1}={\overrightarrow{F}}_k+{\overrightarrow{q}}_k;\forall k=0,\cdots N-1\hfill \\ & & Z\left({\widehat{O}}_1,\cdots {\widehat{O}}_{N-1},{\overrightarrow{w}}_1,\cdots {\overrightarrow{w}}_{N-1}\right)=\hfill \end{array}$$

(59b)

$$\begin{array}{ccc}& & \int d^{2d}{\overrightarrow{F}}_0{\delta }^{2d}\left(\overrightarrow{Q}\right)\hfill \end{array}$$

(59c)

$$\begin{array}{ccc}& & \overrightarrow{Q}={\overrightarrow{F}}_N-{\overrightarrow{F}}_0=\sum {\overrightarrow{q}}_k;\hfill \end{array}$$

(59d)

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

The complex vector’s integration and delta function is understood as a product of its real and imaginary parts. This integration leads to a Jacobian factor in the measure, depending on the relative angles between all rotations.

$$\begin{array}{ccc}\hfill & & Z\left({\widehat{O}}_1,\cdots {\widehat{O}}_{N-1},{\overrightarrow{w}}_1,\cdots {\overrightarrow{w}}_{N-1}\right)=\hfill \end{array}$$

$$\begin{array}{ccc}& & \frac{1}{\left|det\frac{\partial \overrightarrow{Q}}{\partial {\overrightarrow{F}}_0}\right|};\hfill \end{array}$$

(60)

$$\begin{array}{ccc}& & det\frac{\partial \overrightarrow{Q}}{\partial {\overrightarrow{F}}_0}\equiv \frac{\partial \left(\mathbf{Re}\overrightarrow{Q},\mathbf{Im}\overrightarrow{Q}\right)}{\partial \left(\mathbf{Re}{\overrightarrow{F}}_0,\mathbf{Im}{\overrightarrow{F}}_0\right)};\hfill \end{array}$$

(61)

This closure condition is sufficient by parameter count to eliminate a complex vector ${\overrightarrow{F}}_0$; however, doing this at an arbitrary number of rotation matrices is a nontrivial technical problem.

We will address this and other technical problems in three dimensions using Mathematica^®.

The remarkable universality of this fractal curve in complex space makes these problems worth an effort.

As for the computation of the Jacobian in (60), it is easier than it looks.

There is a recurrent equation (in matrix notations, with complex vector treated as 6 dimensional real vector)

$$\begin{array}{ccc}& & \frac{\partial {\overrightarrow{F}}_{k+1}}{\partial {\overrightarrow{F}}_k}=\widehat{I}+\frac{\partial {\overrightarrow{q}}_k}{\partial {\overrightarrow{F}}_k}=\widehat{I}+{\widehat{O}}_k·\frac{\partial \overrightarrow{u}}{\partial {\overrightarrow{f}}_k}·{\widehat{O}}_k^T\hfill \end{array}$$

(62)

Finally, using the chain rule, we compute the Jacobian (the orthogonal matrices ${\widehat{O}}_k,{\widehat{O}}_k^T$ all cancel)

$$\begin{array}{ccc}& & det\left(-\widehat{I}+\frac{\partial {\overrightarrow{F}}_N}{\partial {\overrightarrow{F}}_0}\right)=det\left(-\widehat{I}+\prod _{k=0}^{N-1}\frac{\partial {\overrightarrow{F}}_{k+1}}{\partial {\overrightarrow{F}}_k}\right)=\hfill \\ & & det\left(-\widehat{I}+\prod _{k=1}^{N-1}\left(\widehat{I}+\frac{\partial \overrightarrow{u}}{\partial {\overrightarrow{f}}_k}\right)\right)\hfill \end{array}$$

(63)

8. Reflected solution and Wilson loop convergence

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 (18). Since the equation for $\overrightarrow{P}$ is nonlinear, we cannot add a complex conjugate solution. Besides, the complex conjugate solution would not satisfy the same nonlinear equation (Section 6) with complex coefficients.

As we shall see later, the distribution of the random variables ${\overrightarrow{F}}_k$ decays only as a power. This makes the exponential in the Wilson loop grow too fast for the power distribution to converge.

However, another solution corresponds to solving the recurrent equations (Section 6) in the opposite direction on the loop. In the continuum limit, it corresponds to the orientation reversal $\theta \to 2\pi -\theta$.

A discrete version corresponds to counting points in the opposite direction on a loop. In the equations (Section 6), only the sign in front of the gap $\Delta \overrightarrow{P}$ will change, equivalent to a complex conjugation of the equation. In other words, the solution with the opposite orientation is the complex conjugate ${\overrightarrow{P}}^{*}\left(\theta \right)$ of the original solution $\overrightarrow{P}\left(\theta \right)$.

One of these two solutions would correspond to the exponential growth of the Wilson loop, and another one - to exponential decrease

$$\begin{array}{ccc}& & \left|exp\left(\frac{ı}{\nu }\oint d{C}_{\alpha }\left(\theta \right)P_{\alpha }\left(\theta \right)\right)\right|\hfill \end{array}$$

$$\begin{array}{ccc}& & =exp\left(\frac{-1}{\nu }\oint d{C}_{\alpha }\left(\theta \right)\mathbf{Im}{P}_{\alpha }\left(\theta \right)\right);\hfill \\ & & \left|exp\left(\frac{ı}{\nu }\oint d{C}_{\alpha }\left(\theta \right)P_{\alpha }^{*}\left(\theta \right)\right)\right|\hfill \end{array}$$

(64)

$$\begin{array}{ccc}& & =exp\left(\frac{1}{\nu }\oint d{C}_{\alpha }\left(\theta \right)\mathbf{Im}{P}_{\alpha }\left(\theta \right)\right);\hfill \end{array}$$

(65)

The increment of exponential growth/decrease generally depends on the loop C.

We conclude that only one of these two solutions is acceptable in various regions of the loop space C.

Another problem is that we expect the correlation functions of vorticity (see the next section) to be real. Presumably, the ambiguities of the probability measure can be used to cancel imaginary parts of observable expectation values.

This issue remains unsolved and requires further investigation.

9. 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}{ccc}& & \overrightarrow{\omega }=\frac{-ı\overrightarrow{F}\left(\theta \right)\times \Delta \overrightarrow{F}\left(\overrightarrow{\theta }\right)}{2(t+t_0)};\hfill \end{array}$$

(66)

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}{ccc}& & {\overrightarrow{\omega }}_k=\frac{-ı{\overrightarrow{G}}_k}{2(t+t_0)};\hfill \end{array}$$

(67)

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

(68)

The time derivative of energy density in our theory is

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

(69)

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

(70)

$$\begin{array}{ccc}& & {〈X\left[F\right]〉}_F=\frac{\int d{\Omega }_3\left[F\right]X\left[F\right]}{\int d{\Omega }_3\left[F\right]}\hfill \end{array}$$

(71)

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

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

(72)

$$\begin{array}{ccc}& & t_0=\frac{\nu 〈\kappa 〉}{4〈E_0〉}\hfill \end{array}$$

(73)

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

(74)

As we shall see later, in Section 11, this is what 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}{ccc}& & \Psi {\left[C\right]}_{t\to \infty }\to \frac{1}{2}{〈exp\left(\frac{ı\oint d\theta {\overrightarrow{C}}^{\prime }\left(\theta \right)·\overrightarrow{F}\left(\theta \right)}{\sqrt{2\nu (t+t_0)}}\right)〉}_F\hfill \\ & & +\frac{1}{2}\{F\leftrightarrow F^{*}\};\hfill \end{array}$$

(75)

The finite steps approximation we considered above

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

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

(76)

$$\begin{array}{ccc}& & \Delta {\overrightarrow{C}}_k=\overrightarrow{C}\left(\frac{(2k+1)\pi }{N}\right)-\overrightarrow{C}\left(\frac{(2k-1)\pi }{N}\right);\hfill \end{array}$$

(77)

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

$$\begin{array}{ccc}& & \overrightarrow{C}\left(\overrightarrow{\theta }\right)=R\{cos\theta ,sin\theta ,0\};\hfill \\ & & \Delta {\overrightarrow{C}}_k=2Rsin\left(\frac{\pi }{N}\right)\hfill \end{array}$$

(78)

$$\begin{array}{ccc}& & \left\{-sin\left(\frac{2\pi k}{N}\right),cos\left(\frac{2\pi k}{N}\right),0\right\}\hfill \end{array}$$

(79)

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

(80)

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

(81)

We can compute our prediction for this function by numerically solving the recurrent equations 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}{ccc}& & 〈\overrightarrow{\omega }\left(\overrightarrow{0}\right)\otimes \overrightarrow{\omega }\left(\overrightarrow{r}\right)〉=\frac{1}{4{(t+t_0)}^2}\hfill \\ & & {〈\sum _{n,m}\frac{{\overrightarrow{G}}_m\otimes {\overrightarrow{G}}_n}{N^2}exp\left(\frac{ı\overrightarrow{r}·\left({\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right)}{2\sqrt{\nu (t+t_0)}}\right)〉}_F\hfill \end{array}$$

(82)

$$\begin{array}{ccc}& & +\{c.c\};\hfill \end{array}$$

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

(83)

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

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

(84)

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

We wrote a Mathematica^® program [32] 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.

For 1000 steps, it takes a few seconds on a laptop to compute the whole path.

The closure of the path requires special efforts. We create a frozen set of random orthogonal matrices ${\widehat{O}}_k$ for the path with N steps.

At each step ${\overrightarrow{F}}_{k+1}={\overrightarrow{F}}_k+{\widehat{O}}_k·{\overrightarrow{q}}_k$ we use the frozen random matrix ${\widehat{O}}_k$ and choose from four solutions for ${\overrightarrow{F}}_{k+1}$ the one closest to the initial point ${\overrightarrow{F}}_0$ by Euclidean distance in six dimensions ${\left|\overrightarrow{X}\right|}^2={\left(\mathbf{Re}\overrightarrow{X}\right)}^2+{\left(\mathbf{Im}\overrightarrow{X}\right)}^2$.

Then we generate a parallel table of 10000 paths with various random initial vectors ${\overrightarrow{F}}_0$.

We use a Monte-Carlo minimization of the absolute distance between initial and final vectors $\left|{\overrightarrow{F}}_N-{\overrightarrow{F}}_0\right|$. Altogether, it took $10000\times 1000={10}^7$ steps of a random walk.

The convergence to zero of the sequence of these absolute distances is shown in Figure 2 in logarithmic scale.

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

(85)

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

The statistical data for parameters

$$\begin{array}{ccc}& \mathrm{Estimate}& \mathrm{Standard}\mathrm{Error}\\ 1& -2.13468& 0.0400099\\ \xi & 0.938229& 0.00758189\end{array}$$

(86)

This data is compatible with $d_f=1.066\pm 0.007$.

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

The statistical table for the parameters of this fit

$$\begin{array}{ccc}\hfill & \mathrm{Estimate}\hfill & \mathrm{Standard}\mathrm{Error}\hfill \\ 1\hfill & 11.1085\hfill & 0.00231018\hfill \\ log\left(s\right)\hfill & -2.15408\hfill & 0.00132643\hfill \end{array}$$

(87)

The mean value and error of a step is

$$\begin{array}{c}\hfill s=1.07048\pm 0.559112\end{array}$$

(88)

Another interesting distribution is $\kappa$ in (69).

The CDF is shown in Figure 6. The tail is compatible with $1/\kappa$ decay, corresponding to the $1/{\kappa }^2$ 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}{ccc}\hfill & \mathrm{Estimate}\hfill & \mathrm{Standard}\mathrm{Error}\hfill \\ 1\hfill & 15.0125\hfill & 0.000236326\hfill \\ log\left(\kappa \right)\hfill & -0.950783\hfill & 0.0000443829\hfill \end{array}$$

(89)

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

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

An expert in the conventional theory of Turbulence may wonder why the solutions of the Loop equation are relevant to the statistical distribution of the vorticity 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 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 described by Wylde diagrams.

As we have seen, the Anzatz (13) satisfies the loop equation and initial condition at $C\to 0$. It has a finite area derivative, which makes it a Stokes-type functional.

Moreover, we have shown in [1,2] (and also here, in Section 4) how the velocity distribution for the random uniform vorticity in the fluid was reproduced by our Anzatz (13) a singular momentum loop $\overrightarrow{P}\left(\theta \right)$.

The solution for $\overrightarrow{P}\left(\theta \right)$ was random complex and had slowly decreasing Fourier coefficients, leading to a discontinuity $sign(\theta -{\theta }^{\prime })$ in a pair correlation function (35).

The exact solution for $\overrightarrow{P}\left(\theta \right)$ in decaying Turbulence which we have found in this paper, also belongs to the Stokes functionals and satisfies the initial boundary value at the shrunk loop.

Therefore, it represents a statistical distribution in some Navier-Stokes flow, a degenerate fixed point of the Hopf equation for velocity circulation, summing up all the Wylde diagrams plus nonperturbative effects. 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 SO\left(d\right)$ and an arbitrary element ${\overrightarrow{w}}_k\in {\mathbb{S}}^{d-3}$. This step also depends upon the previous position ${\overrightarrow{F}}_k$.

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{Q}}_k,{\overrightarrow{w}}_k$.

We simulated this random walk in three dimensions and studied its statistical properties. The fractal dimension of the random walk without the closure condition approximately equals $d_f\approx 1$. The distribution of lengths of steps in Euclidean six-dimensional space $\mathbf{Re}\overrightarrow{F},\mathbf{Im}\overrightarrow{F}$ has a nontrivial power-like tail $\mathit{PDF}\propto x^{-3.15}$.

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.

There are many things to do next with this conjectured solution to the decaying turbulence problem; the first is to look for unnoticed inconsistencies.

As a first test of this solution, let us compare it with a multitude of experimental data and those from DNS [33].

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 [34], there is a regime change at large Reynolds numbers; the numbers achievable in modern DNS are at such a transitional regime.

Besides, fitting is an unreliable 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 7).

Instead, one should compare a theory with a null hypothesis and estimate the probability of this null hypothesis to be correct.

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