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 \mathsf{\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

$$\mathsf{\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 \end{array}$$

(3)

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

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

(4a)

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

(4b)

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

(4c)

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

(4d)

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

(5a)

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

(5b)

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{\mathsf{\Omega }}·\left\{cos\left({\alpha }_k\right),sin\left({\alpha }_k\right)\overrightarrow{w},icos\left(\frac{\beta }{2}\right)\right\};\hfill \end{array}$$

(6)

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

(7)

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

(8)

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

(9)

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

(10)

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

(11)

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

(12)

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

(13)

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

(14)

Here $p,q,r,N$ are positive integers, $\widehat{\mathsf{\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 (12) solves recurrent equation (5) 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 (12), 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}$$

(15)

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

(16)

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

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

(17)

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

(18)

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

(19)

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

(20)

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

(21)

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

(22)

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

(23)

We have found the method to compute this sum of squares of cotangent of all fractions of $\pi$ with a given denominator analytically. This method is described in Appendix of [4].

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

(24)

The CDF of q is proportional to the totient summatory function

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

(25)

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

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

(26)

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

(27)

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

(28)

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

(29)

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

(30)

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

(31)

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

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

(32)

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

(33)

$$\begin{array}{ccc}& & \alpha ={\int }_0^{\infty }\pi \mathsf{\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}$$

(34)

where $\mathsf{\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-1,q\rfloor \right)\to \frac{{\pi }^2}{q^4}=0^{+}\end{array}$$

(35)

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

(36)

Our distribution (32) 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)=\mathsf{\Phi }\left(\lfloor Ny\rfloor \right);\hfill \end{array}$$

(37)

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

(38)

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

(39)

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 A. 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},\mathsf{\Delta }c,\overline{s},\mathsf{\Delta }s,\overline{w},\mathsf{\Delta }w\right\}$, where

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

(40)

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

(41)

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

(42)

and $\overline{x},\mathsf{\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}$$

(43)

The exact asymptotic distribution for the same variable has the form (33), 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}$$

(44)

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

(45)

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

(46)

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

(47)

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

(48)

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

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

(49)

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

(50)

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

(51)

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

(52)

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

(53)

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

(54)

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

(55)

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

(56)

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 \mathsf{\Phi }\left(\lfloor \frac{1}{\pi \sqrt{\left|Y\right|}}\rfloor \right)\right);\end{array}$$

(57)

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 \mathsf{\Phi }\left(\lfloor \frac{1}{\pi \sqrt{\left|Y\right|}}\rfloor \right)\right);\end{array}$$

(58)

4. Vorticity correlation

4.1. Exact relation to the 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\mathsf{\Omega }\sum _{0\le n<m<N}{〈{\overrightarrow{\omega }}_m·{\overrightarrow{\omega }}_{n}exp\left(ı\overrightarrow{\rho }·\mathsf{\Omega }·\left({\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right)\right)〉}_{\mathbb{E}};\hfill \end{array}$$

(59)

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

(60)

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

(61)

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

(62)

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

(63)

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\mathsf{\Omega }{\delta }^3\left(\mathsf{\Omega }·\left({\overrightarrow{S}}_{n,m}-{\overrightarrow{S}}_{m,n}\right)-\overrightarrow{k}\sqrt{\nu (t+t_0)}\right)〉}_{\mathbb{E}}\hfill \end{array}$$

(64)

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

(65)

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

(66)

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\mathsf{\Omega }\int dZ{f}_Z(Z\mid y){\delta }^3\left(\mathsf{\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\mathsf{\Phi }(\lfloor q\rfloor )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\mathsf{\Phi }(\lfloor q\rfloor )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}$$

(67)

The distribution $f_y(y\mid 0)$ was found above, in (58) 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|\mathsf{\Phi }\left(\lfloor \frac{1}{\pi \sqrt{\left|y\right|}}\rfloor \right)\right)}{\partial y}=\hfill \\ & & \pi \mathsf{\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}$$

(68)

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\mathsf{\Phi }(\lfloor q\rfloor )\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 }\mathsf{\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}$$

(69)

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

(70)

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

(71)

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

(72)

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

(73)

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

(74)

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

(75)

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

(76)

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

(77)

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

(78)

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

(79)

4.3. Final results for the energy spectrum

We use parametric representation (71) 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}$$

(80)

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

(81)

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

(82)

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

(83)

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

(84)

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

(85)

This work neglects such normalization factors, not affecting the energy spectrum.

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

(86)

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

(87)

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

(88)

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

(89)

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

(90)

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

(91)

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

(92)

The normalization integral ${\int }_0^{K_1}KdKH\left(K\right)$ is a universal number

$$\begin{array}{c}\hfill {\int }_0^{K_1}KdKH\left(K\right)=\frac{1}{2}{\pi }^{-\frac{b+4}{b}}{A}^{-2/b}\left(\frac{{\pi }^2\zeta \left(\frac{4}{b}-1\right)}{\zeta \left(\frac{4}{b}\right)}+\frac{4\zeta \left(2+\frac{4}{b}\right)}{b\zeta \left(3+\frac{4}{b}\right)}\right)=1.05869\times {10}^7\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}{ccc}& & E(k,t)=〈\overrightarrow{\omega }·\overrightarrow{\omega }\left(k\right)〉=\frac{ϵ}{|\overrightarrow{k}|(t+t_0)}H\left(|\overrightarrow{k}|\sqrt{\nu (t+t_0)}\right);\hfill \end{array}$$

(94)

$$\begin{array}{ccc}& & ϵ=\frac{35{\pi }^4\nu }{1728\zeta \left(3\right){\mu }^2{\int }_0^{K_1}xdxH\left(x\right)}\hfill \end{array}$$

(95)

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

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 originates from the exact analogy between decaying turbulence and the quantum mechanics in loop space [1,4,8].

5. Comparison with real and numerical experiments

These quantum effects were never observed in numerical simulations [9]. On the other hand, the power laws observed had exponents distributed in the range $n=1.2$ to $n=1.6$ in the decay law $E\sim t^{-n}$.

Let us analyze these data and compare them with our theory. First, our theory corresponds to the idealized case of infinite initial energy taking infinite time to dissipate.

Integrating our energy spectrum yields infrared infinity

$$\begin{array}{c}\hfill E\left(t\right)\propto {\int }_0^{\infty }\frac{dk}{kt}H\left(k\sqrt{\nu (t+t_0)}\right)\sim {\int }_0^{\infty }dk/k^{1+b}=\infty \end{array}$$

(96)

In the real world, with a finite spacial box, the wavevectors are cut off at some $k_0$. The energy spectrum below this cutoff rapidly decreases as $k^2$ or $k^4$.

The first case (Saffman turbulence) corresponds to the finite linear momentum of the fluid or a conserved integral

$$\begin{array}{ccc}& & L=\int d^3\overrightarrow{r}〈\overrightarrow{v}\left(0\right)·\overrightarrow{v}\left(\overrightarrow{r}\right)〉;\hfill \end{array}$$

(97)

$$\begin{array}{ccc}& & E\left(k\right)\to \frac{L{k}^2}{4{\pi }^2}\hfill \end{array}$$

(98)

Our theory assumes the fluid at rest so that $L=0$.

In that case (Batchelor turbulence), the small k behavior of the energy spectrum is determined by another invariant I, corresponding to conserved angular momentum [10]

$$\begin{array}{ccc}& & \overrightarrow{H}=\int \overrightarrow{r}\times \overrightarrow{v}\left(\overrightarrow{r}\right)d^{3}r;\hfill \end{array}$$

(99)

$$\begin{array}{ccc}& & I={\overrightarrow{H}}^2/V;\hfill \end{array}$$

(100)

$$\begin{array}{ccc}& & E\left(k\right)\to \frac{I{k}^4}{24{\pi }^2}\hfill \end{array}$$

(101)

So, we must cut our spectrum at some wavelength $k_0$ to match the Batchelor energy pumping at smaller k. This matching takes place at $t=0$.

Integrating our spectrum with this cutoff, we find

$$\begin{array}{ccc}& & E\left(t\right)=\frac{ϵ}{t+t_0}{\int }_{k_0}^{\infty }\frac{dk}{k}H\left(k\sqrt{\nu (t+t_0)}\right)=\frac{ϵ}{t+t_0}\underset{k_0\sqrt{\nu (t+t_0)}}{\overset{K_1}{\int }}\frac{H\left(K\right)dK}{K};\hfill \end{array}$$

(102)

$$\begin{array}{ccc}& & \frac{I{k}_0^4}{24{\pi }^2}=E\left(0\right)\approx \frac{ϵ}{t_0}\underset{k_0\sqrt{\nu t_0}}{\overset{K_1}{\int }}\frac{H\left(K\right)dK}{K}\hfill \end{array}$$

(103)

The last equation fixes the unknown parameter $t_0$ in our solution. The integral of the first term on $H\left(K\right)$ converges at small k, and it is a universal number

$$\begin{array}{ccc}& & {\int }_0^{K_1}\frac{dK}{K}\frac{2}{{\pi }^3}\sum _{n=1}^{\infty }\frac{\phi \left(n\right)}{n^5}\delta \left(A{K}^b-\frac{1}{n^2{\pi }^2}\right)=\frac{2}{\pi b}\sum _1^{\infty }\frac{\phi \left(n\right)}{n^3}=\frac{\pi }{3b\zeta \left(3\right)}\hfill \end{array}$$

(104)

The second integral reduces to a finite sum so that

$$\begin{array}{ccc}& & E\left(t\right)=\frac{ϵ}{t+t_0}{H}_1\left(k_0\sqrt{\nu (t+t_0)}\right);\hfill \end{array}$$

(105)

$$\begin{array}{ccc}& & H_1\left(K\right)=\frac{\pi }{3b\zeta \left(3\right)}+\frac{2\pi }{b}\sum _{n=1}^{\lfloor \frac{1}{\sqrt{A{K}^b}\pi }\rfloor }\mathsf{\Phi }\left(n\right)log\left(\frac{n+1}{n}\right);\hfill \end{array}$$

(106)

The plot of this function is shown in Figure 11, and in a higrer resolution in Figure 12.

Now we are in a position to derive our prediction for effective time index $n\left(t\right)$ in [9]

$$\begin{array}{ccc}& & n\left(t\right)=-\frac{\partial logE\left(t\right)}{\partial log(t+t_0)}=1+\frac{1}{2}\frac{{\displaystyle H_0\left(k_0\sqrt{\nu (t+t_0)}\right)}}{H_1\left(k_0\sqrt{\nu (t+t_0)}\right)};\hfill \end{array}$$

(107)

$$\begin{array}{ccc}& & H_0\left(K\right)=\pi \mathsf{\Phi }\left(⌊\frac{1}{\sqrt{A{K}^b}\pi }⌋\right)\hfill \end{array}$$

(108)

This universal function is plotted in Figure 13.

The index jumps down from the level $n_l$ to $n_{l-1}$ at the discrete times $t_l$

$$\begin{array}{ccc}\hfill & & n_l=1+\frac{b\mathsf{\Phi }\left(l\right)}{{\displaystyle 4\sum _{n=2}^l\mathsf{\Phi }\left(n\right)log\left(\frac{n+1}{n}\right)+\frac{2}{3\zeta \left(3\right)}}};\hfill \end{array}$$

(109)

$$\begin{array}{ccc}& & t_l=-t_0+\frac{K_l^2}{\nu k_0^2};\hfill \end{array}$$

(110)

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

(111)

The dissipation begins at some high level

$$\begin{array}{c}\hfill l_0=\lfloor \frac{1}{\sqrt{A{\left(\sqrt{\nu t_0}\right)}^b}\pi }\rfloor \end{array}$$

(112)

and goes on with levels decreasing to $l=2$, where dissipation ends. These levels $n_l,K_l$ are universal numbers

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}l& K_l& n_l\\ 25& 4.1967& 1.5142\\ 24& 4.541& 1.5016\\ 23& 4.93005& 1.52034\\ 22& 5.37204& 1.49581\\ 21& 5.87709& 1.50558\\ 20& 6.45786& 1.50809\\ 19& 7.13043& 1.52651\\ 18& 7.91533& 1.49892\\ 17& 8.83925& 1.52387\\ 16& 9.93733& 1.49262\\ 15& 11.2566& 1.50047\\ 14& 12.8612& 1.50739\\ 13& 14.8404& 1.53085\\ 12& 17.3216& 1.4942\\ 11& 20.4916& 1.53148\\ 10& 24.6334& 1.48926\\ 9& 30.1929& 1.51975\\ 8& 37.9058& 1.51271\\ 7& 49.0588& 1.54133\\ 6& 66.0729& 1.48598\\ 5& 93.964& 1.55954\\ 4& 144.591& 1.51304\\ 3& 252.034& 1.55189\\ 2& 551.542& 1.47576\end{array}\right)\end{array}$$

(113)

In numerical or real experiments, these quantum jumps in the range $1.47$ to $1.55$ would be interpreted as noise and averaged out, producing $n\approx 1.51\pm 0.1$.

This behavior agrees with the effective index plotted in [9] in Fig.6a. The curves at that plot oscillate between $1.4$ and $1.6$, leveling at $1.5$; however, the resolution is not sufficient to find an exact match.

Measuring our discrete levels of effective index $n_l$ poses a new challenge for experiments.

6. Conclusions

There are several new results reported in this paper.

• We found continuum limit of distribution of scaling variables (33),(38) in the small Euler ensemble.
• The scaling law (77) for the conditional probability leads to the discrete energy spectrum (87) at the edge of each step of a devil’s staircase.
• The local slope of each step of energy spectrum is $-1$, not counting jumps. The average slope is $-1-b\approx -2.03$. The K41 spectrum $-5/3=-1.67$ lies in between these two quantum levels.
• This universal energy spectrum depending on $k\sqrt{\nu t}$ corresponds to the effective spatial scale $L\propto \sqrt{\nu t}$, or $m=\frac{1}{2}$ in notations of [9].
• The decay at a given time goes only at small enough wavelengths. At a fixed wavelength, the whole spectrum moves to the left with time, eventually stopping at some critical time, inversely proportional to the wavelength square.
• With the cutoff of our spectrum at the small wavelength corresponding to initial energy pumping ($E\left(k\right)\propto k^4$), we obtain discrete levels of the effective index $n\left(t\right)$. During decay, the index jumps downstairs from the level $n_l$ to $n_{l-1}$, until the lowest level is reached. These levels are universal numbers, given by the table (113).

In summary, we computed the continuum limit of our Euler ensemble and found quantum fractal laws, with decay indexes quantized at certain universal levels.

These levels are calculable in the continuum limit of the Euler ensemble. Our computation involves two numerical parameters $a,b$, which we computed in a large-scale simulation of the periodic Markov chain with $200,000,000$ steps on a supercomputer cluster.

It would be a new challenge to number theory to compute these two parameters analytically in the big Euler ensemble’s $N\to \infty$ limit.