Preprint
Article

Dual Theory of Decaying Turbulence

This version is not peer-reviewed.

Submitted:

07 March 2024

Posted:

07 March 2024

Read the latest preprint version here

A peer-reviewed article of this preprint also exists.

Abstract
This paper investigates the recently found \cite{migdal2023exact} reduction of decaying turbulence in the Navier-Stokes equation 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 quantum statistical theory of free fermions on a ring, with an external field related to the random fractions of $\pi$. We find the solution of this system in the statistical limit $N\to \infty$ in terms of a complex trajectory (instanton) providing a saddle point to the path integral over the charge 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 in quadratures, but it is too complex for an immediate computation; a separate investigation by Number Theory methods is required.
Keywords: 
;  ;  ;  ;  ;  ;  ;  ;  ;  

1. 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. There are 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 [2] and its solution [1] 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 structure with yet-to-be-proven relation to the decaying turbulence. The mathematical analysis of this ensemble does not rely on the physics behind this theory and represents a well-defined problem in Number theory [3].

1.1. 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 the nonlinearity of 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 [4].
The perturbative approach is inadequate for the turbulence theory, which must be built from the first principles by solving the Navier-Stokes equations beyond the perturbation theory. 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 [4] suggest the decay of the dissipation rate E t n with n 1.2 or 1.6 depending on initial conditions (finite total momentum v d 3 r or zero total momentum but finite total angular momentum r × v d 3 r ). So, two universality classes of decaying turbulence were observed in these experiments.
It is not clear whether the data in [4] 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 [1], 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 [4] used such a setup, with forcing turned off. They observed subsequent energy dissipation, apparently reaching one of the two universal regimes (starting as k 2 or 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 v and averages over an infinite time interval, reads:
V d 3 r ν ω 2 = V d 3 r β v β p + 1 2 v α 2 + ν v α ( β v α α v β )
By the Stokes theorem, the right side reduces to the flow through the boundary V of the integration region V. The left side is the dissipation in this volume, so we find:
E V = V d σ · v p + 1 2 v α 2 + ν ω × v
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 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 ω × v term drops as there is no vorticity at infinity.
Furthermore, the velocity in the Biot-Savart law decreases as | r | 3 at infinity, so that only the v p term survives
E V V d σ · v p
This energy flow on the right side will stay finite in the limit of the expanding sphere in case the pressure grows as p f · r , where f is the local force at a given point on a large sphere.
E = f α lim R R 3 S 2 n α n β v β ( R n )
Where did we lose the Kolmogorov energy flow? It is still there for any finite volume surrounding the vortex sheet
E V = V d 3 r v β β p + v α v β β v α =
V d 3 r v α v β β v α V d σ · v p
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 f acting on the boundary, say, with periodic boundary conditions, the boundary integral would be absent, and we would recover the Kolmogorov relation.
v α v β β v α = E V / V ;
This relation, together with space symmetry properties in R d leads to the Kolmogorov three-point correlation
v α ( r 0 ) v β ( r 0 ) v γ ( r + r 0 ) = E ( d 1 ) ( d + 2 ) V δ α γ r β + δ β γ r α 2 d δ α β r γ ;
In the conventional approach, based on the time averaging of the Navier-Stokes equations, the periodic Gaussian random force f ( r ) is added to the right side. In this case, with periodic boundary conditions
E V = V d 3 r v β β p v β f β ( r ) + v α v β β v α = V d 3 r v β f β ( r )
In the limit when the force becomes uniform in space, we recover another definition with E = f · P , where P = V d 3 r 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 [5]. In the dual picture of our recent theory [1], these are the random gaps in the momentum curve P ( θ ) .
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
v ( r ) = × d 3 r ω ( r ) 4 π | r r |
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 V with radius R ,
lim R R 3 v ( R n ) 1 4 π blobs blob d 3 r ω ( r ) × r < r > blob ;
Here < 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 ([6]). 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 [7] 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 [4].
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 ) k 2 and another for E ( k , 0 ) 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 to all 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 [1] of the loop equations [2,8] describe the internal stochastization of the decaying turbulence by a dual discrete system. Measuring these internal stochastic phenomena challenges real experiments and numerical simulations alike.

1.2. Mathematical Introduction. The Loop Equation and Its Solution

In the previous paper, [1], we have found a family of exact solutions of the loop equation for decaying turbulence [2,8]. 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
Ψ [ γ , C ] = exp ι γ ν d C ( θ ) · v C ( θ ) i n i t ;
Ψ [ γ , C ] exp ι γ ν d C ( θ ) P ( θ ) E ;
In the first equation (the definition), the averaging goes over initial data for the solutions of the Navier-Stokes equation for velocity field v ( r ) . In the second one (the solution), the averaging goes over the distribution of the random variable P ( θ ) to be determined from the loop equation [1]. We choose in this paper the parametrization of the loop with ξ = θ 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 P ( θ ) follows from the Navier-Stokes equation for v
t v α = ν β ω β α v β ω β α α p + v β 2 2 ;
α v α = 0 ;
ω β α = β v α α v β
After some transformations, replacing velocity and vorticity with the functional derivatives of the loop functional, we found the following momentum loop equation in [1,8]
ν t P = γ 2 ( Δ P ) 2 P + Δ P γ 2 P · Δ P + ι γ ( P · Δ P ) 2 Δ P 2 P 2 ;
P ( θ ) P ( θ + ) + P ( θ ) 2 ;
Δ P ( θ ) P ( θ + ) P ( θ ) ;
The momentum loop has a discontinuity Δ P ( θ ) at every parameter 0 < θ 1 , making it a fractal curve in complex space C d . The details can be found in [1,8]. 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 P ( θ ) at t = 0 , providing the initial velocity field distribution. Formally, the initial distribution W 0 [ P ] of the momentum field P ( θ ) is given by inverse functional Fourier transform.
W 0 [ P ] = D C δ 3 ( C [ 0 ] ) Ψ [ C , γ ] t = 0 exp ι γ ν d C ( θ ) · P ( θ )
This path integral was computed in [1,8] 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)
Ψ 0 [ C ] Ψ [ C , 1 ] t = 0 = exp 1 2 ν 2 C d C ( θ ) · d C ( θ ) f C ( θ ) C ( θ )
where f ( r ) is the velocity correlation function
v α ( r ) v β ( r ) = δ α β α β μ 2 f ( r r )
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 ) f 0 f 1 ( r r ) 2 +
The first term f 0 is proportional to initial energy density,
1 2 v α 2 = d 1 2 f 0
and the second one is proportional to initial energy dissipation rate E 0
f 1 = E 0 2 d ( d 1 ) ν
where d = 3 is the dimension of space. The constant term in (23) 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
Ψ 0 [ C ] exp f 1 ν 2 d C α ( θ ) C β ( θ ) 2
This distribution is almost Gaussian: it reduces to Gaussian one by extra integration
Ψ 0 [ C ] const ( d ϕ ) exp ϕ α β 2 exp 2 ι f 1 ν ϕ μ ν d C μ ( θ ) C ν ( θ )
The integration here involves all 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 P ( θ ) , defined by the Fourier expansion on the circle
P α ( θ ) = odd n = 1 P α , n exp ι n θ + P ¯ α , n exp ι n θ ;
P α , n = N ( 0 , 1 ) ;
P ¯ α , n = 4 f 1 n ν ϕ α β P β , n ;
ϕ α β = ϕ β α ;
ϕ α β = N ( 0 , 1 ) α < β ;
At fixed tensor ϕ the correlations are
P α , n P β , m t = 0 = 4 f 1 m ν δ n m ϕ α β ;
P α ( θ ) P β ( θ ) t = 0 = 2 ι f 1 ν ϕ α β sign ( θ θ ) ;
Ψ 0 [ C ] = exp ι ν d C ( θ ) P ( θ ) P , ϕ
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 P ( θ ) .
Rather than solving the Cauchy problem, we are looking for an attractor: the fixed trajectory for P ( θ , t ) with some universal probability distribution related to the decaying turbulence statistics.
The following transformation reveals the hidden scaling invariance of decaying turbulence
P = ν 2 ( t + t 0 ) F γ
The new vector function F satisfies an equation
2 τ F = 1 ( Δ F ) 2 F + Δ F γ 2 F · Δ F + ι γ ( F · Δ F ) 2 Δ F 2 F 2 ;
τ = log ( t + t 0 )
This equation is invariant under translations of the new variable τ = log ( t + t 0 ) , corresponding to the rescaling/translation of the original time.
t λ t + ( λ 1 ) t 0
There are two consequences of this invariance.
  • There is a fixed point for F .
  • The approach to this fixed point is exponential in τ , which is power-like in original time.
Both of these properties were used in [1]: 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.

1.3. The Big and Small Euler Ensembles

Let us remember the basic properties of the fixed point for F in [1]. It is defined as a limit N of the polygon F 0 F N = F 0 with the following vertices
F k = cos ( α k ) , sin ( α k ) , i cos β 2 2 sin β 2 ;
θ k = k N ; β = 2 π p q ; N ;
α k + 1 = α k + σ k β ; σ k = ± 1 , β σ k = 2 π p r ;
The parameters Ω ^ , p , q , r , σ 0 σ N = σ 0 are random, making this solution for F ( θ ) a fixed manifold rather than a fixed point.
It is a fixed point of (37) with the discrete version of discontinuity and principal value:
Δ F F k + 1 F k ;
F F k + 1 + F k 2
Both terms of the right side (37) vanish; the term proportional to Δ F and the term proportional to F . Otherwise, we would have F Δ F , leading to zero vorticity [1]. The ensemble of all the different solutions is called the big Euler ensemble. The integer numbers σ k = ± 1 came as the solution of the loop equation, and the requirement of the rational p q came from the periodicity requirement.
We can use integration (summation) by parts to write the circulation as follows (in virtue of periodicity):
d C ( θ ) · P ( θ ) = d P ( θ ) · C ( θ ) ;
k Δ C k P k = Δ P k · C k ;
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 [9]. 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 p q with fixed denominator are counted by Euler totient function φ ( q ) [10]
φ ( q ) = p = 1 ( p , q ) q 1 1 = q p | q 1 1 p ;
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 .

2. 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 [1]. 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 , σ 1 one after another. At each step, there will be M = N , 0 remaining σ . We get a transition n n 1 with probability n M and n 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
M = 1 N Q ( M ) ;
Q ( M ) | n = M n M | n + n M | n 1 ;
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
P [ X ] = n = 0 N + n M = 1 N Q ^ ( M ) · X · N + ;
Q ^ ( M ) · X | n = n M X σ M 1 | n 1 + M n M X σ M 1 | n
Here N + = ( N + σ l ) / 2 = ( N + q r ) / 2 is the number of positive sigmas. The operator Q ^ ( M ) sets in X | n the variable σ M to 1 with probability n M and to 1 with complementary probability. The generalization of the Markov matrix Q ( M ) to the operator 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 [11], 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 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 ] δ n , 0 ):
a i , a j + = δ i j ;
a i , a j + = a i , a j + = 0 ;
a n | σ 1 , , σ N = δ [ σ n + 1 ] | σ 1 , , σ n 1 , σ N ;
a n | σ 1 , , σ N = δ [ σ n 1 ] | σ 1 , , σ n 1 , σ N ;
ν ^ n = a n a n ;
ν ^ n | σ 1 , , σ N = δ [ σ n 1 ] | σ 1 , , σ N
The number n ( M ) of positive sigmas l = 1 M δ [ σ l 1 ] coincides with the occupation number of these Fermi particles.
n ^ ( M ) = l = 1 M ν ^ l ;
This relation leads to the representation
Q ^ ( M ) = ν ^ M n ^ ( M ) M + ( 1 ν ^ M ) M n ^ ( M ) M ;
The variables σ l can also be expressed in terms of this operator algebra by using
σ ^ l = 2 ν ^ l 1 .
The Wilson loop in (12) can now be represented as an average over the small Euler ensemble E ( N ) of a quantum trace expression
Ψ [ γ , C ] = tr Z ^ ( q r ) exp i γ ν l Δ C l · Ω ^ · P l ( t ) M = 1 N Q ^ ( M ) Ω ^ , E ( N ) tr Z ^ ( q r ) M = 1 N Q ^ ( M ) E ( N ) ,
Z ^ ( s ) = d ω 2 π exp ι ω l σ ^ l s ;
Δ C l = C l + 1 N C l N ,
P l ( t ) = ν 2 t + t 0 F l γ , Ω ^ O ( 3 ) ,
F l = cos α ^ l , sin α ^ l , 0 2 sin β 2 ,
E ( N ) : p , q , r Z with 0 < p < q < N , gcd ( p , q ) = 1 , N q r N ,
α ^ l = β k = 1 l 1 ( 2 ν ^ k 1 ) ;
The last component of the vector F l is set to 0 as this component does not depend on l and yields zero in the sum over the loop l Δ C l = 0 .
The proof of equivalence to the combinatorial formula with an average over σ l = ± 1 can be given using the following Lemma (obvious for a physicist).
Lemma 1.
The operators ν ^ l all commute with each other.
Proof. 
Using commutation relations, we can write
ν ^ l ν ^ 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
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, ν ^ l ν ^ n = ν ^ n ν ^ l
Theorem 1.
The trace formula (61) equals the expectation value of the momentum loop ansatz (13), (36), (40) in the big Euler ensemble.
Proof. 
As all the operators ν ^ l commute with each other, the operators Q ^ ( M ) can be applied in arbitrary order to the states Σ = | σ 1 , σ N involved in the trace. The same is true about individual terms in the circulation in the exponential of the Wilson loop. These terms F l involve the operators α ^ l , which commute with each other and with each Q ^ ( M ) . Thus, we can use the ordered product of the operators G ^ l = Q ^ ( l ) exp ι ω σ ^ l + i γ ν Δ C l · P l ( t ) . Each of the operators G ^ l acting in turn on arbitrary state Σ will create two terms with δ [ σ l ± 1 ] . The exponential in G ^ l will involve σ ^ k , k l . As a result of the application of the operator Z ^ l = k = 1 l G ^ k to the state vector Σ we get 2 l terms with Σ k = 1 l δ [ σ k η k ] , η k = ± 1 . The factors Z ^ l will involve only σ ^ k , k l , which are all reduced to η k , k l in virtue of the product of the Kronecker deltas. Multiplying all operators G ^ M will lead to superposition Π ^ N of 2 N terms, each with product M = 1 N δ [ σ M η M ] with various choices of the signs η i = ± 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 to a single term corresponding to a particular history η 1 , η 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 i γ ν l Δ C l · Ω ^ · P l ( t ) with the operators σ ^ in P k ( t ) replaced by numbers η leading to the usual numeric 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 [12]. The integration over ω will produce δ l η ^ l s . This delta function will reduce the trace to the required sum over all histories of the Markov process with a fixed 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 . 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 of this paper.

3. The Continuum Limit

As we shall shortly see, in the continuum limit N , the accumulated numbers of Fermi particles ν k = 1 and Dirac holes ν k = 0 tend to some classical functions of “position” ξ = k N , leading to the exact solution for partition function (61).

3.1. 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 ± n ± ( l ) , Δ n ± = 1 )
Π N = exp N Λ N ;
Λ N = 1 N l Δ n + log n + n + + n + Δ n log n n + + n ;
n + = k l ν k ;
n = k l ( 1 ν k ) ;
These n ± are net numbers of η = ± 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
n + + n = l ; l
which follows from the above definition in terms of the occupation numbers.
We can redefine n ± as N times the piecewise constant functions.
n ± = N f ± ( ξ ) ;
f ± ( ξ ) = k = 1 N ξ ν k N ;
f ± ( ξ ) = l = 1 N δ ξ l N k = 1 l ν k N ;
0 < ξ < 1 ;
The sums can be rewritten as Lebesgue integrals
Λ N = 1 0 d f + ( ξ ) log f + ( ξ ) ξ + d f ( ξ ) log f ( ξ ) ξ
The sum over histories of the Markov process will become a path integral over the difference ϕ = f + ( ξ ) f ( ξ )
η . = ± 1 exp N ( Λ N + ι Λ N ( 1 ) ) D ϕ exp N ( Λ N + ι Λ N ( 1 ) )
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
max ϕ Λ N [ ϕ ] ;
Λ N [ ϕ ] = 1 0 d ξ d f + d ξ log f + ξ + d f d ξ log f ξ ;
f ± ( ξ ) = 1 2 ξ ± ϕ ( ξ ) ;
f ± ( ξ ) 0 ;
This problem is, however, a degenerate one, as the functional reduces to the integral of the total derivative:
δ Λ N [ ϕ ] δ ϕ ( ξ ) = 0 ;
Λ N [ ϕ ] = 1 0 d f + log f + + f log f + 1 + 0 1 d ξ log ξ =
1 2 ( 1 s ) log ( 1 s ) 1 2 ( 1 + s ) log ( 1 + s ) + log ( 2 ) ;
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 for the Euler ensemble [1]
lim N log N N ( 1 + s ) / 2 N = 1 2 ( 1 s ) log ( 1 s ) 1 2 ( 1 + s ) log ( 1 + s ) + log ( 2 )
This Λ ( s ) = Λ ( s ) is a smooth even function of s taking positive values from Λ ( ± 1 ) = 0 to the maximal value Λ ( 0 ) = log ( 2 ) (see Figure 1).
Now, let us add the circulation term to the exponential of the partition function (61). This term can be directly expressed in terms of the difference between our two densities N ϕ ( ξ ) = N f + ( ξ ) N f ( ξ ) :
ι N Λ N ( 1 ) [ ϕ , C Ω ] = ι 2 ν t 0 1 d C Ω ( ξ ) · F ( ξ ) ;
F ( ξ ) = sin ( β N ϕ ( ξ ) ) , cos ( β N ϕ ( ξ ) ) , 0 2 sin ( β / 2 ) ;
C Ω ( θ ) = Ω ^ · C ( θ ) ;
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 α ( ξ ) N Δ α = N β σ = ± N β .
As we demonstrate below, in continuum theory, this discrete distribution can be replaced by a Gaussian distribution with the same mean square
α = ± N β d α exp ( α ) 2 2 N 2 β 2
To demonstrate that, we consider in the critical region β 2 N 1 0 the most general term that arises in the moments of the circulation in (87) (see [3] for some exact computations of these moments)
2 N σ i = ± 1 exp ι β i k i σ i = i cos β k i exp β 2 / 2 i k i 2
where k i are some integers. With a large number N of these integers, the sum in the exponential becomes an integral
exp β 2 / 2 i k i 2 exp N β 2 / 2 0 1 d ξ k ( ξ ) 2
The same result follows from the Gaussian integral
i d σ i 2 π exp σ i 2 / 2 exp ι β k i σ i ;
This representation would lead to a standard path integral measure
[ D α ] = D α ( ξ ) exp 0 1 d ξ ( α ) 2 2 N β 2 ; [ D α ] exp ι N 0 1 d ξ α ( ξ ) K ( ξ ) [ D α ] =
exp 1 2 N 2 d ξ 1 d ξ 2 K ( ξ 1 ) K ( ξ 2 ) G ( ξ 1 , ξ 2 ) ;
G ( ξ 1 , ξ 2 ) = α ( ξ 1 ) α ( ξ 2 ) ;
The next section will compute this Green’s function G ( ξ 1 , ξ 2 ) .
Thus, we arrive at the following path integral in the continuum limit
Ψ [ C ] = p < q ; ( p , q ) Ω O ( 3 ) d Ω [ D α ] exp ι 0 1 d ξ Im C Ω ( ξ ) exp ι α ( ξ ) 2 sin ( π p / q ) 2 ν ( t + t 0 ) p < q ; ( p , q ) | O ( 3 ) | [ D α ] ;
C Ω ( θ ) = Ω ^ · C ( θ ) · { ι , 1 , 0 } ;
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 . The effective potential for this theory is a linear function of the loop slope C ( ξ ) .
This model is yet another representation of the Euler ensemble, suitable for the continuum limit.

3.2. Matching Combinatorial Sums of Big Euler Ensemble

The results of the path integration over α must match the combinatorial calculations with σ l = ± 1 in the limit of large N. Without the interaction provided by the circulation term in (97), this path integral is dominated by a linear trajectory
α c l ( ξ ) = β N ξ s ;
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 m )
σ n σ m σ = N s
The direct calculation using methods of [1,3] leads to
σ n σ m σ = N s = d ω 2 π exp ι ω N s cos ( ( N 2 ) ω ) sin 2 ( ω ) / Z =
( N N 2 s 2 ) 2 N 3 ( N 2 N 2 s 2 ) N 2 1 2 ( N + N s 2 ) / Z ;
Z = d ω 2 π exp ι ω N s cos ( N ω ) = 2 N N 1 2 ( s N + N )
Using Gamma function properties, this ratio simplifies to
N s 2 1 N 1
This result can be derived from symmetry without any integration.
σ n σ m σ = N s = A ( N , s ) ( 1 δ n m ) + δ n m ;
n , m σ n σ m σ = N s = s 2 N 2 = A ( N , s ) N ( N 1 ) + N
A ( N , s ) = N s 2 1 N 1
The same limit A ( , s ) = s 2 follows from the classical trajectory
σ n σ m σ = N s α c l ( ξ ) β N α c l ( ξ ) β N = s 2
Let us consider less trivial example [1,3]
U n , m k = n m 1 exp ι α k , n ;
α k , n = β l = 0 l n k σ l ;
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
U n , m = N ξ 1 ξ 2 d ξ exp ι α ( ξ ) = N ξ 1 ξ 2 d ξ exp 1 2 G ( ξ , ξ )
Here G ( ξ 1 , ξ 2 ) is the Green’s function corresponding to a 1D particle on a line interval ξ ( 0 , 1 ) , introduced in the previous section. It satisfies the equation, which follows from our Gaussian Action
ξ 2 G ( ξ , ξ ) = β 2 N δ ( ξ ξ ) ;
G ( 0 , ξ ) = G ( ξ , 0 ) = 0 ;
The solution is
G ( ξ , ξ ) = 1 2 β 2 N ξ + ξ | ξ ξ |
Thus, we find
U n , m = N ξ 1 ξ 2 d ξ exp 1 2 β 2 N ξ = 2 β 2 exp y / 2 exp x / 2 ;
x = β 2 N ξ 1 ;
y = β 2 N ξ 2 ;
in agreement with [1,3] in the critical region N , β 2 1 / N . Finally, the expectation value
U n , m U ¯ n , m = l = n m 1 k = n m 1 exp ι α k n ι α l n
Here, the Gaussian path integration yields
U n , m U ¯ n , m N 2 ξ 1 ξ 2 d ξ ξ 1 ξ 2 d ξ exp 1 2 G ( ξ , ξ ) + G ( ξ , ξ ) 2 G ( ξ , ξ ) = N 2 ξ 1 ξ 2 d ξ ξ 1 ξ 2 d ξ exp 1 2 β 2 N | ξ ξ | = 4 β 4 2 exp ( y x ) / 2 + x y 2
This result also agrees with combinatorial computations in [1,3].

3.3. 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 = q r N is distributed between 1 , 1 with the binomial weight [1,3] 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
W ( r ) = 1 if r = 0 ; 2 π N q exp ( q r ) 2 2 N otherwise ;
As it was conjectured in [1] and supported by rigorous estimates in [3], 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
f p ( p q ) = p = 1 ( p , q ) q 1 δ ( p n ) φ ( q )
As we shall see, rather than p, we would need an asymptotic distribution of a scaling variable
X ( p , q ) = 1 q 2 cot 2 π p q
This distribution for X ( p , q ) at fixed q can be found analytically, using newly established relations for the cotangent sums (see Appendix in [1], and Mathematica® notebook [13]). Asymptotically, at large q, these relations read
X n lim q 1 q p = 1 ( p , q ) q 1 X ( p , q ) n = δ n , 0 + 2 π 2 n ζ ( 2 n ) ( 2 n + 1 ) ζ ( 2 n + 1 )
This relation can be transformed further as
X n = 1 if n = 0 2 k = 1 φ ( k ) k ( 2 n + 1 ) ( 2 n + 1 ) π 2 n if n > 0
The Mellin transform of these moments leads to the following singular distribution
X n = 0 f X ( X ) d X X n ;
f X ( X ) = ( 1 α ) δ ( X ) + π X X Φ 1 π X ;
α = π 0 X X d X Φ 1 π X =
2 π 5 1 Φ ( k ) 1 ( π k ) 5 1 ( π ( k + 1 ) ) 5 = 2 5 π 4 1 φ ( k ) k 5 = 1 225 ζ ( 5 )
where Φ ( n ) is the totient summatory function
Φ ( q ) = n = 1 q φ ( n )
The distribution can also be rewritten as an infinite sum
d x f X ( x ) F ( x ) = ( 1 α ) F ( 0 ) + π n = 1 φ ( n ) 0 1 π 2 n 2 x 3 2 d x F ( x )
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
X m a x = X q 1 , q 1 π 2
Our distribution (124) is consistent with this upper limit, as the argument 1 π X becomes zero at X π 2 > 1 . It is plotted in Figure 2.
Once we are zooming into the tails of the p , q distribution, we also must recall that
P ( q < y N ) = Φ N y Φ N y 2 ;
f y ( y ) = q = 2 N δ y q N φ ( q ) Φ N

3.4. Complex Classical Trajectory in the Path Integral

This classical equation for our path integral reads:
α = ι κ C Ω exp ι α + ( C Ω ) exp ι α ;
κ = 1 2 π y X 2 ν ( t + t 0 )
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 ι α :
z = ( z ) 2 z + κ C Ω z 2 + ( C Ω ) ;
z ( 0 ) = z ( 1 ) = 1
This equation cannot be analytically solved for arbitrary periodic function C Ω ( ξ ) .
The weak and strong coupling expansions by κ are straightforward.
At small κ
z ( ξ ) 1 + 2 κ A ξ + 0 ξ Re C Ω ( ξ ) d ξ + O ( κ 2 ) ;
A = 0 1 Re C Ω ( ξ ) d ξ
At large κ
z ( ξ ) ι exp ι arg C Ω ( ξ ) = ι C Ω ( ξ ) C Ω ( ξ )
This solution is valid at intermediate ξ , not too close to the boundaries ξ = ( 0 , 1 ) . In the region near the boundaries ξ ( 1 ξ ) 1 κ , the following asymptotic agrees with the classical equation
z 1 ± ι ξ 2 κ Re C Ω ( 0 ) + O ( ξ 2 ) ;
z 1 ± ι ( 1 ξ ) 2 κ Re C Ω ( 1 ) + O ( ( 1 ξ ) 2 ) ;
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 [1], the viscosity ν = ν ˜ N 2 0 in our theory. This limit makes κ N , 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.

4. Vorticity Correlation

The simplest observable quantity we can extract from the loop functional is the vorticity correlation function [8], corresponding to the loop C backtracking between two points in space r 1 = 0 , r 2 = r , (see [1] 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 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 [1]
ω ( 0 ) · ω ( r ) = 1 4 t 2 O ( 3 ) d Ω | O ( 3 ) | 0 n < m < N ω m · ω n exp ι ρ · Ω ^ · S n , m S m , n E ;
S n , m = k = n m 1 F k m n ( mod N ) ;
ω k = 0 , 0 , ι σ k 2 cot β 2 ;
ρ = r 2 ν t ;
X [ σ . . . ] p , q , r E X [ σ . . . ] δ [ q r σ ] E δ [ q r σ ] ;
Integrating the global rotation matrix O ( 3 ) is part of the ensemble averaging.

4.1. Correlation Function and Path Integral

Let us apply our path integral to the expectation value over spins σ = ± 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:
ω ( 0 ) · ω ( r ) 1 t 2 Φ ( N ) even q < N p ; ( p | q ) cot 2 ( β / 2 ) β 2 0 < ξ 1 < ξ 2 < 1 d ξ 1 d ξ 2
O ( 3 ) d Ω [ D α ] α ( ξ 1 ) α ( ξ 2 ) exp ι r · Ω ^ · Im V ( ξ 1 , ξ 2 ) ν t | O ( 3 ) | [ D α ] ;
V ( ξ 1 , ξ 2 ) = q X ι , 1 , 0 S ( ξ 1 , ξ 2 ) S ( ξ 2 , 1 + ξ 1 ) ;
S ( a , b ) = a b d ξ exp ι α ( ξ ) b a ;
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 r = 0 computed in previous work [1].
The computations significantly simplify in Fourier space.
ω ( 0 ) · ω ( k ) = d 3 r ω ( 0 ) · ω ( r ) exp ι k · r
1 t 2 Φ ( N ) even q < N p ; ( p | q ) cot 2 ( β / 2 ) β 2 0 < ξ 1 < ξ 2 < 1 d ξ 1 d ξ 2
O ( 3 ) d Ω [ D α ] α ( ξ 1 ) α ( ξ 2 ) δ r · Ω ^ · Im V ( ξ 1 , ξ 2 ) ν t k | O ( 3 ) | [ D α ] ;
The angular integration d Ω reduces the three-dimensional delta function to one-dimensional
O ( 3 ) d Ω δ Ω ^ · Im V ( ξ 1 , ξ 2 ) ν t k ν t k 2 δ Im V ( ξ 1 , ξ 2 ) | k | ν t
Now, using the Lagrange multiplier λ for this condition, we have to minimize effective Action
A [ α , λ ] = π y 2 X 2 ξ 1 1 + ξ 1 ( α ) 2 + λ y X ξ 1 ξ 2 d ξ e ι α ξ 2 ξ 1 ξ 2 1 + ξ 1 d ξ e ι α 1 + ξ 1 ξ 2 ;
λ A [ α , λ ] = | k | ν t ;
This variational problem reduces to two pendulum equations
α 1 + r ξ 2 ξ 1 sin α 1 = 0 ; ξ 1 < ξ < ξ 2
α 2 + r ξ 2 ξ 1 1 sin α 2 = 0 ; ξ 2 < ξ < 1 + ξ 1
r = λ π y X I ( r ) ;
I ( r ) = ξ 1 ξ 2 d ξ e ι α 1 ξ 2 ξ 1 ξ 2 1 + ξ 1 d ξ e ι α 2 1 + ξ 1 ξ 2 ;
The well-known solution is Jacobi amplitude am ( x u ) ,
α 1 ( ξ ) = 2 am ξ α 2 2 a 1 r a 1 2 ( ξ 2 ξ 1 ) ;
α 2 ( ξ ) = 2 am ξ α 2 2 a 2 r a 2 2 ( 1 ξ 2 + ξ 1 ) ;
The free parameters a 1 , a 2 , α 1 , α 2 satisfy four equations
α 1 ( ξ 2 ) = α 2 ( ξ 2 ) ;
α 1 ( ξ 2 ) = α 2 ( ξ 2 ) ;
α 1 ( ξ 1 ) = α 2 ( 1 + ξ 1 ) ;
α 1 ( ξ 1 ) = α 2 ( 1 + ξ 1 ) ;
together with the constraint following from the variation of the Lagrange multiplier λ :
I ( r ) = | k | ν t y X

4.2. 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 , the remaining effective action at the extremum
N A [ α c , λ c ] = π N y 2 X 2 ξ 1 ξ 2 d ξ ( α 1 ( ξ ) ) 2 + ξ 2 1 + ξ 1 d ξ ( α 2 ( ξ ) ) 2
grows as N, unless both α 1 ( ξ ) α 2 ( ξ ) N 1 2 0 . In this case, the above constraint can be expanded in α 1 , α 2 . As we show in [14], the leading constant and linear terms both cancel so that the quadratic term remains
2 | k | ν t y X = ξ 1 ξ 2 d ξ α 1 2 ( ξ ) ξ 2 ξ 1 ξ 2 1 + ξ 1 d ξ α 2 2 ( ξ ) 1 + ξ 1 ξ 2 1 N
This estimate then requires vanishing viscosity in the local limit, at fixed turbulent viscosity
ν ˜ = ν N 2 const .
This phenomenon of renormalization of viscosity by a factor of N 2 was already observed in our first paper [1]. Our Euler ensemble in the local limit N can only solve the inviscid limit of the Navier-Stokes decaying turbulence, with finite ν ˜ 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 a 2 0 . This linearization is equivalent to replacing sin ( α ) α in the differential equation and studying the resulting linear ODE as a boundary problem. We choose different parametrizations in this linear case
α 1 ( ξ ) = a cos K 1 ( ξ ξ 2 ) + b K 1 sin K 1 ( ξ ξ 2 ) ;
α 2 ( ξ ) = a cos K 2 ( ξ ξ 2 ) + b K 2 sin K 2 ( ξ ξ 2 ) ;
K 1 = r Δ ;
K 2 = r Δ 1 ;
Δ = ξ 2 ξ 1 ;
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
b = r Δ 1 sin ( 1 Δ ) r Δ 1 + r Δ sin Δ r Δ cos ( Δ 1 ) r Δ 1 cos Δ r Δ
The remaining matching condition α 1 ( ξ 1 ) = α 2 ( 1 + ξ 1 ) reduces to the root of the function
g ( r , Δ ) = ( 2 Δ 1 ) sin ( Δ 1 ) r sin Δ r ( Δ 1 ) Δ + 2 cos ( Δ 1 ) r cos Δ r 2 r 2
This function has multiple roots, but we are looking for the root r 0 ( Δ ) with minimal value of the action
A c ( r , Δ ) = ξ 1 ξ 2 d ξ α 1 ( ξ ) 2 + ξ 2 1 + ξ 1 d ξ α 2 ( ξ ) 2
This integral is elementary, but the expression is too long to be presented here. It can be found in the Mathematica® notebook[14], where it is used to select the roots r 0 ( Δ ) of g ( r , Δ ) , minimizing A c r 0 ( Δ ) , Δ for a given value of Δ ( 0 , 1 ) .
This lowest action root is plotted in Figure 3. The corresponding value of minimal action L ( Δ ) = A c r 0 ( Δ ) , Δ is plotted in Figure 4.
The constraint (165) is also reduced to elementary functions, too lengthy to quote here (see [14]).
This constraint yields the quadratic relation for the last unknown parameter a in our solution
a 2 = | k | ν t | S ( Δ ) | y X ;
with universal function S ( Δ ) presented in [14] and shown in Figure 5).
The resulting integral (up to the pre-exponential factor Q) is equal to
[ D α ] α ( ξ 1 ) α ( ξ 2 ) δ Im V ( ξ 1 , ξ 2 ) | k | ν t
Q α ( ξ 1 ) α ( ξ 2 ) exp y X | k | ν ˜ t π L ( Δ ) 2 | S ( Δ ) |
The factor α ( ξ 1 ) α ( ξ 2 ) contains two terms :
α ( ξ 1 ) α ( ξ 2 ) = α 1 ( ξ 1 ) α 1 ( ξ 2 ) + δ α ( ξ 1 ) δ α ( ξ 2 )
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 [14])
α 1 ( ξ 1 ) α 1 ( ξ 2 ) = | k | ν ˜ t N y X J Δ | S ( Δ ) | ;
J Δ = r A ( r , Δ ) B ( r , Δ ) Δ ( Δ 1 ) C ( r , Δ ) 2 r = r 0 ( Δ ) ;
A ( r , Δ ) = Δ sin ( Δ 1 ) r ( Δ 1 ) Δ sin Δ r ;
B ( r , Δ ) = Δ sin ( Δ 1 ) r cos Δ r
( Δ 1 ) Δ sin Δ r cos ( Δ 1 ) r ;
C ( r , Δ ) = cos ( Δ 1 ) r cos Δ r ;
This function J ( Δ ) has an essential singularity at Δ = 1 due to oscillations of B ( r 0 ( Δ ) , Δ ) at a large real argument ( Δ 1 ) r 0 ( Δ ) at large negative r 0 ( Δ 1 ) . These oscillations lead to a sequence of zeros, condensing to Δ = 1 . The denominator of J ( Δ ) does not have zeros inside the ( 0 , 1 ) interval, so there are no poles. The function J ( Δ ) is real in physical domain because both factors A ( r 0 ( Δ ) , Δ ) , B ( r 0 ( Δ ) , Δ ) are purely imaginary and the factor C ( r 0 ( Δ ) , Δ ) is real (see Figure 6).
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 L ^ corresponding to linearized effective action (151) in the vicinity of the saddle point λ c , α c ( ξ ) .
A [ α c + δ α , λ c + δ λ ] A [ α c , λ c ] + 1 2 tr V L ^ V ;
V = δ λ , δ α ;
Q ( Δ , τ ) = exp 1 2 tr log L ^ L ^ α α 0 ;
τ = y X ;
The fluctuation correction reduces to the inverse operator L ^ , which we compute in the next section.
Now, we can reduce multiple sum/integral in (149) to the following
ω ( 0 ) · ω ( k ) = ν ˜ 3 2 H k ν ˜ t t ;
H ( κ ) = 1 Z n = 1 φ ( n ) 0 1 / n d τ τ 5 / n 5 τ 10 0 1 d Δ ( 1 Δ ) G Δ , τ , κ ; G Δ , τ , κ = Q Δ , τ
κ N τ J Δ | S ( Δ ) | + δ α ( ξ 1 ) δ α ( ξ 2 ) exp τ κ L ( Δ ) 2 π | S ( Δ ) |
where Z is the normalization constant to be determined later.

4.3. Functional Determinant in the Path Integral

As we have discussed in the previous section, in the limit a 0 the classical solution α 1 , 2 ( ξ ) a 0 .
This observation simplifies the linearized theory corresponding to this quadratic form V L ^ V . First, Let us integrate the fluctuations δ λ of λ around the saddle point solution.
The Lagrange multiplier at the saddle point vanishes, as we show in [14]
λ 0 = τ r 0 ( Δ ) I ( r 0 ( Δ ) ) = 0
The quadratic term comes from the first derivatives I λ = λ I , I r = r I , λ r = r λ , which can be simplified by switching to λ ( r ) = τ t I ( r )
A λ λ = τ I λ = τ I r λ r = τ I r τ r I r = 1 r
The bilinear term λ δ α also simplifies
A α λ ( δ α ) = ι τ δ λ F [ δ α ] ;
F [ δ α ] = ξ 1 ξ 2 d ξ δ α ( ξ ) ξ 2 ξ 1 ξ 2 1 + ξ 1 d ξ δ α ( ξ ) 1 + ξ 1 ξ 2
We can integrate out λ , producing the extra pre-exponential factor Q λ = | r 0 ( Δ ) | / N .
The bilinear term in the exponential after λ integration leads to the following effective quadratic Action for δ α
A e f f [ δ α ] = τ 2 2 δ α 2 + 1 2 r 0 ( Δ ) τ 2 F [ δ α ] 2 ;
The ( δ α ) 2 term corresponds to the linear eigenvalue equation with f 1 , 2 = δ α 1 , 2
f 1 , 2 ( ξ ) + μ 1 , 2 F [ f ] = ω f 1 , 2 ( ξ ) ;
μ 1 , 2 = r Δ , r Δ 1 ;
ϵ = ω τ 2 ;
r = r 0 ( Δ ) ;
The solution matching with first derivative at ξ = ξ 2 , ξ = ( ξ 1 , 1 + ξ 1 ) is built the same way as in (172). Equations for f 1 , 2 being linear homogeneous, we can fix the normalization as F [ f ] = 1 ,
f 1 ( x ) = a sin ( ω ( ξ ξ 2 ) ) + B 1 cos ( ω ( ξ ξ 2 ) ) + r ω Δ ;
f 2 ( ξ ) = a sin ( ω ( ξ ξ 2 ) ) + B 1 cos ( ω ( ξ ξ 2 ) ) + r ω ( Δ 1 ) ;
The spectrum ω = ω n is defined by the transcendental equation (the discriminant of this linear system of equations), which we found in [14]
f ( ω n , Δ ) = 0 ;
f ( ω , Δ ) = ( Δ 1 ) Δ ω ( ( Δ 1 ) Δ ω + r ) sin r ω 2 r ω +
cos 1 2 ( 1 2 Δ ) r ω cos r ω 2 ;
r = r 0 ( Δ ) ;
The spectrum is positive in the whole interval 0 < Δ < 1 except for some regions where the real solution for ω n ( Δ ) ends.
The functional determinant, resulting from the WKB approximation to the α path integral, would be related to the infinite product of eigenvalues ϵ n = τ 2 ω n , which can be written using a contour integral
Q α Δ , τ = ω n 0 τ 2 ω n 1 2 = exp 1 2 Im Γ f ( ω ) f ( ω ) log ( ω τ 2 ) d ω 2 π ω α α 0 ;
and the integration contour Γ encircles anticlockwise the nonzero poles of the meromorphic function f ( ω ) / f ( ω ) . As these poles are all located in the positive real axis, the contour Γ can be redirected around the negative real axis clockwise, where the integrand decreases as log ( ω ) ω 1 2 α . The integral converges at α > 1 2 and should be analytically continued to α = 0 .
Here is how this integral can be analytically continued to Re α < 1 / 2 . Consider the asymptotic of f ( ω ) f ( ω ) at Re ω
f ( ω ) f ( ω ) r 0 ( Δ ) 4 ω + 3 2 ω + O ω 3 2 ;
This asymptotic term represents a holomorphic function inside the contour Γ (which surrounds the negative real axis).
Therefore, this term yields zero at Re α > 1 / 2 and we can subtract it
Γ f ( ω ) f ( ω ) log ( ω τ ) d ω 2 π ω α = Γ f ( ω ) f ( ω ) r 0 ( Δ ) 4 ω 3 2 ω log ( ω τ ) d ω 2 π ω α
This integral already converges at Re α > 1 2 , so that we can set α = 0 :
Q α Δ , τ = exp 1 2 Im Γ f ( ω ) f ( ω ) r 0 ( Δ ) 4 ω 3 2 ω log ( ω τ 2 ) d ω 2 π
The remarkable property of this functional determinant is the factorization of the τ dependence
Q α ( Δ , τ ) = τ μ ( Δ ) Q α ( Δ , 1 ) ;
μ ( Δ ) = Im Γ f ( ω ) f ( ω ) r 0 ( Δ ) 4 ω 3 2 ω d ω 2 π
Our result for the correlation function is given by (188) with
Q Δ , τ = Q α Δ , 1 τ μ ( Δ ) r 0 ( Δ )
and Q α ( Δ , τ ) , μ ( Δ ) given by (208),(210). 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.

4.4. The Fluctuation Term in α ( ξ 1 ) α ( ξ 2 )

The last missing term is the fluctuation contribution to α ( ξ 1 ) α ( ξ 2 ) . In the Gaussian approximation, valid at N , this term equals
δ α ( ξ 1 ) δ α ( ξ 2 ) = 1 N τ 2 ξ ξ G ( ξ , ξ ) ξ = ξ 1 , ξ = ξ 2
where G ( ξ , ξ ) is a resolvent for the effective quadratic Action (195). This resolvent satisfies the equation
ξ 2 G ( ξ , ξ ) + μ ( ξ ) F [ G ] = δ ( ξ ξ ) ;
G ( ξ 1 , ξ ) = G ( ξ 1 + 1 , ξ ) = 0 ;
F [ G ] = ξ 1 ξ 2 d ξ G ( ξ , ξ ) ξ 2 ξ 1 ξ 2 1 + ξ 1 d ξ G ( ξ , ξ ) 1 + ξ 1 ξ 2 ;
μ ( ξ ) = r Δ if ξ 1 ξ < ξ 2 r Δ 1 if ξ 2 ξ < 1 + ξ 1
The solution of this equation, matching with the first derivative at ξ = ξ 2 is
G ( ξ , ξ ) = A + | ξ ξ | 2 + B ( ξ ξ 2 ) F [ G ] r ( ξ ξ 2 ) 2 Δ if ξ 1 ξ < ξ 2 A + | ξ ξ | 2 + B ( ξ ξ 2 ) F [ G ] r ( ξ ξ 2 ) 2 Δ 1 if ξ 2 ξ < 1 + ξ 1
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 [14] is too lengthy to present here. The desired quantity (212) is quite simple
1 N τ 2 ξ ξ G ( ξ , ξ ) ξ = ξ 1 , ξ = ξ 2 = 2 ( r + 3 ) N τ 2 ( 6 r )
Finally, we get the function G ( Δ , τ , κ ) in (188) (absorbing the constant factors in Z )
G Δ , τ , κ = Q α Δ , 1 τ μ ( Δ ) 2 r 0 ( Δ ) τ κ J Δ | S ( Δ ) | + 2 ( r 0 ( Δ ) + 3 ) ( 6 r 0 ( Δ ) ) exp τ κ L ( Δ ) 2 π | S ( Δ ) |

5. The Decaying Energy in Finite System

The vorticity correlation in Fourier space doubles as an energy spectrum
E ( k , t ) = 4 π k 2 v · v k = 4 π ω · ω k
The energy spectrum in a finite system with size L is bounded from below. At low | k | π / 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 [4] 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)
E ( k , t 0 ) P k 2 E ( k , t 0 ) M k 4
At 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 would apply to this regime. This solution corresponds to zero net momentum, which leaves the second regime with k 4 spectrum at small k and some universal decay at large k, reflecting these distributed vortex structures.
Therefore, the decaying energy, given by the part of the spectrum k > k 0 1 / L , has the following structure
E ( t ) = k 0 d k E ( k , t ) = 4 π ν ˜ t k 0 ν ˜ t H ( κ ) d κ
On top of the trivial decrease ν ˜ t , as prescribed by dimensional counting in an infinite system, there is some extra decrease related to the increase of the lower limit.
The detailed investigation of this complex singular representation deserves a separate publication. The asymptotic behavior of the above integral at t is determined by singularities of the functions
r 0 ( Δ ) , L ( Δ ) , J ( Δ ) , Q ( Δ ) , μ ( Δ )
at Δ = 1 . The powerlike singularities would result in a scaling law E ( t ) t n , n > 1 .
The energy decay can be investigated in more detail by considering the energy dissipation rate
E ( t ) = 4 π ν d k k 2 ω ( 0 ) · ω ( k ) = 4 π ν ˜ t 2 k 0 ν ˜ t κ 2 H ( κ ) d κ
Unlike the net remaining energy E ( t ) , the energy dissipation rate has a finite limit in an infinite system in our theory ( k 0 = 0 ). This limit was computed in [1] in a slightly different grand canonical ensemble, where N was fluctuating with the weight exp μ N , μ 0 .
With our current ensemble of fixed even N the results of [1] read:
E ( t ) = ν ˜ t 2 π 2 864 ζ ( 3 ) ;
In our present theory, the same quantity is given by the above integral at k 0 = 0
E ( t ) = 4 π ν ˜ t 2 0 κ 2 H ( κ ) d κ
Comparing these two expressions, we get the normalization of H ( κ )
4 π 0 κ 2 H ( κ ) d κ = π 2 864 ζ ( 3 )
The integral on the left can be further simplified [14] to the following normalization condition:
Z = 69120 ζ ( 3 ) 0 1 d Δ ( 1 Δ ) | r 0 ( Δ ) | Q α ( Δ , 1 ) | S ( Δ ) | ζ ( μ ( Δ ) + 7 ) π J ( Δ ) 6 r 0 ( Δ ) + L ( Δ ) r 0 ( Δ ) + 3 ( μ ( Δ ) + 3 ) ( μ ( Δ ) + 8 ) L ( Δ ) 2 6 r 0 ( Δ ) ζ ( μ ( Δ ) + 8 )
This normalization constant Z can be used in equation (222) for the energy decay in a finite system. All the functions of Δ were defined above.
Asymptotically, at ν ˜ t , the spectrum reduces to the integral of the power term (with Q ( Δ ) = μ ( Δ ) + 4 )
E ( k , t ) 1 t 0 1 d Δ P ( Δ ) k ν ˜ t Q ( Δ ) ;
E ( t ) 1 t 0 1 d Δ P ( Δ ) Q ( Δ ) 1 k 0 ν ˜ t 1 Q ( Δ )
The spectrum will, in general, decay as
E ( k , t ) 1 t log ( k ν ˜ t ) k ν ˜ t min Q ( Δ ) ;
E ( t ) 1 t log ( k ν ˜ t ) k 0 ν ˜ t 1 min Q ( Δ ) ;
There is also a possibility of a lack of the minimum min Q ( Δ ) inside the interval ( 0 , 1 ) , in which case the asymptotic of the integral will come from one of the endpoints.
This problem is left for future investigation.

6. Summary

Here are the new results reported in this paper.
  • We found the continuum limit of distribution of scaling variables (125) in the small Euler ensemble.
  • We reduced the Markov process for the Euler ensemble in its fermionic representation (61) to the path integral.
  • This path integral in the continuum limit (97) is dominated by a complex classical trajectory (instanton), satisfying a nonlinear ODE (132).
  • The continuum limit of this solution, N , corresponds to the inviscid limit of the decaying turbulence in the Navier-Stokes equation. Effective turbulent viscosity is ν ˜ = ν N 2 const .
  • 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 contribution of the instanton to the vorticity correlation function by using zeta regularization of the functional determinant of this linearized operator.
  • We derived an analytic formula for energy decay in finite system (222), (223). The singular integrals involved would require new analysis to obtain critical behavior at large times.
  • The turbulent viscosity ν ˜ is the only internal parameter of our theory, not counting the minimal wavelength k 0 , parametrizing the external forcing at the boundary. The shape of the energy spectrum and decay is universal and, in principle, calculable in quadratures.
  • This amazing simplification of the turbulence problem became possible after reducing the Hopf equation for the decaying turbulence to the quasiclassical limit of the Schrödinger equation in the loop space [2]. The solution of this equation was reduced to the quantum statistical ensemble of Fermi particles on a ring and solved in this paper.

7. Remaining Problems

  • We performed all the calculations up to numerical factors, which we recovered in the end from the known ω 2 . It would be useful to compute all the functions, including normalization, and thus double-check the solution.
  • The spectrum of the fluctuations around the classical solution is determined by transcendental equation (202), see Figure 7. Detailed investigation of this spectrum and corresponding complex integral (208) is needed.
  • The next remaining problem would be the asymptotic behavior of the effective index n ( t ) = t t log E ( t ) . Is there a calculable limit n ( ) in our solution (222)?
  • The spectrum of decay indexes for deviations from our fixed trajectory [1] can be evaluated in the scaling limit, with finite ν ˜ .
  • Another solvable problem is the PDF of the velocity circulation
    P [ Γ , C ] = d γ Ψ [ γ , C ] exp ι γ Γ ν
    The γ dependence is absent in the leading term of the solution, so the correction term must be computed, like in [1], with the corresponding decay index at a large time. The PDF of the velocity circulation was measured in [15,16] and explained as an instanton in [8]. Computing it from first principles would be an important achievement.

Acknowledgments

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 helpful. They helped me understand the physics of decaying turbulence in a finite system and match my solution with the DNS data.
This research was supported by a Simons Foundation award ID 686282 at NYU Abu Dhabi. The computations were done on the High-Performance Computing resources at New York University Abu Dhabi.

Funding

The author is supported by a Simons Foundation award ID 686282.

Data Availability Statement

The Mathematica®notebooks used to verify equations and compute some functions are available for download in [13,14].

References

  1. Migdal, A. To the Theory of Decaying Turbulence. Fractal and Fractional 2023, 7, 754. [Google Scholar] [CrossRef]
  2. Migdal, A. Loop Equation and Area Law in Turbulence. In Quantum Field Theory and String Theory; Baulieu, L., Dotsenko, V., Kazakov, V., Windey, P., Eds.; Springer US, 1995; pp. 193–231. [Google Scholar] [CrossRef]
  3. Basak, D.; Zaharescu, A. Connections between Number Theory and the theory of Turbulence. 2024; To be published. [Google Scholar]
  4. Panickacheril John, J.; Donzis, D.A.; Sreenivasan, K.R. Laws of turbulence decay from direct numerical simulations. Philos. Trans. A Math. Phys. Eng. Sci. 2022, 380, 20210089. [Google Scholar] [CrossRef] [PubMed]
  5. Migdal, A. Topological Vortexes, Asymptotic Freedom, and Multifractals. Fractal Fract. 2023, 7, 351. [Google Scholar] [CrossRef]
  6. Matsuzawa, T.; Irvine, W. Realization of Confined Turbulence Through Multiple Vortex Ring Collisions, 03/12/2019. “Talk at the Flatiron Conference Universality Turbulence Across Vast Scales”. https://www.quantamagazine.org/an-unexpected-twist-lights-up-the-secrets-of-turbulence-20200903/.
  7. Matsuzawa, T.; Mitchell, N.P.; Perrard, S.; Irvine, W.T. Creation of an isolated turbulent blob fed by vortex rings. Nature Physics 2023, 19, 1193–1200. [Google Scholar] [CrossRef]
  8. Migdal, A. Statistical Equilibrium of Circulating Fluids. Physics Reports 2023, 1011C, 1–117. [Google Scholar] [CrossRef]
  9. AGISHTEIN, M.; MIGDAL, A. SIMULATIONS OF FOUR-DIMENSIONAL SIMPLICIAL QUANTUM GRAVITY AS DYNAMICAL TRIANGULATION. Modern Physics Letters A 1992, 7, 1039–1061. [Google Scholar] [CrossRef]
  10. Hardy, G.H.; Wright, E.M. An Introduction to the Theory of Numbers, 6th ed.; Oxford University Press: Oxford, UK, 2008; pp. xxii+621, Revised by D. R. Heath-Brown and J. H. Silverman, With a foreword by Andrew Wiles. [Google Scholar]
  11. Bulatov, M.; Migdal, A. Dual Theory of Decaying Turbulence. 2. Numerical Simulations.
  12. Norris, J.R. Markov Chains; Cambridge Univ. Press, 2007. [Google Scholar]
  13. Migdal, A. “BernSum”. 2024. https://www.wolframcloud.com/obj/sasha.migdal/Published/BernSum.nb.
  14. Migdal, A. “BernSum”. 2024. https://www.wolframcloud.com/obj/sasha.migdal/Published/JacobiEllipticForDecayingTurbulence.nb.
  15. Iyer, K.P.; Sreenivasan, K.R.; Yeung, P.K. Circulation in High Reynolds Number Isotropic Turbulence is a Bifractal. Phys. Rev. X 2019, 9, 041006. [Google Scholar] [CrossRef]
  16. Iyer, K.P.; Bharadwaj, S.S.; Sreenivasan, K.R. Area rule for circulation and minimal surfaces in three-dimensional turbulence. 2020; arXiv:physics.flu-dyn/2007.06723. [Google Scholar]
Figure 1. The plot of the function Λ ( s ) . As required, it is positive, takes a maximal value log ( 2 ) at s = 0 , and vanishes at both ends s = ± 1 of the physical region.
Figure 1. The plot of the function Λ ( s ) . As required, it is positive, takes a maximal value log ( 2 ) at s = 0 , and vanishes at both ends s = ± 1 of the physical region.
Preprints 100819 g001
Figure 2. The log-log-plot of the distribution f X ( X ) . It is equals π Φ ( k ) X X at each internal 1 π 2 ( k + 1 ) 2 < X < 1 π 2 k 2 with positive integer k. Asymptotically f X ( X ) 3 X π 3 at X 0 .
Figure 2. The log-log-plot of the distribution f X ( X ) . It is equals π Φ ( k ) X X at each internal 1 π 2 ( k + 1 ) 2 < X < 1 π 2 k 2 with positive integer k. Asymptotically f X ( X ) 3 X π 3 at X 0 .
Preprints 100819 g002
Figure 3. Log plot of r 0 ( Δ ) .
Figure 3. Log plot of r 0 ( Δ ) .
Preprints 100819 g003
Figure 4. Log plot of L ( Δ ) = A c r 0 ( Δ ) , Δ .
Figure 4. Log plot of L ( Δ ) = A c r 0 ( Δ ) , Δ .
Preprints 100819 g004
Figure 5. Plot of S Δ .
Figure 5. Plot of S Δ .
Preprints 100819 g005
Figure 6. Plot of universal function J ( Δ ) .There are no singularities inside the physical region 0 < Δ < 1 , but there are zeros, condensing at Δ = 1 .
Figure 6. Plot of universal function J ( Δ ) .There are no singularities inside the physical region 0 < Δ < 1 , but there are zeros, condensing at Δ = 1 .
Preprints 100819 g006
Figure 7. The first levels of the spectrum satisfying equation f ( ω , Δ ) = 0 . The broken lines correspond to the end of zones of real solutions.
Figure 7. The first levels of the spectrum satisfying equation f ( ω , Δ ) = 0 . The broken lines correspond to the end of zones of real solutions.
Preprints 100819 g007
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.
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.

Downloads

870

Views

421

Comments

0

Subscription

Notify me about updates to this article or when a peer-reviewed version is published.

Email

Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2025 MDPI (Basel, Switzerland) unless otherwise stated