Dual Theory of Decaying Turbulence. 2. Numerical Simulation and Continuum Limit

1. Introduction

In 1993, we suggested [1] nonperturbative approach to Navier-Stokes turbulence based on the loop equations previously introduced in the gauge theory [2,3].

The Wilson loop average for the turbulence

$$\begin{array}{c}\hfill \Psi [\gamma ,C]={〈exp\left(\frac{ı\gamma }{\nu }\oint d\theta {\overrightarrow{C}}^{\prime }\left(\theta \right)·\overrightarrow{v}\left(\overrightarrow{C}\left(\theta \right)\right)\right)〉}_{init}\end{array}$$

(1)

treated as a function of time and a functional of the periodic function $C:\overrightarrow{r}=\overrightarrow{C}\left(\theta \right);\theta \in (0,2\pi )$ (not necessarily a single closed loop), satisfies certain functional equation which was defined and studied in these papers.

The averaging ${\langle \rangle }_{init}$ corresponds to averaging over an ensemble of solutions with different initial data. One does not need to specify this ensemble to obtain the Hopf or loop equation.

We are not going into details of the loop technology but shall investigate the solution found in the recent paper [4].

2. The theory

This solution has the following form

$$\Psi [\gamma ,C]={〈exp\left(\frac{ı\gamma }{\nu }\oint d\theta {\overrightarrow{C}}^{\prime }\left(\theta \right)·\overrightarrow{P}\left(\theta \right)\right)〉}_F$$

(2)

The averaging ${\langle \rangle }_F$ corresponds to averaging over an ensemble of solutions for $\overrightarrow{F}\left(\theta \right)$ to be specified below.

This periodic vector function $\overrightarrow{P}\left(\theta \right)$ also depends of time as follows

$$\begin{array}{ccc}\hfill & & \overrightarrow{P}(\theta ,t)=\sqrt{\frac{\nu }{2(t+t_0)}}\frac{\overrightarrow{F}\left(\theta \right)}{\gamma };\hfill \\ & & \hfill (4)\end{array}$$

(3)

The time-independent function $\overrightarrow{F}\left(\theta \right)$ satisfies two singular equations

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

(5a)

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

(5b)

$$\begin{array}{ccc}& & \Delta \overrightarrow{F}=\overrightarrow{F}(\theta +0)-\overrightarrow{F}(\theta -0);\hfill \end{array}$$

(5c)

$$\begin{array}{ccc}& & \overrightarrow{F}=\frac{\overrightarrow{F}(\theta +0)+\overrightarrow{F}(\theta -0)}{2};\hfill \end{array}$$

(5d)

This equation describes a fixed point regime of the loop equations for decaying turbulence. No approximations have been made so far.

This $\overrightarrow{F}\left(\theta \right)$ is a fractal curve in complex 3D space, with discontinuity $\Delta \overrightarrow{F}$ at every angle $\theta$.

The solution in [4] builds such a curve as a limit of the polygon with $N\to \infty$ vertexes. The positions ${\overrightarrow{F}}_k=\overrightarrow{F}(2\pi k/N)$ of these vertices satisfy the recurrent equation

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

(6a)

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

(6b)

There could be many solutions of such equations in $d>2$ dimensions, as these are two complex equations relating d complex parameters in ${\overrightarrow{F}}_{k+1}$ to those in ${\overrightarrow{F}}_k$.

Thus, there are $d-2$ complex parameters undetermined in every step ${\overrightarrow{q}}_k={\overrightarrow{F}}_{k+1}-{\overrightarrow{F}}_k$. This degeneration makes this a Markov chain or a random walk where each step depends only on a current position plus some arbitrary (random) parameters.

2.1. The Euler ensemble and its continuum limit

In particular, the following solution, found in [4], has some acceptable physical properties. We call it the Euler ensemble.

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

(7)

$$\begin{array}{ccc}& & \widehat{\Omega }\in O\left(d\right);\hfill \end{array}$$

(8)

$$\begin{array}{ccc}& & \beta =\frac{2\pi p}{q};\hfill \end{array}$$

(9)

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

(10)

$$\begin{array}{ccc}& & {\alpha }_k=\beta \sum _0^{k-1}{\sigma }_l;\hfill \end{array}$$

(11)

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

(12)

$$\begin{array}{ccc}\hfill & & \sum _0^{N-1}{\sigma }_l=qr;\hfill \end{array}$$

(13)

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

(14)

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

(15)

Here $p,q,r,N$ are positive integers, $\widehat{\Omega }$ is a random rotation matrix in $O\left(d\right)$, and $\overrightarrow{w}\in {\mathbb{S}}^{d-3}$ is a unit vector in $d-2$ dimensional space ${\mathbb{R}}^{d-2}$.

This paper only considers the three-dimensional space $d=3$ where $\overrightarrow{w}=\pm 1$. This sign factor can be absorbed into a redefinition of ${\sigma }_l$.

This sequence with arbitrary signs ${\sigma }_k=\pm 1$ subject to periodicity constraint (13) solves recurrent equation (6) independently of $\gamma$.

This ensemble consists of random $p,q,r,N,{\sigma }_0\cdots {\sigma }_{N-1}$ subject to periodicity conditions; we called it an Euler ensemble.

The ergodic hypothesis of [4] states that each element of the Euler ensemble enters with equal weight. This includes every set of ${\sigma }_l$ satisfying periodic boundary condition (13), every pair of coprime numbers $0<p<q<N,\gcd (p,q)=1$ and every $q,r,N$ with $(N-qr)\left(\mathrm{mod}2\right)=0$.

The number theory counts these states using Euler totients.

A unique property of this solution is that the velocity circulation is real in this solution, as it should, despite the complex values of ${\overrightarrow{F}}_k$.

This cancellation of the imaginary part happens because of the closure (i.e., periodicity) of both loops $\overrightarrow{C}\left(\theta \right),\overrightarrow{F}\left(\theta \right)$

$$\begin{array}{ccc}& & \oint d\theta {\overrightarrow{C}}^{\prime }\left(\theta \right)·\overrightarrow{F}\left(\theta \right)=-\oint d\theta \overrightarrow{C}\left(\theta \right)·{\overrightarrow{F}}^{\prime }\left(\theta \right);;\hfill \end{array}$$

(16)

$$\begin{array}{ccc}& & \sum _k\Delta {\overrightarrow{C}}_k·{\overrightarrow{F}}_k=-\sum _k{\overrightarrow{C}}_k·\Delta {\overrightarrow{F}}_k;\hfill \end{array}$$

(17)

The discontinuity (i.e., the polygon side) is a real unit vector :

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

(18)

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

(19)

The sum over all different configurations of the Euler ensemble with exponential weight $exp\left(-\mu N\right)$ is called a grand canonical ensemble, using the language of statistical physics. The number theory counting using Euler totients [5] produces the following general formula [4]

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

(20)

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

(21)

At the critical point $\mu \to 0$, the even N dominate, with the following result

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

(22)

The enstrophy, defined as an expectation value of the square of local vorticity, reduces in this theory to the following time-dependent formula [4]

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

(23)

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

(24)

We have found the method to compute this sum of squares of cotangent of all fractions of $\pi$ with a given denominator analytically. This result was quoted in [4] and then generalized to arbitrary powers of cotangent.

This method is described in Appendix A

2.2. Small Euler ensemble at large N

At a fixed large value of q, the CDF of $x=\beta /\left(2\pi \right)$ is given by the ratio of the Euler totient functions, asymptotically equal to x at large q

$$\begin{array}{ccc}& & \mathbb{P}(\beta <2\pi x\mid q)=\frac{\phi \left(xq\right)}{\phi \left(q\right)}\to x;\hfill \end{array}$$

(25)

The CDF of q is proportional to the totient summatory function

$$\begin{array}{c}\hfill \Phi \left(n\right)=\sum _{k=1}^n\phi \left(k\right)\to \frac{6}{{\pi }^2}{n}^2\end{array}$$

(26)

Thus, the normalized CDF of $y=q/N$

$$\begin{array}{ccc}& & \mathbb{P}(q<yN)\to y^2;\hfill \end{array}$$

(27)

$$\begin{array}{ccc}& & f_y\left(y\right)\to 2y\theta \left(y\right)\theta (1-y);\hfill \end{array}$$

(28)

As we shall see, rather than $\beta$, we would need an asymptotic distribution of a scaling variable

$$\begin{array}{c}\hfill X(p,q)=\frac{1}{q^2}{cot}^2\left(\frac{\pi p}{q}\right)\to \frac{1}{{\pi }^2{p}^2}\end{array}$$

(29)

The $\beta N$ variable can be related to y and $X(p,q)$

$$\begin{array}{ccc}& & \beta N=2NarcCot\left(q\sqrt{X}\right)\to \frac{2}{y\sqrt{X}};\hfill \end{array}$$

(30)

This distribution for $X(p,q)$ at fixed $q\sim N\to \infty$ can be found analytically, using newly established relations for the cotangent sums (see Appendix in [4]). Asymptotically, at large q these relations read

$$\begin{array}{c}\hfill 〈X^n〉\equiv \underset{q\to \infty }{lim}\frac{1}{q}\sum _{\begin{array}{c}p=1\\ (p,q)\end{array}}^{q-1}X{(p,q)}^n={\delta }_{n,0}+\frac{2{\pi }^{-2n}\zeta \left(2n\right)}{(2n+1)\zeta (2n+1)}\end{array}$$

(31)

This relation can be transformed further as

$$\begin{array}{c}\hfill 〈X^n〉=\left\{\begin{array}{cc}1\hfill & \mathrm{if}n=0\hfill \\ \frac{{\displaystyle 2\sum _{k=1}^{\infty }\phi \left(k\right)k^{-(2n+1)}}}{(2n+1){\pi }^{2n}}\hfill & \mathrm{if}n>0\hfill \end{array}\right.\end{array}$$

(32)

The Mellin transform of these moments leads to the following singular distribution

$$\begin{array}{ccc}\hfill & & 〈X^n〉={\int }_0^{\infty }{W}^{\prime }\left(X\right)dX{X}^n;\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill & & W^{\prime }\left(X\right)=(1-\alpha )\delta (X-0^{+})+\pi \Phi \left(\lfloor \frac{1}{\pi \sqrt{X}}\rfloor \right);\hfill \end{array}$$

(34)

$$\begin{array}{ccc}& & \alpha ={\int }_0^{\infty }\pi \Phi \left(\lfloor \frac{1}{\pi \sqrt{X}}\rfloor \right)dX=\pi \sum _{k=1}^{\infty }\phi \left(k\right)\frac{2k+1}{k^2{(k+1)}^2}=0.387153\hfill \end{array}$$

(35)

where $\Phi \left(n\right)$ is the totient summatory function

The first term $(1-\alpha )\delta (X-0^{+})$ accumulates the contribution from the vicinity of the finite lower limit of X

$$\begin{array}{c}\hfill X_{min}=X\left(\lfloor q/2\rfloor -1,q\right)\to \frac{{\pi }^2}{q^4}=0^{+}\end{array}$$

(36)

This lower limit tends to zero, which leaves some part of the delta function out of the support of probability distribution. Normalizing the base moment $〈X^0〉=1$ yields the correction factor $(1-\alpha )$ in front of the delta function. Therefore, the $\alpha$ part of the delta function in the normalized distribution falls below the lower limit, with the remaining part $(1-\alpha )$ staying in the support of the distribution.

The upper limit of X

$$\begin{array}{c}\hfill X_{max}=X\left(q-1,q\right)\to \frac{1}{{\pi }^2}\end{array}$$

(37)

Our distribution (33) is consistent with this upper limit, as the argument $\lfloor \frac{1}{\pi \sqrt{X}}\rfloor$ becomes zero at $X{\pi }^2>1$. It is plotted in Figure 1, except for the delta function part at infinitesimal X.

Once we are zooming into the tails of the $p,q$ distribution, we also must recall that

$$\begin{array}{ccc}& & \mathbb{P}(q<yN)=\Phi \left(\lfloor Ny\rfloor \right);\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill & & f_y\left(y\right)=\sum _{q=2}^{\infty }\delta \left(y-\frac{q}{N}\right)\phi \left(q\right)\hfill \end{array}$$

(39)

The variable s also has a nontrivial distribution in the small Euler ensemble with even $N\to \infty$, with a dominant contribution at $s=0$ as conjectured in [4] and proven in [6]

$$\begin{array}{c}\hfill f_s\left(s\right)=\delta \left(s\right)+O\left(N^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\end{array}$$

(40)

3. The big Euler ensemble as a Markov process

3.1. Numerical Simulations

We took $N=200,000,000$ and generated $T=167,772,160$ random data samples for $\beta ,{\overrightarrow{\omega }}_n·{\overrightarrow{\omega }}_m,\left|{\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right|$ with randomly chosen $0\le n<m<N$.

We devised a fast Python/CPP algorithm for this simulation, described in Appendix B. We ran it on an NYUAD cluster Jubail, which took about a week altogether, for three GPU runs and three CPU runs, which we all pooled into one dataset and binned by $2^{20}$ intervals by the ranking of $log{cot}^2(\beta /2)$. Two datasets were created, one for each sign of ${\overrightarrow{\omega }}_n·{\overrightarrow{\omega }}_m$.

We used a recursive rank-binning algorithm based on the standard library function $\mathit{nth}-\mathit{element}$, which finds the median of unordered data by partial sorting in place in $O\left(N\right)$ time. We recursively applied this function in parallel to each of the two halves of the array of data. The CPU time of this recursive parallel code is linear, as the geometric series $N+N/2+N/4+\cdots$ adds up to $2N$.

Each interval contained only $n=80\pm 1$ data points $\left\{n,\overline{c},\Delta c,\overline{s},\Delta s,\overline{w},\Delta w\right\}$, where

$$\begin{array}{ccc}& & c=2log\left|cot(\beta /2)\right|;\hfill \end{array}$$

(41)

$$\begin{array}{ccc}& & s=log\left|{\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right|;\hfill \end{array}$$

(42)

$$\begin{array}{ccc}& & w=log|{\overrightarrow{\omega }}_n·{\overrightarrow{\omega }}_m|;\hfill \end{array}$$

(43)

and $\overline{x},\Delta x$ denotes mean and standard deviation. Later, we used these data with Gaussian errors for the statistical analysis. We used logarithms because the values varied by several orders of magnitude, preventing statistical analysis.

The code was optimized to take $O\left(N^0\right)$ RAM, so the CPU and GPU resources of all 200 nodes were used to collect statistics. This optimization allows us to simulate astronomical ensembles of random fractions without hitting memory limits.

The first thing we measured was the distribution of $x=\frac{\beta }{2\pi }=\frac{p}{q}$, which we know in advance to be uniform except for small regions of $x\sim 1/N,1-x\sim 1/N$, where quantization plays the role.

The measured distribution supports this hypothesis (see Figure 2).

The CDF is linear, supporting the uniform distribution of $\frac{p}{q}$ except for the endpoints $(0,1)$.

However, the critical phenomena in our system are related to these two endpoints where the scaling factor ${cot}^2(\beta /2)$ grows to $\sim q^2$. The uniform distribution of $\beta$ corresponds to the variable $X=q^{-2}{cot}^2(\beta /2)$ distributed at large q as

$$\begin{array}{c}\hfill f_X(X\mid q)\propto q^{-1}{X}^{-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\end{array}$$

(44)

The exact asymptotic distribution for the same variable has the form (34), found in the previous section.

3.2. Scaling variables in continuum limit

This quantum distribution is too complicated to reproduce in numerical simulation, but we measured distributions of other variables that behave smoother.

We studied the distribution of two variables involved in the correlation function (see next section)

$$\begin{array}{ccc}& & DS=q^{-1}\left|{\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right|;\hfill \end{array}$$

(45)

$$\begin{array}{ccc}& & WW=q^{-2}{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n;\hfill \end{array}$$

(46)

We multiplied both variables by powers of q to make them finite in the statistical limit $q\sim N\to \infty$.

The distributions of these variables are shown in Figure 3 and Figure 4.

The variable $DS$ measures the mean Euclidean distance between two random points on a curve in the embedding 3D complex space.

Let us study the moments of this distribution. Integrating by part in the standard definition, we have

$$\begin{array}{ccc}& & 〈{\left(DS\right)}^n〉=\int W^{\prime }\left(DS\right){\left(DS\right)}^{n}d\left(DS\right)=n\int (1-W\left(DS\right)){\left(DS\right)}^{n-1}d\left(DS\right)=\hfill \end{array}$$

$$\begin{array}{ccc}& & \int exp\left(n\xi +T\left(\xi \right)+logn\right)d\xi ;\hfill \end{array}$$

(47)

$$\begin{array}{ccc}& & T\left(\xi \right)=log\left((1-W\left(exp\left(\xi \right)\right)\right);\hfill \end{array}$$

(48)

Our data was fit perfectly by 20 Chebyshev polynomials for this $T\left(\xi \right)$ (see [7] and Figure 5).

The definition of effective index can be given as

$$\begin{array}{c}\hfill \eta \left(n\right)=\frac{log〈{\left(DS\right)}^n〉}{log〈DS〉};\end{array}$$

(49)

The saddle point computation of this integral gives the following parametric equations

$$\begin{array}{ccc}& & n=-T^{\prime }\left(\xi \right);\hfill \end{array}$$

(50)

$$\begin{array}{ccc}& & log〈{\left(DS\right)}^n〉\to L_n\left(\xi \right)={\left(n\xi +T\left(\xi \right)+logn\right)}_{n=-T^{\prime }\left(\xi \right)};\hfill \end{array}$$

(51)

$$\begin{array}{ccc}& & L_1={\left({\xi }_1+T\left({\xi }_1\right)\right)}_{1=-T^{\prime }\left({\xi }_1\right)};\hfill \end{array}$$

(52)

$$\begin{array}{ccc}& & \eta (-T^{\prime }\left(\xi \right))=\frac{L_n\left(\xi \right)}{L_1}\hfill \end{array}$$

(53)

This function was computed in [7], and here is the plot in Figure 6. It is numerically close to a linear function , though it is analytically nonlinear.

We also study some quantum fractal phenomena in the next section. The observable variables are distributed with quantum jumps, not fitting the multifractal framework.

There is another variable that depends on two spins at two different points ${\sigma }_n,{\sigma }_m,n<m$

$$\begin{array}{c}\hfill Y(p,q,\sigma )=-{\sigma }_n{\sigma }_{m}X(p,q)\end{array}$$

(54)

The even moments $〈Y^{2l}〉$ of the distribution of Y are the same as those for X, but the odd moments are different:

$$\begin{array}{c}\hfill 〈Y^{2l+1}〉=-{〈X^{2l+1}{〈{\sigma }_n{\sigma }_m〉}_{\sigma }〉}_{\mathcal{E}}\end{array}$$

(55)

The analytic computation of the constrained average over sigmas by Mathematica^® yields

$$\begin{array}{c}\hfill {〈{\sigma }_n{\sigma }_m〉}_{\sigma }=-\frac{N-q^2{r}^2}{(N-1)N}\end{array}$$

(56)

We find the following result for the distribution of Y on the whole real axis

$$\begin{array}{ccc}& & f_Y(Y\mid N,q,r)=\left(\frac{1}{2}+signY\frac{N-q^2{r}^2}{2(N-1)N}\right)f_X\left(\right|Y\left|\right)\hfill \end{array}$$

(57)

Asymptotically, at $s=qr/N\sim 1$ we have

$$\begin{array}{c}\hfill f_Y(Y\mid s)\to \frac{1-sign\left(Y\right)s^2}{2}\left((1-\alpha )\delta \left(\right|Y|-0^{+})+\pi \Phi \left(\lfloor \frac{1}{\pi \sqrt{\left|Y\right|}}\rfloor \right)\right);\end{array}$$

(58)

Due to kinematical restrictions $\left|y\right|<1$, our probably stays positive. At $r=0$, we have a special case

$$\begin{array}{c}\hfill f_Y(Y\mid s=0)\to \frac{N+sign\left(Y\right)}{2N}\left((1-\alpha )\delta \left(\right|Y|-0^{+})+\pi \Phi \left(\lfloor \frac{1}{\pi \sqrt{\left|Y\right|}}\rfloor \right)\right);\end{array}$$

(59)

4. Vorticity correlation

4.1. Exact relation with conditional probability density

The simplest observable quantities we can extract from the loop functional are the vorticity correlation functions [8], corresponding to the loop C backtracking between two points in space ${\overrightarrow{r}}_1=0,{\overrightarrow{r}}_2=\overrightarrow{r}$, (see [4] for details and the justification). The vorticity operators are inserted at these two points.

The correlation function reduces to the following average over the big Euler ensemble $\mathbb{E}$ of our random curves in complex space [4]

$$\begin{array}{ccc}\hfill & & 〈\overrightarrow{\omega }\left(\overrightarrow{0}\right)·\overrightarrow{\omega }\left(\overrightarrow{r}\right)〉=\hfill \end{array}$$

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

(60)

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

(61)

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

(62)

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

(63)

$$\begin{array}{ccc}& & {〈X[\sigma ...]〉}_{p,q,r}\equiv \frac{{\sum }_{\mathbb{E}}X[\sigma ...]\delta [qr-\sum \sigma ]}{{\sum }_{\mathbb{E}}\delta [qr-\sum \sigma ]};\hfill \end{array}$$

(64)

Integrating the global rotation matrix $O\left(3\right)$ is part of the ensemble averaging.

These formulas simplify in the Fourier space, corresponding to the correlation of $\overrightarrow{\omega }\left(\overrightarrow{k}\right)·\overrightarrow{\omega }(-\overrightarrow{k})$.

$$\begin{array}{ccc}& & 〈\overrightarrow{\omega }\left(\overrightarrow{k}\right)·\overrightarrow{\omega }(-\overrightarrow{k})〉\propto {\nu }^{\frac{3}{2}}{(t+t_0)}^{-\frac{1}{2}}\sum _{0\le n<m<N}\hfill \\ & & {〈{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n{\int }_{O\left(3\right)}d\Omega {\delta }^3\left(\Omega ·\left({\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right)-\overrightarrow{k}\sqrt{\nu (t+t_0)}\right)〉}_{\mathbb{E}}\hfill \end{array}$$

(65)

This three-dimensional delta function disappears if we introduce the above distribution $f_Y\left(y\right)$ for $Y=q^{-2}{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_n$ and conditional distribution $f_Z(Z\mid y)$ for the second variable $Z=\frac{\left|{\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right|}{q}$ given the first one y. We rely on $r=0$ and even N to dominate the partition function, as it was conjectured in [4] and proven in [6]. We also assume the distribution $f_Z(z\mid y)$ for the second scaling variable Z, conditional by $\left|Y\right|=y$, to be independent of the sign of Y. The heuristic argument is that the sign of Y depends only on ${\sigma }_n{\sigma }_m$, only two out of N variables ${\sigma }_l$, which becomes statistically insignificant at large N.Our large-scale numerical simulation confirms this heuristic argument for $N=2\ast {10}^8$.

We also use the conditional probability, which is an integral of the conditional distribution

$$\begin{array}{ccc}& & f_Z(z\mid Y>y)={\int }_y^{0}d{y}^{\prime }{f}_Z(z\mid y^{\prime });;\hfill \end{array}$$

(66)

$$\begin{array}{ccc}& & f_Z(z\mid y)=-{\partial }_y{f}_Z(z\mid Y>y)\hfill \end{array}$$

(67)

In this case, after switching to $f_Z(z\mid Y>y)$ and integrating by parts we find

$$\begin{array}{ccc}& & 〈\overrightarrow{\omega }\left(\overrightarrow{k}\right)·\overrightarrow{\omega }(-\overrightarrow{k})〉\propto {\nu }^{\frac{3}{2}}{(t+t_0)}^{-\frac{1}{2}}\int dW\left(q\right)q^4\hfill \\ & & \int f_Y\left(y\right)ydy\int d\Omega \int dZ{f}_Z(Z\mid y){\delta }^3\left(\Omega ·\left\{0,0,Zq\right\}-\overrightarrow{k}\sqrt{\nu (t+t_0)}\right)=\hfill \\ & & \frac{\sqrt{\nu }}{{\overrightarrow{k}}^2{(t+t_0)}^{\frac{3}{2}}}\int dq\Phi \left(q\right)q\int dy{f}_Y\left(y\right)y{f}_Z\left(\frac{|\overrightarrow{k}|\sqrt{\nu (t+t_0)}}{q}\mid y\right)=\hfill \\ & & \frac{\sqrt{\nu }}{{\overrightarrow{k}}^2{(t+t_0)}^{\frac{3}{2}}}\int dq\Phi \left(q\right)q\int dy{f}_Z\left(\frac{|\overrightarrow{k}|\sqrt{\nu (t+t_0)}}{q}\mid Y>y\right){\partial }_y\left(f_Y(y\mid 0)y\right);\hfill \end{array}$$

(68)

The distribution $f_y(y\mid 0)$ was found above, in (59) by the number theory methods. We are using the derivative

$$\begin{array}{ccc}& & {\partial }_y\left(f_Y(y\mid 0)y\right)=\pi \frac{\partial \left(\left|y\right|\Phi \left(\lfloor \frac{1}{\pi \sqrt{\left|y\right|}}\rfloor \right)\right)}{\partial y}=\hfill \\ & & \pi \Phi \left(\lfloor \frac{1}{\pi \sqrt{\left|y\right|}}\rfloor \right)+\frac{2}{{\pi }^3}\sum _{n=1}^{\infty }\frac{\phi \left(n\right)}{n^5}\delta \left(y-\frac{1}{{\pi }^2{n}^2}\right)\hfill \end{array}$$

(69)

The resulting formula for the correlation function

$$\begin{array}{ccc}\hfill & & 〈\overrightarrow{\omega }\left(\overrightarrow{k}\right)·\overrightarrow{\omega }(-\overrightarrow{k})〉\propto \frac{\sqrt{\nu }}{{\overrightarrow{k}}^2{(t+t_0)}^{\frac{3}{2}}}\int dqq\Phi \left(q\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & \left(\frac{2}{{\pi }^3}\sum _{n=1}^{\infty }\frac{\phi \left(n\right)}{n^5}{f}_Z\left(z\mid Y<\frac{1}{{\pi }^2{n}^2}\right)+\pi \sum _{n=1}^{\infty }\Phi (n+1){\int }_{\frac{1}{{\pi }^2{(n+1)}^2}}^{\frac{1}{{\pi }^2{n}^2}}dy{f}_Z(z\mid Y<y)\right);\hfill \end{array}$$

(70)

$$\begin{array}{ccc}& & z=\frac{|\overrightarrow{k}|\sqrt{\nu (t+t_0)}}{q};\hfill \end{array}$$

(71)

4.2. Extracting conditional distribution from numerical data

We measured conditional probability in our data set of $167,772,160$ samples for $N=200,000,000$ (see Figure 7).

The conditional probability $\mathbb{P}\left(DS\mid WW<y\right)$ was defined as a rate of events with given $DS$ and $WW<y$. It is shown as a 3D plot on the log-log-log scale in Figure 7.

In a classical fractal system, one would expect to observe the two-dimensional surface $\mathbb{P}\left(logZ\mid logWW<Y\right)$ in a three-dimensional space $Y,Z,\mathbb{P}$.

Instead, we found a thin, smooth line in three dimensions, which can be fitted as a parametric equation

$$\begin{array}{ccc}\hfill & & logY=a+blogZ;\hfill \end{array}$$

(72)

$$\begin{array}{ccc}& & logZ=f\left(L\right);\hfill \end{array}$$

(73)

$$\begin{array}{ccc}& & log\mathbb{P}\left(Z\mid WW>Y\right)=L;\hfill \end{array}$$

(74)

$$\begin{array}{ccc}& & f\left(L\right)\approx f_0+f_{1}L+f_2{L}^2+f_3{L}^3;\hfill \end{array}$$

(75)

Projections of this line on three planes reveal its simplicity.

The relation between $logZ,log\mathbb{P}$ is shown in Figure 8.

The fit parameter table for $f\left(x\right)$ is

$$\begin{array}{c}\hfill \begin{array}{ccccc} \hfill & \mathrm{Estimate}\hfill & \mathrm{Standard}\mathrm{Error}\hfill & \mathrm{t}-\mathrm{Statistic}\hfill & \mathrm{P}-\mathrm{Value}\hfill \\ f_0\hfill & -10.6123\hfill & 0.0030473\hfill & -3482.54\hfill & 0.\hfill \\ f_1\hfill & 5.43603\hfill & 0.00542859\hfill & 1001.37\hfill & 0.\hfill \\ f_2\hfill & 0.849258\hfill & 0.00224021\hfill & 379.097\hfill & 0.\hfill \\ f_3\hfill & 0.0500626\hfill & 0.000234038\hfill & 213.908\hfill & 0.\hfill \end{array}\end{array}$$

(76)

If we directly plot $Y\left(Z\right)$, we get a straight line corresponding to a scaling law (see Figure 9). The fitting table for this linear relation $logY=a+blogZ$ in the log-log scale is

$$\begin{array}{c}\hfill \begin{array}{ccccc} \hfill & \mathrm{Estimate}\hfill & \mathrm{Standard}\mathrm{Error}\hfill & \mathrm{t}-\mathrm{Statistic}\hfill & \mathrm{P}-\mathrm{Value}\hfill \\ a\hfill & -10.2124\hfill & 0.0039544\hfill & -2582.54\hfill & 0.\hfill \\ b\hfill & 1.03547\hfill & 0.00025789\hfill & 4015.16\hfill & 0.\hfill \end{array}\end{array}$$

(77)

So, a scaling law is hidden inside the discontinuous distributions induced by Euler totients.

$$\begin{array}{ccc}\hfill & & Y=A{Z}^b;\hfill \end{array}$$

(78)

$$\begin{array}{ccc}& & A=exp\left(a\right)=0.0000367128;\hfill \end{array}$$

(79)

The probability density $f_Z(z\mid Y>y)$, which we need in correlation function, is related to this number of events $\mathbb{P}\left(Z\mid WW>y\right)$ as follows

$$\begin{array}{c}\hfill f_Z(z\mid Y>y)=\frac{d\mathbb{P}\left(Z\mid WW>y\right)}{dZ}=\frac{{\partial }_L\mathbb{P}\left(Z\mid WW>y\right)}{{\partial }_{L}Z}=\frac{exp\left(L-f\left(L\right)\right)}{|f^{\prime }\left(L\right)|}\end{array}$$

(80)

4.3. Final results for the energy spectrum

We use parametric representation (72) and transform the integral, inverting relation between L and $logK/q$ (see details in [7]).

We arrive at the following correlation

$$\begin{array}{ccc}\hfill & & 〈\overrightarrow{\omega }·\overrightarrow{\omega }\left(k\right)〉=\frac{{\nu }^2}{\sqrt{\nu (t+t_0)}}{G}_N\left(K\right)H\left(K\right)/K;\hfill \end{array}$$

(81)

$$\begin{array}{ccc}& & K=|\overrightarrow{k}|\sqrt{\nu (t+t_0)};\hfill \end{array}$$

(82)

$$\begin{array}{ccc}& & G_N\left(K\right)={\int }_K^N\frac{q^2\Phi \left(q\right)exp\left(L(q/K)\right)}{|f^{\prime }\left(L(q/K)\right)|}dq;\hfill \end{array}$$

(83)

$$\begin{array}{ccc}& & L(q/k)=f^{-1}(logK/q);\hfill \end{array}$$

(84)

$$\begin{array}{ccc}& & H\left(K\right)=\frac{2}{{\pi }^3}\sum _{n=1}^{\infty }\frac{\varphi \left(n\right)}{n^5}\delta \left(A{K}^b-\frac{1}{n^2{\pi }^2}\right)+\pi \Phi \left(⌊\frac{1}{\sqrt{A{K}^b}\pi }⌋\right);\hfill \end{array}$$

(85)

Here $w=f^{-1}\left(\xi \right)$ is the real root of equation $f\left(w\right)=\xi$. In our case $\xi =log(K/q)\to -\infty$, and $f^{-1}\left(\xi \right)\propto {\xi }^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\to -{(logq)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$ These estimates show that function $G_N\left(K\right)$ in continuum limit $N\to \infty$ tends to a constant, independent of K, but growing with N

$$\begin{array}{c}\hfill G_N\left(K\right)\propto F\left(N\right)\end{array}$$

(86)

This work neglects such normalization factors, not affecting the energy spectrum. We restore normalization, comparing our dissipation with the exact result from [4].

The summation ${\sum }_n$ in H(K) at any given K selects only one term with a quantized spectrum in the energy

$$\begin{array}{c}\hfill E(k,t)\propto k^2〈\overrightarrow{v}·\overrightarrow{v}\left(k\right)〉\propto 〈\overrightarrow{\omega }·\overrightarrow{\omega }\left(k\right)〉\propto \frac{{\nu }^2}{\sqrt{\nu (t+t_0)}}H\left(K\right)/K;\end{array}$$

(87)

The spectrum has delta function peaks at

$$\begin{array}{ccc}\hfill & & K_n=\frac{1}{{\left(A{n}^2{\pi }^2\right)}^{\frac{1}{b}}},n=1,2,\cdots ;\hfill \end{array}$$

(88)

$$\begin{array}{ccc}& & |\overrightarrow{k}{|}_n=\frac{K_n}{\sqrt{\nu (t+t_0)}};\hfill \end{array}$$

(89)

The whole spectrum support lies at $K<K_1$; otherwise, both terms in $H\left(K\right)$ vanish. The critical region is $K\to 0$, which means the smaller and smaller wavevectors at larger and larger times. In the limit $K\to 0$

$$\begin{array}{ccc}& & H\left(K\right)/K\propto \frac{1}{K^{b+1}};\hfill \end{array}$$

(90)

$$\begin{array}{ccc}& & E(k,t)\propto \frac{F\left(N\right){\nu }^2}{|\overrightarrow{k}{|}^{b+1}{\left(\nu (t+t_0)\right)}^{1+b/2}};\hfill \end{array}$$

(91)

This asymptotic power spectrum is only part of the story. There is a maximal time for any wavevector

$$\begin{array}{c}\hfill t_{max}\left(k\right)=-t_0+\frac{K_1^2}{\nu {\overrightarrow{k}}^2}\end{array}$$

(92)

The decay stops at this wavelength after this moment, as $H(K>K_1)=0$.

The quantum effects lead to peaks and jumps at finite $K=k\sqrt{\nu (t+t_0)}$. Only the average over a large range of the spectrum approaches the scaling law.

The integral over $\overrightarrow{k}$ leads to a simple power law for the energy dissipation rate

$$\begin{array}{ccc}& & \mathcal{E}\propto \frac{4\pi }{{\left(2\pi \right)}^3}{\int }_0^{\infty }dk{k}^2〈\overrightarrow{\omega }·\overrightarrow{\omega }\left(k\right)〉\propto \frac{{\nu }^2}{2{\pi }^2\nu (t+t_0)}{\int }_0^{\infty }dkkH\left(k\sqrt{\nu (t+t_0)}\right)=\hfill \\ & & \frac{{\int }_0^{K_1}KdKH\left(K\right)}{2{\pi }^2{(t+t_0)}^2}\hfill \end{array}$$

(93)

This law agrees with the exact computation of the same quantity made in [4]. Comparing coefficients in front of ${(t+t_0)}^{-2}$, one can restore the missing normalization factor in our energy spectrum.

$$\begin{array}{c}\hfill E(k,t)=〈\overrightarrow{\omega }·\overrightarrow{\omega }\left(k\right)〉=\frac{35{\pi }^4}{1728\zeta \left(3\right){\mu }^2}\frac{\nu }{|\overrightarrow{k}|(t+t_0)}\frac{H\left(|\overrightarrow{k}|\sqrt{\nu (t+t_0)}\right)}{{\int }_0^{K_1}xdxH\left(x\right)}\end{array}$$

(94)

where $\mu \sim \sqrt{\nu }\to 0$ is the chemical potential introduced in [4].

Finally, let us plot the universal function $H\left(K\right)/K$, using a width $\sigma ={10}^{-3}$ for normal distribution instead of a pure delta function with $\sigma =0$ (see Figure 10). The narrow peaks at the $K=K_n$ are not visible because of the very small heights of these peaks. The curve has tilted steps with diminishing width as $K\to 0$. The local slope is $K^{-1}$, but the average slope over wide range is $K^{-1-b}\approx K^{-2.03}$.

5. Conclusions

There are several new results reported in this paper.

• Continuum limit of distribution of scaling variables (34),(39) in the small Euler ensemble.
• Relation between the vorticity correlation and conditional probability distribution $f_Z(z\mid Y)$.
• Fast algorithm with $O\left(1\right)$ memory to simulate the big Euler ensemble as a Markov chain.
• The scaling law (78) for the conditional probability leads to the discrete energy spectrum (88) on top of a continuous background in (81).
• The decay at a given time goes only at small enough wavelengths. At a fixed wavelength, the decay stops after some critical time, inversely proportional to the wavelength.

Nothing like these quantum effects was observed in numerical simulations [9]. There could be several explanations (see discussion in the paper [10]). It is very difficult to numerically distinguish the devil’s staircase like that in Figure 10 from the good old power law.

Discrete energy spectrum in classical turbulence sounds like heresy, though it does not contradict any known fluid dynamics principles. Let us stress that we did not assume a discrete spectrum: it came from large-scale simulations on a supercomputer with very small statistical errors and some exact distributions derived from the number theory.

This discrete spectrum stems from the exact analogy between decaying turbulence and the quantum mechanics in loop space [1,4,8].

We would like to hear the response from leading experts in fluid dynamics before finally publishing this paper in a journal.