Quantum Solution of Classical Turbulence. Decaying Energy Spectrum

Preprint

Article

Quantum Solution of Classical Turbulence. Decaying Energy Spectrum

Altmetrics

Downloads

846

Views

405

Comments

A peer-reviewed article of this preprint also exists.

Alexander Migdal^*

This version is not peer-reviewed

Submitted:

11 April 2024

Posted:

12 April 2024

Read the latest preprint version here

Alerts

Abstract

This paper brings to fruition recently found \cite{migdal2023exact} exact reduction of decaying turbulence in the Navier-Stokes equation in $3 + 1$ dimensions to a Number Theory problem of finding the statistical limit of the Euler ensemble. We reformulate the Euler ensemble as a Markov chain and show the equivalence of this formulation to the \textbf{quantum} statistical theory of $N$ fermions on a ring, with an external field related to the random fractions of $\pi$. We find the solution of this system in the turbulent limit $N\to \infty, \nu \to 0$ in terms of a complex trajectory (instanton) providing a saddle point to the path integral over the density of these fermions. This results in an analytic formula for the observable correlation function of vorticity in wavevector space. This is a full solution of decaying turbulence from the first principle without approximations or fitted parameters. The energy spectrum decays as $k^{-\frac{7}{2}}$. We compute resulting integrals in \Mathematica{} and present effective index $n(t) = -t\partial_t \log E $ for the energy decay as a function of time Fig.\ref{fig::NPlot}. The asymptotic value of the effective index in energy decay $n(\infty) = \frac{5}{4}$, but the universal function $n(t)$ is neither constant nor linear.

Keywords:

Subject: Physical Sciences - Fluids and Plasmas Physics

I. Prologue

Turbulence seems a hopelessly complex problem. Here is a picture from a recent simulation [37] that produced vortex lines of various shapes and sizes tangled like spaghetti.

Figure 1. The three-edimensional vortices in turbulent flow, from [37].

How could something so complex be described analytically?

The sheer complexity of this picture raises hopes of a simple statistical description.

Molecular dynamics with a large number of bodies is also hopelessly complex. Still, there is a simple statistical description that becomes increasingly accurate as the complexity of the original dynamical problem increases.

As Maxwell, Boltzmann, and Gibbs discovered, Newton’s mechanics cover the energy surface uniformly, given enough time. The rest of the statistical mechanics is a well-defined theory, though sometimes computations are difficult, like critical phenomena.

So, why not use the same approach in case of the Navier-Stokes turbulence?

Unfortunately, there is no known analog of the Gibbs distribution in the case of the turbulent flow.

So, the cornerstone of the turbulence theory must be a replacement for the Gibbs distribution.

Hopf made the first step in this direction in 1952 (see the recent review in [34]). He wrote down a functional equation satisfied by the probability distribution of the turbulent velocity field. Iterations of this equation in the nonlinear term in the Navier-Stokes equation generate expansion in inverse powers of viscosity, which is the opposite of our goal.

This equation is beautiful but mathematically as complex as the above picture of the vortex spaghetti. This apparent complexity placed the turbulence high in a hierarchy of physics theories, somewhere between the critical phenomena and the quark confinement problem.

Our theory, which started in 1993 (see Figure 2 for the history), says there is a much simpler version of the Hopf equation (the loop equation), still sufficient to determine the statistics.

The loop equation is equivalent to the Schrödinger equation in loop space. This deep analogy is not a poetic metaphor but a precise mathematical equivalence with far-reaching consequences, such as quantum oscillations resulting in regime changes of the effective indexes in scaling laws of classical turbulence.

We will present here and investigate the one-dimensional quantum statistical system satisfying the loop equation of the classical three-dimensional Navier-Stokes theory.

We have found a new example of duality between the strong coupling phase of one theory (a fluctuating velocity field in three dimensions) and the weak coupling of another (a one-dimensional quantum ring of Fermi particles).

This weak coupling limit is analytically solvable, yielding explicit formulas for observables in decaying turbulence, such as the energy spectrum. The experimental data for the decaying turbulence (see the compilation [35]) do not agree with each other, though there is a qualitative agreement with our predictions. We hope more precise future measurements will verify our theory.

II. Summary

Here are the main results reported in this paper.

We review the theory of Navier-Stokes loop equation, its relation to the Hopf functional equation, and the representation of the loop functional in terms of momentum loop.
We present the solution of the loop equation in the inviscid limit of the three-dimensional Navier-Stokes theory in terms of the Euler ensemble. This ensemble consists of a one-dimensional ring of Ising spins in an external field related to random fractions of $π$ .
We reduced the Markov process for the Euler ensemble in its fermionic representation (60) to the path integral.
This path integral in the continuum limit is dominated by a complex classical trajectory (instanton), satisfying a nonlinear ODE (130).
We solved this classical equation corresponding to the vorticity correlation function and found the spectrum of the linear operator for small fluctuations around this solution.
We computed the instanton’s contribution to the vorticity correlation function by using zeta regularization of the functional determinant of this linearized operator.
The continuum limit of this solution, $N \to \infty$ , corresponds to the inviscid limit of the decaying turbulence in the Navier-Stokes equation. Effective turbulent viscosity is $\tilde{ν} = ν N^{2} \to const$ . There are no quantum corrections to this instanton contribution in the turbulent limit.
We derived an analytic formula for energy spectrum and dissipation in finite system (233a), (233c), (223) and investigated it in Appendix A.
The energy spectrum decays asymptotically as

$\frac{{(\tilde{ν})}^{\frac{3}{2}}}{\sqrt{t}} κ^{- \frac{7}{2}}$

where $κ = k \sqrt{\tilde{ν} t}$ .
The total remaining energy decays as $E (t) \propto t^{- \frac{5}{4}}$
Both effective indexes $n (t) = - \frac{d log E}{d log t}, s (κ) = - \frac{\partial log E (k, t)}{\partial log k}$ are nontrivial functions of the logarithm of scale and time, approaching $n (\infty) = \frac{5}{4}, s (\infty) = \frac{7}{2}$ , see Figure 5), Figure 4
This oscillation of the effective decay indexes in classical turbulence exemplifies the quantum effects in our solution (singularities in the Mellin transform of the energy spectrum, related to complex zeros of the Riemann $ζ$ function).

Figure 3. Data for effective index

n (t) = - t \partial_{t} log E - 1

as compiled in [35]. Various curves correspond to various experiments, some for incompressible, others for compressible turbulence, at various Reynolds numbers. These data contradict each other, which indicates a lack of the true universal turbulent limit in these experiments.

Figure 3. Data for effective index

n (t) = - t \partial_{t} log E - 1

Figure 4. Orange curve is the effective index

s (t) = - κ \partial_{κ} log H (κ)

as a function of

κ = k \sqrt{\tilde{ν} t}

. It starts at zero and slowly rises to the limit

\frac{7}{2}

Figure 4. Orange curve is the effective index

s (t) = - κ \partial_{κ} log H (κ)

as a function of

κ = k \sqrt{\tilde{ν} t}

. It starts at zero and slowly rises to the limit

\frac{7}{2}

Figure 5. Green curve is the effective index

n (t) = - t \partial_{t} log E - 1

as a function of time. It starts at 1, and slowly rises to its limit

\frac{5}{4}

Figure 5. Green curve is the effective index

n (t) = - t \partial_{t} log E - 1

as a function of time. It starts at 1, and slowly rises to its limit

\frac{5}{4}

III. Introduction

Our paper will have two introductions: one for physicists and another for mathematicians.

The physics introduction discusses the possible correspondence between our theory and decaying turbulence as observed in real or numerical experiments. For a physicist, this theory provides a solution of the Hopf functional equation for the statistical distribution of velocity field in the unforced Navier-Stokes equation. This distribution provides a missing analog of the Gibbs distribution in decaying turbulence. There are strong indications that our theory applies to one of the two universality classes observed in these experiments.

The second introduction summarizes the mathematics behind the loop equation [20] and its solution [22] in terms of the Euler ensemble. This introduction is addressed to mathematicians, who can skip the first introduction and study this ensemble as a new Number Theory set with a conjectured relation to the decaying turbulence.

There is a third category of readers: applied mathematicians and engineers who are not interested in abstract theories but want to have data or formulas to compute the energy decay curve. These readers could go directly to the last Section VII where they will find final formulas (233c), (233a) and the Mathematica^® code [24,26], which allows them to compute both the spectrum and the dissipation.

The fourth category (my favorite) is curious but skeptical readers who doubt that quantum tricks can solve the hardest problem of classical physics. These readers would still have questions, even after reviewing this paper’s lengthy computations. There is the last Section IX, where some of their questions are answered.

A. Physical Introduction. The Energy Flow and Random Vorticity Structures

The decaying turbulence is an old subject addressed before within a weak turbulence framework (truncated perturbative expansion in inverse powers of viscosity

ν

in the forced Navier-Stokes equation). Some phenomenological models were also fitted to the experiments: see the recent review of these models and the experimental data in [35].

The perturbative approach is inadequate for the turbulence theory, which must be built by solving the Navier-Stokes equations in the strong coupling limit

ν \to 0

, opposite to the weak coupling. The problem of universality of the strong turbulence with and without random forcing is the first question to ask when building such a theory.

The experimental data for the energy decay in turbulent flows, fitted in [35] suggest the decay of the dissipation rate

E \sim t^{- n}

with

n \approx 1.2

1.6

depending on initial conditions (finite total momentum

\oint \vec{v} d^{3} r

or zero total momentum but finite total angular momentum

\int \vec{r} \times \vec{v} d^{3} r

). So, two universality classes of decaying turbulence were observed in these experiments.

It is unclear which data in [35] (if any) reached the strong turbulence limit corresponding to our regime. It is also possible that the stochastic forces added to the Navier-Stokes equation in simulations pollute the decaying turbulence. By design, these forces were supposed to trigger and amplify the spontaneous stochasticity of the turbulent flow. However, in our theory [22], this natural stochasticity is related to a dual quantum system and is discrete.

The Gaussian forcing with a continuous wavevector spectrum can distort these quantum stochastic phenomena as these forces stir the flow "with a large dirty spoon" all over the space at every moment. With forcing turned off at some moment and the turbulence reaching the universal stage in its subsequent decay, the energy dissipation would occur in vorticity structures deep inside the volume by the pure turbulent dynamics we are studying.

Some recent DNS reviewed in [35] used such a setup, with forcing turned off. They observed subsequent energy dissipation, apparently reaching one of the two universal regimes (with spectrum starting as

k 2

k 4

), presumably due to these random vorticity structures.

The following calculation supports this scenario. The energy balance in the pure Navier-Stokes reduces to the energy dissipation by the enstrophy in the bulk, compensated by the energy pumping by forces from the boundary (say, the large sphere around the flow).

The general identity, which follows from the Navier-Stokes if one multiplies both sides by

\vec{v}

and averages over an infinite time interval, reads:

\begin{matrix} \int_{V} d^{3} r 〈ν {\vec{ω}}^{2}〉 = - \int_{V} d^{3} r \partial_{β} 〈v_{β} (p + \frac{1}{2} v_{α}^{2}) + ν v_{α} (\partial_{β} v_{α} - \partial_{α} v_{β})〉 \end{matrix}

(1)

By the Stokes theorem, the right side reduces to the flow through the boundary

\partial V

of the integration region V. The left side is the dissipation in this volume, so we find:

\begin{matrix} E_{V} = - \int_{\partial V} d \vec{σ} \cdot 〈\vec{v} (p + \frac{1}{2} v_{α}^{2}) + ν \vec{ω} \times \vec{v}〉 \end{matrix}

(2)

This identity holds for an arbitrary volume. The left side represents the viscous dissipation inside V, while the right represents the energy flow through the boundary

\partial V

If there is a finite collection of vortex structures in the bulk, we can expand this volume to an infinite sphere; in this case, the

\vec{ω} \times \vec{v}

term drops as there is no vorticity at infinity.

Furthermore, the velocity in the Biot-Savart law decreases as

| \vec{r} |^{- 3}

at infinity, so that only the

\vec{v} p

term survives

\begin{matrix} 〈E_{V}〉 \to - \int_{\partial V} d \vec{σ} \cdot 〈\vec{v} p〉 \end{matrix}

(3)

This energy flow on the right side will stay finite in the limit of the expanding sphere in case the pressure grows as

p \to - \vec{f} \cdot \vec{r}

, where

\vec{f}

is the local force at a given point on a large sphere.

\begin{matrix} 〈E〉 = f_{α} lim_{R \to \infty} R^{3} \int_{S_{2}} n_{α} n_{β} 〈v_{β} (R \vec{n})〉 \end{matrix}

(4)

Where did we lose the Kolmogorov energy flow? It is still there for any finite volume surrounding the vortex sheet

\begin{matrix} 〈E〉 = - \int_{V} d^{3} r 〈v_{β} \partial_{β} p + v_{α} v_{β} \partial_{β} v_{α}〉 = \\ - \int_{V} d^{3} r 〈v_{α} v_{β} \partial_{β} v_{α}〉 - \int_{\partial V} d \vec{σ} \cdot 〈\vec{v} p〉 \end{matrix}

(5)

The first term is the Kolmogorov energy flow inside the volume V, and the second is the energy flow through the boundary.

Without finite force

\vec{f}

acting on the boundary, say, with periodic boundary conditions, the boundary integral would be absent, and we would recover the Kolmogorov relation.

\begin{matrix} 〈v_{α} v_{β} \partial_{β} v_{α}〉 = - 〈E〉 / V; \end{matrix}

(6)

This relation, together with space symmetry properties in

R_{d}

leads to the Kolmogorov three-point correlation

\begin{matrix} 〈v_{α} ({\vec{r}}_{0}) v_{β} ({\vec{r}}_{0}) v_{γ} (\vec{r} + {\vec{r}}_{0})〉 = \frac{E}{(d - 1) (d + 2) V} (δ_{α γ} r_{β} + δ_{β γ} r_{α} - \frac{2}{d} δ_{α β} r_{γ}); \end{matrix}

(7)

In the conventional approach, based on the time averaging of the Navier-Stokes equations, the periodic Gaussian random force

\vec{f} (\vec{r}, t)

is added to the right side. In this case, with periodic boundary conditions

\begin{matrix} 〈E_{V}〉 = - \int_{V} d^{3} r 〈v_{β} \partial_{β} p - v_{β} f_{β} (\vec{r}, t) + v_{α} v_{β} \partial_{β} v_{α}〉 = \int_{V} d^{3} r 〈v_{β} f_{β} (\vec{r}, t)〉 \end{matrix}

(8)

In the limit when the force becomes uniform in space, we recover another definition with

E = \vec{f} \cdot \vec{P}

, where

\vec{P} = \int_{V} d^{3} r \vec{v}

is a total momentum.

The turbulence phenomenon we study is a universal spontaneous stochasticity independent of the boundary conditions.

As long as there is an energy flow from the boundaries, the confined turbulence in the middle would dissipate this flow in singular vortex tubes. The spontaneous stochasticity results from the random distribution of these singular tubes inside the volume in the velocity flow picture [23]. In the dual picture of our recent theory [22], these are the random gaps in the momentum curve

\vec{P} (θ)

(see below).

The relation between the energy pumping on the large sphere and the distribution of the vortex blobs in bulk follows from the Biot-Savart integral

\begin{matrix} \vec{v} (r) = - \vec{\nabla} \times \int d^{3} r^{'} \frac{\vec{ω} (r^{'})}{4 π | r - r^{'} |} \end{matrix}

(9)

In case there is a finite total momentum of the fluid, there will also be a constant velocity added to the Biot-Savart integral. We are considering the fluid at rest, where no such term exists. Therefore, our theory belongs to the

k^{4}

universality class.

On a large sphere

\partial V

with radius

R \to \infty

\begin{matrix} lim_{R \to \infty} R^{3} \vec{v} (R \vec{n}) \propto \\ \frac{1}{4 π} \sum_{blobs} \int_{blob} d^{3} r^{'} \vec{ω} (r^{'}) \times ({\vec{r}}^{'} - < \vec{r} >_{blob}); \end{matrix}

(10)

Here

< \vec{r} >_{blob}

is the geometric center of each blob. Substituting this into the identity (4), we directly relate the energy pumping with the forces at the boundary and the blob’s dipole moments of vorticity.

No forcing inside the flow is needed for this energy pumping; the energy flow starts at the boundary and propagates to numerous singular vorticity blobs, where it is finally dissipated. The distribution of these vorticity blobs is all we need for the turbulence theory. The forcing is required only as a boundary condition at infinity.

These assumptions about confined turbulence as stochastic dynamics of isolated vortex structures were confirmed in a beautiful experimental work by William Irvine and collaborators at Chicago University ([15]). The energy was pumped in by vortex rings flying from eight corners of a large glass cube and colliding in the center, making a turbulent blob.

They measured the (approximately) Kolmogorov energy spectrum, proving that periodic boundary conditions were unnecessary.

The latest paper [16] also observed how the singular vortex structures move and reconnect inside this confined turbulence.

As for the decaying turbulence, these authors observed (William Irvine, private communication) two distinct decay regimes, not just a single power law like the old works [35].

Another critical comment: with the velocity correlations growing with distance by the approximate K41 law, even the forcing at the remote boundary would influence the potential part of velocity in bulk. This boundary influence makes the energy cascade picture non-universal; it may depend upon the statistics of the random forcing.

Two asymptotic regimes manifesting this non-universality were observed for the energy spectrum

E (k, t)

: one for initial spectrum

E (k, 0) \propto k^{2}

and another for

E (k, 0) \propto k^{4}

. The potential velocity part differs for these regimes, as the first adds a constant velocity to the Biot-Savart integral. In the general case, it will be a harmonic potential flow with certain boundary conditions at infinity, with explicit continuous dependence of the boundary forces.

Only the statistics of the rotational part of velocity, i.e., vorticity, could reach some universal regime independent of the boundary conditions at infinity. Certain discrete universality classes could exist as it is common in critical phenomena.

Unlike the potential part of velocity, the vorticity is localized in singular regions – tubes and sheets, sparsely filling the space, as observed in numerical simulations. The potential part of velocity drops in the loop equations, and the remaining stochastic motion of the velocity circulation is equivalent to the vorticity statistics. Therefore, our solutions [22] of the loop equations [20,21] describe the internal stochastization of the decaying turbulence by a dual discrete system.

B. Mathematical Introduction. The Loop Equation and Its Solution

We derived a functional equation for the so-called loop average or Wilson loop in turbulence a while ago. All the references to our previous works can be found in a recent review paper [21].

The path to an exact solution by a dimensional reduction in this equation was proposed in the 1993 paper but has just been explored (see Figure 2). 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, 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 [18,19]. Recently, some numerical methods were found to solve loop equations beyond perturbation theory [2,13]. The loop dynamics was extended to quantum gravity, where it was used to study nonperturbative phenomena [4,38].

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 [3,11,12,31] 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 [36].

Theoretically, we studied the loop equation in the confinement region (large circulation over large loop C) and justified the Area law suggested in ’93 on heuristic arguments. This law says that the tails of velocity circulation PDF in the confinement region are functions of the minimal area inside this loop. It was verified in DNS a few years ago [12], which triggered the further development of the geometric theory of turbulence [3,21,31]. In particular, the Area law was justified for flat and quadratic minimal surfaces, and an exact scaling law in confinement region

Γ \propto \sqrt{A r e a}

was derived [21]. The area law was verified with better precision in [11].

In the previous paper, [22], we have found a family of exact solutions of the loop equation for decaying turbulence [20,21]. This family describes a fixed trajectory of solutions with the universal time decay factor. The solutions are formulated in terms of the Wilson loop or loop average

\begin{matrix} Ψ [γ, C] = {〈exp (\frac{ı γ}{ν} \oint d \vec{C} (θ) \cdot \vec{v} (\vec{C} (θ)))〉}_{i n i t}; \end{matrix}

(11)

\begin{matrix} Ψ [γ, C] \Rightarrow {〈exp (\frac{ı γ}{ν} \oint d \vec{C} (θ) \vec{P} (θ))〉}_{E}; \end{matrix}

(12)

In the first equation (the definition), the averaging

〈\dots〉

goes over initial data for the solutions of the Navier-Stokes equation for velocity field

\vec{v} (\vec{r})

. In the second one (the solution), the averaging goes over the space of solutions

\vec{P} (θ)

of the loop equation [22]. We choose in this paper the parametrization of the loop with

ξ = \frac{θ}{2 π}

to match with the fermionic coordinates below (the parametrization is arbitrary, in virtue of parametric invariance of the loop dynamics).

The loop equation for the momentum loop

\vec{P} (θ)

follows from the Navier-Stokes equation for

\vec{v}

\begin{matrix} \partial_{t} v_{α} = ν \partial_{β} ω_{β α} - v_{β} ω_{β α} - \partial_{α} (p + \frac{v_{β}^{2}}{2}); \end{matrix}

(13)

\begin{matrix} \partial_{α} v_{α} = 0; \end{matrix}

(14)

\begin{matrix} ω_{β α} = \partial_{β} v_{α} - \partial_{α} v_{β} \end{matrix}

(15)

After some transformations, replacing velocity and vorticity with the functional derivatives of the loop functional, we found the following momentum loop equation in [21,22]

\begin{matrix} ν \partial_{t} \vec{P} = - γ^{2} {(Δ \vec{P})}^{2} \vec{P} + Δ \vec{P} (γ^{2} \vec{P} \cdot Δ \vec{P} + ı γ (\frac{{(\vec{P} \cdot Δ \vec{P})}^{2}}{Δ {\vec{P}}^{2}} - {\vec{P}}^{2})); \end{matrix}

(16)

\begin{matrix} \vec{P} (θ) \equiv \frac{\vec{P} (θ^{+}) + \vec{P} (θ^{-})}{2}; \end{matrix}

(17)

\begin{matrix} Δ \vec{P} (θ) \equiv \vec{P} (θ^{+}) - \vec{P} (θ^{-}); \end{matrix}

(18)

The momentum loop has a discontinuity

Δ \vec{P} (θ)

at every parameter

0 < θ \leq 1

, making it a fractal curve in complex space

C_{d}

. The details can be found in [21,22]. We will skip the arguments

t, θ

in these loop equations, as there is no explicit dependence of these equations on either of these variables.

This Anzatz represents a plane wave in loop space, solving the loop equation for the Wilson loop due to the lack of direct dependence of the loop operator on the shape of the loop.

The superposition of these plane wave solutions would solve the Cauchy problem in loop space: find the stochastic function

\vec{P} (θ)

t = 0

, providing the initial velocity field distribution. Formally, the initial distribution

W_{0} [P]

of the momentum field

\vec{P} (θ)

is given by inverse functional Fourier transform.

\begin{matrix} W_{0} [P] = \int D C δ^{3} (\vec{C} [0]) Ψ {[C, γ]}_{t = 0} exp (- \frac{ı γ}{ν} \int d \vec{C} (θ) \cdot \vec{P} (θ)) \end{matrix}

(19)

This path integral was computed in [21,22] for a special stochastic solution of the Navier-Stokes equation: the global rotation with Gaussian random rotation matrix. The initial velocity distribution is Gaussian, with a slowly varying correlation function. The corresponding loop field reads (we set

γ = 1

for simplicity in this section)

\begin{matrix} Ψ_{0} [C] \equiv Ψ {[C, 1]}_{t = 0} = exp (- \frac{1}{2 ν^{2}} \int_{C} d \vec{C} (θ) \cdot d \vec{C} (θ^{'}) f (\vec{C} (θ) - \vec{C} (θ^{'}))) \end{matrix}

(20)

where

f (\vec{r})

is the velocity correlation function

〈v_{α} (r) v_{β} (r^{'})〉 = (δ_{α β} - \partial_{α} \partial_{β} \partial_{μ}^{- 2}) f (r - r^{'})

(21)

The potential part drops out in the closed loop integral. The correlation function varies at the macroscopic scale, which means that one could expand it in the Taylor series

f (r - r^{'}) \to f_{0} - f_{1} {(r - r^{'})}^{2} + \dots

(22)

The first term

f_{0}

is proportional to initial energy density,

\frac{1}{2} 〈v_{α}^{2}〉 = \frac{d - 1}{2} f_{0}

(23)

and the second one is proportional to initial energy dissipation rate

E_{0}

f_{1} = \frac{E_{0}}{2 d (d - 1) ν}

(24)

where

d = 3

is the dimension of space. The constant term in (22) as well as

r^{2} + r^{' 2}

terms drop from the closed loop integral, so we are left with the cross-term

r r^{'}

, which reduces to a full square

\begin{matrix} Ψ_{0} [C] \to exp (- \frac{f_{1}}{ν^{2}} {(\oint d C_{α} (θ) C_{β} (θ))}^{2}) \end{matrix}

(25)

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

\begin{matrix} Ψ_{0} [C] \to const \int (d ϕ) exp (- ϕ_{α β}^{2}) exp (2 ı \frac{\sqrt{f_{1}}}{ν} ϕ_{μ ν} \oint d C_{μ} (θ) C_{ν} (θ)) \end{matrix}

(26)

The integration here involves all

\frac{d (d - 1)}{2} = 3

independent

α < β

components of the antisymmetric tensor

ϕ_{α β}

. Note that this is ordinary integration, not the functional one.

This distribution can be translated into the momentum loop space. Here is the resulting stochastic function

\vec{P} (θ)

, defined by the Fourier expansion on the circle

\begin{matrix} P_{α} (θ) = \sum_{odd n = 1}^{\infty} P_{α, n} e^{ı n θ} + {\bar{P}}_{α, n} e^{- ı n θ}; \end{matrix}

(27)

\begin{matrix} P_{α, n} = N (0, 1); \end{matrix}

(28)

\begin{matrix} {\bar{P}}_{α, n} = \frac{4 \sqrt{f_{1}}}{n ν} ϕ_{α β} P_{β, n}; \end{matrix}

(29)

\begin{matrix} ϕ_{α β} = - ϕ_{β α}; \end{matrix}

(30)

\begin{matrix} ϕ_{α β} = N (0, 1) \forall α < β; \end{matrix}

(31)

At fixed tensor

ϕ

the correlations are

\begin{matrix} {〈P_{α, n} P_{β, m}〉}_{t = 0} = \frac{4 \sqrt{f_{1}}}{m ν} δ_{n m} ϕ_{α β}; \end{matrix}

(32)

\begin{matrix} {〈P_{α} (θ) P_{β} (θ^{'})〉}_{t = 0} = 2 ı \frac{\sqrt{f_{1}}}{ν} ϕ_{α β} sign (θ^{'} - θ); \end{matrix}

(33)

\begin{matrix} Ψ_{0} [C] = {〈exp (\frac{ı}{ν} \oint d \vec{C} (θ) \vec{P} (θ))〉}_{P, ϕ} \end{matrix}

(34)

Though this special solution does not describe isotropic turbulence, it helps understand the mathematical properties of the loop technology. In particular, it shows the significance of the discontinuities of the momentum loop

\vec{P} (θ)

Rather than solving the Cauchy problem, we are looking for an attractor: the fixed trajectory for

\vec{P} (θ, t)

with some universal probability distribution related to the decaying turbulence statistics.

The following transformation reveals the hidden scaling invariance of decaying turbulence

\begin{matrix} \vec{P} = \sqrt{\frac{ν}{2 (t + t_{0})}} \frac{\vec{F}}{γ} \end{matrix}

(35)

The new vector function

\vec{F}

satisfies an equation

\begin{matrix} 2 \partial_{τ} \vec{F} = (1 - {(Δ \vec{F})}^{2}) \vec{F} + Δ \vec{F} (γ^{2} \vec{F} \cdot Δ \vec{F} + ı γ (\frac{{(\vec{F} \cdot Δ \vec{F})}^{2}}{Δ {\vec{F}}^{2}} - {\vec{F}}^{2})); \end{matrix}

(36)

\begin{matrix} τ = log (t + t_{0}) \end{matrix}

(37)

This equation is invariant under translations of the new variable

τ = log (t + t_{0})

, corresponding to the rescaling/translation of the original time.

\begin{matrix} t \Rightarrow λ t + (λ - 1) t_{0} \end{matrix}

(38)

There are two consequences of this invariance.

There is a fixed point for $\vec{F}$ .
The approach to this fixed point is exponential in $τ$ , which is power-like in original time.

Both of these properties were used in [22]: the first one was used to find a fixed point, and the second one was used to derive the spectral equation for the anomalous dimensions

λ_{i}

of decay

t^{- λ_{i}}

of the small deviations from these fixed points. In this paper, we only consider the fixed point, leaving the exciting problem of the spectrum of anomalous dimensions for future research.

C. The Big and Small Euler Ensembles

Let us remember the basic properties of the fixed point for

\vec{F}

in [22]. It is defined as a limit

N \to \infty

of the polygon

{\vec{F}}_{0} \dots {\vec{F}}_{N} = {\vec{F}}_{0}

with the following vertices

\begin{matrix} {\vec{F}}_{k} = \frac{\{cos (α_{k}), sin (α_{k}), i cos (\frac{β}{2})\}}{2 sin (\frac{β}{2})}; \end{matrix}

(39)

\begin{matrix} θ_{k} = \frac{k}{N}; β = \frac{2 π p}{q}; N \to \infty; \end{matrix}

(40)

\begin{matrix} α_{k + 1} = α_{k} + σ_{k} β; σ_{k} = \pm 1, β \sum σ_{k} = 2 π p r; \end{matrix}

(41)

The parameters

\hat{Ω}, p, q, r, σ_{0} \dots σ_{N} = σ_{0}

are random, making this solution for

\vec{F} (θ)

a fixed manifold rather than a fixed point. We suggested calling this manifold the big Euler ensemble of just the Euler ensemble.

It is a fixed point of (36) with the discrete version of discontinuity and principal value:

\begin{matrix} Δ \vec{F} \equiv {\vec{F}}_{k + 1} - {\vec{F}}_{k}; \end{matrix}

(42)

\begin{matrix} \vec{F} \equiv \frac{{\vec{F}}_{k + 1} + {\vec{F}}_{k}}{2} \end{matrix}

(43)

Both terms of the right side (36) vanish; the term proportional to

Δ \vec{F}

and the term proportional to

\vec{F}

. Otherwise, we would have

\vec{F} ‖ Δ \vec{F}

, leading to zero vorticity [22]. The ensemble of all the different solutions is called the big Euler ensemble. The integer numbers

σ_{k} = \pm 1

came as the solution of the loop equation, and the requirement of the rational

\frac{p}{q}

came from the periodicity requirement.

We can use integration (summation) by parts to write the circulation as follows (in virtue of periodicity):

\begin{matrix} \oint d \vec{C} (θ) \cdot \vec{P} (θ) = - \oint d \vec{P} (θ) \cdot \vec{C} (θ); \end{matrix}

(44)

\begin{matrix} \sum_{k} Δ C_{k} {\vec{P}}_{k} = - \oint Δ {\vec{P}}_{k} \cdot {\vec{C}}_{k}; \end{matrix}

(45)

A remarkable property of this solution

\vec{P} (θ, t)

of the loop equation is that even though it satisfies the complex equation and has an imaginary part, the resulting circulation (44) is real! The imaginary part of the

\vec{P} (θ, t)

does not depend on

θ

and thus drops from the integral

\oint d \vec{C} (θ) \cdot \vec{P} (θ, t)

There is, in general, a larger manifold of periodic solutions to the discrete loop equation, which has all three components of

{\vec{F}}_{k}

complex and varying along the polygon.

We could not find a global parametrization of such a solution1. Instead, we generated it numerically by taking a planar closed polygon and evolving its vertices

{\vec{F}}_{k}

by a stochastic process in the local tangent plane to the surface of the equations in multi-dimensional complex space.

We could not submit such a solution to an extra restriction

Im {\vec{F}}_{k} = const

needed to make circulation real. We cannot prove that such a general solution does not exist but rather take the Euler ensemble as a working hypothesis and investigate its properties.

This ensemble can be solved analytically in the statistical limit and has nice physical properties, matching the expected behavior of the decaying turbulence solution.

We assign equal weights to all elements of this set; we call this conjecture the ergodic hypothesis. This prescription is similar to assigning equal weights to each triangulation of curved space with the same topology in dynamically triangulated quantum gravity [1]. Mathematically, this is the most symmetric weight assignment, and there are general expectations that various discrete theories converge into the same symmetry classes of continuum theories in the statistical limit. This method works remarkably well in two dimensions, providing the same correlation functions as continuum gravity (Liouville theory).

The fractions

\frac{p}{q}

with fixed denominator are counted by Euler totient function

φ (q)

[10]

\begin{matrix} φ (q) = \sum_{\begin{matrix} p = 1 \\ (p, q) \end{matrix}}^{q - 1} 1 = q \prod_{p | q} (1 - \frac{1}{p}); \end{matrix}

(46)

For example

φ (16) = 8

and

φ (17) = 16

In some cases, one can analytically average over spins

σ

in the big Euler ensemble, reducing the problem to computations of averages over the small Euler ensemble

E (N) : N, p, q, r

with the measure induced by averaging over the spins in the big Euler ensemble.

In this paper, we perform this averaging over

σ

analytically, without any approximations, reducing it to a partition function of a certain quantum mechanical system with Fermi particles. This partition function is calculable using a WKB approximation in the statistical limit

N \to \infty

IV. The Markov Chain and Its Fermionic Representation

Here is a new representation of the Euler ensemble, leading us to the exact analytic solution below.

We start by replacing independent random variables

σ

with fixed sum by a Markov process, as suggested in [22]. We start with n random values of

σ_{i} = 1

and remaining

N - n

values of

σ_{i} = - 1

. Instead of averaging over all of these values simultaneously, we follow a Markov process of picking

σ_{N}, \dots σ_{1}

one after another. At each step, there will be

M = N, \dots 0

remaining

σ

. We get a transition

n \Rightarrow n - 1

with probability

\frac{n}{M}

and

n \Rightarrow n

with complementary probability.

Multiplying these probabilities and summing all histories of the Markov process is equivalent to the computation of the product of the Markov matrices

\begin{matrix} \prod_{M = 1}^{N} Q (M); \end{matrix}

(47)

\begin{matrix} Q (M) | n 〉 = \frac{M - n}{M} | n 〉 + \frac{n}{M} | n - 1 〉; \end{matrix}

(48)

This Markov process will be random until

n = 0

. After that, all remaining

σ_{k}

will have negative signs and be taken with probability 1, keeping

n = 0

The expectation value of some function

X (\{σ\})

reduces to the matrix product

\begin{matrix} P [X] = \sum_{n = 0}^{N_{+}} 〈n| (\prod_{M = 1}^{N} \hat{Q} (M)) \cdot X \cdot |N_{+}〉; \end{matrix}

(49)

\begin{matrix} \hat{Q} (M) \cdot X | n 〉 = \frac{n}{M} X (σ_{M} \to 1) | n - 1 〉 + \frac{M - n}{M} X (σ_{M} \to - 1) | n 〉 \end{matrix}

(50)

Here

N_{+} = (N + \sum σ_{l}) / 2 = (N + q r) / 2

is the number of positive sigmas. The operator

\hat{Q} (M)

sets in

X | n 〉

the variable

σ_{M}

to 1 with probability

\frac{n}{M}

and to

- 1

with complementary probability. The generalization of the Markov matrix

Q (M)

to the operator

\hat{Q} (M)

will be presented shortly.

Once the whole product is applied to X, all the sigma variables in all terms will be specified so that the result will be a number.

This Markov process is implemented as a computer code in [6], leading to a fast simulation with

O (N^{0})

memory requirement.

Now, we observe that quantum Fermi statistics can represent the Markov chain of Ising variables. Let us construct the operator

\hat{Q} (M)

with Fermionic creation and annihilation operators, with occupation numbers

ν_{k} = (1 + σ_{k}) / 2 = (0, 1)

. These operators obey (anti)commutation relations, and they create/annihilate

σ = 1

state as follows (with Kronecker delta

δ [n] \equiv δ_{n, 0}

\begin{matrix} {[a_{i}, a_{j}^{†}]}_{+} = δ_{i j}; \end{matrix}

(51)

\begin{matrix} {[a_{i}, a_{j}]}_{+} = {[a_{i}^{†}, a_{j}^{†}]}_{+} = 0; \end{matrix}

(52)

\begin{matrix} a_{n}^{†} | σ_{1}, \dots, σ_{N} 〉 = δ [σ_{n} + 1] | σ_{1}, \dots, σ_{n} \to 1, \dots σ_{N} 〉; \end{matrix}

(53)

\begin{matrix} a_{n} | σ_{1}, \dots, σ_{N} 〉 = δ [σ_{n} - 1] | σ_{1}, \dots, σ_{n} \to - 1, \dots σ_{N} 〉; \end{matrix}

(54)

\begin{matrix} {\hat{ν}}_{n} = a_{n}^{†} a_{n}; \end{matrix}

(55)

\begin{matrix} {\hat{ν}}_{n} | σ_{1}, \dots, σ_{N} 〉 = δ [σ_{n} - 1] | σ_{1}, \dots, σ_{N} 〉 \end{matrix}

(56)

The number

n (M)

of positive sigmas

\sum_{l = 1}^{M} δ [σ_{l} - 1]

coincides with the occupation number of these Fermi particles.

\begin{matrix} \hat{n} (M) = \sum_{l = 1}^{M} {\hat{ν}}_{l}; \end{matrix}

(57)

This relation leads to the representation

\begin{matrix} \hat{Q} (M) = {\hat{ν}}_{M} \frac{\hat{n} (M)}{M} + (1 - {\hat{ν}}_{M}) \frac{M - \hat{n} (M)}{M}; \end{matrix}

(58)

The variables

σ_{l}

can also be expressed in terms of this operator algebra by using

\begin{matrix} {\hat{σ}}_{l} = 2 {\hat{ν}}_{l} - 1 . \end{matrix}

(59)

The Wilson loop in (11) can now be represented as an average over the small Euler ensemble

E (N)

of a quantum trace expression

\begin{matrix} Ψ [γ, C] = & \frac{{〈tr (\hat{Z} (q r) exp (\frac{i γ}{ν} \sum_{l} Δ {\vec{C}}_{l} \cdot \hat{Ω} \cdot {\vec{P}}_{l} (t)) \prod_{M = 1}^{N} \hat{Q} (M))〉}_{\hat{Ω}, E (N)}}{{〈tr (\hat{Z} (q r) \prod_{M = 1}^{N} \hat{Q} (M))〉}_{E (N)}}, \end{matrix}

(60)

\begin{matrix} \hat{Z} (s) = & \oint \frac{d ω}{2 π} exp (ı ω (\sum_{l} {\hat{σ}}_{l} - s)); \end{matrix}

(61)

\begin{matrix} Δ {\vec{C}}_{l} = & \vec{C} (\frac{l + 1}{N}) - \vec{C} (\frac{l}{N}), \end{matrix}

(62)

\begin{matrix} {\vec{P}}_{l} (t) = & \sqrt{\frac{ν}{2 (t + t_{0})}} \frac{{\vec{F}}_{l}}{γ}, \hat{Ω} \in O (3), \end{matrix}

(63)

\begin{matrix} {\vec{F}}_{l} = & \frac{\{cos ({\hat{α}}_{l}), sin ({\hat{α}}_{l}), 0\}}{2 sin (\frac{β}{2})}, \end{matrix}

(64)

\begin{matrix} E (N) : & p, q, r \in Z with 0 < p < q < N, \gcd (p, q) = 1, - N \leq q r \leq N, \end{matrix}

(65)

\begin{matrix} {\hat{α}}_{l} = & β \sum_{k = 1}^{l - 1} (2 {\hat{ν}}_{k} - 1); \end{matrix}

(66)

The last component of the vector

{\vec{F}}_{l}

is set to 0 as this component does not depend on l and yields zero in the sum over the loop

\sum_{l} Δ {\vec{C}}_{l} = 0

The proof of equivalence to the combinatorial formula with an average over

σ_{l} = \pm 1

can be given using the following Lemma (obvious for a physicist).

Lemma 1.

The operators

{\hat{ν}}_{l}

all commute with each other.

Proof.

Using commutation relations, we can write

\begin{matrix} {\hat{ν}}_{l} {\hat{ν}}_{n} = a_{l}^{†} (δ_{l n} - a_{n}^{†} a_{l}) a_{n} = a_{l}^{†} a_{n} δ_{l n} - a_{l}^{†} a_{n}^{†} a_{l} a_{n} \end{matrix}

(67)

Interchanging indexes

l, n

in this relation, we see that the first term does not change due to Kronecker delta, and the second term does not change because

a_{l}^{†}, a_{n}^{†}

anti-commute, as well as

a_{l}, a_{n}

, so the second term is symmetric as well. Therefore,

{\hat{ν}}_{l} {\hat{ν}}_{n} = {\hat{ν}}_{n} {\hat{ν}}_{l}

□

Quantum Trace Theorem.The trace formula (60) equals the expectation value of the momentum loop ansatz (12), (35), (39) in the big Euler ensemble.

Proof.

As all the operators

{\hat{ν}}_{l}

commute with each other, the operators

\hat{Q} (M)

can be applied in arbitrary order to the states

Σ = | σ_{1}, \dots σ_{N} 〉

involved in the trace. The same is true about individual terms in the circulation in the exponential of the Wilson loop. These terms

{\vec{F}}_{l}

involve the operators

{\hat{α}}_{l}

, which commute with each other and with each

\hat{Q} (M)

. Thus, we can use the ordered product of the operators

{\hat{G}}_{l} = \hat{Q} (l) exp (ı ω {\hat{σ}}_{l} + \frac{i γ}{ν} Δ {\vec{C}}_{l} \cdot {\vec{P}}_{l} (t))

. Each of the operators

{\hat{G}}_{l}

acting in turn on arbitrary state

Σ

will create two terms with

δ [σ_{l} \pm 1]

. The exponential in

{\hat{G}}_{l}

will involve

{\hat{σ}}_{k}, k \leq l

. As a result of the application of the operator

{\hat{Z}}_{l} = \prod_{k = 1}^{l} {\hat{G}}_{k}

to the state vector

Σ

we get

2^{l}

terms with

Σ \prod_{k = 1}^{l} δ [σ_{k} - η_{k}], η_{k} = \pm 1

. The factors

{\hat{Z}}_{l}

will involve only

{\hat{σ}}_{k}, k \leq l

, which are all reduced to

η_{k}, k \leq l

in virtue of the product of the Kronecker deltas. Multiplying all operators

{\hat{G}}_{M}

will lead to superposition

{\hat{Π}}_{N}

2^{N}

terms, each with product

\prod_{M = 1}^{N} δ [σ_{M} - η_{M}]

with various choices of the signs

η_{i} = \pm 1

for each i. Furthermore, the product of Kronecker deltas will project the total sum of

2^{N}

combinations of the states

Σ

in the trace

tr \dots

to a single term corresponding to a particular history

η_{1}, \dots η_{N}

of the Markov process. The product of Kronecker deltas in each history will be multiplied by the same state vector

Σ

, by the product of Markov transition probabilities, and by the exponential

exp (\frac{i γ}{ν} \sum_{l} Δ {\vec{C}}_{l} \cdot \hat{Ω} \cdot {\vec{P}}_{l} (t))

with the operators

\hat{σ}

{\vec{P}}_{k} (t)

replaced by numbers

η

leading to the usual numeric

{\vec{P}}_{l} (t)

. The transition probabilities of the Markov process are designed to reproduce combinatorial probabilities of random sigmas, adding up to one after summation over histories [32]. The integration over

ω

will produce

δ [\sum_{l} {\hat{η}}_{l} - s]

. This delta function will reduce the trace to the required sum over all histories of the Markov process with a fixed

\sum_{l} η_{l}

. □

We have found a third vertex of the triangle of equivalent theories: the decaying turbulence in three-dimensional space, the fractal curve in complex space, and Fermi particles on a ring. By degrees of freedom, this is a one-dimensional Fermi-gas in the statistical limit

N \to \infty

. However, there is no local Hamiltonian in this quantum partition function, just a trace of certain products of operators in Fock space. So, an algebraic (or quantum statistical) problem remains to find the continuum limit of this theory of the fermion ring.

This problem is addressed in the next section.

V. The Continuum Limit

As we shall shortly see, in the continuum limit

N \to \infty

, the accumulated numbers of Fermi particles

ν_{k} = 1

and Dirac holes

ν_{k} = 0

tend to some classical functions of "position"

ξ = \frac{k}{N}

, leading to the exact solution.

A. Path Integral Over Markov Histories

Let us represent the product

Π_{N}

of the transitional probabilities of the particular history of the Markov processes as follows (with

n_{\pm} \equiv n_{\pm} (l), Δ n_{\pm} = - 1

)

\begin{matrix} Π_{N} = exp (N Λ_{N}); \end{matrix}

(68)

\begin{matrix} Λ_{N} = \frac{1}{N} \sum_{l} (Δ n_{+} log (\frac{n_{+}}{n_{+} + n_{-}}) + Δ n_{-} log (\frac{n_{-}}{n_{+} + n_{-}})); \end{matrix}

(69)

\begin{matrix} n_{+} = \sum_{k \leq l} ν_{k}; \end{matrix}

(70)

\begin{matrix} n_{-} = \sum_{k \leq l} (1 - ν_{k}); \end{matrix}

(71)

These

n_{\pm}

are net numbers of

η = \pm 1

in terms of Ising spins or occupation numbers

ν_{k} = (1, 0)

in the Fermi representation. There is an extra constraint on the Markov process

\begin{matrix} n_{+} + n_{-} = l; \forall l \end{matrix}

(72)

which follows from the above definition in terms of the occupation numbers. We can redefine

n_{\pm}

as N times the piecewise constant functions.

\begin{matrix} n_{\pm} = N f_{\pm} (ξ); \end{matrix}

(73)

\begin{matrix} f_{\pm} (ξ) = \sum_{k = 1}^{N ξ} \frac{ν_{k}}{N}; \end{matrix}

(74)

\begin{matrix} f_{\pm}^{'} (ξ) = \sum_{l = 1}^{N} δ (ξ - \frac{l}{N}) \sum_{k = 1}^{l} \frac{ν_{k}}{N}; \end{matrix}

(75)

\begin{matrix} 0 < ξ < 1; \end{matrix}

(76)

The sums can be rewritten as Lebesgue integrals

\begin{matrix} Λ_{N} = \int_{1}^{0} (d f_{+} (ξ) log (\frac{f_{+} (ξ)}{ξ}) + d f_{-} (ξ) log (\frac{f_{-} (ξ)}{ξ})) \end{matrix}

(77)

The sum over histories of the Markov process will become a path integral over the difference

ϕ = f_{+} (ξ) - f_{-} (ξ)

\begin{matrix} \sum_{η_{.} = \pm 1} exp (N (Λ_{N} + ı Λ_{N}^{(1)})) \to \int D ϕ exp (N (Λ_{N} + ı Λ_{N}^{(1)})) \end{matrix}

(78)

This path integral will be dominated by the "classical history," maximizing the product of transitional probabilities if such a classical trajectory exists.

The first term (without the circulation term) brings the variational problem

\begin{matrix} max_{ϕ} Λ_{N} [ϕ]; \end{matrix}

(79)

\begin{matrix} Λ_{N} [ϕ] = \int_{1}^{0} d ξ (\frac{d f_{+}}{d ξ} log (\frac{f_{+}}{ξ}) + \frac{d f_{-}}{d ξ} log (\frac{f_{-}}{ξ})); \end{matrix}

(80)

\begin{matrix} f_{\pm} (ξ) = \frac{1}{2} (ξ \pm ϕ (ξ)); \end{matrix}

(81)

\begin{matrix} f_{\pm} (ξ) \geq 0; \end{matrix}

(82)

This problem is, however, a degenerate one, as the functional reduces to the integral of the total derivative:

\begin{matrix} \frac{δ Λ_{N} [ϕ]}{δ ϕ (ξ)} = 0; \\ Λ_{N} [ϕ] = \int_{1}^{0} d (f_{+} log f_{+} + f_{-} log f_{-}) + 1 + \int_{0}^{1} d ξ log ξ = \end{matrix}

(83)

\begin{matrix} - \frac{1}{2} (1 - s) log (1 - s) - \frac{1}{2} (1 + s) log (1 + s) + log (2); \end{matrix}

(84)

It depends on the boundary condition

ϕ (0) = 0, ϕ (1) = s

but not on the shape of

ϕ (ξ)

This expression matches the Stirling formula for the logarithm of the binomial coefficient in the combinatorial solution [22] for the Euler ensemble

\begin{matrix} lim_{N \to \infty} \frac{log (\binom{N}{N (1 + s) / 2})}{N} = - \frac{1}{2} (1 - s) log (1 - s) - \frac{1}{2} (1 + s) log (1 + s) + log (2) \end{matrix}

(85)

This

Λ (s) = Λ_{\infty} (s)

is a smooth even function of s taking positive values from

Λ (\pm 1) = 0

to the maximal value

Λ (0) = log (2)

(see Figure 6).

Now, let us add the circulation term to the exponential of the partition function (60). This term can be directly expressed in terms of the difference between our two densities

N ϕ (ξ) = N f_{+} (ξ) - N f_{-} (ξ)

\begin{matrix} ı N Λ_{N}^{(1)} [ϕ, C_{Ω}] = \frac{ı}{\sqrt{2 ν t}} \int_{0}^{1} d {\vec{C}}_{Ω} (ξ) \cdot \vec{F} (ξ); \end{matrix}

(86)

\begin{matrix} \vec{F} (ξ) = \frac{\{sin (β N ϕ (ξ)), cos (β N ϕ (ξ)), 0\}}{2 sin (β / 2)}; \end{matrix}

(87)

\begin{matrix} {\vec{C}}_{Ω} (θ) = \hat{Ω} \cdot \vec{C} (θ); \end{matrix}

(88)

The key assumption is, of course, the existence of the smooth limit of the charge density

ϕ (ξ)

of these fermions when they are densely covering this loop.

We are working with

α (ξ) = β N ϕ (ξ)

in the following.

The measure for paths

[D α]

is undetermined. The derivatives of these alphas were quantized in the original Fermi theory: each step

α^{'} (ξ) \approx N Δ α = N β σ = \pm N β

As we demonstrate below, in continuum theory, this discrete distribution can be replaced by a Gaussian distribution with the same mean square

\begin{matrix} \sum_{α^{'} = \pm N β} \leftrightarrow \int d α^{'} exp (- \frac{{(α^{'})}^{2}}{2 N^{2} β^{2}}) \end{matrix}

(89)

To demonstrate that, we consider in the critical region

β^{2} \sim N^{- 1} \to 0

the most general term that arises in the moments of the circulation in (86) (see [5] for some exact computations of these moments)

\begin{matrix} 2^{- N} \sum_{σ_{i} = \pm 1} exp (ı β \sum_{i} k_{i} σ_{i}) = \prod_{i} cos β k_{i} \to exp (- β^{2} / 2 \sum_{i} k_{i}^{2}) \end{matrix}

(90)

where

k_{i}

are some integers. With a large number N of these integers, the sum in the exponential becomes an integral, which is equivalent to a Gaussian integral

\begin{matrix} exp (- β^{2} / 2 \sum_{i} k_{i}^{2}) = \prod_{i} \int_{- \infty}^{\infty} \frac{d σ_{i}}{\sqrt{2 π}} exp (- σ_{i}^{2} / 2) exp (ı β k_{i} σ_{i}) \to exp (- N β^{2} / 2 \int_{0}^{1} d ξ k {(ξ)}^{2}) \end{matrix}

(91)

We arrive at the standard path integral measure

\begin{matrix} \int [D α] = \int D α (ξ) exp (- \int_{0}^{1} d ξ \frac{{(α^{'})}^{2}}{2 N β^{2}}); \end{matrix}

(92)

\begin{matrix} \frac{\int [D α] exp (ı N \int_{0}^{1} d ξ α (ξ) K (ξ))}{\int [D α]} = exp (- \frac{1}{2} N^{2} \int \int d ξ_{1} d ξ_{2} K (ξ_{1}) K (ξ_{2}) G (ξ_{1}, ξ_{2})); \end{matrix}

(93)

\begin{matrix} G (ξ_{1}, ξ_{2}) = 〈α (ξ_{1}) α (ξ_{2})〉; \end{matrix}

(94)

The next section will compute this Green’s function

G (ξ_{1}, ξ_{2})

Thus, we arrive at the following path integral in the continuum limit

\begin{matrix} Ψ [C] = \frac{\sum_{p < q; (p, q)} \int_{Ω \in O (3)} d Ω \int [D α] exp (ı \frac{\int_{0}^{1} d ξ Im (C_{Ω}^{'} (ξ) exp (ı α (ξ)))}{2 sin (π p / q) \sqrt{2 ν (t + t_{0})}})}{\sum_{p < q; (p, q)} | O (3) | \int [D α]}; \end{matrix}

(95)

\begin{matrix} C_{Ω} (θ) = \vec{C} (θ) \cdot \hat{Ω} \cdot {ı, 1, 0}; \end{matrix}

(96)

We get the

U (1)

statistical model with the boundary condition

α (1) = α (0) + β N s

. The period

β N s = 2 π p r

is a multiple of

2 π

, which is irrelevant at

N \to \infty

. The effective potential for this theory is a linear function of the loop slope

{\vec{C}}^{'} (ξ)

This model is yet another representation of the Euler ensemble, suitable for the continuum limit.

B. Matching Path Integral with Combinatorial Sums in Big Euler Ensemble

The results of the path integration over

α

must match the combinatorial calculations with

σ_{l} = \pm 1

in the limit of large N. Without the interaction provided by the circulation term in (95), this path integral is dominated by a linear trajectory

\begin{matrix} α_{c l} (ξ) = β N ξ s; \end{matrix}

(97)

We already saw the match between the classical Action

Λ_{N} [ϕ (ξ) = ξ s]

and the asymptotic value of the logarithm of the Binomial coefficient of the combinatorial solution for the sum over

σ

variables.

Let us verify some examples of the expectation values over

σ

. The simplest is (with

n \neq m

)

\begin{matrix} {〈σ_{n} σ_{m}〉}_{\sum σ = N s} \end{matrix}

(98)

The direct calculation using methods of [5,22] leads to

\begin{matrix} {〈σ_{n} σ_{m}〉}_{\sum σ = N s} = - \oint \frac{d ω}{2 π Z} e^{ı ω N s} cos ((N - 2) ω) {sin}^{2} (ω) = \frac{(1 - N s^{2})}{2^{N - 3} Z N (1 - s^{2})} (\binom{N - 2}{\frac{1}{2} (N + N s - 2)}); \end{matrix}

(99)

\begin{matrix} Z = \oint \frac{d ω}{2 π} e^{ı ω N s} cos (N ω) = 2^{- N} (\binom{N}{\frac{1}{2} (s N + N)}) \end{matrix}

(100)

Using Gamma function properties, this ratio simplifies to

\begin{matrix} \frac{N s^{2} - 1}{N - 1} \end{matrix}

(101)

This result can be derived from symmetry without any integration.

\begin{matrix} {〈σ_{n} σ_{m}〉}_{\sum σ = N s} = A (N, s) (1 - δ_{n m}) + δ_{n m}; \end{matrix}

(102)

\begin{matrix} \sum_{n, m} {〈σ_{n} σ_{m}〉}_{\sum σ = N s} = s^{2} N^{2} = A (N, s) N (N - 1) + N \end{matrix}

(103)

\begin{matrix} A (N, s) = \frac{N s^{2} - 1}{N - 1} \end{matrix}

(104)

The same limit

A (\infty, s) = s^{2}

follows from the classical trajectory

\begin{matrix} {〈σ_{n} σ_{m}〉}_{\sum σ = N s} \to \frac{α_{c l}^{'} (ξ)}{β N} \frac{α_{c l}^{'} (ξ^{'})}{β N} = s^{2} \end{matrix}

(105)

Let us consider less trivial example [5,22]

\begin{matrix} U_{n, m} \to \sum_{k = n}^{m - 1} exp (ı α_{k, n}); \end{matrix}

(106)

\begin{matrix} α_{k, n} = β \sum_{\begin{matrix} l = 0 \\ l \neq n \end{matrix}}^{k} σ_{l}; \end{matrix}

(107)

We shall set

s = 0

, as this is the leading contribution to the partition function. The expectation value of

U_{n, m}

in our continuum limit becomes

\begin{matrix} 〈U_{n, m}〉 = N \int_{ξ_{1}}^{ξ_{2}} d ξ 〈exp (ı α (ξ))〉 = N \int_{ξ_{1}}^{ξ_{2}} d ξ exp (- \frac{1}{2} G (ξ, ξ)) \end{matrix}

(108)

Here

G (ξ_{1}, ξ_{2})

is the Green’s function corresponding to a 1D particle on a line interval

ξ \in (0, 1)

, introduced in the previous section. It satisfies the equation, which follows from our Gaussian Action

\begin{matrix} \partial_{ξ}^{2} G (ξ, ξ^{'}) = - β^{2} N δ (ξ - ξ^{'}); \end{matrix}

(109)

\begin{matrix} G (0, ξ^{'}) = G (ξ, 0) = 0; \end{matrix}

(110)

The solution is

\begin{matrix} G (ξ, ξ^{'}) = \frac{1}{2} β^{2} N (ξ + ξ^{'} - | ξ - ξ^{'} |) \end{matrix}

(111)

Thus, we find

\begin{matrix} 〈U_{n, m}〉 = N \int_{ξ_{1}}^{ξ_{2}} d ξ exp (- \frac{1}{2} β^{2} N ξ) = \frac{2}{β^{2}} (exp (- y / 2) - exp (- x / 2)); \end{matrix}

(112)

\begin{matrix} x = β^{2} N ξ_{1}; \end{matrix}

(113)

\begin{matrix} y = β^{2} N ξ_{2}; \end{matrix}

(114)

in agreement with [5,22] in the critical region

N \to \infty, β^{2} \sim 1 / N

. Finally, the expectation value

\begin{matrix} 〈U_{n, m} {\bar{U}}_{n, m}〉 = \sum_{l = n}^{m - 1} \sum_{k = n}^{m - 1} 〈exp (ı α_{k n} - ı α_{l n})〉 \end{matrix}

(115)

Here, the Gaussian path integration yields

\begin{matrix} 〈U_{n, m} {\bar{U}}_{n, m}〉 \to N^{2} \int_{ξ_{1}}^{ξ_{2}} d ξ \int_{ξ_{1}}^{ξ_{2}} d ξ^{'} exp (- \frac{1}{2} (G (ξ, ξ) + G (ξ^{'}, ξ^{'}) - 2 G (ξ, ξ^{'}))) = \\ N^{2} \int_{ξ_{1}}^{ξ_{2}} d ξ \int_{ξ_{1}}^{ξ_{2}} d ξ^{'} exp (\frac{1}{2} β^{2} N | ξ - ξ^{'} |) = \frac{4}{β^{4}} (2 exp ((y - x) / 2) + x - y - 2) \end{matrix}

(116)

This result also agrees with combinatorial computations in [5,22].

C. Small Euler Ensemble in Statistical Limit

The remaining problem is averaging over the variables

N, p, q, r

of the small Euler ensemble.

The variable

s = \frac{q r}{N}

is distributed between

- 1, 1

with the binomial weight [5,22]

(\binom{N}{N (1 + s) / 2})

peaked at

s = 0

. There is a finite term coming from

r = 0

plus a continuum spectrum coming from large r

\begin{matrix} W (r) = \{\begin{matrix} 1 & if r = 0; \\ \frac{\sqrt{2 π N}}{q} exp (- \frac{{(q r)}^{2}}{2 N}) & otherwise; \end{matrix} \end{matrix}

(117)

As it was conjectured in [22] and supported by rigorous estimates in [5], the

r = 0

term dominates the sums, after which the variables

y, x

can be treated as continuous variables.

The variable p at fixed q has a discrete distribution

\begin{matrix} f_{p} (p ∣ q) = \frac{\sum_{\begin{matrix} p = 1 \\ (p, q) \end{matrix}}^{q - 1} δ (p - n)}{φ (q)} \end{matrix}

(118)

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

\begin{matrix} X (p, q) = \frac{1}{q^{2}} {cot}^{2} (\frac{π p}{q}) \end{matrix}

(119)

This distribution for

X (p, q)

at fixed

q \to \infty

can be found analytically, using newly established relations for the cotangent sums (see Appendix in [22], and Mathematica^® notebook [25]). Asymptotically, at large q, these relations read

\begin{matrix} 〈X^{n}〉 \equiv lim_{q \to \infty} \frac{\sum_{\begin{matrix} p = 1 \\ (p, q) \end{matrix}}^{q - 1} X {(p, q)}^{n}}{q} = δ_{n, 0} + \frac{2 π^{- 2 n} ζ (2 n)}{(2 n + 1) ζ (2 n + 1)} \end{matrix}

(120)

This relation can be transformed further as

\begin{matrix} 〈X^{n}〉 = \{\begin{matrix} 1 & if n = 0 \\ \frac{2 \sum_{k = 1}^{\infty} φ (k) k^{- (2 n + 1)}}{(2 n + 1) π^{2 n}} & if n > 0 \end{matrix} \end{matrix}

(121)

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

\begin{matrix} 〈X^{n}〉 = \int_{0}^{\infty} f_{X} (X) d X X^{n}; \end{matrix}

(122)

\begin{matrix} f_{X} (X) = (1 - α) δ (X) + π X \sqrt{X} Φ (⌊ \frac{1}{π \sqrt{X}} ⌋); \\ α = π \int_{0}^{\infty} X \sqrt{X} d X Φ (⌊ \frac{1}{π \sqrt{X}} ⌋) = \end{matrix}

(123)

\begin{matrix} \frac{2 π}{5} \sum_{1}^{\infty} Φ (k) (\frac{1}{{(π k)}^{5}} - \frac{1}{{(π (k + 1))}^{5}}) = \frac{2}{5 π^{4}} \sum_{1}^{\infty} \frac{φ (k)}{k^{5}} = \frac{1}{225 ζ (5)} \end{matrix}

(124)

where

Φ (n)

is the totient summatory function

\begin{matrix} Φ (q) = \sum_{n = 1}^{q} φ (n) \end{matrix}

(125)

The distribution can also be rewritten as an infinite sum

\begin{matrix} \int d x f_{X} (x) F (x) = (1 - α) F (0) + π \sum_{n = 1}^{\infty} φ (n) \int_{0}^{\frac{1}{π^{2} n^{2}}} x^{\frac{3}{2}} d x F (x) \end{matrix}

(126)

The normalization of this distribution comes out 1 as it should, with factor

1 - α

in front of the delta function.

The upper limit of X

\begin{matrix} X_{m a x} = X (q - 1, q) \to \frac{1}{π^{2}} \end{matrix}

(127)

Our distribution (122) is consistent with this upper limit, as the argument

⌊ \frac{1}{π \sqrt{X}} ⌋

becomes zero at

X π^{2} > 1

. It is plotted in Figure 7.

Once we are zooming into the tails of the

p, q

distribution, we also must recall that

\begin{matrix} P (q < y N) = \frac{Φ (⌊ N y ⌋)}{Φ (N)} \to y^{2}; \end{matrix}

(128)

\begin{matrix} f_{y} (y) = \frac{\sum_{q = 2}^{N} δ (y - \frac{q}{N}) φ (q)}{Φ (N)} \end{matrix}

(129)

D. Complex Classical Trajectory in the Path Integral

This classical equation for our path integral reads:

\begin{matrix} α^{''} = - ı κ (C_{Ω}^{'} exp (ı α) + {(C_{Ω}^{'})}^{★} exp (- ı α)); \end{matrix}

(130)

\begin{matrix} κ = \frac{1}{2 π y \sqrt{X} \sqrt{2 ν (t + t_{0})}} \end{matrix}

(131)

The parameter

κ

is distributed according to the above distributions for

y, X

in a small Euler ensemble in the statistical limit.

This complex equation leads to a complex classical solution (instanton). It simplifies for

z = exp (ı α)

\begin{matrix} z^{''} = \frac{{(z^{'})}^{2}}{z} + κ (C_{Ω}^{'} z^{2} + {(C_{Ω}^{'})}^{★}); \end{matrix}

(132)

\begin{matrix} z (0) = z (1) = 1 \end{matrix}

(133)

This equation cannot be analytically solved for arbitrary periodic function

C_{Ω}^{'} (ξ)

The weak and strong coupling expansions by

κ

are straightforward.

At small

κ

\begin{matrix} z (ξ) \to 1 + 2 κ (- A ξ + \int_{0}^{ξ} Re C_{Ω} (ξ^{'}) d ξ^{'}) + O (κ^{2}); \end{matrix}

(134)

\begin{matrix} A = \int_{0}^{1} Re C_{Ω} (ξ^{'}) d ξ^{'} \end{matrix}

(135)

At large

κ

\begin{matrix} z (ξ) \to ı exp (- ı arg C_{Ω}^{'} (ξ)) = ı \frac{|C_{Ω}^{'} (ξ)|}{C_{Ω}^{'} (ξ)} \end{matrix}

(136)

This solution is valid at intermediate

ξ

, not too close to the boundaries

ξ = (0, 1)

. In the region near the boundaries

ξ (1 - ξ) ≪ \frac{1}{\sqrt{κ}}

, the following asymptotic agrees with the classical equation

\begin{matrix} z \to 1 \pm ı ξ \sqrt{2 κ Re C_{Ω}^{'} (0)} + O (ξ^{2}); \end{matrix}

(137)

\begin{matrix} z \to 1 \pm ı (1 - ξ) \sqrt{2 κ Re C_{Ω}^{'} (1)} + O ({(1 - ξ)}^{2}); \end{matrix}

(138)

One can expand in small or large values of

κ

and use the above distributions for

X, y

term by term.

As it was noticed in the previous paper [22], the viscosity

ν = \frac{\tilde{ν}}{N^{2}} \to 0

in our theory. This limit makes

κ \sim N \to \infty

, justifying the strong coupling limit for the Wilson loop solution. In the next section we are considering an important calculable case of the vorticity correlation function, where the full solution in quadratures is available.

VI. Dual Theory of Vorticity Correlation

The simplest observable quantity we can extract from the loop functional is the vorticity correlation function [21], corresponding to the loop C backtracking between two points in space

{\vec{r}}_{1} = 0, {\vec{r}}_{2} = \vec{r}

, (see [22] for details and the justification).

The vorticity operators are inserted at these two points. Let us outline an analytical solution. We shift the time variable by

t + t_{0} \Rightarrow t

to simplify the formulas.

The correlation function reduces to the following average over the big Euler ensemble

E

of our random curves in complex space [22]

\begin{matrix} 〈\vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{r})〉 = \int_{O (3)} \frac{d Ω}{4 t^{2} | O (3) |} \sum_{0 \leq n < m < N} {〈{\vec{ω}}_{m} \cdot {\vec{ω}}_{n} e^{ı \vec{ρ} \cdot \hat{Ω} \cdot ({\vec{S}}_{n, m} - {\vec{S}}_{m, n})}〉}_{E}; \end{matrix}

(139)

\begin{matrix} {\vec{S}}_{n, m} = \frac{\sum_{k = n}^{m - 1} {\vec{F}}_{k}}{m - n (mod N)}; \end{matrix}

(140)

\begin{matrix} {\vec{ω}}_{k} = \{0, 0, \frac{ı σ_{k}}{2} cot (\frac{β}{2})\}; \end{matrix}

(141)

\begin{matrix} \vec{ρ} = \frac{\vec{r}}{2 \sqrt{ν t}}; \end{matrix}

(142)

\begin{matrix} {〈X [σ . . .]〉}_{p, q, r} \equiv \frac{\sum_{E} X [σ . . .] δ [q r - \sum σ]}{\sum_{E} δ [q r - \sum σ]}; \end{matrix}

(143)

Integrating the global rotation matrix

O (3)

is part of the ensemble averaging.

A. Correlation Function and Path Integral

Let us apply our path integral to the expectation value over spins

σ = \pm 1

in the big Euler ensemble, with the distribution of

q, X

established in the previous section. In the continuum limit, we replace summation with integration. We arrive at the following expression for the correlation function:

\begin{matrix} 〈\vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{r})〉 \propto \frac{\sum_{even q < N} \sum_{p; (p | q)} \frac{{cot}^{2} (π p / q)}{{(p / q)}^{2}} \int_{0 < ξ_{1} < ξ_{2} < 1} d ξ_{1} d ξ_{2} \int_{O (3)} d Ω \int [D α] α^{'} (ξ_{1}) α^{'} (ξ_{2}) e^{ı \frac{\vec{r} \cdot \hat{Ω} \cdot Im \vec{V} (ξ_{1}, ξ_{2})}{\sqrt{ν t}}}}{t^{2} Φ (N) | O (3) | \int [D α]}; \end{matrix}

(144)

\begin{matrix} \vec{V} (ξ_{1}, ξ_{2}) = q \sqrt{X} \{ı, 1, 0\} (S (ξ_{1}, ξ_{2}) - S (ξ_{2}, 1 + ξ_{1})); \end{matrix}

(145)

\begin{matrix} S (a, b) = \frac{\int_{a}^{b} d ξ e^{ı α (ξ)}}{b - a}; \end{matrix}

(146)

Here and in the following, we skip all positive constant factors, including powers of N. Ultimately, we restore the correct normalization of the vorticity correlation using its value at

\vec{r} = 0

computed in previous work [22].

The computations significantly simplify in Fourier space.

\begin{matrix} 〈\vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k})〉 = \int d^{3} \vec{r} 〈\vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{r})〉 exp (- ı \vec{k} \cdot \vec{r}) \propto \\ \frac{\sum_{even q < N} \sum_{p; (p | q)} \frac{{cot}^{2} (π p / q)}{{(p / q)}^{2}} \int_{0 < ξ_{1} < ξ_{2} < 1} d ξ_{1} d ξ_{2} \int_{O (3)} d Ω \int [D α] α^{'} (ξ_{1}) α^{'} (ξ_{2}) δ (\frac{\vec{r} \cdot \hat{Ω} \cdot Im \vec{V} (ξ_{1}, ξ_{2})}{\sqrt{ν t}} - \vec{k})}{t^{2} Φ (N) | O (3) | \int [D α]}; \end{matrix}

(147)

The angular integration

\int d Ω

yields

\begin{matrix} \int_{O (3)} d Ω δ (\frac{\hat{Ω} \cdot Im \vec{V} (ξ_{1}, ξ_{2})}{\sqrt{ν t}} - \vec{k}) \propto \frac{\sqrt{ν t}}{{\vec{k}}^{2}} δ (|Im \vec{V} (ξ_{1}, ξ_{2})| - | \vec{k} | \sqrt{ν t}) \end{matrix}

(148)

Now, using the Lagrange multiplier

λ

for this condition, and shifting integration over

λ

to the real axis, we have to minimize effective action

\begin{matrix} A [α, λ] = \frac{π y^{2} X}{2} \int_{ξ_{1}}^{1 + ξ_{1}} {(α^{'})}^{2} + λ y \sqrt{X} |\frac{\int_{ξ_{1}}^{ξ_{2}} d ξ e^{ı α}}{ξ_{2} - ξ_{1}} - \frac{\int_{ξ_{2}}^{1 + ξ_{1}} d ξ e^{ı α}}{1 + ξ_{1} - ξ_{2}}|; \end{matrix}

(149)

\begin{matrix} \partial_{λ} A [α, λ] = | \vec{k} | \sqrt{ν t}; \end{matrix}

(150)

This variational problem reduces to two pendulum equations

\begin{matrix} α_{1}^{''} + \frac{r}{ξ_{2} - ξ_{1}} sin α_{1} = 0; \forall ξ_{1} < ξ < ξ_{2} \end{matrix}

(151)

\begin{matrix} α_{2}^{''} + \frac{r}{ξ_{2} - ξ_{1} - 1} sin α_{2} = 0; \forall ξ_{2} < ξ < 1 + ξ_{1} \end{matrix}

(152)

\begin{matrix} r = \frac{λ}{π y \sqrt{X} I (r)}; \end{matrix}

(153)

\begin{matrix} I (r) = |\frac{\int_{ξ_{1}}^{ξ_{2}} d ξ e^{ı α_{1}}}{ξ_{2} - ξ_{1}} - \frac{\int_{ξ_{2}}^{1 + ξ_{1}} d ξ e^{ı α_{2}}}{1 + ξ_{1} - ξ_{2}}|; \end{matrix}

(154)

The well-known solution is Jacobi amplitude

am (x ∣ u)

\begin{matrix} α_{1} (ξ) = 2 am (\frac{ξ - α_{2}}{2} a_{1} ∣ \frac{r}{a_{1}^{2} (ξ_{2} - ξ_{1})}); \end{matrix}

(155a)

\begin{matrix} α_{2} (ξ) = 2 am (\frac{ξ - α_{2}}{2} a_{2} ∣ - \frac{r}{a_{2}^{2} (1 - ξ_{2} + ξ_{1})}); \end{matrix}

(155b)

The free parameters

a_{1}, a_{2}, α_{1}, α_{2}

satisfy four equations

\begin{matrix} α_{1} (ξ_{2}) = α_{2} (ξ_{2}); \end{matrix}

(156a)

\begin{matrix} α_{1}^{'} (ξ_{2}) = α_{2}^{'} (ξ_{2}); \end{matrix}

(156b)

\begin{matrix} α_{1} (ξ_{1}) = α_{2} (1 + ξ_{1}); \end{matrix}

(156c)

\begin{matrix} α_{1}^{'} (ξ_{1}) = α_{2}^{'} (1 + ξ_{1}); \end{matrix}

(156d)

together with the constraint following from the variation of the Lagrange multiplier

λ

\begin{matrix} I (r) = \frac{| \vec{k} | \sqrt{ν t}}{y \sqrt{X}} \end{matrix}

(157)

B. Turbulent Viscosity and the Local Limit

These five equations, in general, are quite complex, but there is one simplifying property. In the local limit

N \to \infty

, the remaining effective action at the extremum

\begin{matrix} N A [α_{c}, λ_{c}] = \frac{π N y^{2} X}{2} (\int_{ξ_{1}}^{ξ_{2}} d ξ {(α_{1}^{'} (ξ))}^{2} + \int_{ξ_{2}}^{1 + ξ_{1}} d ξ {(α_{2}^{'} (ξ))}^{2}) \end{matrix}

(158)

grows as N, unless both

α_{1} (ξ) \sim α_{2} (ξ) \sim N^{- \frac{1}{2}} \to 0

. In this case, the above constraint can be expanded in

α_{1}, α_{2}

. As we show in [27], the leading constant and linear terms both cancel so that the quadratic term remains

\begin{matrix} \frac{2 | \vec{k} | \sqrt{ν t}}{y \sqrt{X}} = |\int_{ξ_{1}}^{ξ_{2}} \frac{d ξ α_{1}^{2} (ξ)}{ξ_{2} - ξ_{1}} - \int_{ξ_{2}}^{1 + ξ_{1}} \frac{d ξ α_{2}^{2} (ξ)}{1 + ξ_{1} - ξ_{2}}| \sim \frac{1}{N} \end{matrix}

(159)

This estimate then requires vanishing viscosity in the local limit, at fixed turbulent viscosity

\begin{matrix} \tilde{ν} = ν N^{2} \to const . \end{matrix}

(160)

This phenomenon of renormalization of viscosity by a factor of

N^{2}

was already observed in our first paper [22]. Our Euler ensemble in the local limit

N \to \infty

can only solve the inviscid limit of the Navier-Stokes decaying turbulence, with finite

\tilde{ν}

acting as a turbulent viscosity.

The desired anomalous dissipation phenomenon takes place in this limit of our theory.

Returning to the elliptic function solution, we rewrite it in the linearized case at

a_{1} \sim a_{2} \to 0

. This linearization is equivalent to replacing

sin (α) \to α

in the differential equation and studying the resulting linear ODE as a boundary problem. We choose different parametrizations in this linear case

\begin{matrix} α_{1} (ξ) = a (cos (K_{1} (ξ - ξ_{2})) + \frac{b}{K_{1}} sin (K_{1} (ξ - ξ_{2}))); \end{matrix}

(161)

\begin{matrix} α_{2} (ξ) = a (cos (K_{2} (ξ - ξ_{2})) + \frac{b}{K_{2}} sin (K_{2} (ξ - ξ_{2}))); \end{matrix}

(162)

\begin{matrix} K_{1} = \sqrt{\frac{r}{Δ}}; \end{matrix}

(163)

\begin{matrix} K_{2} = \sqrt{\frac{r}{Δ - 1}}; \end{matrix}

(164)

\begin{matrix} Δ = ξ_{2} - ξ_{1}; \end{matrix}

(165)

In the physical region

0 < Δ < 1, r < 0

K_{2}

is real, and

K_{1}

imaginary, but the solution stays real. The matching conditions at

α_{1} (ξ_{2}) = α_{2} (ξ_{2}), α_{1}^{'} (ξ_{2}) = α_{2}^{'} (ξ_{2})

are identically satisfied with this Anzatz. The derivative match

α_{1}^{'} (ξ_{1}) = α_{2}^{'} (1 + ξ_{1})

can be solved exactly for b

\begin{matrix} b = \frac{\sqrt{\frac{r}{Δ - 1}} sin ((1 - Δ) \sqrt{\frac{r}{Δ - 1}}) + \sqrt{\frac{r}{Δ}} sin (Δ \sqrt{\frac{r}{Δ}})}{cos ((Δ - 1) \sqrt{\frac{r}{Δ - 1}}) - cos (Δ \sqrt{\frac{r}{Δ}})} \end{matrix}

(166)

The remaining matching condition

α_{1} (ξ_{1}) = α_{2} (1 + ξ_{1})

reduces to the root of the function

\begin{matrix} g (r, Δ) = \frac{(2 Δ - 1) sin (\sqrt{(Δ - 1) r}) sin (\sqrt{Δ r})}{\sqrt{(Δ - 1) Δ}} + 2 cos (\sqrt{(Δ - 1) r}) cos (\sqrt{Δ r}) - 2 \end{matrix}

(167)

This function has multiple roots, but we are looking for the real root

r_{0} (Δ)

with minimal value of the action at given

Δ

\begin{matrix} A_{c} (r, Δ) = \int_{ξ_{1}}^{ξ_{2}} d ξ α_{1}^{'} {(ξ)}^{2} + \int_{ξ_{2}}^{1 + ξ_{1}} d ξ α_{2}^{'} {(ξ)}^{2} \end{matrix}

(168)

This integral is elementary, but the expression is too long to be presented here. It can be found in the Mathematica^® notebook [27], where it is used to select the roots

r_{0} (Δ)

g (r, Δ)

, minimizing

A_{c} (r_{0} (Δ), Δ)

for a given value of

Δ \in (0, 1)

This lowest action root is plotted in Figure 8. The corresponding value of minimal action

L (Δ) = A_{c} (r_{0} (Δ), Δ)

is plotted in Figure 12.

There are phase transitions at

\begin{matrix} Δ = \{Δ_{1}, Δ_{2}, Δ_{3}\}; \end{matrix}

(169)

\begin{matrix} Δ_{1} = 0.157143; \end{matrix}

(170)

\begin{matrix} Δ_{2} = 0.43015; \end{matrix}

(171)

\begin{matrix} Δ_{3} = \frac{1}{2}; \end{matrix}

(172)

These branch points in

Δ

correspond to the switch of the lowest action solution. At small positive

Δ - \frac{1}{2}

\begin{matrix} r_{0} (Δ) \to - 48 (Δ - \frac{1}{2}) - \frac{1536 {(Δ - \frac{1}{2})}^{3}}{7} + \dots; \end{matrix}

(173)

\begin{matrix} A_{0} (r_{0} (Δ), Δ) \to 256 {(Δ - \frac{1}{2})}^{2} + \dots \end{matrix}

(174)

Δ \to 1

all solutions go to

- \infty

\begin{matrix} r_{n} (Δ) \to - \frac{6}{1 - Δ} + O (1) \end{matrix}

(175)

This behavior matches numerical computations in Mathematica^® [27].

The constraint (159) is also reduced to elementary functions, too lengthy to quote here (see [27]).

This constraint yields the quadratic relation for the last unknown parameter a in our solution

\begin{matrix} a^{2} = \frac{| \vec{k} | \sqrt{ν t}}{S (Δ) y \sqrt{X}}; \end{matrix}

(176)

with universal function

S (Δ)

presented in [27] and shown in Figure 9).

The resulting integral (up to the pre-exponential factor Q) is equal to

\begin{matrix} \int [D α] α^{'} (ξ_{1}) α^{'} (ξ_{2}) δ (|Im \vec{V} (ξ_{1}, ξ_{2})| - | \vec{k} | \sqrt{ν t}) \propto \\ Q 〈α^{'} (ξ_{1}) α^{'} (ξ_{2})〉 exp (- y \sqrt{X} | \vec{k} | \sqrt{\tilde{ν} t} \frac{π L (Δ)}{2 | S (Δ) |}) \end{matrix}

(177)

The factor

〈α^{'} (ξ_{1}) α^{'} (ξ_{2})〉

contains two terms :

\begin{matrix} 〈α^{'} (ξ_{1}) α^{'} (ξ_{2})〉 = α_{1}^{'} (ξ_{1}) α_{1}^{'} (ξ_{2}) + 〈δ α^{'} (ξ_{1}) δ α^{'} (ξ_{2})〉 \end{matrix}

(178)

The first term is the contribution of the classical solution we have just found, and the second term comes from Gaussian fluctuations

δ α (ξ)

around this solution.

The classical term is calculable (see [27])

\begin{matrix} α_{1}^{'} (ξ_{1}) α_{1}^{'} (ξ_{2}) = \frac{| \vec{k} | \sqrt{\tilde{ν} t}}{N y \sqrt{X}} \frac{J (Δ)}{| S (Δ) |}; \end{matrix}

(179)

\begin{matrix} J (Δ) = {\frac{r A (r, Δ) B (r, Δ)}{Δ (Δ - 1) C {(r, Δ)}^{2}}|}_{r = r_{0} (Δ)}; \end{matrix}

(180)

\begin{matrix} A (r, Δ) = Δ sin (\sqrt{(Δ - 1) r}) - \sqrt{(Δ - 1) Δ} sin (\sqrt{Δ r}); \\ B (r, Δ) = Δ sin (\sqrt{(Δ - 1) r}) cos (\sqrt{Δ r}) - \end{matrix}

(181)

\begin{matrix} \sqrt{(Δ - 1) Δ} sin (\sqrt{Δ r}) cos (\sqrt{(Δ - 1) r}); \end{matrix}

(182)

\begin{matrix} C (r, Δ) = cos (\sqrt{(Δ - 1) r}) - cos (\sqrt{Δ r}); \end{matrix}

(183)

(see Figure 11).

The fluctuation term

〈δ α^{'} (ξ_{1}) δ α^{'} (ξ_{2})〉

is also proportional to

1 / N

, therefore we must keep this term as well.

As for the pre-exponential factor Q in the saddle point integral, it is given by the functional determinant of the operator

\hat{L}

corresponding to linearized effective action (168) in the vicinity of the saddle point

λ_{c}, α_{c} (ξ)

\begin{matrix} A [α_{c} + δ α, λ_{c} + δ λ] \to A [α_{c}, λ_{c}] + \frac{1}{2} tr (V^{†} ∣ \hat{L} ∣ V); \end{matrix}

(184)

\begin{matrix} V = \{δ λ, δ α\}; \end{matrix}

(185)

\begin{matrix} Q (Δ, τ) = exp {(- \frac{1}{2} tr (\frac{log \hat{L}}{{\hat{L}}^{α}}))}_{α \to 0}; \end{matrix}

(186)

\begin{matrix} τ = y \sqrt{X}; \end{matrix}

(187)

The fluctuation correction reduces to the inverse operator

\hat{L}

, which we compute in the next section.

Figure 12. Log plot of

L (Δ) = A_{c} (r_{0} (Δ), Δ)

Figure 12. Log plot of

L (Δ) = A_{c} (r_{0} (Δ), Δ)

Now, we can reduce multiple sum/integral in (147) to the following

\begin{matrix} 〈\vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k})〉 = \frac{{\tilde{ν}}^{\frac{3}{2}} H (k \sqrt{\tilde{ν} t})}{\sqrt{t}}; \end{matrix}

(188)

\begin{matrix} H (κ) = \frac{1}{Z} \sum_{n = 1}^{\infty} φ (n) \int_{0}^{1 / n} d τ (τ^{5} / n^{5} - τ^{10}) \int_{0}^{1} d Δ (1 - Δ) G (Δ, τ, κ); \end{matrix}

(189)

\begin{matrix} G (Δ, τ, κ) = Q (Δ, τ) (\frac{κ}{N τ} \frac{J (Δ)}{| S (Δ) |} + 〈δ α^{'} (ξ_{1}) δ α^{'} (ξ_{2})〉) exp (- \frac{τ κ L (Δ)}{2 π | S (Δ) |}) \end{matrix}

(180)

where

Z

is the normalization constant to be determined later.

C. Functional Determinant in the Path Integral

As we have discussed in the previous section, in the limit

a \to 0

the classical solution

α_{1, 2} (ξ) \propto a \to 0

This observation simplifies the linearized theory corresponding to this quadratic form

〈V ∣ \hat{L} ∣ V〉

. First, integrate the fluctuations

δ λ

λ

around the saddle point solution.

The Lagrange multiplier at the saddle point vanishes, as we show in [27]

\begin{matrix} λ_{0} = τ r_{0} (Δ) I (r_{0} (Δ)) = 0 \end{matrix}

(191)

The quadratic term comes from the first derivatives

I_{λ} = \partial_{λ} I, I_{r} = \partial_{r} I, λ_{r} = \partial_{r} λ

, which can be simplified by switching to

λ (r) = τ t I (r)

\begin{matrix} A_{λ λ} = τ I_{λ} = \frac{τ I_{r}}{λ_{r}} = \frac{τ I_{r}}{τ r I_{r}} = \frac{1}{r} \end{matrix}

(192)

The bilinear term

λ δ α

also simplifies

\begin{matrix} A_{α λ} (δ α) = ı τ δ λ F [δ α]; \end{matrix}

(193)

\begin{matrix} F [δ α] = \int_{ξ_{1}}^{ξ_{2}} \frac{d ξ δ α (ξ)}{ξ_{2} - ξ_{1}} - \int_{ξ_{2}}^{1 + ξ_{1}} \frac{d ξ δ α (ξ)}{1 + ξ_{1} - ξ_{2}} \end{matrix}

(194)

We can integrate out

λ

, producing the extra pre-exponential factor

Q_{λ} = \sqrt{| r_{0} (Δ) |} / \sqrt{N}

The bilinear term in the exponential after

λ

integration leads to the following effective quadratic Action for

δ α

\begin{matrix} A_{e f f} [δ α] = \frac{τ^{2}}{2} \int_{ξ_{1}}^{1 + ξ_{1}} δ α^{' 2} + \frac{r_{0} (Δ) τ^{2}}{2} F {[δ α]}^{2}; \end{matrix}

(195)

There is a zero-mode

δ α (ξ) = const

, related to translational invariance of

A_{e f f} [δ α]

. Naturally, this zero-mode must be eliminated from the spectrum when we compute the functional determinant and the resolvent below.

After discarding the zero-mode, this effective action becomes a positive definite functional of

δ α

only in the region of

Δ

where

r_{0} (Δ) > 0

, i.e., for

Δ_{1} < Δ < Δ_{2}

As we shall see below, the spectrum of fluctuations is positive only in this region. Therefore, we restrict our integration to this region.

The

{(δ α)}^{2}

term corresponds to the linear eigenvalue equation with

f_{1, 2} = δ α_{1, 2}

\begin{matrix} f_{1, 2}^{''} (ξ) - μ_{1, 2} F [f] = - ω f_{1, 2} (ξ); \end{matrix}

(196)

\begin{matrix} μ_{1, 2} = \{\frac{r}{Δ}, \frac{r}{Δ - 1}\}; \end{matrix}

(197)

\begin{matrix} ϵ = ω τ^{2}; \end{matrix}

(198)

\begin{matrix} r = r_{0} (Δ); \end{matrix}

(199)

The solution matching with first derivative at

ξ = ξ_{2}, ξ = (ξ_{1}, 1 + ξ_{1})

is built the same way as in (166). Equations for

f_{1, 2}

being linear homogeneous, we can fix the normalization as

F [f] = 1

\begin{matrix} f_{1} (x) = a sin (\sqrt{ω} (ξ_{-} ξ_{2})) + B_{1} cos (\sqrt{ω} (ξ_{-} ξ_{2})) + \frac{r}{ω Δ}; \end{matrix}

(200)

\begin{matrix} f_{2} (ξ) = a sin (\sqrt{ω} (ξ_{-} ξ_{2})) + B_{2} cos (\sqrt{ω} (ξ_{-} ξ_{2})) + \frac{r}{ω (Δ - 1)}; \end{matrix}

(201)

The spectrum

ω = ω_{n}

is defined by the transcendental equation (the discriminant of this linear system of equations), which we found in [27]

\begin{matrix} f (ω_{n}, Δ) = 0; \\ f (ω, Δ) = (Δ - 1) Δ \sqrt{ω} sin (\frac{\sqrt{ω}}{2}) ((Δ - 1) Δ ω + r) + \end{matrix}

(202)

\begin{matrix} r cos (\frac{1}{2} (1 - 2 Δ) \sqrt{ω}) - r cos (\frac{\sqrt{ω}}{2}); \end{matrix}

(203)

\begin{matrix} r = r_{0} (Δ); \end{matrix}

(204)

The spectrum is positive in the interval

Δ_{1} < Δ < Δ_{2}

where

r_{0} (Δ) > 0

so that the solution for

ω_{n} (Δ)

is stable. In the following, we only select the stable region with positive

r_{0} (Δ)

Figure 13. The first levels of the spectrum satisfying equation

f (ω, Δ) = 0

. The colored lines correspond to four phases. Red:

0 < Δ < Δ_{1}

, Green:

Δ_{1} < Δ < Δ_{2}

, Blue:

Δ_{2} < Δ < \frac{1}{2}

, Brown:

\frac{1}{2} < Δ < 1

. The green zone is left as stable, and others are eliminated because

r_{0} (Δ) < 0

in these zones. Naturally, we eliminate the zero-mode

ϵ_{0} = 0

corresponding to translational invariance of the effective Action.

Figure 13. The first levels of the spectrum satisfying equation

f (ω, Δ) = 0

. The colored lines correspond to four phases. Red:

0 < Δ < Δ_{1}

, Green:

Δ_{1} < Δ < Δ_{2}

, Blue:

Δ_{2} < Δ < \frac{1}{2}

, Brown:

\frac{1}{2} < Δ < 1

. The green zone is left as stable, and others are eliminated because

r_{0} (Δ) < 0

in these zones. Naturally, we eliminate the zero-mode

ϵ_{0} = 0

corresponding to translational invariance of the effective Action.

The functional determinant, resulting from the WKB approximation to the

α

path integral, would be related to the infinite product of positive eigenvalues

ϵ_{n} = τ^{2} ω_{n}

, which can be written using a contour integral

\begin{matrix} Q_{α} (Δ, τ) = \prod_{ω_{n} > 0} {(τ^{2} ω_{n})}^{- \frac{1}{2}} = {exp (\frac{1}{2} \partial_{α} Im \oint_{Γ} \frac{f^{'} (ω)}{f (ω)} \frac{d ω}{2 π {(ω τ^{2})}^{α}})|}_{α \to 0}; \end{matrix}

(205)

and the integration contour

Γ

encircles anticlockwise the positive real poles of the meromorphic function

f^{'} (ω) / f (ω)

. The integral converges at

α > \frac{1}{2}

and should be analytically continued to

α = 0

For this purpose, let us introduce another function

\begin{matrix} Φ (ω) = \frac{f (ω)}{cos (\sqrt{ω} / 2)} ω^{- 3 / 2} \end{matrix}

(206)

We show in [27] that at large

ω = ı y

this function reaches finite limits

\begin{matrix} Φ (ı y) \to ı sign y {(Δ - 1)}^{2} Δ^{2} + \frac{(Δ - 1) Δ r}{| y |} \end{matrix}

(207)

The logarithmic derivative of the original function differs from

\frac{Φ^{'} (ω)}{Φ (ω)}

by the following meromorphic function

\begin{matrix} \frac{f^{'} (ω)}{f (ω)} - \frac{Φ^{'} (ω)}{Φ (ω)} = - \frac{tan (\frac{\sqrt{ω}}{2})}{4 \sqrt{ω}} + \frac{3}{2 ω} \end{matrix}

(208)

This difference produces a calculable contribution to our integral. By summing residues of the poles of the tangent, we get

\begin{matrix} \oint_{Γ} (- \frac{tan (\frac{\sqrt{ω}}{2})}{4 \sqrt{ω}} + \frac{3}{2 ω}) \frac{d ω}{2 π {(ω τ^{2})}^{α}} = ı (1 - 2^{2 α}) {(2 π τ)}^{- 2 α} ζ (2 α) \end{matrix}

(209)

The derivative at

α = 0

yields a constant

\begin{matrix} \frac{1}{2} \partial_{α} Im ı (1 - 2^{2 α}) {(2 π τ)}^{- 2 α} ζ (2 α) \to \frac{log (2)}{2} \end{matrix}

(210)

leading to an irrelevant renormalization of

Q (α, τ)

by a factor

\sqrt{2}

The remaining integral with

f (ω) \Rightarrow Φ (ω)

already converges at

Re α > - 1

, so that we can set

α = 0

and rotate the integration contour

Γ

parallel to the imaginary axis at

Re Γ = ϵ > 0

\begin{matrix} Q_{α} (Δ, τ) = \\ exp (\frac{1}{2} Im \int_{ϵ - ı \infty}^{ϵ + ı \infty} \frac{Φ^{'} (ω)}{Φ (ω)} \frac{log (ω τ^{2}) d ω}{2 π}) \end{matrix}

(211)

The remarkable property of this functional determinant is the factorization of the

τ

dependence

\begin{matrix} Q_{α} (Δ, τ) = τ^{μ (Δ)} Q_{α} (Δ, 1); \end{matrix}

(212)

\begin{matrix} μ (Δ) = Im \int_{ϵ - ı \infty}^{ϵ + ı \infty} \frac{Φ^{'} (ω)}{Φ (ω)} \frac{d ω}{2 π} \end{matrix}

(213)

The index

μ (Δ)

has a topological origin and can be computed analytically.

\begin{matrix} μ (Δ) = \frac{arg Φ (ı \infty) - arg Φ (- ı \infty)}{2 π} = \frac{1}{2} \end{matrix}

(214)

Our result for the correlation function is given by (188) with

\begin{matrix} Q (Δ, τ) = Q_{α} (Δ, 1) τ^{\frac{1}{2}} \sqrt{r_{0} (Δ)} \end{matrix}

(215)

and

Q_{α} (Δ, 1)

given by (211). All the constant factors we have omitted here are absorbed by the normalization factor

Z

, which we determine at the end of the next section.

D. The Fluctuation Term in $α^{'} (ξ_{1}) α^{'} (ξ_{2})$

The last missing term is the fluctuation contribution to

α^{'} (ξ_{1}) α^{'} (ξ_{2})

. In the Gaussian approximation, valid at

N \to \infty

, this term equals

\begin{matrix} 〈δ α^{'} (ξ_{1}) δ α^{'} (ξ_{2})〉 = \frac{- 1}{N τ^{2}} {[\partial_{ξ} \partial_{ξ^{'}} G (ξ, ξ^{'})]}_{ξ = ξ_{1}, ξ^{'} = ξ_{2}} \end{matrix}

(216)

where

G (ξ, ξ^{'})

is a resolvent for the effective quadratic Action (195). This resolvent satisfies the equation

\begin{matrix} \partial_{ξ}^{2} G (ξ, ξ^{'}) - μ (ξ) F [G] = δ (ξ - ξ^{'}); \end{matrix}

(217)

\begin{matrix} G (ξ_{1}, ξ^{'}) = G (ξ_{1} + 1, ξ^{'}) = 0; \end{matrix}

(218)

\begin{matrix} F [G] = \int_{ξ_{1}}^{ξ_{2}} \frac{d ξ G (ξ, ξ^{'})}{ξ_{2} - ξ_{1}} - \int_{ξ_{2}}^{1 + ξ_{1}} \frac{d ξ G (ξ, ξ^{'})}{1 + ξ_{1} - ξ_{2}}; \end{matrix}

(219)

\begin{matrix} μ (ξ) = \{\begin{matrix} \frac{r}{Δ} & if ξ_{1} \leq ξ < ξ_{2} \\ \frac{r}{Δ - 1} & if ξ_{2} \leq ξ < 1 + ξ_{1} \end{matrix} \end{matrix}

(220)

The solution of this equation, matching with the first derivative at

ξ = ξ_{2}

\begin{matrix} G (ξ, ξ^{'}) = \{\begin{matrix} A + \frac{| ξ - ξ^{'} |}{2} + B (ξ - ξ_{2}) + \frac{F [G] r {(ξ - ξ_{2})}^{2}}{2 Δ} & if ξ_{1} \leq ξ < ξ_{2} \\ A + \frac{| ξ - ξ^{'} |}{2} + B (ξ - ξ_{2}) + \frac{F [G] r {(ξ - ξ_{2})}^{2}}{2 (Δ - 1)} & if ξ_{2} \leq ξ < 1 + ξ_{1} \end{matrix} \end{matrix}

(221)

The linear functional

F [G]

on this solution becomes a linear function of these unknown parameters

A, B

. Two boundary conditions

G (ξ_{1}, ξ^{'}) = G (ξ_{1} + 1, ξ^{'}) = 0

fix these parameters as functions of

ξ_{1}, ξ_{2}, ξ^{'}

The result derived in [27] is too lengthy to present here. The desired quantity (216) is quite simple

\begin{matrix} \frac{1}{N τ^{2}} {[\partial_{ξ} \partial_{ξ^{'}} G (ξ, ξ^{'})]}_{ξ = ξ_{1}, ξ^{'} = ξ_{2}} = \frac{2 (r - 6)}{(r + 12) N τ^{2}} \end{matrix}

(222)

Finally, we get the following correlation (188) (absorbing the constant factors in

Z

)

\begin{matrix} 〈\vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k})〉 = \frac{{\tilde{ν}}^{\frac{3}{2}} H (k \sqrt{\tilde{ν} t})}{\sqrt{t}}; \end{matrix}

(223a)

\begin{matrix} H (κ) = \int_{Δ_{1}}^{Δ_{2}} d Δ (1 - Δ) \sum_{n = 1}^{\infty} φ (n) \int_{0}^{1 / n} \frac{d τ}{τ^{\frac{5}{2}}} (τ^{5} / n^{5} - τ^{10}) G (Δ, τ κ); \end{matrix}

(223b)

\begin{matrix} G (Δ, x) = \frac{Q_{α} (Δ, 1) \sqrt{r_{0} (Δ)}}{Z} (x \frac{J (Δ)}{S (Δ)} + \frac{2 (r_{0} (Δ) - 6)}{(12 + r_{0} (Δ))}) exp (- \frac{x L (Δ)}{2 π S (Δ)}) \end{matrix}

(223c)

VII. The Decaying Energy in Finite System

The vorticity correlation in Fourier space doubles as an energy spectrum

\begin{matrix} E (k, t) = 4 π {\vec{k}}^{2} 〈\vec{v} \cdot {\vec{v}}_{k}〉 = 4 π 〈\vec{ω} \cdot {\vec{ω}}_{k}〉 \end{matrix}

(224)

The energy spectrum in a finite system with size L is bounded from below. At low

| \vec{k} | \leq π / L

, the spectrum is no longer related to the turbulence but is given by the energy pumping by external forces at the boundaries.

This energy pumping [35] takes place at

t < t_{0}

, after which the pumping stops. At this moment, the energy spectrum is growing with wavevector by one of two possible laws (with P being the net momentum of the fluid and M being the rotation moment)

\begin{matrix} \{\begin{matrix} E (k, t_{0}) \propto P k^{2} \\ E (k, t_{0}) \propto M k^{4} \end{matrix} \end{matrix}

(225)

t > t_{0}

, without the forcing, the pumped energy dissipates at large k corresponding to smaller spatial scales of the hierarchy of vortex structures of all scales, ending with dissipative scales, or wavevectors

k ≫ π / L

. After sufficient time, the universal regime kicks in, corresponding to the decaying turbulence. It is implied that a large amount of energy was pumped in, so it takes a long time to reach this decaying regime, corresponding to some fixed trajectory.

Our solution presumably applies to the second regime with

k^{4}

spectrum at small k and some universal decay at large k, reflecting these distributed vortex structures. This issue requires more investigation.

The decaying energy, given by the part of the spectrum

k > k_{0} \sim 1 / L

, has the following form

\begin{matrix} E (t) = \int_{k_{0}}^{\infty} d k E (k, t) = \frac{4 π \tilde{ν}}{t} \int_{k_{0} \sqrt{\tilde{ν} t}}^{\infty} H (κ) d κ \end{matrix}

(226)

On top of the trivial decrease

\frac{\tilde{ν}}{t}

, as prescribed by dimensional counting in an infinite system, there is some extra decrease related to the increase of the lower limit.

The energy in our theory does not have a finite statistical limit as the integral in (226) diverges at the lower limit when

k_{0} \to 0

. Thus, we compute the energy as

\begin{matrix} E (t) = \int_{t}^{\infty} E (t^{'}) d t^{'} \end{matrix}

(227)

This energy dissipation rate

E (t^{'})

is calculable

\begin{matrix} E (t) = 4 π ν \int d k k^{2} 〈\vec{ω} (\vec{0}) \cdot \vec{ω} (\vec{k})〉 = 4 π \frac{\tilde{ν}}{t^{2}} \int_{k_{0} \sqrt{\tilde{ν} t}}^{\infty} κ^{2} H (κ) d κ \end{matrix}

(228)

In our theory, this integral has a finite limit in an infinite system (

k_{0} = 0

This limit was computed in [22] in a slightly different grand canonical ensemble, where N was fluctuating with the weight

exp (- μ N), μ \to 0

With our current ensemble of fixed even

N \to \infty

the results of [22] read:

\begin{matrix} E_{\infty} (t) = \frac{\tilde{ν}}{t^{2}} \frac{π^{2}}{864 ζ (3)}; \end{matrix}

(229)

In our present theory, the same quantity is given by the above integral at

k_{0} = 0

\begin{matrix} E_{\infty} (t) = 4 π \frac{\tilde{ν}}{t^{2}} \int_{0}^{\infty} κ^{2} H (κ) d κ \end{matrix}

(230)

Comparing these two expressions, we get the normalization of

H (κ)

\begin{matrix} 4 π \int_{0}^{\infty} κ^{2} H (κ) d κ = \frac{π^{2}}{864 ζ (3)} \end{matrix}

(231)

The integral on the left can be further reduced [27] to the following normalization condition:

\begin{matrix} Z = \frac{276480 ζ (3) ζ (\frac{15}{2})}{119 ζ (\frac{17}{2})} \int_{Δ_{1}}^{Δ_{2}} d Δ \frac{(1 - Δ) \sqrt{r_{0} (Δ)} S (Δ) Q_{α} (Δ, 1) (π J (Δ) (r_{0} (Δ) + 12) + L (Δ) (r_{0} (Δ) - 6))}{L {(Δ)}^{2} (r_{0} (Δ) + 12)} \end{matrix}

(232)

This normalization constant

Z

can be used in equation (226) for the energy decay in a finite system. All the functions of

Δ

were defined above.

As for the energy spectrum, this is not an independent function in our theory. Comparing the two expressions (224) and (228), we arrive at the following relation

\begin{matrix} t^{2} E (t) = 4 π \tilde{ν} F (k_{0} \sqrt{\tilde{ν} t}); \end{matrix}

(233a)

\begin{matrix} F (κ) = \int_{κ}^{\infty} H (x) x^{2} d x, \end{matrix}

(233b)

\begin{matrix} \sqrt{t} E (k, t) = 4 π \tilde{ν} \sqrt{\tilde{ν}} H (k \sqrt{\tilde{ν} t}); \end{matrix}

(233c)

Both the energy dissipation and the energy spectrum are related to the same function

H (κ)

, but the energy spectrum related to this function at large argument

κ = k \sqrt{\tilde{ν} t}

, whereas the energy dissipation is related to integral of this function from small argument

κ = k_{0} \sqrt{\tilde{ν} t}

to infinity.

We computed this function in the Appendix A and numerically integrated it to obtain

F (κ)

. Asymptotically, at large

κ = k \sqrt{\tilde{ν} t}

\begin{matrix} E (k, t) \propto \frac{\tilde{ν} \sqrt{\tilde{ν}}}{\sqrt{t}} {(k \sqrt{\tilde{ν} t})}^{- \frac{7}{2}}; \end{matrix}

(234)

\begin{matrix} E (t) \propto t^{- \frac{9}{4}} \end{matrix}

(235)

We computed the effective decay indexes

\begin{matrix} n (t) = - t \partial_{t} log E - 1; \end{matrix}

(236)

\begin{matrix} s (κ) = - κ \partial_{κ} log H (κ) \end{matrix}

(237)

numerically in the Appendix A, using Mathematica^®. The accuracy is just 4-5 digits, but it can be easily improved by taking more CPU time once experimental data gets more precise. These curves are universal, and they change the regime before approaching their limits from below

n (\infty) = \frac{5}{4}, s (\infty) = \frac{7}{2}

. These regime changes are due to the quantum effects (complex zeros of the zeta function contributing to the energy spectrum’s Mellin transform, as shown in Appendix A).

Current data for the energy decay index n [35] are inconsistent. There are large discrepancies between various experiments, presumably explained by attempts to fit the energy decay to a single power. In our theory, the energy decay is an oscillating function on the log-log scale.

The best approach to experimental data is to present it as an effective curve

n (t)

, as it was done in [35] and presented here in Figure 3 with permission from Sreeni.

Some of these curves match our theoretical prediction in Figure 5, but comparing the theory with inconsistent experimental data is premature.

Our universal curves were computed directly from the analytic solution of the loop equation in the turbulent limit without any fitting parameters. It will be very interesting to compare these curves with more precise experiments (real or numerical).

VIII. The Spectrum of Decay Indexes in the Energy Spectrum

In addition to numerically computing the energy spectrum and plotting effective decay indexes, we can go one step further in the mathematical analysis of the concept of the scaling laws in turbulence.

The notion of power decay can be formulated in precise mathematical terms by studying the Mellin transform.

Namely, the spectrum of decay indexes is given by the positions of poles of the Mellin transform in the negative semiplane. We change the sign of p in conventional definition to better describe our functions.

In our case, there is a universal function

H (κ)

in (223). The Mellin transform is the following transformation ( assuming our function decreases at infinity)

\begin{matrix} H (κ) = \int_{+ 0 - ı \infty}^{+ 0 + ı \infty} \frac{d p}{2 π ı} h (p) κ^{p}; \end{matrix}

(238)

\begin{matrix} h (p) = \int_{0}^{\infty} \frac{d κ}{κ} κ^{- p} H (κ) \end{matrix}

(239)

The pure power law of decay would correspond to the Mellin transform

h (p)

having a single pole in the left semi-plane. The position

p = - a

of this pole becomes an index of the power law.

\begin{matrix} H (κ) \propto κ^{- a} \end{matrix}

(240)

The next level of complexity would be a function, depending on an extra parameter n, such that the Mellin transform has a simple pole, moving with this parameter

p = - a (n)

. This will produce multifractal scaling laws

\begin{matrix} κ^{- a (n)} \end{matrix}

(241)

Our energy spectrum has a more complex singularity structure in its Mellin transform

\begin{matrix} h (p) = \frac{f (p) ζ (p + \frac{15}{2}) Γ (- p)}{(2 p + 7) (2 p + 17) ζ (p + \frac{17}{2})}; \end{matrix}

(242)

\begin{matrix} f (p) = 20 \int_{Δ_{1}}^{Δ_{2}} d Δ (1 - Δ) C^{p - 1} (A C - B p) \end{matrix}

(243)

where

A, B, C

are some smooth positive functions of

Δ

varying in finite limits (see Appendix A). Given these properties, it is simple to prove that the Taylor series of

f (p)

at the origin converges as the expansion coefficients decrease as a factorial of the expansion order.

This makes this function an entire function without any finite singularities. The values of

C (Δ)

are bounded by two positive limits

\begin{matrix} 0.0541984 < C < 0.0630755, \forall Δ_{1} < Δ < Δ_{2} \end{matrix}

(244)

Therefore, the entire function

f (p)

decreases as

e^{- 2.76342 p}

in the right semiplane and grows as

e^{- 2.9151 p}

in the left semiplane.

The singularities of the Mellin transform

h (p)

are given by the following table of simple poles

\begin{matrix} |\begin{matrix} p o s i t i o n & m u l t i p l i c i t y \\ - \frac{17}{2} & 1 \\ - \frac{13}{2} & 1 \\ - \frac{7}{2} & 1 \\ \frac{1}{2} (- 17 - 4 c_{1}) if c_{1} \in Z \land c_{1} \geq 1 & 1 \\ 2 c_{1} if c_{1} \in Z \land c_{1} \geq 0 & 1 \\ 1 + 2 c_{1} if c_{1} \in Z \land c_{1} \geq 0 & 1 \\ \frac{1}{2} (- 17 + 2 t_{c_{1}}) if c_{1} \in Z & 1 ? \end{matrix}| \end{matrix}

(245)

Here

\pm t_{n}

are imaginary parts of zeros of the

ζ

function on the critical line

z = \frac{1}{2} + ı t

. About

10^{13}

zeros are already known, though the Riemann hypothesis (that there are no other complex zeros) is still unproven. These poles are shown in Figure 14

The real negative poles yield decaying power terms, but the infinite series of complex poles at Riemann zeros adds the oscillations in log scale

\begin{matrix} | W_{n} | κ^{- 8} cos (t_{n} log κ + arg W_{n}) \end{matrix}

(246)

These slow oscillations are visible as regime change in the log-log plots of effective indexes in Figure 5, Figure 4.

IX. Discussion

Some parts of the following discussion were already published in my paper [22]. We use them with significant modifications, reflecting our increased understanding of this theory.

A. Stochastic Solution of the Navier-Stokes Equation and Quantum Hypothesis

Statistical "analysis of circulating or turbulent fluids" was defined by Feynman [7] as the last unsolved problem of classical physics. We are pursuing this problem by finding a stochastic solution of the Navier-Stokes equation, asymptotically covering a certain manifold at infinite time. Our singularities arise in correlation functions after averaging over this manifold in the statistical limit when the manifold’s dimension goes to infinity. However, the origin of singularities is an infinite time needed to cover this fixed manifold uniformly.

We find this manifold (Euler ensemble) by solving the loop equation (a subset of the Hopf functional equation for generating functional for the velocity field probability). None of the particular solutions in this manifold has a finite-time blow-up. We encounter a singularity of the fixed point of the loop equation rather than its solution for a finite time.

We take the most natural invariant measure from the point of the number theory: each element of the Euler ensemble enters with equal weight. We call it our ergodic hypothesis. This hypothesis is unnecessary to solve the loop (i.e., Hopf) equation, as any linear superposition of the found solutions would satisfy the loop equation. The singularities of our Euler ensemble at

N \to \infty

remain in local variables such as enstrophy and its PDF, indicating singularities of the Navier-Stokes equation in the inviscid limit.

With this particular invariant measure, the Euler ensemble is singled out by its equivalence to the trace of the evolution operator of a quantum system: the Fermi particles on a ring in the external field related to random fractions of

π

. Every distinct state in the discrete quantum system enters a quantum trace with unit weight. In our case, this is equivalent to counting every element in the Euler ensemble with equal weight.

In other words, our ergodic hypothesis is equivalent to the requirement of exact equivalence of the loop functional to the quantum trace of the one-dimensional ring of Fermi particles. This quantum representation opened the way to the analytic solution in the turbulent limit, equivalent to the quasiclassical limit of this Fermi system. The viscosity acts as Planck’s constant, and the inverse Reynolds number acts as a weak coupling constant of the dual theory.

There is an unexpected gift of dimensional reduction. This quantum system is one-dimensional despite its equivalence to the three-dimensional decaying turbulence. Such equivalence would be impossible given the local relation between the dynamical variables of these two systems; however, with our nonlocal relation through the loop average, the duality becomes possible.

A mathematical transformation relates the quantum amplitude to the turbulent probability, but this is not an analytic continuation to imaginary time, like in ordinary quantum statistics. It is a Fourier transformation from the probability distribution of velocity circulation to the loop functional. This functional is dual to the trace of the evolution operator for a one-dimensional quantum system.

B. The Physical Meaning of the Loop Equation and Dimensional Reduction

The long-term evolution of Newton’s dynamical system with many particles eventually covers the energy surface (microcanonical ensemble). The ergodic hypothesis (accepted in Physics but still not proven mathematically) states that this energy surface is covered uniformly. The Turbulence theory aims to find a replacement of the microcanonical ensemble for the Navier-Stokes equation. This surface would also take part in the decay in the pure Navier-Stokes equation without artificial forcing. In both cases, Newton and Navier-Stokes, the probability distribution must satisfy the Hopf equation, which follows from the dynamics without specification of the mechanism of the stochastization. Indeed, the Gibbs, as well as the microcanonical distributions in Newton’s dynamics, satisfy the Hopf equation in a rather trivial way: It reduces to the conservation of the probability measure (Liouville theorem), which suggests the energy surface as the only additive integral of motion to use as a fixed point manifold.

In the case of the decaying turbulence, the loop equations represent a closed subset of the Hopf equations, which subset is still sufficient to generate the dynamics of vorticity. In this case, the exact solution we have found for the Hopf functional also follows from the integrals of motion, this time, the conservation laws in the loop space. The loop space Hamiltonian we have derived from the unforced Navier-Stokes equation does not have any potential terms (those with explicit dependence upon the shape of the loop). The Schrödinger equation with only kinetic energy in the Hamiltonian conserves the momentum. The corresponding wave function is the plane wave

exp (ı \vec{p} \cdot \vec{x})

. This plane wave is the solution we have found, except the dot product

\vec{p} \cdot \vec{x}

becomes a symplectic form

\oint \vec{P} (θ) \cdot d \vec{C} (θ)

in the loop space. Our momentum

\vec{P} (θ, t)

is not an integral of motion, but simple scaling properties of the pure Navier-Stokes equation lead to the solution with

\vec{P} (θ, t) \propto \vec{F} (θ) / \sqrt{t}

, with

\vec{F} (θ)

being the integral of motion. The rest is a purely technical task: substituting this scaling solution into the Navier-Stokes equation and solving the resulting universal equation for a fixed point

\vec{F} (θ)

This equation led us to the Fermi ring in the quasiclassical limit. This limit resulted in the energy spectrum and dissipation in a finite system found in the last Section VII.

C. Classical Flow and Quantum Mechanics

Our computations heavily rely on the number theory, particularly Jordan’s multitotients

φ_{l} (q)

, generalizing [30] the Euler totient function. What could the number theory have in common with the turbulent flow? The quantization of parameters of the fixed manifold of decaying turbulence stems from the deep quantum correspondence we have discovered. The statistical distribution of a nonlinear classical Navier-Stokes PDE is exactly related to the wave functional of a linear Schrödinger equation in the loop space.

This relation is indirect: the loop functional, a Fourier transform of the classical probability distribution for circulation, equals the complex quantum amplitude of the loop space theory, as represented by a Fermion ring. Probability is real and positive, but its Fourier transform is complex, so there is no contradiction between classical and quantum mechanics.

The quantization mechanism is the same as in ordinary quantum mechanics; this requires the solution’s periodicity.

An expert in the traditional approach to turbulence may wonder why the loop equation’s solutions are related to the velocity field’s statistics in a decaying turbulent flow. How could this fractal curve in complex momentum space

\vec{P} (θ) \in C_{3}

describe velocity field in a physical space

\vec{v} (\vec{r}) \in R_{3}

? Is this a miracle or a mirage?

The equivalence of a strong coupling phase of the fluctuating vector field to quantum geometry is a well-known duality phenomenon in gauge theory (the ADS/CFT duality), ringing a bell to the modern theoretical physicist. In our case, this is a simpler quantum geometry: a fractal curve in complex space.

Miracles are known to happen in nonlinear PDE: some lucky nonlinear equations are related to linear problems in higher dimensions [8]. In our case, the Navier-Stokes equation is related to the linear problem in infinite-dimensional loop space. The miracle already happened when Hopf derived his linear functional equation. The loop equation is a reduced version of the Hopf equation, luckily solvable by a plane wave in loop space due to translational invariance.

The loop technology was thoroughly discussed in the last few decades in the gauge theories, including QCD [2,13,18,29] where the loop equations were derived first [14,17]. The short answer is that duality only applies to the correlation functions of two theories with different dynamical variables; there is no correspondence between these variables, but the correlation functions are identical.

Mathematical physics sometimes has alternative languages for the same phenomena; examples are the duality between Schrödinger wave equation and Heisenberg’s matrix mechanics, between dynamical triangulation and Liouville theory in 2D quantum gravity.

Extra complications in the gauge theory are the short-distance singularities related to the infinite number of fluctuating degrees of freedom in quantum field theory. The Wilson loop functionals in coordinate space are singular in the gauge field theory and cannot be multiplicatively renormalized.

Perturbatively, there is no short-distance divergence in the Navier-Stokes equations nor the Navier-Stokes loop dynamics. The Euler equations represent the singular limit, which, as we argued, should be resolved using singular topological solitons regularized by the Burgers vortex. In the present theory, we do not encounter any singularities in coordinate space. The anomalous dissipation is achieved in our solution via a completely different mechanism: large fluctuations of the fractal curve at

p ≪ q

. However, these singularities only happen in the inviscid limit

ν \propto 1 / N^{2} \to 0

Therefore, these singularities correspond to the Euler singularities, such as line vortexes. Those vortexes are regularized by finite viscosity, just like our singularities.

D. Stokes-Type Functionals and Vorticity Correlations

The loop equation describes the gauge invariant sector of the gauge field theory. Therefore, the gauge degrees of freedom are lost in the loop functional. However, the gauge-invariant correlations of the field strength are recoverable from the solutions of the loop equation [14,17].

There is no gauge invariance regarding the velocity field in fluid dynamics (though there is such invariance in the Clebsch variables [21]). The longitudinal, i.e., a potential part of the velocity, has a physical meaning – it is responsible for pressure and energy pumping. This part is lost in the loop functional but is recoverable from the rotational part (the vorticity) using the Biot-Savart integral.

In the Fourier space, the correlation functions of the velocity field are algebraically related to those of vorticity

{\vec{v}}_{k} = \frac{ı \vec{k} \times {\vec{ω}}_{k}}{{\vec{k}}^{2}}

. Thus, the general solution for the Wilson loop functional

Ψ [γ, C]

allows computing both vorticity and velocity correlation functions.

We demonstrated that by computing the moments of the enstrophy and resulting anomalous dissipation. This computation is nonperturbative: it corresponds to the extreme turbulent limit and cannot be expanded in inverse powers of viscosity.

E. Renormalizability of the Inviscid Limit of the Loop Equation

Let us discuss the relation

ν N^{2} = const

between the vanishing viscosity and growing number of discontinuities on the momentum loop

\vec{P} (θ)

The Navier-Stokes equation itself is an idealization of the molecular dynamics: the nonlocal theory is approximated by truncated expansion in powers of gradients.

The laminar flow does not show any problems with this truncated expansion. However, the velocity field in our solution for the loop equation becomes singular in the local limit.

The mathematical statement is that the velocity and vorticity are not the ordinary functions in

R_{3}

but rather stochastic variables with

Δ \vec{v} \sim {(Δ \vec{r})}^{α}, α \approx \frac{1}{3}

The Fractal Calculus [9] was introduced for a general description of such fields. That, however, by itself does not explain turbulence: in particular, the more general power laws with multifractal dimensions

〈{(Δ \vec{v})}^{n}〉 \sim {(Δ \vec{r})}^{ζ_{n}}

cannot be explained this way.

In our theory, there is a universal function for the expectation of Fourier components of vorticity as functions of the wavevector. Still, this relation is not a single power: the Mellin transform of this function has infinitely many poles in the complex plane.

Our solution describes a fractal vorticity field, but this is not a known fractal type. The fractal curve

\vec{P} (θ)

is not a simple random walk on a circle: dynamic restrictions (such as periodicity) lead to nontrivial critical behavior not describable by any fractal power laws.

Our discretization of the loop equation (replacement of continuum loop by a polygon with a growing number of vertices) defines the turbulent limit of the velocity field distribution.

The relation of the vanishing viscosity with the growing number of degrees of freedom is of the same nature as the RG relation between the bare coupling constant

g_{0}

and the lattice spacing a in QCD:

a \propto g_{0}^{α} exp (- β / g_{0}^{2})

The naive local limit

a \to 0, g_{0} = c o n s t

does not exist in

Q C D

, but the RG limit describes the strong interaction of hadrons very well.

In the same way, as in QCD, there is renormalizability in our solution. The dissipation rate stays finite when

\tilde{ν} = ν N^{2} \to c o n s t

. Moreover, the whole energy spectrum is expressed in

\tilde{ν}

and stays finite in the local limit. So, like in the renormalizable QFT, there is a "dimensional transmutation": infinities can be hidden into the dimensional parameter

\tilde{ν}

defining the time scale.

F. Relation of Our Solution to the Weak Turbulence

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

The Navier-Stokes Wilson loop is a case of the Abelian loop functional, with commuting components of the vector field

\vec{v}

. As we discussed in detail in [14,17,21], any Stokes-type functional

Ψ [γ, C]

satisfying boundary condition at shrunk loop

Ψ [0] = 1

, and solving the loop equation can be iterated in the nonlinear term in the Navier-Stokes equations (which iterations would apply at large viscosity).

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

We have demonstrated in [20,21] (and also here, in Section B) how the velocity distribution for the random uniform vorticity in the fluid was reproduced by a singular momentum loop

\vec{P} (θ)

The solution for

\vec{P} (θ)

in this special fixed point of the loop equation was random complex and had slowly decreasing Fourier coefficients, leading to a discontinuity

sign (θ - θ^{'})

in a pair correlation function (33). The corresponding Wilson loop was equal to the Stokes-type functional (25).

Using this example, we demonstrated how a discontinuous momentum loop describes the vorticity distribution in the stochastic Navier-Stokes flow. In this example, the vorticity is a global random variable corresponding to a random uniform fluid rotation: a well-known exact solution of the Navier-Stokes equation.

This example corresponds to a special fixed point for the loop equation, not general enough to describe the turbulent flow but mathematically ideal as a toy model for the loop technology. It demonstrates how the momentum loop solution sums up all the terms of the

1 / ν

expansion in the Navier-Stokes equation.

In our general solution, with the Euler ensemble, the summation of a divergent perturbation expansion occurs at an extreme level, leading to a universal fixed point in the limit of vanishing viscosity.

At a given initial condition, after a finite time, the solution will still depend on viscosity and initial condition. There will be no singularities for a smooth initial velocity field.

At the large time, though, it will approach our universal fixed manifold and (supposedly, for random initial data) cover it uniformly, according to the Euler ensemble measure.

So, the

ν \to 0, N \to \infty

limit of our solution for the fixed point corresponds to the limit

ν \to 0, t \to \infty

of the solution for the time-dependent loop equation.

The vorticity will become a random variable with a singular distribution in the local limit, corresponding to the inviscid limit of the Navier-Stokes equation.

X. Remaining Problems

We performed all the calculations up to numerical factors in the vorticity correlation function, which we recovered from the previously computed $〈{\vec{ω}}^{2}〉$ (see [22] ). It would be useful to compute all the normalization factors and thus double-check the solution.
The loop functional for the circular loop is the simplest object in this theory. It can be computed using the methods developed in this paper, with even simpler results. In this case, the classical equation is trivial, so the computations reduce to the functional determinant and the resolvent. On the other hand, this is an observable quantity, and one could measure it in DNS. It would be an interesting problem to solve and compare with the DNS.
The spectrum of indexes for deviations from our fixed trajectory [22] can be evaluated in the scaling limit, with finite $\tilde{ν}$ . This computation will produce results for the vorticity correlation functions in the Navier-Stokes equation, perturbed by an infinitesimal random force.

Data Availability Statement

The Mathematica^®notebooks used to verify equations and compute some functions are available for download in [24,25,26,27,28].

Acknowledgments

Maxim Bulatov helped me with mathematical suggestions and discussions. He also participated in an unpublished paper [6] where we simulated the Euler ensemble as a Markov chain. The strong cancellations in the vorticity correlation function prevented this simulation from producing reliable numerical results. I filled the gaps in my knowledge of Number Theory in discussions of the Euler ensemble with Peter Sarnak, Alexandru Zaharesku, Konstantin Khanin, and Debmalya Basak. Yang-Hui and I discussed Euler’s totients at the Cambridge University workshop, where this theory was first reported in November 2023. His comments helped me derive the asymptotic distribution for scaling variables. I am also grateful to the organizers and participants of the "Field Theory and Turbulence " workshop at ICTS in Bengaluru, India, where this work was advanced in December 2023. Discussions with Katepalli Sreenivasan, Rahul Pandit, and Gregory Falkovich were especially useful. They helped me understand the physics of decaying turbulence in a finite system and match my solution with the DNS data. Recently, this work was discussed at the "Conformal Field Theory, Integrability, and Geometry" conference in Stony Brook in March 11-15, 2024. I am very grateful to Nikita Nekrasov, Sasha Polyakov, Sasha Zamolodchikov, Dennis Sullivan and other participants for deep and inspiring discussions. Finally, this work was thoroughly discussed for a week at the Perimeter Institute for Theoretical Physics, where I spoke at a colloquium on April 3. I received many interesting questions and comments, which are reflected in the last section, "Discussion," of this paper. I thank Davide Gaiotto, Luis Lehner, Sabrina Pastersky, Sergey Sibiryakov, Lee Smolin, and Andrey Shkerin for inspiring discussions. The help from Zechuan Zheng in Mathematica computations of the energy spectrum was also very useful. This research was supported by a Simons Foundation award ID 686282 at NYU Abu Dhabi. Small computations were done on the laptop, the big ones- on the High-Performance Computing resources at NYU Abu Dhabi.

Appendix A Mellin Integral for the Energy Spectrum

The spectral function

H (κ)

can be computed as follows. Let us represent the theta function as inverse Mellin transform (at

n > 0, τ > 0

)

\begin{matrix} θ (1 - n τ) = \int_{+ 0 - ı \infty}^{+ 0 + ı \infty} \frac{d p}{2 π ı p} n^{- p} τ^{- p} \end{matrix}

(A1)

Substituting this representation into our integral in (223), and interchanging summations/integrations, we find

\begin{matrix} H (κ) = \int_{Δ_{1}}^{Δ_{2}} d Δ (1 - Δ) \int_{+ 0 - ı \infty}^{+ 0 + ı \infty} \frac{d p}{2 π ı p} \int_{0}^{\infty} d t \sum_{n = 1}^{\infty} \frac{φ (n)}{n^{p}} (\frac{t^{5 / 2}}{n^{5}} - t^{15 / 2}) t^{- p} (A + B t κ) e^{- C t κ}; \end{matrix}

(A2)

\begin{matrix} A = \frac{Q_{α} (Δ, 1) \sqrt{r_{0} (Δ)}}{Z} \frac{2 (r_{0} (Δ) - 6)}{(12 + r_{0} (Δ))}; \end{matrix}

(A3)

\begin{matrix} B = \frac{Q_{α} (Δ, 1) \sqrt{r_{0} (Δ)}}{Z} \frac{J (Δ)}{S (Δ)}; \end{matrix}

(A4)

\begin{matrix} C = \frac{L (Δ)}{2 π S (Δ)}; \end{matrix}

(A5)

The last integral reduces to Gamma functions and powers of n, after which the sum over n reduces to zeta functions

\begin{matrix} H (κ) = \int_{Δ_{1}}^{Δ_{2}} d Δ (1 - Δ) \int_{- 2 - ı \infty}^{- 2 + ı \infty} \frac{d p κ^{p}}{2 π ı} \frac{20 C^{p - 1} (A C - B p) ζ (p + \frac{15}{2}) Γ (- p)}{(2 p + 7) (2 p + 17) ζ (p + \frac{17}{2})} \end{matrix}

(A6)

According to the Riemann hypothesis, the function

ζ (z), z = p + \frac{17}{2}

in the denominator has trivial zeros at negative even points

z = - 2, - 4, \dots

and nontrivial zeros at the critical line

Re z = \frac{1}{2}

There is also a contribution from the poles at

p + \frac{15}{2} = 1, (2 p + 7) = 0, (2 p + 17) = 0

We re not going to sum the residues in these poles, but rather integrate along the line

p = - 2 + ı y

, where the integrand exponentially decreases in both directions.

There are oscillations on top of this exponential decrease. The integrals can be computed with arbitrary accuracy in Mathematica^®using the "DoubleExponentialOscillatory" integration method.

Figure A1. Log-log plot of the universal function

H (κ)

. It starts as a constant, then turns down asymptotically decaying as

H (κ) \propto κ^{- \frac{7}{2}}

Figure A1. Log-log plot of the universal function

H (κ)

. It starts as a constant, then turns down asymptotically decaying as

H (κ) \propto κ^{- \frac{7}{2}}

Each integration takes a fraction of a second. The remaining integral over

Δ

was performed by computing the table of (positive) values of integrand

I (Δ)

on a grid with step

δ Δ = \frac{Δ_{2} - Δ_{1}}{100}

, and interpolating

log (I (Δ))

between these values by fifth order polynomial

P_{5} (Δ)

. Integral of the exponential of this polynomial

\int d Δ (1 - Δ) exp (P_{5} (Δ))

yields good enough accuracy ( floating precision

10^{- 8}

This process was repeated for all M =1000 values of

κ = 0.1 {(1000 / 0.1)}^{n / M}, n = 0, \dots, M

, correspobnding to equidistant log grid.

Here is the resulting function (see Figure A1).

Figure A2. The log-log plot energy dissipation

t^{2} E (t)

as a function of t has a regime change due to quantum effects. It starts as constant at small t and asymptotically decays as

t^{2} E \propto t^{- \frac{1}{4}}

Figure A2. The log-log plot energy dissipation

t^{2} E (t)

as a function of t has a regime change due to quantum effects. It starts as constant at small t and asymptotically decays as

t^{2} E \propto t^{- \frac{1}{4}}

The asymptotic behavior of this function is

\begin{matrix} H (κ) \propto κ^{- \frac{7}{2}} \end{matrix}

(A7)

The effective index

s (t) = - \frac{d log H (κ)}{d log κ}

slowly increases reaching

s (\infty) = \frac{7}{2}

The function

F (κ) = \int_{κ}^{\infty} x^{2} H (x) d x

was obtained by numerical integration of exponential of interpolated function as before, with asymptotic tail

H (κ) \to c o n s t κ^{- \frac{7}{2}}

integrated analytically.

The log-log plot of the energy dissipation as a function of time is shown in Figure A2, and it is curved in the log-log scale.

References

M.E. AGISHTEIN and A.A. MIGDAL. Simulations of four-dimensional simplicial quantum gravity as dynamical triangulation. Modern Physics Letters A, 07(12):1039–1061, 1992.
Peter D. Anderson and Martin Kruczenski. Loop equations and bootstrap methods in the lattice. Nuclear Physics B, 921:702–726, 2017.
Gabriel Apolinario, Luca Moriconi, Rodrigo Pereira, and Victor valadão. Vortex gas modeling of turbulent circulation statistics. PHYSICAL REVIEW E, 102:041102, 10 2020.
Abhay Ashtekar. New variables for classical and quantum gravity. Physical Review Letters, 57:2244–2247, 11 1986.
Debmalya Basak and Alexandru Zaharescu. Connections between number theory and the theory of turbulence, 2024. To be published.
Maxim Bulatov and Alexander Migdal. Numerical simulation of euler ensemble. "unpublished", February 2024.
Richard P Feynman, Robert B Leighton, and Matthew Sands. 3-7 how did it get that way? In Michael A. Gottlieb and Rudolf Pfeiffer, editors, The Feynman lectures on physics, Vol. I: Mainly mechanics, radiation, and heat, volume 1. Basic books, 2011.
Clifford S Gardner, John M Greene, Martin D Kruskal, and Robert M Miura. Method for solving the korteweg-devries equation. Physical review letters, 19(19):1095, 1967.
Alireza Khalili Golmankhaneh. Fractal Calculus and its Applications. WORLD SCIENTIFIC, 2022.
G. H. Hardy and E. M. Wright. An introduction to the theory of numbers. Oxford University Press, Oxford, sixth edition, 2008. Revised by D. R. Heath-Brown and J. H. Silverman, With a foreword by Andrew Wiles.
Kartik P Iyer, Sachin S Bharadwaj, and Katepalli R Sreenivasan. The area rule for circulation in three-dimensional turbulence. Proceedings of the National Academy of Sciences of the United States of America, 118(43):e2114679118, 10 2021.
Kartik P. Iyer, Katepalli R. Sreenivasan, and P. K. Yeung. Circulation in high reynolds number isotropic turbulence is a bifractal. Phys. Rev. X, 9:041006, 10 2019.
Vladimir Kazakov and Zechuan Zheng. Bootstrap for lattice Yang-Mills theory. Phys. Rev. D, 107(5):L051501, 2023.
Yu.M. Makeenko and A.A. Migdal. Exact equation for the loop average in multicolor qcd. Physics Letters B, 88(1):135–137, 1979.
T. Matsuzawa and W. Irvine. Realization of confined turbulence through multiple vortex ring collisions, 03/12/2019. "Talk at the Flatiron Conference Universality Turbulence Across Vast Scales".
Takumi Matsuzawa, Noah P. Mitchell, Stéphane Perrard, and William T. Irvine. Creation of an isolated turbulent blob fed by vortex rings. Nature Physics, 19(8):1193–1200, 2023.
A. Migdal. Loop equations and $\frac{1}{N}$ expansion. Physics Reports, 201, 1983.
A.A. Migdal. Momentum loop dynamics and random surfaces in qcd. Nuclear Physics B, 265(4):594–614, 1986.
A.A. Migdal. Second quantization of the wilson loop. Nuclear Physics B - Proceedings Supplements, 41(1):151–183, 1995.
Alexander Migdal. Loop equation and area law in turbulence. In Laurent Baulieu, Vladimir Dotsenko, Vladimir Kazakov, and Paul Windey, editors, Quantum Field Theory and String Theory, pages 193–231. Springer US, 1995.
Alexander Migdal. Statistical equilibrium of circulating fluids. Physics Reports, 1011C:1–117, 2023.
Alexander Migdal. To the theory of decaying turbulence. Fractal and Fractional, 7(10):754, Oct 2023.
Alexander Migdal. Topological vortexes, asymptotic freedom, and multifractals. MDPI Fractals and Fractional, Special Issue, 2023.
Alexander Migdal. "abc". "https://www.wolframcloud.com/obj/sasha.migdal/Published/ABCInterpolator.m.nb", 04 2024.
Alexander Migdal. "bernsum". "https://www.wolframcloud.com/obj/sasha.migdal/Published/BernSum.nb", 02 2024.
Alexander Migdal. "instantoncomputations". "https://www.wolframcloud.com/obj/sasha.migdal/Published/DualTheory3.nb", 04 2024.
Alexander Migdal. "instantoncomputations". "https://www.wolframcloud.com/obj/sasha.migdal/Published/DualTheory1.nb", 03 2024.
Alexander Migdal. "instantoncomputations". "https://www.wolframcloud.com/obj/sasha.migdal/Published/DualTheory2.nb", 03 2024.
Alexander A. Migdal. Hidden symmetries of large N QCD. Prog. Theor. Phys. Suppl., 131:269–307, 1998.
Pieter Moree, Sumaia Saad Eddin, Alisa Sedunova, and Yuta Suzuki. Jordan totient quotients. Journal of Number Theory, 209:147–166, 2020.
Nicolás P. Müller, Juan Ignacio Polanco, and Giorgio Krstulovic. Intermittency of velocity circulation in quantum turbulence. Phys. Rev. X, 11:011053, 03 2021.
James R. Norris. Markov chains. Cambridge Univ. Press, 2007.
Nikita Nekrasov (private communication) suggested to me an algorithm of generating this solution as a set of adjacent triangles in complex 3-space and pointed out an invariant measure in phase space, made of lengths of shared sides and angles between them. Unfortunately, this beautiful construction does not guarantee real circulation, requiring further study.
Koji Ohkitani. Study of the hopf functional equation for turbulence: Duhamel principle and dynamical scaling. Physical Review E, 101(1):013104, 2020.
John Panickacheril John, Diego A Donzis, and Katepalli R Sreenivasan. Laws of turbulence decay from direct numerical simulations. Philos. Trans. A Math. Phys. Eng. Sci., 380(2218):20210089, March 2022.
G Parisi and U Frisch. On the singularity structure of fully developed turbulence turbulence and predictability. In R Benzi M Ghil and G Parisi, editors, Geophysical Fluid Dynamics: Proc. Intl School of Physics E. Fermi, pages 84–88. Amsterdam: North-Holland, 1985.
Juan Ignacio Polanco, Nicolás P. Müller, and Giorgio Krstulovic. Vortex clustering, polarisation and circulation intermittency in classical and quantum turbulence. Nature Communications, 12(1):7090, Dec 2021.
Carlo Rovelli and Lee Smolin. Knot theory and quantum gravity. Phys. Rev. Lett., 61:1155–1158, 09 1988.

Nikita Nekrasov (private communication) suggested to me an algorithm of generating this solution as a set of adjacent triangles in complex 3-space and pointed out an invariant measure in phase space, made of lengths of shared sides and angles between them. Unfortunately, this beautiful construction does not guarantee real circulation, requiring further study.

Figure 2. The path to the microscopic theory.

Figure 6. The plot of the function

Λ (s)

. As required, it is positive, takes a maximal value

log (2)

s = 0

, and vanishes at both ends

s = \pm 1

of the physical region.

Figure 6. The plot of the function

Λ (s)

. As required, it is positive, takes a maximal value

log (2)

s = 0

, and vanishes at both ends

s = \pm 1

of the physical region.

Figure 7. The log-log-plot of the distribution

f_{X} (X)

. It is equals

π Φ (k) X \sqrt{X}

at each internal

\frac{1}{π^{2} {(k + 1)}^{2}} < X < \frac{1}{π^{2} k^{2}}

with positive integer k. Asymptotically

f_{X} (X) \to \frac{3 \sqrt{X}}{π^{3}}

X \to 0

Figure 7. The log-log-plot of the distribution

f_{X} (X)

. It is equals

π Φ (k) X \sqrt{X}

at each internal

\frac{1}{π^{2} {(k + 1)}^{2}} < X < \frac{1}{π^{2} k^{2}}

with positive integer k. Asymptotically

f_{X} (X) \to \frac{3 \sqrt{X}}{π^{3}}

X \to 0

Figure 8. Log plot of

r_{0} (Δ) .

Figure 8. Log plot of

r_{0} (Δ) .

Figure 9. Plot of

S (Δ)

Figure 9. Plot of

S (Δ)

Figure 10. Plot of

A_{c} (r_{0} (Δ), Δ)

for the three real solutions of the equation

g (r, Δ) = 0

. At

Δ = Δ_{1}

and

Δ = Δ_{2}

, the action curves intersect; at

Δ = \frac{1}{2}

, there is a gap between the lowest action (

= 0

) and the lowest of the other two. So, there are second-order phase transitions at

Δ_{1}, Δ_{2}

and the first-order phase transition at

Δ = \frac{1}{2}

Figure 10. Plot of

A_{c} (r_{0} (Δ), Δ)

for the three real solutions of the equation

g (r, Δ) = 0

. At

Δ = Δ_{1}

and

Δ = Δ_{2}

, the action curves intersect; at

Δ = \frac{1}{2}

, there is a gap between the lowest action (

= 0

) and the lowest of the other two. So, there are second-order phase transitions at

Δ_{1}, Δ_{2}

and the first-order phase transition at

Δ = \frac{1}{2}

Figure 11. Plot of universal function

J (Δ)

.the four corves correspond to four phases (solutions for

r_{0} (Δ)

Figure 11. Plot of universal function

J (Δ)

.the four corves correspond to four phases (solutions for

r_{0} (Δ)

Figure 14. Complex poles of the Mellin transform

h (p)

in (242) defining the spectrum of decay indexes.

Figure 14. Complex poles of the Mellin transform

h (p)

in (242) defining the spectrum of decay indexes.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Quantum Solution of Classical Turbulence. Decaying Energy Spectrum

Abstract

I. Prologue

II. Summary

III. Introduction

A. Physical Introduction. The Energy Flow and Random Vorticity Structures

B. Mathematical Introduction. The Loop Equation and Its Solution

C. The Big and Small Euler Ensembles

IV. The Markov Chain and Its Fermionic Representation

V. The Continuum Limit

A. Path Integral Over Markov Histories

B. Matching Path Integral with Combinatorial Sums in Big Euler Ensemble

C. Small Euler Ensemble in Statistical Limit

D. Complex Classical Trajectory in the Path Integral

VI. Dual Theory of Vorticity Correlation

A. Correlation Function and Path Integral

B. Turbulent Viscosity and the Local Limit

C. Functional Determinant in the Path Integral

D. The Fluctuation Term in α ′ ( ξ 1 ) α ′ ( ξ 2 )

VII. The Decaying Energy in Finite System

VIII. The Spectrum of Decay Indexes in the Energy Spectrum

IX. Discussion

A. Stochastic Solution of the Navier-Stokes Equation and Quantum Hypothesis

B. The Physical Meaning of the Loop Equation and Dimensional Reduction

C. Classical Flow and Quantum Mechanics

D. Stokes-Type Functionals and Vorticity Correlations

E. Renormalizability of the Inviscid Limit of the Loop Equation

F. Relation of Our Solution to the Weak Turbulence

X. Remaining Problems

Data Availability Statement

Acknowledgments

Appendix A Mellin Integral for the Energy Spectrum

References

MDPI Initiatives

Important Links

Subscribe

D. The Fluctuation Term in $α^{'} (ξ_{1}) α^{'} (ξ_{2})$