
\(\nu ^k_{t,x} \in \mathcal {P}(U^k)\)  for a.e. \((t,x)\in [0,T]\times D^k\) , and the map \(x \mapsto \bigl \langle \nu ^k_{t,x}, f\bigr \rangle \)  is measurable for every \(f\in C_b(U^k)\)  and almost every \(t\in [0,T]\) . (In other words, \(\nu ^k_t\)  is a Young measure from \(D^k\)  to \(U^k\) .)

\(L^p\)  integrability:
\(\operatornamewithlimits{ess\, sup}_{t\in [0,T)}\left(\int _{D} \bigl \langle \nu ^1_{t,x}, |\xi |^p\bigr \rangle \,dx\right)^{1/p} \leqslant c< +\infty \) 

Diagonal continuity (DC):
\(\int _0^{T^{\prime }}\omega _r^p\big (\nu ^2_t\big )\,dt \rightarrow 0 \qquad \text{ as } r\rightarrow 0 \text{ for all } T^{\prime }\in (0,T),\) 
where
\(\omega _r^p(\nu ^2_t) \int _D \mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







Br(x) 2t,x,y, |1-2|p dy dx
\(is called the \emph {modulus of continuity} of \) 
\(.\) 
We denote the set of all time-dependent correlation measures by \({L}^{p}([0,T),D;U)\) .
In [1]} (and see [2]} for a time-independent version), the following equivalence between time-dependent correlation measures and parametrized probability measures on \(L^p(D)\)  was proved:
Theorem 2.2 
For every time-dependent correlation measure \(\mathbf {\nu }\in {L}^p([0,T),D;U)\)  there is a unique (up to subsets of \([0,T)\)  of Lebesgue measure 0) map \(\mu \colon [0,T)\rightarrow \mathcal {P}(L^p(D;U))\)  such that

the map
\(t\mapsto \bigl \langle \mu _t, L_g\bigr \rangle  = \int _{L^p} \int _{D^k}g(x,u(x))\,dx\,d\mu _t(u)\) 
is measurable for all \(g\in {H}^k_0(D;U)\) ,

\(\mu \)  is \(L^p\) -bounded:
\(\operatornamewithlimits{ess\, sup}_{t\in [0,T)}\int _{L^p}\Vert u\Vert _{L^p}^p\,d\mu _t(u) \leqslant c^p < \infty \) 

\(\mu \)  is dual to \(\mathbf {\nu }\) : the identity
\(\int _{D^k}\bigl \langle \nu _t^k, g(x)\bigr \rangle \,dx = \int _{L^p} \int _{D^k}g(x,u(x))\,dx\,d\mu _t(u)\) 
holds for a.e. \(t\in [0,T)\) , every \(g\in {H}^k_0(D;U)\)  and all \(k\in \mathbb {N}\) .

Conversely, for every \(\mu \colon [0,T)\rightarrow \mathcal {P}(L^p(D;U))\) 
satisfying (i) and (ii), there is a unique correlation measure \(\mathbf {\nu }\in {L}^p([0,T),D;U)\)  satisfying (iii).

We also have the following “Compactness” Theorem for time-dependent correlation measures [1]}
Theorem 2.3 
Let \(\mathbf {\nu }_n \in {L}^{p}([0,T),D;U)\)  for \(n=1,2,\dots \)  be a sequence of correlation measures such that
\(\sup _{n\in \mathbb {N}}\operatornamewithlimits{ess\, sup}_{t\in [0,T)}\left(\int _{D} \bigl \langle \nu ^1_{n;t,x}, |\xi |^p\bigr \rangle \,dx\right)^{1/p} \leqslant c< +\infty  \\\lim _{r\rightarrow 0}\limsup _{n\rightarrow \infty }\int _0^{T^{\prime }}\omega _r^p\bigl (\nu ^2_{n,t}\bigr )\,dt = 0 \) 

for some \(c>0\)  and all \(T^{\prime }\in [0,T)\) . Then there exists a subsequence \((n_j)_{j=1}^\infty \)  and some \(\mathbf {\nu }\in {L}^{p}([0,T),D;U)\)  such that

 \(\mathbf {\nu }_{n_j} \overset{*}{\rightharpoonup }\mathbf {\nu }\)  as \(j\rightarrow \infty \) , that is, \(\bigl \langle \nu ^k_{n_j}, g\bigr \rangle _{{H}^k} \rightarrow \bigl \langle \nu ^k, g\bigr \rangle _{{H}^k}\)  for every \(g\in {H}^k_0([0,T),D;U)\)  and every \(k\in \mathbb {N}\)

 \(\bigl \langle \nu ^1_t, |\xi |^p\bigr \rangle _{{H}^1} \leqslant c^p\)  for a.e. \(t\in [0,T)\)

 \(\int _0^{T^{\prime }}\omega _r^p\big (\nu ^2_t\big )\,dt \leqslant \liminf _{n\rightarrow \infty } \int _0^{T^{\prime }}\omega _r^p\big (\nu ^2_{n,t}\big )\,dt\)  for every \(r>0\)  and \(T^{\prime }\in [0,T)\)

 for \(k\in \mathbb {N}\) , let \(\varphi \in L^1_{\textrm {loc}}([0,T)\times D^k)\)  and \(\kappa \in C(U^k)\)  be nonnegative, and let \(g(t,x,\xi ) \varphi (t,x)\kappa (\xi )\) . Then
\(\bigl \langle \nu ^k, g\bigr \rangle _{{H}^k} \leqslant \liminf _{j\rightarrow \infty } \bigl \langle \nu _{n_j}^k, g\bigr \rangle _{{H}^k}.\) 


Assume moreover that \(D\subset \mathbb {R}^d\)  is compact, \(T<\infty \)  and that \(\mathbf {\nu }_n\)  have uniformly bounded support, in the sense that
\(\Vert u\Vert _{L^p} \leqslant R \qquad \text{for $\mu ^n_t$-a.e.\ } u\in L^p(D;U) \text{ for every } n\in \mathbb {N}, ~ a.e~t \in (0,T), \) 
with \(\mu _t^n\in \mathcal {P}_T(L^p(D;U))\)  being dual to \(\mathbf {\nu }_n\) ,
then the following observables converge strongly:
\(\lim _{j\rightarrow \infty }\int _{D^k} \left|\int _0^T\left(\bigl \langle \nu ^k_{n_j;t,x}, g(t,x)\bigr \rangle  - \bigl \langle \nu ^k_{t,x}, g(t,x)\bigr \rangle \right)\,dt\right| \,dx = 0\) 
for every \(g \in {H}^{k,p}_1([0,T],D;U)\) .


Statistical solutions
The goal of this section is to show that the statistical solutions of Navier–Stokes as introduced by Foiaş and Prodi [4]}, [5]}, [6]}, [7]}, [8]} are equivalent to families of correlation measures as introduced in [2]} that satisfy the Friedman–Keller system of moment equations. For the sake of simplicity we will assume that the support of the initial measure \(\mu _0\)  lies in a bounded set \(B\subset L^2(D;U)\) , that is,
\(\text{supp}(\mu _0)\subset B\subset L^2(D;U).\) 
The Leray projector
We recall first that the Helmholtz–Leray projector, or simply Leray projector, is the linear map \(\mathbb {P}\colon L^2(D;U)\rightarrow L^2_{\operatorname{div}}(D;U)\lbrace v\in L^2(D;U): \operatorname{div}v = 0,\,\int v dx =0\rbrace \)  that projects a vector field \(f\in L^2(D;U)\)  to its divergence free component, that is \(f=\mathbb {P}f +\nabla \psi _f\)  with \(\operatorname{div}(\mathbb {P}f)=0\)  and \(\mathbb {P}f, \nabla \psi _f\in L^2(D;U)\) . One can show that \(\nabla \psi _f\)  is orthogonal in \(L^2(D;U)\)  to any function \(u\in L^2_{\operatorname{div}}(D;U)\) ,
\(\int _{D}u \cdot \nabla \psi _f\,dx = 0.\) 
For functions in the tensor product space \(L^2(D^k;U^k)\)  we let \(\mathbb {P}_{x_i}\)  denote the Leray projector in the \(i\) th component, i.e., \(\varphi = \mathbb {P}_{x_i}\varphi +\nabla _{x_i}\psi _{\varphi ,i}\)  where \(\operatorname{div}_{x_i}(\mathbb {P}_{x_i}\varphi ) = 0\) , and \(\psi _{\varphi ,i}\in L^2(D^k; U^{k-1})\)  with \(\nabla _{x_i}\psi _{\varphi ,i}\in L^2(D^k;U^k)\) .

Definitions
We will start by recalling the different definitions of statistical solutions introduced in [6]}, [11]}, [2]}.
Definition 3.1 (Definition 3.2 in [2]})  Let \(\varepsilon \geqslant 0\) . The Friedman–Keller system of moment equations, defined for time-dependent correlation measures \(\mathbf {\nu }\in {L}^2([0,T),D;U)\) , is the hierarchy of equations
\(\begin{split}&\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\frac{\partial \varphi }{\partial t}(t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt \\&+\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\varphi (0,x)\,d\nu _{0,x}^k(\xi ) \,dx \\&+\sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}(\xi _1\otimes \cdots \otimes (\xi _i\otimes \xi _i)\otimes \cdots \xi _k):\nabla _{x_i}\varphi (t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt\\=&\ -\varepsilon \sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\Delta _{x_i}\varphi (t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt\end{split}\) 
for all \(k\in \mathbb {N}\) , for all \(\varphi \in C^2_c([0,T)\times D^k;U^k)\)  with \(\operatorname{div}_{x_i}\varphi =0\)  for all \(i=1,\dots , k\) ,
along with the divergence constraint
\(\int _{D^{k}}\int _{U^{k}} \xi _1\otimes \dots \otimes \xi _{\ell }\otimes \alpha _{\ell +1}(\xi _{\ell +1})\otimes \dots \otimes \alpha _k(\xi _k) \,d\nu _{t,x}^k(\xi )\cdot \nabla _{x_1,\dots , x_\ell }\varphi (x) \,dx = 0,\) 
where \(\nabla _{x_1,\dots ,x_\ell }=(\nabla _{x_{1}},\dots ,\nabla _{x_\ell })^\top \) , \(1\leqslant \ell \leqslant k\in \mathbb {N}\) , for all \(\varphi \in H^1(D^k;U^{k-\ell })\) , \(\alpha _j\in C(U;U)\) , with \(\alpha _j(v)\leqslant C(1+|v|^2)\)  for all \(j=1,\dots , k\) .
If \(\mathbf {\nu }\)  solves the Friedman–Keller system of moments equations and in addition satisfies the energy inequality
\(\begin{aligned}&\sum _{k=0}^K a_k\int _{D^k}\int _{U^k}|\xi _1|^2\dots |\xi _k|^2 \,d\nu _{t,x}^k(\xi ) \,dx\\&+ 2\varepsilon \sum _{k=0}^K a_k\sum _{i=1}^k\sum _{j=1}^d \lim _{h\rightarrow 0}\frac{1}{h^2}\int _0^t\!\!\int _{D^{k}}\!\int _{U^{k+1}} \!\!\!|\xi _1|^2\dots |\xi _i - \xi _{k+1}|^2\dots |\xi _k|^2\,d\nu _{t,(x,x_i+h\mathbf {e}_j)}^{k+1}(\xi ,\xi _{k+1})\,dx \,ds\\&\leqslant \sum _{k=0}^K a_k\int _{D^k}\!\!\int _{U^{k}}|\xi _1|^2\dots |\xi _k|^2 \,d\nu _{0,x}^k(\xi )\,dx\end{aligned}\) 
for all \(K\in \mathbb {N}\)  and \(a_k\in \mathbb {R}\) , \(k=0,\dots , K\)  such that \(p_K(s)=\sum _{k=0}^K a_k s^k\)  is a nonnegative, nondecreasing polynomial for \(s\in [0,R]\)  for \(R\)  sufficiently large related to the support of the correlation measure (see (REF )),
then we call \(\mathbf {\nu }\)  a Friedman–Keller statistical solution of the Navier–Stokes (when \(\varepsilon >0\) ) or Euler (when \(\varepsilon =0\) ) equations.

Remark 3.2 Friedman–Keller statistical solutions are analogous to the definition of statistical solutions for hyperbolic systems of conservation laws as introduced in [2]}, [1]} for compressible flows.

Remark 3.3 
By a standard argument for weak solutions to continuity equations, the map \(t \mapsto \langle \nu ^k_{t,\cdot },\, \xi _1\otimes \cdots \otimes \xi _k\rangle \)  is weakly continuous for every \(k\in \mathbb {N}\) ; see e.g. [16]}.

Remark 3.4 (Formulation of (REF ) with non-divergence free test functions) 
Denoting \(\mathbb {P}_k\mathbb {P}_{x_1}\dots \mathbb {P}_{x_k}\)  (cf. Section REF ), we can replace the divergence-free test function \(\varphi \in C^2_c([0,T]\times D^k;U^k)\)  in (REF ) by \(\mathbb {P}_k\varphi \)  for an arbitrary \(\varphi \in C^2_c([0,T]\times D^k;U^k)\) . Using the fact that the Leray projection is self-adjoint and that, by (REF ),
\(\operatorname{div}_{x_i}\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k) \,d\nu _{t,x}^k(\xi ) = 0 \qquad \forall \ i=1,\dots , k,\quad \text{a.e. }\, (t,x)\in [0,T]\times D^k\) 
we observe that
\(\mathbb {P}_k \int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k) \,d\nu _{t,x}^k(\xi ) = \int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k) \,d\nu _{t,x}^k(\xi ), \quad \text{a.e. }\, t\in [0,T],\,x\in D^k.\) 
Therefore,
\(&\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\frac{\partial \mathbb {P}_k\varphi }{\partial t}(t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt\\= &\int _0^T\!\!\int _{D^k}\!\!\mathbb {P}_k\left(\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k)\,d\nu _{t,x}^k(\xi )\right):\frac{\partial \varphi }{\partial t}(t,x) \,dx \,dt \\= &\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k)\,d\nu _{t,x}^k(\xi ):\frac{\partial \varphi }{\partial t}(t,x) \,dx \,dt.\) 

Similarly,
\(\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\mathbb {P}_k\varphi (0,x)\,d\nu _{0,x}^k(\xi ) \,dx = \int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\varphi (0,x)\,d\nu _{0,x}^k(\xi ) \,dx,\) 
and
\(&\sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\Delta _{x_i}\mathbb {P}_k\varphi (t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt\\= &\sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\Delta _{x_i}\varphi (t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt,\) 

the last one being true due to the fact that the Laplacian and the Leray projection commute on the torus.
Hence, the weak formulation (REF ) can be rewritten as
\(\begin{split}&\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\frac{\partial \varphi }{\partial t}(t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt+\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\varphi (0,x)\,d\nu _{0,x}^k(\xi ) \,dx \\&+\sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}(\xi _1\otimes \cdots \otimes (\xi _i\otimes \xi _i)\otimes \cdots \xi _k):\nabla _{x_i}\varphi (t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt\\=&-\varepsilon \sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\Delta _{x_i}\varphi (t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt\\&+\sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}(\xi _1\otimes \cdots \otimes (\xi _i\otimes \xi _i)\otimes \cdots \xi _k):\nabla _{x_i}\nabla _{x_i}\psi _{\varphi ,i}(t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt\end{split}\) 
where \(\varphi \in C^2_c([0,T]\times D^k;U^k)\) . The terms
\(\sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}(\xi _1\otimes \cdots \otimes (\xi _i\otimes \xi _i)\otimes \cdots \xi _k):\nabla _{x_i}\nabla _{x_i}\psi _{\varphi ,i}(t,x) \,d\nu _{t,x}^k(\xi ) \,dx \,dt\) 
correspond to the pressure in the deterministic setting.

To define statistical solutions in the sense of Foiaş and Prodi, we need to introduce some notation.
We denote by \(L^2_{\operatorname{div}}(D;U)\)  the space of divergence free \(L^2(D;U)\) -vector fields and by \( the space of divergence free functions in \) H1(D;U)\( (these can be obtained as the closures of \) C(D;U){divu=0}\( in \) L2(D;U)\( and \) H1(D;U)\(, respectively, with suitable integral conditions:\begin{equation*}\begin{split}L^2_{\operatorname{div}}(D;U)& = \left\lbrace u\in L^2(\mathbb {T}^d)\ :\ \operatorname{div}u = 0,\ \int _{\mathbb {T}^d} u(x)\,dx=0\right\rbrace ,\\=\left\lbrace u\in H^1(\mathbb {T}^d)\ :\ \operatorname{div}u = 0,\ \int _{\mathbb {T}^d} u(x)\,dx=0\right\rbrace ,\end{split}\end{equation*}for periodic boundary conditions.We denote the \) L2\(-inner product by\begin{equation*}(u,v) = \int _{D} u(x) v(x) \,dx,\end{equation*}and for \) >0\( the \) H1\(-inner product by\begin{equation*}a(u,v)= \varepsilon \sum _{i=1}^d\int _{D}\frac{\partial u}{\partial x^i}\cdot \frac{\partial u}{\partial x^i} \,dx,\end{equation*}Define the Stokes operator \) A\( by\begin{equation*}\begin{split}Au &= -\mathbb {P}\Delta u, \quad \text{for all }u\in D(A)= H^2(D;U),\\\varepsilon (Au,v)&= a(u,v),\quad \text{for all } u,v\in D(A^{1/2}),\end{split}\end{equation*}where \) P\( is the Leray projector, and the skew-symmetric trilinear form \) b\( by\begin{equation}\begin{split}b(u,v,w) &\int _{D} (u\cdot \nabla ) v\cdot w \,dx = (B(u,v),w), \quad u,v,w\in D(A^{1/2})\\B(u)&B(u,u).\end{split}\end{equation}We can then write the Navier--Stokes equations in the functional formulation: Let \) T>0\(, \) u0L2div(D;U)\(, find \) uL([0,T];L2div(D;U))L2([0,T];\( with \) u'ddtuL1([0,T];D(A-1/2))\( such that\begin{equation}u^{\prime }+ \varepsilon Au+B(u)= 0\end{equation}and \) u(0)= u0\( in a suitable sense. This corresponds to the weak formulation\begin{equation}\frac{d}{dt}(u,v) +a(u,v)+b(u,u,v) =0 \qquad \text{for all } v\in \end{equation}If we denote\begin{equation}F(t,u)-\varepsilon A u -B(u),\end{equation}the functional formulation becomes\begin{equation}u^{\prime }(t)=F(t,u(t)).\end{equation}We need the following class of test functions:\)
Notation 3.5 [6]}
Let \(\mathcal {T}_{\text{cyl}}\)  denote the class of cylindrical test functions consisting of the real-valued functionals \(\Phi =\Phi (u)\)  that depend on a finite number \(k\in \mathbb {N}\)  of components of \(u\) , that is,
\(\Phi (u)= \varphi \bigl ((u,g_1),\dots , (u,g_k)\bigr ),\) 
where \(\varphi \in C^1_c(\mathbb {R}^k)\)  and \(g_1,\dots , g_k\in H^1(D;U)\) .
Let \(\mathcal {T}_{\text{cyl}}^0\)  denote the subset of such functions which satisfy \(g_1,\dots , g_k\in .We denote by \) '\( the differential of \)\( in \) L2div(D;U)\(, which can be expressed as\begin{equation*}\Phi ^{\prime }(u)= \sum _{j=1}^k \partial _j \varphi \bigl ((u,g_1),\dots ,(u,g_k)\bigr ) g_j,\end{equation*}where \) j \( is the derivative of \)\( with respect to its \) j\(th component.\) 
We can now define statistical solutions in the sense of Foiaş and Prodi. We will use the definition as it stated in their newer work [11]}:
Definition 3.6 (Foiaş–Prodi [6]}, [4]}, [7]}, [11]})
A family of probability measures \((\mu _t)_{0\leqslant t\leqslant T}\)  on \(L^2_{\operatorname{div}}(D;U)\)  is a Foiaş–Prodi statistical solution of the Navier–Stokes equations on \(L^2_{\operatorname{div}}(D;U)\)  with initial data \(\mu _0\) 
if

The function
\(t\mapsto \int _{L^2_{\operatorname{div}}} \varphi (u)\,d\mu _t(u),\) 
is measurable on \([0,T]\)  for every \(\varphi \in C_b(L^2_{\operatorname{div}}(D;U))\) ;

\(\mu \)  satisfies the weak formulation
\(\int _{L^2_{\operatorname{div}}} \Phi (u) \,d\mu _t(u)=\int _{L^2_{\operatorname{div}}}\Phi (u) \,d\mu _0(u) + \int _0^t\int _{L^2_{\operatorname{div}}} (F(s,u),\Phi ^{\prime }(u)) \,d\mu _s(u) \,ds\) 
for all \(t\in [0,T]\)  and all cylindrical test functions \(\Phi \in \mathcal {T}_{\text{cyl}}^0\) , where \(F\)  is given in ().

\(\mu \)  satisfies the strengthened mean energy inequality:
For any \(\psi \in C^1(\mathbb {R},\mathbb {R})\)  nonnegative, nondecreasing with bounded derivative and \(t\in [0,T]\) , the inequality
\(\int _{L^2_{\operatorname{div}}} \psi (\Vert u\Vert _{L^2(D)}^2) \,d\mu _t(u) + 2\varepsilon \int _0^t\int _{L^2_{\operatorname{div}}} \psi ^{\prime }(\Vert u\Vert _{L^2(D)}^2) |u|_{H^1(D)}^2 \,d\mu _s(u)\\\leqslant \int _{L^2_{\operatorname{div}}} \psi (\Vert u\Vert _{L^2(D)}^2 )\,d\mu _{0}(u)\) 
holds.

The function
\(t\mapsto \int _{L^2_{\operatorname{div}}}\psi (\left\Vert  u \right\Vert _{L^2(D)}^2)\, d\mu _t(u)\) 
is continuous at \(t=0\)  from the right, for any function \(\psi \in C^1(\mathbb {R},\mathbb {R})\)  nonnegative, nondecreasing with bounded derivative.


Remark 3.7 
Note that, as a consequence of the energy inequality (REF ) for \(\psi (s)=s\) , the function
\(t\mapsto \int _{L^2_{\operatorname{div}}}\Vert u\Vert _{L^2(D)}^2 \,d\mu _t(u),\) 
belongs to \(L^\infty ([0,T])\)  and the function
\(t\mapsto \int _{L^2_{\operatorname{div}}}|u|_{H^1(D)}^2 \,d\mu _t(u),\) 
belongs to \(L^1([0,T])\) . Notice also that (REF ) implies that
\(t\mapsto \int _{L^2_{\operatorname{div}}}\Phi (u)d\mu _t(u)\) 
for \(\Phi (u)\)  a cylindrical test function, is continuous since
\(\int _{L^2_{\operatorname{div}}} (F(s,u),\Phi ^{\prime }(u)) \,d\mu _t(u)\) 
is locally integrable. Combining this fact with condition (c), conditionn (d) follows directly.

Equivalence between the solution concepts
Next, we show that the Friedman–Keller statistical solutions in Definition REF  and the Foiaş–Prodi statistical solutions in Definition REF  are in fact the same.
Theorem 3.8 (Foiaş–Temam statistical solutions satisfy the Friedman–Keller system) 
Let \(\mu \)  be a Foiaş–Prodi statistical solution such that the initial condition \(\mu _0\)  has bounded support,
\(\text{supp}(\mu _0)\subset B\subset L^2_{\operatorname{div}}(D;U), \qquad B\subset \big \lbrace u\in L^2(D;U)\ :\ \Vert u\Vert _{L^2(D)}\leqslant R\big \rbrace \) 
for some \(R>0\) . Then \(\mu \)  corresponds (cf. Theorem REF ) to a correlation measure \(\mathbf {\nu }_t=(\nu ^1,\nu ^2,\dots )\)  that is a statistical solution in the Friedman–Keller sense (cf. Definition REF ).

Conversely, we have:
Theorem 3.9 (Friedman–Keller solutions are Foiaş–Prodi statistical solutions) 
Let \(\mathbf {\nu }=(\mathbf {\nu }_t)_{0\leqslant t\leqslant T}\)  be a Friedman–Keller statistical solution of Navier–Stokes (cf. Definition REF ) with bounded support, i.e.,
\(\int _{D^k}\int _{U^{k}} |\xi _1|^2\dots |\xi _k|^2 \,d\nu _{t,x}^k(\xi ) \,dx\leqslant R^k<\infty ,\) 
for some \(0<R<\infty \) , every \(k\in \mathbb {N}\) , and almost every \(t\in [0,T]\) . Then \(\mathbf {\nu }\)  corresponds to a probability measure \(\mu =(\mu _t)_{0\leqslant t\leqslant T}\)  on a bounded set of \(L^2_{\operatorname{div}}(D;U)\)  which is a Foiaş–Prodi statistical solution of the Navier–Stokes equations (cf. Definition REF ).

The proofs of these two results are given in Appendix .
Foiaş et al. have shown existence of Foiaş–Prodi statistical solutions for the (forced) Navier–Stokes equations, see e.g. [4]}, [5]}, [6]}. Using these equivalence theorems, this implies existence of statistical solutions via correlation measures as in Definition REF .
Remark 3.10 The equivalence theorems REF  and REF  are restricted to probability measures with bounded support. It should be possible to extend these results to probability measures having sufficiently fast decay near infinity; however, the proofs would become significantly more technical. We have therefore decided to restrict ourselves to probability measures with bounded support.


Vanishing viscosity limit of statistical solutions of Navier–Stokes
The goal of this section is to pass to the inviscid limit \(\varepsilon \rightarrow 0\)  under the assumption of weak statistical scaling (c.f. Section REF , Assumption REF ). We will first prove a rigorous result on the longitudinal third order structure function
\(S_\Vert ^3(\tau ,r)\int _0^\tau \int _{L^2_{\operatorname{div}}}\!\!_{\mathbb {S}^2}\int _{D}\big ((u(x+r n)-u(x))\cdot n\big )^3 \,dxdS(n) \,d\mu _t(u)\,dt,\) 
and then relate it to the similarly defined second order structure function using the weak scaling assumption. Together with weak statistical anisotropy, this yields diagonal continuity of the correlation measures \(\mathbf {\nu }^\varepsilon \)  that is needed to apply the compactness theorem REF  and pass to the limit. The proof of the scaling estimate for the third order structure function (REF ) in Lemma REF  and  REF  largely follows the proof of a similar result for martingale solutions of stochastic Navier–Stokes equations in [26]}. To simplify notation, we will omit writing the dependence of \(\mathbf {\nu }\)  and \(\mu \)  on \(\varepsilon \)  in the following sections.
Kármán–Howarth–Monin relation
The key to deriving an estimate on the behavior of the third order structure function (REF ) is the so-called Kármán–Howarth–Monin (KHM) relation [27]} that describes the evolution of the second correlation marginal. Similar relations have been derived before for various settings (stochastic, forced, etc.), see [28]}, [29]}, [27]}, [31]}, [26]}, [33]}, [34]}. For statistical solutions we derive:
Proposition 4.1 
Let \(\mathbf {\nu }\)  be a Friedman–Keller statistical solution of the Navier–Stokes equations. Then the second correlation marginal \(\nu ^{2}\)  satisfies the Kármán–Howarth–Monin relation for correlation measures:
\( \begin{split}&\quad \sum _{ij}\int _D\!\int _D\!\int _{U^2}\!\!\xi _1^i \xi _2^j\,d\nu _{\tau ,x,x+h}^2(\xi ) \,dx\, \sigma ^{ij}(h) \,dh-\sum _{ij}\int _{D}\!\int _D\!\int _{U^2}\!\!\xi _1^i \xi _2^j\,d\nu _{0,x,x+h}^2(\xi ) \,dx \,\sigma ^{ij}(h) \,dh \\&\quad +\frac{1}{2}\sum _{ijk}\int _0^\tau \!\!\int _D\!\int _D\!\int _{U^2}(\xi _2^i-\xi ^i_1)(\xi _2^j-\xi _1^j)(\xi _2^k-\xi ^k_1)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,\partial _{h^k}\sigma ^{ij}(h)\,dh \,dt\\&=-\varepsilon \sum _{ij}\int _0^\tau \!\!\int _D\!\int _D\!\int _{U^2}\!\!(\xi _1^i-\xi _2^i)(\xi _1^j- \xi _2^j)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,\Delta _h\sigma ^{ij}(h) \,dh \,dt.\end{split}\) 
for any \(\tau >0\) , where \(\sigma =(\sigma _{ij})_{i,j=1}^3\)  is any smooth, compactly supported, isotropic rank 2 tensor – that is, any \(\sigma \in C^2_c(\mathbb {R}^d, \mathbb {R}^{d\times d})\)  of the form
\(\sigma (h)= \left(\omega _1(|h|) \mathbf {I}+ \omega _2(|h|)\hat{h}\otimes \hat{h}\right) \qquad \text{where } \hat{h}={\left\lbrace \begin{array}{ll}h/|h| & h\ne 0 \\ 0 & h=0\end{array}\right.}\) 
for \(\omega _1,\omega _2\in C^2_c(\mathbb {R})\) .

We consider equation (REF ) for \(k=2\)  with the test function \(\varphi (t,x,y)= \eta (t,y-x)\)  (for simplicity replacing \(x_1\)  and \(x_2\)  by \(x\)  and \(y\)  and writing in component form \(\eta = (\eta ^{ij})_{ij}\) ),
\(\begin{split}&\sum _{ij}\int _0^T\!\!\int _D\!\int _D\int _{U^2}\xi _1^i \xi _2^j\frac{\partial \eta ^{ij}}{\partial t}(t,y-x)\,d\nu _{t,x,y}^2(\xi ) \,dx\,dy \,dt\\&\quad +\sum _{ij}\int _D\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\eta ^{ij}(0,y-x)\,d\nu _{0,x,y}^2(\xi ) \,dx \,dy \\&\quad +\sum _{ijk}\int _0^T\!\!\int _D\!\int _D\int _{U^2}\xi _1^i\xi _1^k\xi _2^j\partial _{x^k}\eta ^{ij}(t,y-x)\,d\nu _{t,x,y}^2(\xi ) \,dx\,dy \,dt\\&\quad +\sum _{ijk}\int _0^T\!\!\int _D\!\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\partial _{y^k}\eta ^{ij}(t,y-x)\,d\nu _{t,x,y}^2(\xi ) \,dx\,dy \,dt\\&=-\varepsilon \sum _{ij}\int _0^T\!\!\int _D\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j(\Delta _{x}+\Delta _y)\eta ^{ij}(t,y-x)\,d\nu _{t,x,y}^2(\xi ) \,dx \,dy \,dt\\&\quad +\sum _{ijk}\int _0^T\!\!\int _D\!\int _D\int _{U^2}\xi _1^i\xi _1^k\xi _2^j\partial _{x^k}\partial _{x^i}\psi ^j_{\eta ,1}(t,y-x)\,d\nu _{t,x,y}^2(\xi ) \,dx \,dy \,dt\\&\quad +\sum _{ijk}\int _0^T\!\!\int _D\!\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\partial _{y^k}\partial _{y^j}\psi ^i_{\eta ,2}(t,y-x)\,d\nu _{t,x,y}^2(\xi ) \,dx \,dy \,dt\end{split}\) 
Since \(\psi _{\eta ,1}\)  solves \(\Delta _x\psi _{\eta ,1}=\operatorname{div}_x\eta \)  and \(\psi _{\eta ,2}\)  solves \(\Delta _y\psi _{\eta ,2}=\operatorname{div}_y\eta =-\operatorname{div}_x\eta \) , we have \(\psi _{\eta ,1}= (\Delta _x)^{-1}\operatorname{div}_x\eta \)  and \(\psi _{\eta ,2}=-(\Delta _y)^{-1}\operatorname{div}_x\eta =-(\Delta _x)^{-1}\operatorname{div}_x\eta \)  and so \(\psi _\eta \psi _{\eta ,1}=-\psi _{\eta ,2}\)  (up to additive constants). Using this and changing the integration variables to \(x\)  and \(hy-x\) , we obtain
\(\begin{split}&\sum _{ij}\int _0^T\int _D\!\!\int _D\int _{U^2}\xi _1^i \xi _2^j\frac{\partial \eta ^{ij}}{\partial t}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&\quad +\sum _{ij}\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\eta ^{ij}(0,h)\,d\nu _{0,x,x+h}^2(\xi ) \,dx \,dh \\&\quad -\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _1^k\xi _2^j\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&\quad +\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&=-2\varepsilon \sum _{ij}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\Delta _h\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt\\&\quad +\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _1^k\xi _2^j\partial _{h^k}\partial _{h^i}\psi ^j_{\eta }(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt\\&\quad -\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\partial _{h^k}\partial _{h^j}\psi ^i_{\eta }(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt\end{split}\) 
where \(\psi _\eta = \psi _{\eta ,1}\) .
The cubic terms can be rewritten using the following simple fact (which can also be found in Frisch [28]}, and in a similar, weak form in [26]}):
\(\begin{split}&\quad -2\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _2^j(\xi _1^k-\xi ^k_2)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&=\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _1^i-\xi ^i_2)(\xi _1^j-\xi _2^j)(\xi _1^k-\xi ^k_2)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt.\end{split}\) 
The proof of this is postponed to the end of this proof. Using this, we can rewrite (REF ) as
\(\begin{split}&\quad \sum _{ij}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\frac{\partial \eta ^{ij}}{\partial t}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&\quad +\sum _{ij}\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\eta ^{ij}(0,h)\,d\nu _{0,x,x+h}^2(\xi ) \,dx \,dh \\&\quad -\frac{1}{2}\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _1^i-\xi ^i_2)(\xi _1^j-\xi _2^j)(\xi _1^k-\xi ^k_2)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&=-2\varepsilon \sum _{ij}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\Delta _h\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt\\&\quad +\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _1^k\xi _2^j\partial _{h^k}\partial _{h^i}\psi ^j_{\eta }(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt\\&\quad -\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\partial _{h^k}\partial _{h^j}\psi ^i_{\eta }(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt.\end{split}\) 
Since
\(&\quad \sum _{ij}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _1^j\Delta _h\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt\\&= \sum _{ij}\int _0^T\!\!\int _{D}\!\int _{U}\!\!\xi _1^i \xi _1^j \,d\nu _{t,x}^1(\xi ) \,dx \int _{D} \Delta _h\eta ^{ij}(t,h)\,dh \,dt=0,\) 

we can rewrite
\(&\quad -2\varepsilon \sum _{ij}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\Delta _h\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt\\&=\varepsilon \sum _{ij}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\!\!(\xi _1^i-\xi _2^i)(\xi _1^j- \xi _2^j)\Delta _h\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt.\) 

Moreover, we have for symmetric, smooth and compactly supported rank 2 tensors \(\eta \)  of the form
\(\eta (t,h)= \omega _1(t,|h|) \mathbf {I}+ \omega _2(t,|h|)\hat{h}\otimes \hat{h},\) 
with \(\psi ^i_\eta = -(\Delta _h)^{-1}\operatorname{div}_h\eta ^{i,\cdot }\) ,
\(\begin{split}&\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _1^k\xi _2^j\partial _{h^k}\partial _{h^i}\psi ^j_{\eta }(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt = 0,\\&\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\partial _{h^k}\partial _{h^j}\psi ^i_{\eta }(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt=0,\end{split}\) 
whose proof is postponed to the end of this proof. Using this, (REF ) becomes
\(\begin{split}&\sum _{ij}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\frac{\partial \eta ^{ij}}{\partial t}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&\quad +\sum _{ij}\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\eta ^{ij}(0,h)\,d\nu _{0,x,x+h}^2(\xi ) \,dx \,dh \\&\quad -\frac{1}{2}\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _1^i-\xi ^i_2)(\xi _1^j-\xi _2^j)(\xi _1^k-\xi ^k_2)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&=\varepsilon \sum _{ij}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\!\!(\xi _1^i-\xi _2^i)(\xi _1^j- \xi _2^j)\Delta _h\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt.\end{split}\) 
Let \(\theta _\delta \)  be a sequence of smooth, uniformly bounded functions with the property that \(\theta _\delta \rightarrow \mathbf {1}_{(0,\tau ]}(t)\)  for every \(t\)  as \(\delta \rightarrow 0\) . If we now use a test function
\(\eta (t,h)= \sigma (h)\theta _\delta (t),\) 
where \(\sigma \)  is of the form (REF ), then we can use the weak continuity in time of the moments \(\int _{U^k}\xi _1\otimes \dots \otimes \xi _k \,d\nu ^k_{t,x}(\xi )\)  to obtain for any \(\tau >0\) , as \(\delta \rightarrow 0\) ,
\(\begin{split}&\quad -\sum _{ij}\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\,d\nu _{\tau ,x,x+h}^2(\xi ) \,dx \sigma ^{ij}(h) \,dh +\sum _{ij}\int _D\!\!\int _D\int _{U^2}\!\!\xi _1^i \xi _2^j\,d\nu _{0,x,x+h}^2(\xi ) \,dx \sigma ^{ij}(h) \,dh \\&\quad -\frac{1}{2}\sum _{ijk}\int _0^\tau \!\!\int _D\!\!\int _D\int _{U^2}(\xi _1^i-\xi ^i_2)(\xi _1^j-\xi _2^j)(\xi _1^k-\xi ^k_2)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \partial _{h^k}\sigma ^{ij}(h)\,dh \,dt\\&=\varepsilon \sum _{ij}\int _0^\tau \!\!\int _D\!\!\int _D\int _{U^2}\!\!(\xi _1^i-\xi _2^i)(\xi _1^j- \xi _2^j)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\Delta _h\sigma ^{ij}(h) \,dh \,dt.\end{split}\) 
We expand the right hand side:
\(&\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _1^i-\xi ^i_2)(\xi _1^j-\xi _2^j)(\xi _1^k-\xi ^k_2)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&= \sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _1^i\xi _1^j\xi _1^k-\xi _2^i\xi _2^j\xi _2^k)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,\partial _{h^k}\eta ^{ij}(t,h)\,dh \,dt\\&\quad +\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _2^i\xi _2^j\xi _1^k-\xi _1^i\xi _1^j\xi _2^k)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&\quad -\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _1^i\xi _2^j\xi _1^k-\xi _2^i\xi _1^j\xi _2^k)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\&\quad - \sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _2^i\xi _1^j\xi _1^k-\xi _1^i\xi _2^j\xi _2^k)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt.\) 

The first term on the right hand side is zero since \(\eta \)  is compactly supported (after changing the integration variable from \(x\)  to \(x-h\)  in one of the terms). The second term on the right hand side vanishes using the divergence constraint (REF ). Using that \(\eta \)  and \(\nu \)  are symmetric, the last two terms are identical and so
\(\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _1^i-\xi ^i_2)(\xi _1^j-\xi _2^j)(\xi _1^k-\xi ^k_2)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt\\= -2\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}(\xi _2^i\xi _1^j\xi _1^k-\xi _1^i\xi _2^j\xi _2^k)\partial _{h^k}\eta ^{ij}(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,dh \,dt,\) 
which proves the claim.
We consider the second expression,
assume \(\eta \)  is of the form
\(\eta (t,h)= \omega _1(t,|h|) \mathbf {I}+ \omega _2(t,|h|)\hat{h}\otimes \hat{h},\) 
where \(\omega _i\) , \(i=1,2\)  are compactly supported in the torus. Then using that \(\psi ^i_\eta = -(\Delta _h)^{-1}\operatorname{div}_h\eta ^{i,\cdot }\)  (the first term is treated in a similar way)
\(E &-\sum _{ijk}\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\partial _{h^k}\partial _{h^j}\psi ^i_{\eta }(t,h)\,d\nu _{t,x,x+h}^2(\xi ) \,dx \,dh \,dt\\&= \sum _{ijk\ell }\int _0^T\!\!\int _D\!\!\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\,d\nu _{t,x,x+h}^2(\xi ) \,dx\,\partial _{h^k}\partial _{h^j}\Delta ^{-1}\left(\partial _{h^\ell }\eta ^{i\ell }(t,h)\right) \,dh \,dt\\&= \sum _{ijk\ell }\int _0^T\!\!\int _{D}\!\partial _{h^k}\partial _{h^j}\Delta ^{-1}\left(\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\,d\nu _{t,x,x+h}^2(\xi ) \,dx\right)\partial _{h^\ell }\eta ^{i\ell }(t,h)dh \,dt.\) 

We note that since \(\omega _i\)  have compact support, we can write in polar coordinates
\(\sum _\ell \partial _{h^\ell } \eta ^{i\ell }(t,h) = \biggl (\underbrace{\omega ^{\prime }_1(t,|h|)+\omega _2^{\prime }(t,|h|)+2\frac{\omega _2(t,|h|)}{|h|}}_{\,G(t,|h|)}\biggr )\hat{h}^i\) 
and so
\(\begin{split}E&=\sum _{ijk\ell }\int _0^T\!\!\int _{D}\!\partial _{h^k}\partial _{h^j}\Delta ^{-1}\left(\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\,d\nu _{t,x,x+h}^2(\xi ) \partial _{h^\ell }\eta ^{i\ell }(t,h)\,dx\right)dh \,dt\\&= \sum _{ijk}\int _0^T\!\!\int _{D}\!\partial _{h^k}\partial _{h^j}\Delta ^{-1}\left(\int _D\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\,d\nu _{t,x,x+h}^2(\xi ) \,dx\right) G(t,|h|)\hat{h}^i \,dh \,dt\\&= \sum _{jk}\int _0^T\!\!\int _0^\infty \int _{|h|=r}\!\partial _{h^k}\partial _{h^j}\Delta ^{-1}\left(\int _D\!\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\,d\nu _{t,x,x+h}^2(\xi ) \,dx\right) \hat{h}^i dS(h)\, G(t,r)\, dr \,dt\\&= \sum _{ijk}\int _0^T\!\!\int _0^\infty \int _{|h|\leqslant r}\!\partial _{h^i}\partial _{h^k}\partial _{h^j}\Delta ^{-1}\left(\int _D\!\int _{U^2}\xi _1^i\xi _2^j\xi _2^k\,d\nu _{t,x,x+h}^2(\xi ) \,dx\right)dh\, G(t,r)\, dr \,dt\\&= \sum _{jk}\int _0^T\!\!\int _0^\infty \int _{|h|\leqslant r}\!\partial _{h^k}\partial _{h^j}\Delta ^{-1}\operatorname{div}_h\left(\int _D\!\int _{U^2}\xi _1\xi _2^j\xi _2^k\,d\nu _{t,x,x+h}^2(\xi ) \,dx\right)dh\, G(t,r)\, dr \,dt=0,\end{split}\) 
where we used the divergence theorem in the second to last identity and the divergence constraint (REF ) for the last identity.

Scaling of third order structure functions
Next, we use the KHM-relation (REF ) to derive a scaling relation for the averaged third order structure function in terms of the measure \(\mu _t\) ,
\(S^{3}_{0}(\tau ,r)=\int _0^\tau \int _{L^2_{\operatorname{div}}}\!_{\mathbb {S}^2}\int _{D}\big |u(x)-u(x+r n)\big |^2\big (u(x+r n)-u(x)\big )\cdot n \,dx\, dS(n) \,d\mu _t(u)\,dt,\) 
which will be more convenient to work with for this purpose. We have:
Lemma 4.2 
Let \(\mu _t\)  be a Foiaş–Prodi statistical solution of the Navier–Stokes equations (cf. Definition REF ). Then
\(\left|\frac{S^{3}_{0}(\tau ,r)}{r}\right|\leqslant 2 E_0,\) 
where \(E_0\)  is the initial energy,
\(E_0\int _{L^2_{\operatorname{div}}(D;U)}\left\Vert  u(x) \right\Vert ^2_{L^2(D)} \,d\mu _0(u).\) 

We take a test function of the form \(\sigma (h)=\omega (|h|)\mathbf {I}\)  in the KHM-relation (REF ) with \(\omega \)  having compact support in \([0,0.5)\) . A little bit of algebra yields (denoting \(\hat{h}_{k}=h^k/|h|\) )
\(\partial _{h^k}\omega (|h|)=\omega ^{\prime }(|h|)\hat{h}^k,\) 
and so (REF ) for this particular test function reads
\(\begin{split}&\int _D\int _D\!\int _{U^2}\!\!\xi _1\cdot \xi _2 \,d\nu _{\tau ,x,x+h}^2(\xi ) \,dx \omega (|h|) \,dh-\int _D\int _D\!\int _{U^2}\!\!\xi _1\cdot \xi _2\,d\nu _{0,x,x+h}^2(\xi ) \,dx \omega (|h|) \,dh \\&+\frac{1}{2}\int _0^\tau \!\!\int _D\int _D\!\int _{U^2}|\xi _1-\xi _2|^2(\xi _2-\xi _1)\cdot \hat{h}\,d\nu _{t,x,x+h}^2(\xi ) \,dx \omega ^{\prime }(|h|)\,dh \,dt\\=&-\varepsilon \int _0^\tau \!\!\int _D\int _D\!\int _{U^2}\!\!|\xi _1-\xi _2|^2 \,d\nu _{t,x,x+h}^2(\xi ) \,dx\Delta _h\omega (|h|) \,dh \,dt.\end{split}\) 
In terms of the statistical solution \((\mu _t)_{t>0}\)  this is
\( \begin{split}&\quad \int _{L^2_{\operatorname{div}}}\!\int _D\int _D\!\! u(x)\cdot u(x+h) \,dx \omega (|h|) \,dh \,d\mu _\tau (u)-\int _{L^2_{\operatorname{div}}}\int _D\int _D\! u(x)u(x+h) \,dx\, \omega (|h|) \,dh \,d\mu _0(u) \\&\quad +\frac{1}{2}\int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}\int _D\int _D\!|u(x)-u(x+h)|^2(u(x+h)-u(x))\cdot \hat{h}\, \,dx \omega ^{\prime }(|h|)\,dh \,d\mu _t(u) \,dt\\&=-\varepsilon \int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}\int _D\int _D\!|u(x)-u(x+h)|^2 \,dx\Delta _h\omega (|h|) \,dh \,d\mu _t(u) \,dt.\end{split}\) 
The last term can also be written as
\(\begin{split}&\varepsilon \int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}\int _D\int _D\!|u(x)-u(x+h)|^2 \,dx\Delta _h\omega (|h|) \,dh \,d\mu _t(u) \,dt\\& = 2\varepsilon \int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}\int _D\int _D\!\nabla _h u(x+h)\cdot (u(x)-u(x+h)) \,dx\nabla _h\omega (|h|) \,dh \,d\mu _t(u) \,dt\\& = 2\varepsilon \int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}\int _D\int _D\!\nabla _x u(x+h)\cdot (u(x)-u(x+h)) \,dx\nabla _h\omega (|h|) \,dh \,d\mu _t(u) \,dt\\& = 2\varepsilon \int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}\int _D\int _D\!\nabla _x u(x)\cdot (u(x-h)-u(x)) \,dx\nabla _h\omega (|h|) \,dh \,d\mu _t(u) \,dt\\& = -2\varepsilon \int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}\int _D\int _D\!\nabla _x u(x): \nabla _h u(x-h) \,dx \omega (|h|) \,dh \,d\mu _t(u) \,dt\\& = 2\varepsilon \int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}\int _D\int _D\!\nabla _x u(x+h): \nabla _x u(x) \,dx \omega (|h|) \,dh \,d\mu _t(u) \,dt\end{split}\) 
Changing to spherical coordinates and using the definition of \(S^{3}_0\) , (REF ), we obtain
\(&\frac{1}{2}\int _0^\infty S_0^3(\tau ,r) r^2\omega ^{\prime }(r)dr\\&=-2\varepsilon \int _0^\infty \int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}_{\mathbb {S}^2}\int _{D}\!\nabla _x u(x):\nabla _x u(x+r n) \,dxdS(n)\,d\mu _t(u) \,dt r^2\omega (r) dr \\&\quad -\int _0^\infty \int _{L^2_{\operatorname{div}}}\!_{\mathbb {S}^2}\int _{D}\! u(x)\cdot u(x+r n) \,dxdS(n) \,d\mu _\tau (u)r^2\omega (r) dr \\&\quad +\int _0^\infty \int _{L^2_{\operatorname{div}}}_{\mathbb {S}^2}\int _{D}\! u(x)\cdot u(x+r n) \,dx dS(n) \,d\mu _0(u) r^2\omega (r) dr.\) 

We denote
\(m_2(\tau ,r)&\int _{L^2_{\operatorname{div}}}\!_{\mathbb {S}^2}\int _{D}\! u(x)\cdot u(x+r n) \,dxdS(n) \,d\mu _\tau (u)\\v(\tau ,r)&\int _0^\tau \!\!\int _{L^2_{\operatorname{div}}}_{\mathbb {S}^2}\int _{D}\!\nabla _x u(x):\nabla _x u(x+r n) \,dxdS(n)\,d\mu _t(u) \,dt.\) 

Since \((\mu _t)_{t>0}\)  is supported on functions in \(L^\infty ([0,\infty );L^2_{\operatorname{div}}(D;U))\cap L^2([0,\infty );\) , \(S_0^3\)  is a continuous function. Moreover, notice that due to the a priori bounds following from the energy inequality (REF ), both \(m_2\)  and \(v\)  are uniformly bounded and continuous in \(r\)  and \(\tau \)  (for the continuity in \(\tau \)  of the first quantity, one needs weak time continuity of the moments which follows the fact that they satisfy the equations (REF ) where all the terms are integrable). We obtain
\(\frac{1}{2}\int _0^\infty S_0^3(\tau ,r) r^2\omega ^{\prime }(r)dr &=-2\varepsilon \int _0^\infty v(\tau ,r)r^2\omega (r) dr \\&\quad -\int _0^\infty m_2(\tau ,r)r^2\omega (r) dr +\int _0^\infty m_2(0,r) r^2\omega (r) dr,\) 

which is an ODE in the sense of distributions for \(S_0^3(\tau ,\cdot )\) , and because the right hand side is uniformly bounded and continuous, we can consider it in the strong sense (note that boundary terms when integrating the \(S_0^3\)  term by parts vanish):
\(\frac{1}{r^2}\partial _r(r^2 S_0^3(r)) = 4\varepsilon v(\tau ,r)+2 m_2(\tau ,r)-2 m_2(0,r),\) 
or
\(\frac{S_0^3(r)}{r} = \frac{2}{r^3}\int _0^r s^2\big (2\varepsilon v(\tau ,s)+m_2(\tau ,s)-m_2(0,s)\big ) \,ds.\) 
The energy inequality (REF ) and the Cauchy–Schwarz inequality imply that \(2\varepsilon v(\tau ,s)\)  and \(m_2(\tau ,s)\)  are both bounded by \(E_0\)  (defined in (REF )), uniformly in \(\tau ,s\) . Hence,
\(\left|\frac{S_0^3(r)}{r}\right|\leqslant \frac{2}{r}\int _0^r \left(\frac{s}{r}\right)^2\big (2\varepsilon |v(\tau ,s)|+|m_2(\tau ,s)|+|m_2(0,s)|\big ) \,ds\leqslant 2E_0.\) 
(See also [26]} for a related result.)
Using this lemma, we can derive a scaling relation for the averaged longitudinal structure function \(S^3_\Vert \) , where
\(S_\Vert ^p(\tau ,r)=\int _0^\tau \int _{L^2_{\operatorname{div}}}\!\!_{\mathbb {S}^2}\int _{D}\big ((u(x+r n)-u(x))\cdot n\big )^p \,dxdS(n) \,d\mu _t(u)\,dt.\) 
Lemma 4.3 
Let \(\mu _t\)  be a Foiaş–Prodi statistical solution of the Navier–Stokes equations (cf. Definition REF ). Then
\(\left|\frac{S^{3}_\Vert (\tau ,r)}{r}\right|\leqslant 2 E_0,\) 
where \(C>0\)  is some constant independent of \(\varepsilon \)  and \(E_0\)  is the initial energy,
\(E_0\int _{L^2_{\operatorname{div}}(D;U)}\left\Vert  u(x) \right\Vert ^2_{L^2(D)} \,d\mu _0(u)\) 

Again, we start with the KHM relation (REF ). This time we use the test function \(\sigma (h) = \omega (|h|)\hat{h}\otimes \hat{h}\)  where \(\omega \in C_c^\infty (\mathbb {R})\)  is an even function. We have
\(\begin{split}\partial _{h^k} (\omega (|h|)\hat{h}^i \hat{h}^j)& = \omega ^{\prime }(|h|)\hat{h}^i\hat{h}^j\hat{h}^k + \frac{1}{|h|}\left(\delta _{ik}\hat{h}^j+\delta _{jk}\hat{h}^i-2\hat{h}^i\hat{h}^j\hat{h}^k\right)\omega (|h|)\\& =\left(\omega ^{\prime }(|h|)-2\frac{\omega (|h|)}{|h|}\right)\hat{h}^i\hat{h}^j\hat{h}^k + \frac{1}{|h|}\left(\delta _{ik}\hat{h}^j+\delta _{jk}\hat{h}^i\right)\omega (|h|)\end{split}\) 
Therefore, (REF ) becomes
\(\begin{split}&-\int _D\int _D\!\int _{U^2}\!\!\xi _1\cdot \hat{h}\, \xi _2\cdot \hat{h} \,d\nu _{\tau ,x,x+h}^2(\xi ) \,dx \omega (|h|) \,dh\\& +\int _D\int _D\!\int _{U^2}\!\!\xi _1\cdot \hat{h} \,\xi _2 \cdot \hat{h} \, \,d\nu _{0,x,x+h}^2(\xi ) \,dx \omega (|h|) \,dh \\&-\frac{1}{2}\int _0^\tau \!\!\int _D\int _D\!\int _{U^2}\big ((\xi _2-\xi _1)\cdot \hat{h}\big )^3\,d\nu _{t,x,x+h}^2(\xi ) \,dx\left(\omega ^{\prime }(|h|)-2|h|^{-1}\omega (|h|)\right) \,dh \,dt\\&+\int _0^\tau \!\!\int _D\int _D\!\int _{U^2}|\xi _1-\xi _2|^2(\xi _1-\xi _2)\cdot \hat{h} \,d\nu _{t,x,x+h}^2(\xi ) \,dx |h|^{-1}\omega (|h|) \,dh \,dt\\=&\ \varepsilon \sum _{ij}\int _0^\tau \!\!\int _D\int _D\!\int _{U^2}\!\!(\xi _1^i-\xi _2^i)(\xi _1^j- \xi _2^j)\,d\nu _{t,x,x+h}^2(\xi ) \,dx\Delta _h\sigma ^{ij}(h) \,dh \,dt.\end{split}\) 
Again, in terms of \((\mu _t)_{t>0}\) , this means
\(\begin{split}&\quad -\int _{L^2_{\operatorname{div}}}\int _D\int _D\! u(x)\cdot \hat{h}\, u(x+h)\cdot \hat{h} \, \,dx \omega (|h|) \,dh \,d\mu _\tau (u)\\&\quad +\int _{L^2_{\operatorname{div}}}\int _D\int _D\!u(x)\cdot \hat{h} \,u(x+h) \cdot \hat{h} \,dx \omega (|h|) \,dh \,d\mu _0(u) \\&\quad -\frac{1}{2}\int _0^\tau \int _D\int _{L^2_{\operatorname{div}}}\!\!\int _D\big ((u(x+h)-u(x))\cdot \hat{h}\big )^3 \,dx \,d\mu _t(u)\left(\omega ^{\prime }(|h|)-2|h|^{-1}\omega (|h|)\right) \,dh \,dt\\&\quad +\int _0^\tau \int _D\int _{L^2_{\operatorname{div}}}\!\!\int _D|u(x)-u(x+h)|^2\big (u(x)-u(x+h)\big )\cdot \hat{h} \,dx \,d\mu _t(u) |h|^{-1}\omega (|h|) \,dh \,dt\\&=\varepsilon \sum _{ij}\int _0^\tau \int _D\int _{L^2_{\operatorname{div}}}\!\!\int _D\!\big (u^i(x)-u^i(x+h)\big )\big (u^j(x)- u^j(x+h)\big ) \,dx \,d\mu _t(u)\Delta _h\sigma ^{ij}(h) \,dh \,dt.\end{split}\) 
Similar to the computation in (REF ), we have
\(&\quad \varepsilon \sum _{ij}\int _0^\tau \int _D\int _{L^2_{\operatorname{div}}}\!\!\int _D\!\big (u^i(x)-u^i(x+h)\big )\big (u^j(x)- u^j(x+h)\big ) \,dx \,d\mu _t(u)\Delta _h\sigma ^{ij}(h) \,dh \,dt\\&=2\varepsilon \sum _{ij}\int _0^\tau \int _D\int _{L^2_{\operatorname{div}}}\!\!\int _D\!\nabla _xu^i(x+h)\nabla _x u^j(x) \,dx \,d\mu _t(u)\sigma ^{ij}(h) \,dh \,dt\\&=2\varepsilon \int _0^\tau \int _D\int _{L^2_{\operatorname{div}}}\!\!\int _D\!\big (\nabla _x u(x+h)\cdot \hat{h}\big )\cdot \big (\nabla _x u(x)\cdot \hat{h}\big ) \,dx \,d\mu _t(u)\omega (|h|) \,dh \,dt,\) 

so equation (REF ) becomes (after switching to polar coordinates)
\(\begin{split}&\quad -\int _0^\infty \int _{L^2_{\operatorname{div}}}_{\mathbb {S}^2}\int _{D}\! u(x)\cdot n\, u(x+r n)\cdot n \, \,dx\,dS(n) \,d\mu _\tau (u)r^2\omega (r) \,dr\\&\quad +\int _0^\infty \int _{L^2_{\operatorname{div}}}_{\mathbb {S}^2}\int _{D}\! u(x)\cdot n \,u(x+r n) \cdot n \,dx \,dS(n) \,d\mu _0(u)r^2\omega (r) \,dr\\&\quad -\frac{1}{2}\int _0^\infty \int _0^\tau \int _{L^2_{\operatorname{div}}}\!\!_{\mathbb {S}^2}\int _{D}\big (\big (u(x+r n)-u(x)\big )\cdot n\big )^3 \,dx\,dS(n) \,d\mu _t(u)\,dt\left(r^2\omega ^{\prime }(r)-2r\omega (r)\right) \,dr \\&\quad +\int _0^\infty \int _0^\tau \int _{L^2_{\operatorname{div}}}\!\!_{\mathbb {S}^2}\int _{D}|u(x)-u(x+r n)|^2\big (u(x)-u(x+r n)\big )\cdot n \,dx \,dS(n) \,d\mu _t(u)\,dt r\omega (r) \,dr\\&=2\varepsilon \int _0^\infty \int _0^\tau \int _{L^2_{\operatorname{div}}}\!\!_{\mathbb {S}^2}\int _{D}\!\big (\nabla _x u(x+r n)\cdot n\big )\cdot \big (\nabla _x u(x)\cdot n\big ) \,dx \,dS(n) \,d\mu _t(u)\,dt r^2\omega (r) \,dr\end{split}\) 
Denote
\(\widetilde{m}_2(\tau ,r)&\int _{L^2_{\operatorname{div}}}_{\mathbb {S}^2}\int _{D}\! u(x)\cdot n\, u(x+r n)\cdot n \, \,dx \,dS(n) \,d\mu _\tau (u)\\\widetilde{v}(\tau ,r)&\int _0^\tau \int _{L^2_{\operatorname{div}}}\!\!_{\mathbb {S}^2}\int _{D}\!\big (\nabla _x u(x+r n)\cdot n\big )\cdot \big (\nabla _x u(x)\cdot n\big ) \,dx \,dS(n) \,d\mu _t(u)\,dt.\) 

Writing \(\omega ^{\prime }(r)-2r^{-1}\omega (r) = r^2 \frac{d}{dr}\big (\omega (r)r^{-2}\big )\) , (REF ) becomes
\(&\quad \frac{1}{2}\int _0^\infty r^4 S_\Vert ^3(\tau ,r)\partial _r\big (r^{-2}\omega (r)\big ) \,dr +\int _0^\infty S_0^3(\tau ,r) r\omega (r) \,dr\\&=-\int _0^\infty \widetilde{m}_2(\tau ,r)r^2\omega (r) \,dr+\int _0^\infty \widetilde{m}_2(0,r)r^2\omega (r) \,dr-2\varepsilon \int _0^\infty \widetilde{v}(\tau ,r)r^2\omega (r) \,dr.\) 

Again, we note that due to the estimates from the energy inequality (REF ), \(S^3_\Vert \) , \(\widetilde{m}\)  and \(\widetilde{v}\)  are continuous and bounded quantities in \(\tau \)  and \(r\) . And so we can consider this ODE in the sense of distributions as an ODE in the strong sense,
\(\partial _r \big (r^4S_\Vert ^3(\tau ,r)\big ) = 2r^4\left(\frac{S_0^3(\tau ,r)}{r} +\widetilde{m}_2(\tau ,r)-\widetilde{m}_2(0,r)+2\varepsilon \widetilde{v}(\tau ,r)\right),\) 
or
\(S_\Vert ^3(\tau ,r) = 2\int _0^r \frac{s^4}{r^4}\left(\frac{S_0^3(\tau ,s)}{s} +\widetilde{m}_2(\tau ,s)-\widetilde{m}_2(0,s)+2\varepsilon \widetilde{v}(\tau ,s)\right)\,ds.\) 
By the energy bound and Cauchy–Schwarz inequality, \(\widetilde{m}_2\)  and \(\varepsilon \widetilde{v}\)  are uniformly bounded in \(\varepsilon \)  for all \(s,\tau >0\) . Moreover from Lemma REF , we have that \(s^{-1}S_0^3(s,\tau )\)  is uniformly bounded in \(\varepsilon \)  by \(2 E_0\) . Hence,
\(\left|\frac{S_\Vert ^3(\tau ,r)}{r}\right|\leqslant 2 E_0\) 
for some \(C>0\)  independent of \(\varepsilon \) .
Remark 4.4 Combining (REF ) and (REF ), we also obtain a uniform bound on \(r^{-1} S^3_\perp (\tau ,r)\) , where \(S^3_\perp \)  is the transversal structure function
\(S^3_\perp (\tau ,r)\int _0^\tau \int _{L^2_{\operatorname{div}}}\!\!_{\mathbb {S}^2}\int _{D}|\delta ^\perp _{r n}u(x)|^2\delta _{r n} u(x)\cdot n \,dx \,dS(n) \,d\mu _t(u)\,dt = S^3_0(\tau ,r)-S^3_{\Vert }(\tau ,r),\) 
where
\(\delta _h u(x) =u(x+h)-u(x),\quad \delta ^\perp _{h}u = (\mathbf {I}-\hat{h}\otimes \hat{h})\delta _h u, \quad \hat{h}\frac{h}{|h|}.\) 
None of these quantities has a sign and therefore the previously derived bounds do not imply compactness without further assumptions.


Scaling assumption
In order to pass to the limit \(\varepsilon \rightarrow 0\) , we need to an additional assumption about the behavior of structure functions. Specifically, we need
Assumption 1 (Weak statistical scaling) 
For any \(\varepsilon >0\) , let \(\mu ^\varepsilon \)  be a Foiaş–Prodi statistical solution of the incompressible Navier–Stokes equations. We assume that for \(r\ll 1\) , the second and third order longitudinal structure functions (REF ) are related by
\(\left|S^2_{\Vert }(\tau ,r)\right|\leqslant C \left|S^3_{\Vert }(\tau ,r)\right|^{\alpha },\) 
where \(C\)  is a constant independent of \(\varepsilon \)  and \(\alpha >0\) .

Remark 4.5 (Weak statistical scaling) 
Assumption REF  is inspired by the following stronger scaling assumption often encountered in turbulence theory: For any \(p,q\)  with \(q\geqslant p\geqslant 1\) , the \(p\) -th and \(q\) -th order longitudinal structure functions (REF ) are related by
\(\left|S^p_{\Vert }(\tau ,r)\right|\leqslant C \left|S^q_{\Vert }(\tau ,r)\right|^{\frac{\lambda (p)}{\lambda (q)}},\) 
where \(C\)  is a constant independent of \(\varepsilon \)  and \(\lambda (p)>0\)  for \(p\leqslant p_0\)  where \(3\leqslant p_0\in \mathbb {R}\cup \lbrace \infty \rbrace \) . In Kolmogorov's 1941 (“K41”) theory [38]}, [39]}, [40]}, \(\lambda (p)=p\) . However, this cannot be confirmed with physical experiments [41]}, [42]}. Various physicists therefore suggested intermittency corrections to account for the deviation from Kolmogorov's original theory, among others, Kolmogorov himself in 1962 [43]} in his refined theory of turbulence, Frisch et al. the \(\beta \) -model [44]}, as well as Novikov and Stewart [45]}. Assumption (REF ) can also accommodate the frequently used model by She and Leveque [46]} who suggested
\(\lambda (p)= \frac{p}{9}+2\left(1-\left(\frac{2}{3}\right)^{p/3}\right).\) 

Remark 4.6 Combining the bound on the third order structure function in Lemma REF  with Assumption REF , we obtain
\(\left|S^2_{\Vert }(\tau ,r)\right|\leqslant C r^{\alpha }.\) 

We will combine Assumption REF  with the following lemma, whose proof is given in Appendix :
Lemma 4.7 (Weak anisotropy) 
Let \(\mu _t\)  be a statistical solution of the Navier–Stokes equation. Then \(\mu \)  satisfies
\(\begin{split}&3\int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)}_{\partial B_r(0)} (\delta _{rn} u\cdot n)^2dS(n)\,dx \,d\mu _t(u)\,dt\\&=\int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)}_{B_r(0)} |\delta _\ell u(x)|^2 d\ell \,dx \,d\mu _t(u)\,dt .\end{split}\) 

Under Assumption REF , we obtain
\(\int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)}_{B_r(0)} |\delta _\ell u(x)|^2 d\ell \,dx \,d\mu _t(u)\,dt\leqslant C r^\alpha .\) 
Using the equivalence theorem REF , we can write this as
\(\int _0^T\int _{D}_{B_r(0)} \int _{U^2}|\xi _1-\xi _2|^2 \,d\nu _{x,x+y}^2(\xi ) \,dy \,dx \,dt\leqslant C r^\alpha .\) 
and since this is uniform with respect to the viscosity coefficient \(\varepsilon \)  (by the weak scaling assumption), it implies uniform diagonal continuity of the sequence \(\lbrace \mathbf {\nu }\rbrace _{\varepsilon >0}\) .

Passage to the limit \(\varepsilon \rightarrow 0\)
Now we are in a position to prove our main result. We will keep track of the superscript \(\varepsilon \)  again in order to distinguish between the approximating sequence \(\lbrace \mathbf {\nu }^\varepsilon \rbrace _{\varepsilon >0}\)  and the limiting measure \(\mathbf {\nu }\)  for \(\varepsilon =0\) .
Theorem 4.8 
Let \(\lbrace \mathbf {\nu }^\varepsilon \rbrace _{\varepsilon >0}\)  be a sequence of (either Foiaş–Temam or Friedman–Keller) statistical solutions to the Navier–Stokes equations with initial data \(\mu _0\)  with bounded support (cf. (REF )). Assume that \(\mathbf {\nu }^\varepsilon \)  all satisfy Assumption REF . Then, as \(\varepsilon \rightarrow 0\) , \(\mathbf {\nu }^\varepsilon \)  converges (along a subsequence) to a correlation measure \(\mathbf {\nu }\)  on \(L^2\)  with bounded support (cf. (REF ))
\(\int _{D^k}\int _{U^{k}} |\xi _1|^2\dots |\xi _k|^2 \,d\nu _{t,x}^k(\xi ) \,dx\leqslant R^k<\infty ,\) 
for some \(0<R<\infty \) , any \(k\in \mathbb {N}\)  and that satisfies the “inviscid Friedman–Keller system”:
\(\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\frac{\partial \varphi }{\partial t}(t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt+\int _{D^k}\!\!\int _{U^k}\!\!(\xi _1\otimes \cdots \otimes \xi _k):\varphi (0,x)\,d\nu _{0,x}^k(\xi ) \,dx \\+\sum _{i=1}^k\int _0^T\!\!\int _{D^k}\!\!\int _{U^k}(\xi _1\otimes \cdots \otimes (\xi _i\otimes \xi _i)\otimes \cdots \xi _k):\nabla _{x_i}\varphi (t,x)\,d\nu _{t,x}^k(\xi ) \,dx \,dt=0\) 
for all \(k\in \mathbb {N}\) , for all \(\varphi \in C^2_c([0,T]\times D^k;U^k)\)  with \(\operatorname{div}_{x_i}\varphi (x)=0\) , a.e. \(x\in D^k\)  for all \(i=1,\dots , k\) 
and (corresponding to the divergence constraint)
\(\int _{D^{k}}\int _{U^{k}} \xi _1\otimes \dots \otimes \xi _{\ell }\otimes \alpha _{\ell +1}(\xi _{\ell +1})\otimes \dots \otimes \alpha _k(\xi _k) \,d\nu _{t,x}^k(\xi )\cdot \nabla _{x_1,\dots , x_\ell }\psi (x) \,dx = 0,\) 
where \(\nabla _{x_1,\dots ,x_\ell }=(\nabla _{x_{1}},\dots ,\nabla _{x_\ell })^\top \) , \(1\leqslant \ell \leqslant k\in \mathbb {N}\) , for all \(\psi \in H^1(D^k;U^{k-\ell })\) , \(\alpha _j\in C(U;U)\) , \(\alpha _j(v)\leqslant C(1+|v|^2)\)  and \(j=1,\dots , k\) .

From the condition on the initially bounded support (REF ) and the energy inequality (REF ), we obtain that the sequence \(\mathbf {\nu }^\varepsilon \)  satisfies (REF ) for \(p=2\)  uniformly in \(\varepsilon >0\) . The reasoning of Subsection REF , REF  and REF  resulting in (REF ) imply that \(\mathbf {\nu }^\varepsilon \)  is uniformly diagonal continuous as in (). Hence, using Theorem REF , we obtain, up to subsequence, the existence of a limiting correlation measure \(\mathbf {\nu }\in {L}^{2}([0,T),D;U)\) . So it remains to check whether \(\mathbf {\nu }\)  satisfies the equations (REF ) and (REF ). We note that the functions
\(\begin{split}g_1(t,x,\xi )&(\xi _1\otimes \cdots \otimes \xi _k):\frac{\partial \varphi }{\partial t}(t,x),\\g_2(t,x,\xi )&\sum _{i=1}^k (\xi _1\otimes \cdots \otimes \xi _k):\Delta _{x_i}\varphi (t,x),\\g_3(t,x,\xi )&\sum _{i=1}^k (\xi _1\otimes \cdots \otimes (\xi _i\otimes \xi _i)\otimes \cdots \xi _k):\nabla _{x_i}\varphi (t,x)\end{split}\) 
for \(\varphi \in C^2_c([0,T]\times D^k;U^k)\) , \(1\leqslant \ell \leqslant k\in \mathbb {N}\) , are all functions in \({H}^{k,p}_1([0,T],D;U)\) .
Hence, we can pass to the limit in all the terms in the Friedman–Keller system (REF ). The term that is multiplied by \(\varepsilon \)  vanishes because it is a uniformly bounded in \(\varepsilon >0\)  quantity that is multiplied by \(\varepsilon \) . For the divergence constraint (REF ), we note that the function
\(g_4(t,x,\xi )\theta (t)\xi _1\otimes \dots \otimes \xi _{\ell }\otimes \alpha _{\ell +1}(\xi _{\ell +1})\otimes \dots \otimes \alpha _k(\xi _k)\cdot \nabla _{x_1,\dots , x_\ell }\psi (x),\) 
lies in \({H}^{k,p}_1([0,T],D;U)\)  for any \(\theta \in C^\infty _c((0,T))\) , and \(\psi \in H^1(D^k;U^{k-\ell })\) , and \(\alpha _j\in C(U;U)\)  with \(|\alpha _j(v)|\leqslant C(1+|v|^2)\) . Passing \(\varepsilon \rightarrow 0\)  in \(\langle \nu ^{\varepsilon ,k},\, g_4\rangle \)  and using that \(\nu ^{\varepsilon ,k}\)  satisfy the divergence constraint (REF ), we can conclude, by the arbitrariness of \(\theta \) , that (REF ) holds for a.e. \(t\in [0,T]\) .
Remark 4.9 By the equivalence theorem REF , we know that the limiting correlation measure \(\mathbf {\nu }\)  corresponds to a parametrized measure \(\mu =(\mu _t)_{t>0}:[0,T)\rightarrow \mathcal {P}(L^2_{\operatorname{div}}(D;U))\)  that satisfies
\(\int _{L^2(D)} \Phi (u) \,d\mu _t(u)=\int _{L^2(D)}\Phi (u) \,d\mu _0(u)+\int _0^t\int _{L^2(D)}\int _{D} (u(x)\otimes u(x)):\nabla _x \Phi ^{\prime }(u)(x) \,dx \,d\mu _s(u) \,ds\) 
for all cylindrical test functions \(\Phi \in \mathcal {T}_{\text{cyl}}^0\)  that satisfy \(g_j\in C^2(D)\) ,
and the energy inequality
\(\int _{L^2(D)} \Vert u\Vert _{L^2(D)}^2 \,d\mu _t(u)\leqslant \int _{L^2(D)} \Vert u\Vert _{L^2(D)}^2 \,d\mu _0(u),\quad \text{for all } t\in [0,T].\) 
The proof of this fact follows along the lines of the proof of Theorem REF  while ignoring the terms involving \(\varepsilon \)  and not attempting to recover \(B(u)\)  as it may be unbounded.


Discussion
It is well-known that many incompressible fluid flows of interest are characterized by very-high Reynolds number. Hence, a precise characterization of the vanishing viscosity (\(\varepsilon \rightarrow 0\) ) limit of the Navier-Stokes equations (REF ) is of great interest. Formally, one would expect that the vanishing viscosity limit of Navier-Stokes equations is related to the incompressible Euler equations. However as mentioned in the introduction, rigorous results in this direction are only available in two space dimensions, even in the case of periodic boundary conditions. The key aim of this article was to investigate the vanishing viscosity limit of the Navier-Stokes equations, including in three space dimensions.
It is well known that fluid flows at high Reynolds numbers are characterized by turbulence, loosely speaking, marked by the presence of energy containing eddies at smaller and smaller scales. This phenomenon is clearly linked to the lack of compactness in the Leray-Hopf Navier-Stokes solutions as well as their possible instabilities/non-uniqueness.
Hence, one needs to make further assumptions on the Leray-Hopf solutions that can yield additional information and facilitate passage to the limit. One avenue for making such assumptions, which are realistic and possibly observed in experiments, comes from physical theories of turbulence. In particular, Kolmogorov's well-known K41 theory is based on several verifiable assumptions on the underlying fluid flow and results in a precise characterization of quantities such as structure functions and energy spectra.
In [47]}, Chen and Glimm relate the K41 energy spectra to compactness results on the Leray-Hopf solutions, in appropriate Sobolev and Hölder spaces. Consequently, under the assumption of the K41 energy spectrum, the authors prove that the underlying Leray-Hopf solutions converge to weak solutions of the incompressible Euler equations as \(\varepsilon \rightarrow 0\) . However, Kolmogorov's derivation of the decay of energy spectra is based on a probabilistic charectization of the underlying fluid flow. In particular, assumptions such as (statistical) homogeneity, isotropy and scaling, which form the foundation of Kolmogorov's theory, are too stringent if imposed at the deterministic level, as done in [47]}. Moreover, it is now well-established that the strong scaling assumptions of Kolmogorov might not hold in real fluid flows and intermittent corrections are necessary. Hence, the applicability of the assumptions and results of [47]} can be questioned from this perspective.
Nevertheless, the connection with Kolmogorov's theories of turbulence and their variants forms the basis of our work. We start with the realization that a probabilistic description of the solutions of Navier-Stokes equations is necessary to relate physical theories of turbulence to rigorous mathematical statements. To this end, we focus on statistical solutions of Navier-Stokes equations. Two possible frameworks of such statistical solutions are available, namely the Foiaş-Prodi statistical solutions (see Definition REF ) and the Friedman-Keller statistical solutions (see Definition REF ), which is based on the concept of correlation measures of [2]}. We prove that both these solution concepts are equivalent as long as a statistical version of the energy inequality holds. This also allows us to prove the existence of Friedmann-Keller statistical solutions of the incompressible Navier-Stokes equations.
Then, we investigated the vanishing viscosity limit of the statistical solutions of the incompressible Navier-Stokes equations. To this end, we derived a suitable statistical version of the well-known Kármán-Howarth-Monin relation and used it to prove precise rates for the asymptotic decay of a averaged third-order structure function in Lemma REF . However, these estimates do not suffice to pass to the \(\varepsilon \rightarrow 0\)  limit. To this end, we assumed a weak statistical scaling of the Navier-Stokes statistical solutions (see Assumption REF ). This assumption is a weaker version of Kolmogorov's scaling assumptions in his K41 theory. Moreover, it is consistent with and incorporates different variants of scaling that are proposed in the physics literature to explain intermittent corrections to Kolmogorov's theory. Under this assumption, we proved a weak anisotropy result and invoked compactness results of [1]} to rigorously prove that the statistical solutions of the Navier-Stokes equations converge, in a suitable sense, to a statistical solution of the incompressible Euler equations. Thus, we were able to characterize the vanishing viscosity limit of the Navier-Stokes equations in a relevant regime.
At this juncture, it is essential to point that that no assumption, other than weak statistical scaling, is made in our results and all other estimates are derived rigorously. This should be contrasted with the results of [47]} where the authors directly assume a decay of the energy spectrum for the weak solutions of the Navier-Stokes equations. It is currently unclear if one can relax the weak statistical scaling assumption or even if it holds for all incompressible fluid flows. Experimental evidence strongly supports that this assumption is verified in practice, see e.g., [41]}, [46]}, [42]}.
To the best of our knowledge, the only rigorous study of the vanishing viscosity limit of the (Foiaş-Prodi) statistical solutions was carried out by Chae in [56]} where he proved that these statistical solutions converge to a measure-valued solution of the incompressible Euler equations. In contrast, we prove convergence to statistical solutions of the incompressible Euler equations and recall that statistical solutions are much more informative than measure-valued solutions as they also incorporate knowledge of all multi-point correlations.
Finally, our characterization of the vanishing viscosity limit can be viewed in connection to recent results in [57]} where the authors proved convergence of numerical spectral viscosity approximations to the statistical solutions of the Euler equations under very similar weak scaling assumptions.

appendix

Equivalence of different definitions of statistical solutions for the incompressible Navier–Stokes equations
This appendix is devoted to the proof of Theorems REF  and REF . For convenience, we restate the result:
Theorem 1.1 (Foiaş–Prodi statistical solutions are Friedman–Keller solutions) 
Let \(\mu \)  be a Foiaş–Prodi statistical solution such that the initial condition \(\mu _0\)  has bounded support, \(\text{supp}(\mu _0)\subset B\subset L^2_{\operatorname{div}}(D;U)\) , \(B=\lbrace u\in L^2(D;U)\, :\, \Vert u\Vert _{L^2(D)}\leqslant R\rbrace \)  for some \(0<R\in \mathbb {R}\)  large enough. Then \(\mu \)  is a Friedmann–Keller statistical solution (cf. Definition REF ).

It is shown in [4]} that Foiaş–Prodi statistical solutions with initial measure \(\mu _0\)  having bounded support in \(B_{L^2_{\operatorname{div}}(D;U)}(R)\lbrace u\in L^2_{\operatorname{div}}(D;U): \Vert u\Vert _{L^2(D)}\leqslant R\rbrace \)  have bounded support for all times, i.e., \(\text{supp}(\mu _t)\subset B_{L^2_{\operatorname{div}}(D;U)}(R)\) .
Therefore, we can assume that \(\mu _t\)  has uniformly bounded support. This implies in particular that \(\mu _t\)  have bounded moments:
\(\int _{L^2_{\operatorname{div}}(D;U)}\Vert u\Vert _{L^2(D)}^{2k} \,d\mu _t(u)\leqslant R^{2k}\quad \text{for a.e. } t\in [0,T].\) 
Moreover, by [4]}, we have that statistical solutions of Navier–Stokes satisfy
\(\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\left[-\partial _t\Phi (t,u) + a(u,\partial _u\Phi (t,u))+ b(u,u,\partial _u\Phi (t,u))\right]\,d\mu _t(u)\,dt\\=\int _{L^2_{\operatorname{div}}(D;U)}\Phi (0,u)\,d\mu _0(u)\) 
for any test function \(\Phi (t,u)\)  that is Fréchet differentiable on \([0,T]\times L^2_{\operatorname{div}}(D;U)\)  with \(\Phi (t,\cdot )=0\)  near \(t=T\)  and \(|\partial _u\Phi (t,u)|\leqslant C\)  and \(|\partial _t\Phi (t,u)|\leqslant C_1+C_2\Vert u\Vert _{L^2(D)}\)  for all \(u\)  and \(t\)  and some constants \(C, C_1, C_2\)  (equation (\(3.13_{\mathrm {I}}\) ) and condition (3.8) in [4]}). Therefore, we can choose test functions
\(\Phi (u)=q((u,\varphi _1),\dots ,(u,\varphi _k))\theta (t)\) 
for \(q\)  a polynomial on \(\mathbb {R}^k\)  and \(\varphi _j\in H^1_{\operatorname{div}}(D;U)\cap C^2_c(D,U)\)  and \(\theta \in C^1_c([0,T))\) , \(j=1,\dots , k\)  in (REF ), so that we get
\(\int _0^T\!\!\!\!\int _{L^2_{\operatorname{div}}(D;U)}\!\!-\theta ^{\prime }(t)q((u,\varphi _1),\dots ,(u,\varphi _k))d\mu _t(u)dt\\+ \int _0^T\!\!\!\!\int _{L^2_{\operatorname{div}}(D;U)}\!\!\Big [\theta (t)\sum _{i=1}^k\partial _i q((u,\varphi _1),\dots ,(u,\varphi _k))\left[a(u,\varphi _i)+b(u,u,\varphi _i)\right]\Big ]\,d\mu _t(u)\,dt\\=\int _{L^2_{\operatorname{div}}(D;U)}\!\!\!\!\!\theta (0)q((u,\varphi _1),\dots ,(u,\varphi _k)) \,d\mu _0(u),\) 
Note that we can integrate by parts in the terms involving \(a(u,\varphi _i)\)  and \(b(u,u,\varphi _i)\) , \(i=1,\dots , k\) , so that all the derivatives are on the test functions \(\varphi _i\) , \(i=1,\dots , k\)  and \(\theta \) :
\(0=& \int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\!\!\theta ^{\prime }(t)q((u,\varphi _1),\dots ,(u,\varphi _k))\,d\mu _t(u)\,dt+\int _{L^2_{\operatorname{div}}(D;U)}\!\!\theta (0)q((u,\varphi _1),\dots ,(u,\varphi _k)) \,d\mu _0(u)\\&\quad + \sum _{i=1}^k\int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\!\!\theta (t)\partial _i q((u,\varphi _1),\dots ,(u,\varphi _k))\\&\hphantom{\quad +\sum _{i=1}^k\int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\!\!}\times \int _{D}\left[\varepsilon u(x_i)\Delta _{x_i}\varphi _i(x_i)+(u(x_i)\otimes u(x_i)):\nabla _{x_i}\varphi _i(x_i)\right] \,dx_i\,d\mu _t(u)\,dt\) 

Now take \(q(s_1,\dots ,s_k)=s_1\cdots s_k\) , so that the last identity becomes (denote \(\,dx=\,dx_1\dots \,dx_k\) )
\(0=\int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\int _{D^k}\theta ^{\prime }(t)u(x_1)\cdot \varphi _1(x_1) \cdots u(x_k)\cdot \varphi _k(x_k) \,dx \,d\mu _t(u)\,dt\\+\sum _{i=1}^k\int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\int _{D^k}\theta (t) \big (u(x_1)\cdot \varphi _1(x_1) \cdots \varepsilon u(x_i)\cdot \Delta _{x_i}\varphi _i(x_i) \cdots u(x_k)\cdot \varphi _k(x_k)\big ) \,dx \,d\mu _t(u)\,dt\\+ \sum _{i=1}^k\int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\int _{D^k}\theta (t) \big (u(x_1)\cdot \varphi _1(x_1)\cdots (u(x_i)\otimes u(x_i)):\nabla _{x_i}\varphi _i(x_i)\cdots u(x_k)\cdot \varphi _k(x_k)\big ) \,dx \,d\mu _t(u)\,dt\\+\int _{L^2_{\operatorname{div}}(D;U)}\int _{D^k}\theta (0)u(x_1)\cdot \varphi _1(x_1) \cdots u(x_k)\cdot \varphi _k(x_k) \,dx \,d\mu _0(u), \) 
which is, denoting \(x=(x_1,\dots ,x_k)\)  and \(\varphi (t,x)\theta (t) \varphi _1(x_1)\otimes \dots \otimes \varphi _k(x_k)\) , equivalent to
\(0=&\ \int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\int _{D^k}(u(x_1)\otimes \dots \otimes u(x_k)):\partial _t\varphi (t,x) \,dx \,d\mu _t(u)\,dt\\&+ \varepsilon \sum _{i=1}^k\int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\int _{D^k}(u(x_1)\otimes \dots \otimes u(x_k)):\Delta _{x_i}\varphi (t,x) \,dx \,d\mu _t(u)\,dt\\&+ \sum _{i=1}^k\int _0^T\!\!\int _{L^2_{\operatorname{div}}(D;U)}\int _{D^k} (u(x_1)\otimes \dots \otimes (u(x_i)\otimes u(x_i))\otimes \dots \otimes u(x_k)):\nabla _{x_i}\varphi (t,x) \,dx \,d\mu _t(u)\,dt\\&+\int _{L^2_{\operatorname{div}}(D;U)}\int _{D^k}(u(x_1)\otimes \dots \otimes u(x_k)):\varphi (0,x) \,dx \,d\mu _0(u), \) 

Since we assume that the support of \(\mu _0\)  is bounded (and therefore also the support of \(\mu _t\)  for almost all \(t\in [0,T]\) ), we can use a density argument to conclude that the above identity holds for all \(\varphi \in L^2([0,T];(H^1_{\operatorname{div}}(D;U))^k)\cap C^2_c([0,T]\times D^k;U^k)\) , that is, all \(\varphi \in C^2_c([0,T]\times D^k;U^k)\)  with \(\operatorname{div}_{x_i}\varphi =0\)  for all \(i=1,\dots ,k\) . Observe next that all the terms in the above identity have the form required to apply Theorem REF  (recall that \(\mu _t\)  has bounded support). Therefore, by the same theorem, there exists a unique correlation measure \(\mathbf {\nu }= (\mathbf {\nu }_t)_{0\leqslant t\leqslant T}=\big ((\nu ^1_{t},\nu ^2_t,\dots )\big )_{0\leqslant t\leqslant T}\)  corresponding to \(\lbrace \mu _t\rbrace \)  satisfying
\(0=&\ \int _0^T\!\!\int _{D^k}\int _{U^k}(\xi _1\otimes \dots \otimes \xi _k):\partial _t\varphi (t,x) \,d\nu _{t,x}^k(\xi )\,dx\,dt\\&+ \varepsilon \sum _{i=1}^k\int _0^T\!\!\int _{D^k}\int _{U^k}(\xi _1\otimes \dots \otimes \xi _k):\Delta _{x_i}\varphi (t,x) \,d\nu _{t,x}^k(\xi )\,dx\,dt\\&+ \sum _{i=1}^k\int _0^T\!\!\int _{D^k}\int _{U^k} (\xi _1\otimes \dots \otimes (\xi _i\otimes \xi _i)\otimes \dots \otimes \xi _k):\nabla _{x_i}\varphi (t,x) \,d\nu _{t,x}^k(\xi )\,dx \,dt\\&+\int _{D^k}\int _{U^k}(\xi _1\otimes \dots \otimes \xi _k):\varphi (0,x) \,d\nu _{0,x}^k(\xi )\,dx,\) 

which is (REF ).
By previous arguments, each measure \(\mu _t\)  has bounded support (and in particular bounded moments), so we can take any nonnegative, nondecreasing polynomial \(\psi (s)=p_K(s)\) , of degree \(K\) , for \(K\in \mathbb {N}\)  in the strengthened energy inequality (REF ), so that we get
\(\int _{L^2(D)}p_K( \Vert u\Vert _{L^2(D)}^{2}) \,d\mu _t(u) + 2\varepsilon \int _0^t\int _{L^2(D)}p_K^{\prime }(\Vert u\Vert _{L^2(D)}^{2})|u|_{H^1(D)}^2 \,d\mu _s(u)\,ds\\\leqslant \int _{L^2(D)} p_K(\Vert u\Vert _{L^2(D)}^{2}) \,d\mu _0(u)\qquad \text{for a.e. } t\in [0,T].\) 
Thanks to the bounded support of \(\mu _t\) , all the moment terms have the form required to apply Theorem REF , moreover, with Lemma REF , we can rewrite the gradient term such that we obtain for a.e. \(t\in [0,T]\) , the energy inequality (REF ). Indeed, we have with Theorem  REF
\(\int _{L^2(D)}p_K( \Vert u\Vert _{L^2(D)}^{2}) \,d\mu _t(u) &= \sum _{k=0}^Ka_k \int _{L^2(D)}\Vert u\Vert _{L^2(D)}^{2k} \,d\mu _t(u)\\&= \sum _{k=0}^Ka_k \int _{L^2(D)}\int _{D^k} |u(x_1)|^2\dots |u(x_k)|^2 dx \,d\mu _t(u)\\&= \sum _{k=0}^Ka_k \int _{D^k}\int _{U^k} |\xi _1|^2\dots |\xi _k|^2 d\nu _{t,x}^k(\xi ) dx.\) 

Similarly, with Lemma REF , (via Lemma REF )
\(&\int _0^t\int _{L^2(D)}p_K^{\prime }\big (\Vert u\Vert _{L^2(D)}^{2}\big )|u|_{H^1(D)}^2 \,d\mu _s(u)\,ds\\&\qquad = \sum _{k=0}^Ka_k k\int _0^t\int _{L^2(D)}\Vert u\Vert _{L^2(D)}^{2(k-1)}|u|_{H^1(D)}^2 \,d\mu _s(u) \,ds\\&\qquad = \sum _{k=0}^Ka_k k\int _0^t\int _{L^2(D)}\int _{D^{k}}|u(x_1)|^2\dots |u(x_{k-1})|^2 |\nabla u(x_k)|^2 dx \,d\mu _s(u)\,ds\\&\qquad = \sum _{k=0}^Ka_k\sum _{i=1}^k \int _0^t\int _{L^2(D)}\int _{D^{k}}|u(x_1)|^2\dots |\nabla u(x_i)|^2\dots |u(x_{k})|^2 dx \,d\mu _s(u)\,ds\\&\qquad = \sum _{k=0}^Ka_k\sum _{i=1}^k \sum _{j=1}^d\lim _{h\rightarrow 0}\frac{1}{h^2}\int _0^t\int _{D^{k}}\int _{U^{k+1}}|\xi _1|^2\dots |\xi _i-\xi _{k+1}|^2\dots |\xi ^k|^2 d\nu ^{k+1}_{t,(x,x_i+h\mathbf {e}_j)}(\xi ,\xi _{k+1}) dx\,ds\) 

and a similar computation for the term involving \(\mu _0\)  yields (REF ).
The bounded support of \(\mu _0\)  implies, using the equivalence theorem [2]},
\(\int _{D^k}\!\!\int _{U^{k}}|\xi _1|^2\dots |\xi _k|^2 \,d\nu _{0,x}^k(\xi )\,dx<\infty \) 
for any \(k\in \mathbb {N}\)  and thus the boundedness of all terms in (REF ).
It remains to show that each \(\nu ^k\)  satisfies (REF ). We let \(\varphi \in C^1_c(D^k;U^k)\) , \(1\leqslant \ell \leqslant k\) , \(\alpha _j\in C(U;U)\) , with \(\alpha _j(v)\leqslant C(1+|v|^2)\) , \(j=\ell ,\dots , k\)  and compute using the equivalence theorem REF ,
\(\begin{split}& \int _{D^{k}}\int _{U^{k}} \xi _1\otimes \dots \otimes \xi _{\ell }\otimes \alpha _{\ell +1}(\xi _{\ell +1})\otimes \dots \otimes \alpha _k(\xi _k) \,d\nu _{t,x}^k(\xi )\cdot \nabla _{x_1,\dots , x_\ell }\varphi (x) \,dx\\& \qquad = \int _{L^2(D)}\int _{D^k}\!\!\! \begin{aligned}[t]&u(x_1)\otimes \dots \otimes u(x_\ell )\otimes \alpha _{\ell +1}(u(x_{\ell +1}))\otimes \dots \otimes \alpha _k(u(x_k))\\&\cdot \nabla _{x_1,\dots , x_\ell }\varphi (x) \,dx \,d\mu _t(u)\end{aligned}\\&\qquad = 0,\end{split}\) 
since \(\mu _t\)  is supported on divergence free functions. This concludes the proof.
Let us show the reverse direction now, that is, that any correlation measure that solves the Friedman–Keller system and satisfies in addition an energy inequality, is a statistical solution of the Navier–Stokes equations in the sense of Foiaş–Prodi.
Theorem 1.2 (Friedman–Keller solutions are Foiaş–Prodi solutions) 
Let \(\mathbf {\nu }=(\mathbf {\nu }_t)_{0\leqslant t\leqslant T}\)  be a Friedman–Keller statistical solution of Navier–Stokes (cf. Definition REF ) with bounded support, i.e.,
\(\int _{D^k}\int _{U^{k}} |\xi _1|^2\dots |\xi _k|^2 \,d\nu _{t,x}^k(\xi ) \,dx\leqslant R^k<\infty ,\) 
for some \(0<R<\infty \) , every \(k\in \mathbb {N}\) , and almost every \(t\in [0,T]\) . Then \(\mathbf {\nu }\)  corresponds to a probability measure \(\mu =(\mu _t)_{0\leqslant t\leqslant T}\)  on a bounded set of \(L^2_{\operatorname{div}}(D;U)\)  which is a Foiaş–Prodi statistical solution of the Navier–Stokes equations (cf. Definition REF ).

From the equivalence theorem REF , we obtain that \((\mathbf {\nu }_t)_{0\leqslant t\le T}\)  corresponds to a family of measures \((\mu _t)_{0\leqslant t\leqslant T}\subset \mathcal {P}(L^2(D;U))\)  with bounded support, that is \(\text{supp}(\mu _t)\subset \lbrace u\in L^2(D;U)\,:\, \Vert u\Vert _{L^2(D;U)}\leqslant R\rbrace \)  for almost every \(t\in [0,T]\) . Property (a) of Definition REF  – the fact that (REF ) is measurable for all \(\varphi \in C_b(L^2)\)  – follows from a monotone class argument, which we include here. By property (i) of Theorem REF , the function
\(t\mapsto \int _{L^2(D;U)}L_f(u) \,d\mu _t(u) = \int _{L^2(D;U)}\int _D f(x,u(x)) \,dx\,d\mu _t(u)\) 
is measurable for every \(f\in L^1(D^k,C_0(U^k))\) . Let \(\mathbf {M}\)  be the collection of sets
\(\mathbf {M} = \left\lbrace E\in \mathcal {B}(L^2(D;U)) \text{ such that } t \mapsto \int _{L^2(D;U)} \mathbb {1}_E(u)\,d\mu _t(u) \text{ is measurable} \right\rbrace .\) 
By the monotone convergence theorem, \(\mathbf {M}\)  is a monotone class, that is, it is closed under (countable) unions of increasing sequences of sets and intersections of decreasing sequences of sets. By the same argument as in the proof of [2]}, \(\mathbf {M}\)  contains the collection of cylinder sets \(\mathrm {Cyl}(L^2) = \lbrace \text{all cylinder sets on }L^2(D;U)\rbrace \) , that is, all sets of the form \(E = \big \lbrace u\in L^2\ :\ (\langle \varphi _1,\, u\rangle ,\dots ,\langle \varphi _n,\, u\rangle )\in F\big \rbrace \)  for some \(n\in \mathbb {N}\) , a Borel set \(F\subset \mathbb {R}^n\)  and \(\varphi _1,\dots ,\varphi _n\in L^2(D;U)\) . Since \(\mathrm {Cyl}(L^2)\)  is an algebra which generates \(\mathcal {B}(L^2)\)  (cf. e.g. [2]}), it follows from the monotone class lemma that \(\mathbf {M} = \sigma (\mathrm {Cyl}(L^2)) = \mathcal {B}(L^2)\) . Approximating an arbitrary \(\varphi \in C_b(L^2(D;U))\)  by simple functions now gives the desired conclusion.
We claim that \(\mu _t\)  is supported on \(L^2_{\operatorname{div}}(D;U)\) . This follows from Lemma REF : Since \(\mu \)  has bounded support on \(L^2(D)\)  we may take \(\varphi (\xi )=|\xi |^2\)  in that lemma, which is continuous and bounded on any compact subset of \(\mathbb {R}\)  and satisfies \(\varphi (0)=0\) . Then by the lemma, for any \(g\in H^1(D)\) ,
\(\int _{L^2(D)}\left|\int _{D} u\cdot \nabla g \,dx\right|^2\,d\mu _t(u) =0.\) 
Hence, by Chebychev's inequality,
\(\mu _t\left(\left\lbrace u\in L^2(D): \,\int _{D} u\cdot \nabla g \,dx\ne 0\right\rbrace \right)=0.\) 
Since \(g\in H^1(D;U)\)  was arbitrary and \(H^1(D;U)\)  is separable, this implies that \(\mu _t\)  is supported on \(L^2\) -functions that are weakly divergence free, which is exactly the space \(L^2_{\operatorname{div}}(D;U)\)  (see e.g. [64]}).
Next we claim that \(\mu \)  satisfies condition (c) in Definition REF . As \(\mathbf {\nu }\)  is assumed to satisfy (REF ), we can apply [2]} combined with Lemma REF  to each of the terms and obtain that \(\mu _t\)  satisfies (REF ) for any nonnegative and nondecreasing polynomial \(p_K(s)=\sum _{k=0}^K a_k s^k\)  on \([0,R]\) . This implies in particular that \(\mu _t\)  is supported on \(H^1(D;U)\) , for a.e. \(t\in [0,T)\) , and (REF ). Now any differentiable nondecreasing function \(\psi \)  on a bounded interval can be approximated by nondecreasing polynomials [66]}, from which (REF ) follows after passing to the limit in a suitable polynomial approximation.
Specifically, for a given \(\psi \) , let \(\lbrace p_{K_n}\rbrace _{n\in \mathbb {N}}\)  be a sequence of nonnegative, nondecreasing polynomials (\(K_n\rightarrow \infty \) ) satisfying
\(\left\Vert  p_{K_n}-\psi  \right\Vert _{C^1([0,R])}\leqslant \frac{1}{n}.\) 
Then, for \(t\geqslant 0\) , by the compact support property of \(\mu \)  on \(B\) ,
\(\left|\int _{L^2(D)}\psi (\left\Vert  u \right\Vert _{L^2(D)}^2)d\mu _t(u)-\int _{L^2(D)} p_{K_n}(\left\Vert  u \right\Vert ^2_{L^2(D)})d\mu _t(u)\right|\leqslant \left\Vert  p_{K_n}-\psi  \right\Vert _{C^1([0,R])}\leqslant \frac{1}{n},\) 
and
\(&\left|\int _0^t\int _{L^2(D)}\psi ^{\prime }(\left\Vert  u \right\Vert _{L^2(D)}^2) |u|^2_{H^1(D)}d\mu _t(u)ds-\int _0^t\int _{L^2(D)} p_{K_n}^{\prime }(\left\Vert  u \right\Vert ^2_{L^2(D)})|u|^2_{H^1(D)} d\mu _t(u)ds \right|\\& \qquad \leqslant \int _0^t\int _{L^2(D)}\left|\psi ^{\prime }(\left\Vert  u \right\Vert _{L^2(D)}^2)-p_{K_n}^{\prime }(\left\Vert  u \right\Vert ^2_{L^2(D)})\right| |u|^2_{H^1(D)}d\mu _t(u)ds\\&\qquad \leqslant \left\Vert  p_{K_n}-\psi  \right\Vert _{C^1([0,R])}\int _0^t\int _{L^2(D)}|u|^2_{H^1(D)}d\mu _t(u)ds\\&\qquad \leqslant \frac{C}{n},\) 

where the last inequality follows from (REF ).
It remains to show that the correlation measures satisfy the evolution equation (REF ) for all \(\Phi \in \mathcal {T}_{\text{cyl}}^0\) . From this also part (d) of the definition will follow (see Remark REF ). Let \(\Phi (u) = \varphi ((u,g_1),\dots ,(u,g_k))\)  be an arbitrary function in \(\mathcal {T}_{\text{cyl}}^0\)  with \(g_j\in H^1(D;U)\)  and let \(\widetilde{R}=\max _{1\leqslant j\leqslant k}(\left\Vert  g_j \right\Vert _{L^2})\) . Since \(\varphi \)  is continuously differentiable and bounded on \(B[-R\widetilde{R}-1,R\widetilde{R}+1]^k\) , we can approximate it arbitrarily well by polynomials thanks to the Weierstrass approximation theorem. Let
\(p_n(\zeta )\sum _{|\underline{j}|=0}^{N_n}\beta _{\underline{j}}\zeta _1^{j_1}\dots \zeta _k^{j_k},\quad j_i\geqslant 0,\quad \underline{j}=(j_1,\dots ,j_k),\quad |\underline{j}|j_1+\dots +j_k,\) 
for some \(N_n\in \mathbb {N}\)  large enough, be approximations of \(\varphi \)  that satisfy
\(\left\Vert  \varphi -p_n \right\Vert _{C^1(B)}\leqslant \frac{1}{n},\quad n\in \mathbb {N}.\) 
Let \(\theta \in C_c^1((0,T))\)  be an arbitrary compactly supported test function. Since by the equivalence theorem, for any \(\underline{j}=(j_1,\dots , j_k)\) , \(j_i\geqslant 0\) , \(|\underline{j}|=j_1+\dots +j_k\) ,
\(&\int _0^T\theta ^{\prime }(t)\int _{D^{|\underline{j}|}}\left(\int _{U^{|\underline{j}|}}\xi _{1}\otimes \cdots \otimes \xi _{|\underline{j}|} \,d\nu _{t,x}^{|\underline{j}|}(\xi )\right) \\&: \Big (g_1(x_1)\otimes \cdots \otimes g_1(x_{j_1})\otimes \cdots \otimes g_k(x_{|\underline{j}|-j_k+1})\otimes \cdots \otimes g_k(x_{|\underline{j}|})\Big ) \,dx \,dt\\&= \int _0^T\theta ^{\prime }(t) \int _{L^2(D)}(u,g_1)^{j_1}\dots (u,g_k)^{j_k} \,d\mu _t(u) \,dt,\) 

and if \(\widetilde{j}=j_1+\dots + j_{i-1}\) , then
\(&\sum _{\ell =\widetilde{j}+1}^{\widetilde{j}+j_i}\int _0^T\!\!\!\theta (t)\!\int _{D^{|\underline{j}|}}\left(\int _{U^{|\underline{j}|}} \xi _{1}\otimes \cdots \otimes \xi _{|\underline{j}|} \,d\nu _{t,x}^{|\underline{j}|}(\xi )\right)\\&: \Big (g_1(x_1)\otimes \cdots \otimes g_i(x_{\widetilde{j}+1})\otimes \dots \otimes \Delta _{x_{\ell }}g_i(x_{\ell })\otimes \dots \otimes g_i(x_{\widetilde{j}+j_i})\otimes \dots \otimes g_k(x_{|\underline{j}|})\Big ) \,dx \,dt\\=&\ j_i\int _0^T\!\!\theta (t)\int _{L^2(D;U)}(u,g_1)^{j_1}\cdots (u, g_i)^{j_i-1}(u,\Delta g_i)\cdots (u,g_k)^{j_k} \,d\mu _t(u) \,dt\\=&\ j_i\int _0^T\!\!\theta (t)\int _{L^2(D;U)}(u,g_1)^{j_1}\cdots (u,g_i)^{j_i-1}(Au,g_i)\cdots (u,g_k)^{j_k} \,d\mu _t(u) \,dt,\) 

where the last equality follows because \(g_i\)  is divergence free because \(\mu _t\)  is supported on \(H^1(D;U)\) -functions by the energy inequality (REF ). Furthermore, we have
\(&\sum _{\ell =\widetilde{j}+1}^{\widetilde{j}+j_i}\int _0^T\!\!\!\theta (t)\!\int _{D^{|\underline{j}|}}\left(\int _{U^{|\underline{j}|}}\xi _{1}\otimes \cdots \otimes \xi _{\widetilde{j}+1}\otimes \cdots \otimes (\xi _{\ell }\otimes \xi _{\ell })\otimes \cdots \otimes \xi _{\widetilde{j}+j_i}\otimes \cdots \otimes \xi _{|\underline{j}|} \,d\nu _{t,x}^{|\underline{j}|}(\xi )\right)\\&: \Big ( g_1(x_1)\otimes \cdots \otimes \nabla _{x_{\ell }}g_i(x_{\ell })\otimes \cdots \otimes g_k(x_{|\underline{j}|})\Big ) \,dx \,dt\\=&\ j_i\int _0^T\!\!\theta (t) \int _{L^2(D)}(u,g_1)^{j_1}\cdots (u, g_i)^{j_i-1}(u\otimes u,\nabla g_i)\cdots (u,g_k)^{j_k} \,d\mu _t(u) \,dt\\=&\ -j_i\int _0^T\!\!\theta (t) \int _{L^2(D)}(u,g_1)^{j_1}\cdots (u,g_i)^{j_i-1}(B(u),g_i)\cdots (u,g_k)^{j_k} \,d\mu _t(u) \,dt,\) 

again, the last equality following because \(g_i\)  is divergence free and \(\mu _t\)  is supported on \(H^1(D;U)\) -functions for almost every \(t\in [0,T]\) .
We obtain, by combining the terms to the Friedman–Keller system, that if
\(\Phi _n(u) p_n\big ((u,g_1),\dots ,(u,g_k)\big )\) 
then \(\mu \)  satisfies the equation
\(\int _0^T\int _{L^2(D)}\theta ^{\prime }(s)\Phi _n(u) \,d\mu _s(u)\,ds + \int _0^T\theta (s)\int _{L^2(D)} (F(s,u),\partial _u\Phi _n(u)) \,d\mu _s(u) \,ds=0\) 
for every test function \(\theta \in C_c^1((0,T))\) .
Next, note that
\(\begin{split}&\left|\int _{L^2(D)}\Phi _{n}(u)\,d\mu _t(u) - \int _{L^2(D)}\Phi (u) \,d\mu _t(u)\right|=\left|\int _{L^2(D)}\left(\Phi _{n}(u)-\Phi (u)\right) \,d\mu _t(u)\right|\\&\quad \leqslant \int _{L^2(D)}\big |p_n\big ((u,g_1),\dots ,(u,g_k)\big )-\varphi \big ((u,g_1),\dots ,(u,g_k)\big )\big | \,d\mu _t(u)\\&\quad \leqslant \int _{L^2(D)}\sup _{\zeta \in B}\left|p_n(\zeta )-\varphi (\zeta )\right| \,d\mu _t(u)\\&\quad = \left\Vert  \varphi -p_n \right\Vert _{C^0(B)} \leqslant \frac{1}{n},\end{split}\) 
because \(\mu _t\)  has support on \(\lbrace u\in L^2(D)\,:\, \Vert u\Vert _{L^2(D)}\leqslant R\rbrace \) . Moreover,
\(\begin{split}&\left|\int _0^T\theta (s)\int _{L^2(D)}\big (\partial _u\Phi _n(u),F(s,u)\big )\,d\mu _s(u)\,ds - \int _0^T\theta (s)\int _{L^2(D)}\big (\partial _u\Phi (u),F(s,u)\big ) \,d\mu _s(u)\,ds\right|\\&\quad \leqslant \Vert \theta \Vert _{L^\infty }\int _0^T\int _{L^2(D)}\big |\big (\partial _u(\Phi (u)-\Phi _n(u)),F(s,u)\big )\big | \,d\mu _s(u)\,ds\\&\quad \leqslant \Vert \theta \Vert _{L^\infty }\sum _{j=1}^k\int _0^T\int _{L^2(D)}\big |\partial _j\varphi \big ((u,g_1),\dots ,(u,g_k)\big )-\partial _jp_n\big ((u,g_1),\dots ,(u,g_k)\big )\big |\big |(g_j,F(s,u))\big | \,d\mu _s(u)\,ds\\&\quad \leqslant \Vert \theta \Vert _{L^\infty }\sum _{j=1}^k\sup _{B}\left|\partial _j(\varphi -p_n)\right|\int _0^T\int _{L^2(D)}\big |(g_j,F(s,u))\big | \,d\mu _s(u)\,ds\\&\quad \leqslant \Vert \theta \Vert _{L^\infty }\left\Vert  \varphi -p_n \right\Vert _{C^1(B)}\sum _{j=1}^k\int _0^T\int _{L^2(D)}\big |(g_j,F(s,u))\big | \,d\mu _s(u)\,ds\\&\quad \leqslant \frac{1}{2n}\Vert \theta \Vert _{L^\infty } \sum _{j=1}^k\left(\left\Vert  g_j \right\Vert _{H^1(D)}^2+\int _0^T\int _{L^2(D)}\left\Vert  \nabla u \right\Vert _{L^2(D)}^2 \,d\mu _s(u)\,ds \right)\end{split}\) 
where we used Young's inequality for the last inequality.
Thanks to (REF ) and (REF ), we can pass to the limit \(n\rightarrow \infty \)  in (REF ) and obtain
\(\int _0^T\theta ^{\prime }(s)\int _{L^2(D)}\Phi (u) \,d\mu _s(u)\,ds + \int _0^T\theta (s)\int _{L^2(D)} (F(s,u),\partial _u\Phi (u)) \,d\mu _s(u) \,ds=0\) 
for every \(\theta \in C_c^1((0,T))\) . It follows that the distributional derivative of \(\varphi (s)\int _{L^2(D)}\Phi (u) \,d\mu _s(u)\)  is
\(\frac{d}{ds}\varphi (s) = \int _{L^2(D)}(F(s,u),\partial _u\Phi (s,u))\,d\mu _s(u),\) 
and since the right-hand side lies in \(L^1((0,T))\)  we see that \(\varphi \)  is absolutely continuous. Consequently, we can (through a standard approximation procedure) insert \(\theta =\mathbb {1}_{(0,t)}\)  in (REF ) and conclude that (REF ) is true.

Auxiliary lemmata
Lemma 2.1 (Representation of the gradient) 

Let \((\mu _t)_{0\leqslant t\leqslant T}\)  be a measure on \(L^2(D)\)  satisfying
\(\int _0^T\int _{L^2(D)}\left\Vert  \nabla u \right\Vert _{L^2(D)}^2 \,d\mu _t(u) \,dt\leqslant C<\infty .\) 
Define for \(h\in \mathbb {R}\)  the finite difference gradient \(\nabla _h\)  by
\(\nabla _h f(x) = \left( D^h_1 f(x),\dots , D^h_d f(x)\right), \quad D^h_j f(x) = \frac{f(x+ h \mathbf {e}_j)- f(x)}{h},\) 
where \(\mathbf {e}_j\in \mathbb {R}^d\)  is the \(j\) th unit vector. Then we have for any \(V\subset \subset D\)  (or in the case that \(D\)  is the torus also for \(V=D\) )
\(\int _0^T\int _{L^2(D)}\left\Vert  \nabla u \right\Vert _{L^2(V)}^2 \,d\mu _t(u) \,dt = \lim _{h\rightarrow 0} \int _0^T\int _{L^2(D)}\left\Vert  \nabla _h u \right\Vert _{L^2(V)}^2 \,d\mu _t(u) \,dt.\) 

Consider a compact \(V\subset \subset D\)  and a smooth \(u\in C^\infty (D)\) . Then, since we can write
\(D^h_j u(x) = \int _0^1 \partial _{x^j}u(x+t h \mathbf {e}_j) \,dt,\) 
and hence
\(|\nabla _h u(x)|\leqslant \sum _{j=1}^d\int _{0}^{1}|\partial _{x^j}u(x+t h \mathbf {e}_j)| \,dt\leqslant \sum _{j=1}^{d}\left\Vert  \partial _{x^j}u \right\Vert _{L^\infty (V)}\leqslant C\left\Vert  \nabla u \right\Vert _{L^\infty (V)},\) 
for all \(x\in V\)  and \(\int _{V}C\left\Vert  \nabla u \right\Vert _{L^\infty (V)} \,dx<\infty \) , we can use the dominated convergence theorem to obtain that
\(\int _{V}|\nabla u|^2 \,dx =\lim _{h\rightarrow 0}\int _V |\nabla _h u|^2 \,dx.\) 
Moreover, by [67]}, we have for any \(u\in H^1(D)\) ,
\(\left\Vert  \nabla _h u \right\Vert _{L^2(V)}\leqslant C_V \left\Vert  \nabla u \right\Vert _{L^2(D)},\) 
for some constants \(0\leqslant c_V,C_V<\infty \)  and any \(V\subset \subset D\) .
Since \(C^\infty (D)\)  is dense in \(H^1(D)\) ,
we can find a sequence of functions \((u_n)_{n\in \mathbb {N}}\subset C^\infty (D)\)  such that
\(\lim _{n\rightarrow \infty }\left\Vert  u_n-u \right\Vert _{H^1(V)} = 0,\) 
for any \(V\subset D\)  bounded.
So fix \(V\subset \subset D\)  and an arbitrary \(\varepsilon >0\)  and pick \(N_\varepsilon \in \mathbb {N}\)  large enough such that for all \(n\geqslant N_\varepsilon \) ,
\(\left\Vert  u-u_n \right\Vert _{H^1(V)}\leqslant \frac{1}{2}\frac{1}{1+C_V}\varepsilon ,\) 
where \(C_V\)  is the constant in (REF ).
Then for \(n\geqslant N_\varepsilon \)
\(\left\Vert  \nabla u - \nabla _h u \right\Vert _{L^2(V)}&\leqslant \left\Vert  \nabla (u-u^n) \right\Vert _{L^2(V)}+ \left\Vert  \nabla u^n -\nabla _h u^n \right\Vert _{L^2(V)}+ \left\Vert  \nabla _h (u^n-u) \right\Vert _{L^2(V)}\\&\leqslant \left\Vert  \nabla (u-u^n) \right\Vert _{L^2(V)}+ \left\Vert  \nabla u^n -\nabla _h u^n \right\Vert _{L^2(V)}+ C_V\left\Vert  \nabla (u^n-u) \right\Vert _{L^2(V)}\\&\leqslant \frac{\varepsilon }{2}+\left\Vert  \nabla u^n -\nabla _h u^n \right\Vert _{L^2(V)}.\) 

Now pick \(h_0\)  with \(|h_0|\)  small enough that
\(\left\Vert  \nabla u^n -\nabla _h u^n \right\Vert _{L^2(V)}\leqslant \frac{\varepsilon }{2}\) 
for all \(h\)  with \(|h|\leqslant |h_0|\) .
Then
\(\left\Vert  \nabla u - \nabla _h u \right\Vert _{L^2(V)}\leqslant \varepsilon \) 
for all \(h\)  with \(|h|\leqslant |h_0|\) . Since \(\varepsilon >0\)  was arbitrary, we obtain
\(\lim _{h\rightarrow 0} \int _V |\nabla _h u|^2 \,dx = \int _V |\nabla u|^2 \,dx.\) 
Combining (REF ) with (REF ), we can apply the dominated convergence theorem and use (REF ) to pass to the limit:
\(\begin{split}\int _0^T\int _{L^2(D)}\left\Vert  \nabla u \right\Vert _{L^2(V)}^2 \,d\mu _t(u) \,dt& = \int _0^T\int _{L^2(D)}\lim _{h\rightarrow 0}\left\Vert  \nabla _h u \right\Vert _{L^2(V)}^2 \,d\mu _t(u) \,dt\\&= \lim _{h\rightarrow 0} \int _0^T\int _{L^2(D)}\left\Vert  \nabla _h u \right\Vert _{L^2(V)}^2 \,d\mu _t(u) \,dt.\end{split}\) 
Corollary 2.2 
Let \((\mu _t)_{0\leqslant t\leqslant T}\)  be a measure on \(L^q(D)\) . If \(\mu _t\)  satisfies
\(\int _0^T\int _{L^q(D)}\left\Vert  \nabla u \right\Vert _{L^q(D)}^q \,d\mu _t(u) \,dt\leqslant C<\infty ,\) 
then for any \(V\subset \subset D\)  (or if \(D\)  is the torus also for \(V=D\) )
\(\int _0^T\int _{L^q(D)}\left\Vert  \nabla u \right\Vert _{L^q(V)}^q \,d\mu _t(u) \,dt =\lim _{h\rightarrow 0}\int _0^T\int _{L^q(D)}\left\Vert  \nabla _h u \right\Vert _{L^q(V)}^q \,d\mu _t(u) \,dt.\) 

Lemma 2.3 (More representation of gradients) 
Let \((\mu _t)_{0\leqslant t\leqslant T}\)  be a measure on \(L^q(D)\) , where \(D\)  is the \(d\) -dimensional torus satisfying
\(\int _0^T\int _{L^q(D)} \left\Vert  \nabla u \right\Vert _{L^q(D)}^q \,d\mu _t(u)\,dt\leqslant C<\infty .\) 
Then we can represent the integral of the gradient in terms of correlation measures as
\(\int _0^T\int _{L^q(D)} \left\Vert  \nabla u \right\Vert _{L^q(D)}^q \,d\mu _t(u)\,dt = \lim _{h\rightarrow 0}\frac{1}{h^q}\int _0^T \int _{D}\sum _{j=1}^d\int _{U^2}|\xi _{1}-\xi _2|^q \,d\nu ^2_{x,x+h\mathbf {e}_j}(\xi ) \,dx \,dt.\) 

From Lemma REF , we know that we can represent the left hand side of (REF ) as the limit as \(h\)  goes to zero of
\(\frac{1}{h^q}\int _0^T\int _{L^q(D)}\int _{D}|u(x+h)-u(x)|^q \,dx \,d\mu _t(u) \,dt.\) 
Note that for arbitrary nonnegative functions \(\varphi \in L^1(D)\)  with \(\varphi _\varepsilon (x)\varepsilon ^{-d}\varphi (x/\varepsilon )\) , we have that
\(\int _{D}\int _{D}|f(x)-f(x-y)|^q\varphi _\varepsilon (y) \,dy \,dx\stackrel{\varepsilon \rightarrow 0}{\longrightarrow } 0,\) 
for any \(f\in L^q(D)\) , since shifts are continuous in \(L^p\)  for \(L^p\) -functions.
In particular, choosing \(\varphi (x)=\frac{1}{|B_1(0)|}\mathbb {1}_{B_1(0)}(x)\)  this yields
\(\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x) |f(x)-f(y)|q  dy  dx0 0,

and
\(\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x+hej) |f(x)-f(y)|q  dy  dx0 D |f(x)-f(x+hej)|q  dx,

for any fixed \(h>0\)  and \(\mathbf {e}_j\)  is the \(j\) th unit vector.
We have for any \(\varepsilon >0\) ,
\(&\Biggl |\frac{1}{h}\left(\int _0^T\!\!\!\int _{L^q(D)}\int _{D}\!\!\!|u(x)-u(x+h\mathbf {e}_j)|^q \,dx\,d\mu _t(u) \,dt\right)^{1/q}\\&\quad \qquad -\frac{1}{h}\left(\int _0^T\!\!\!\int _{L^q(D)}\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\right.\) 






B(x+hej) |u(x)-u(y)|q  dy  dx dt(u)  dt)1/q|
   (1hq0TLq(D)D







B(x+hej) |u(x+hej)-u(y)|q  dy  dx dt(u)  dt)1/q
   = (1hq0TLq(D)D







B(x) |u(x)-u(y)|q  dy  dx dt(u)  dt)1/q 0 0,

by (REF ) and Lebesgue dominated convergence theorem.
Hence
\(\frac{1}{h^q}\int _0^T\!\!\!\int _{L^q(D)}\int _{D}\!\!\!|u(x)-u(x+h\mathbf {e}_j)|^q \,dx\,d\mu _t(u) \,dt \\=\frac{1}{h^q}\lim _{\varepsilon \rightarrow 0}\int _0^T\!\!\!\int _{L^q(D)}\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x+hej) |u(x)-u(y)|q  dy  dx dt(u)  dt.

Now we can apply [2]} to rewrite the last expression in terms of the unique correlation measure \(\mathbf {\nu }_t\)  to which \(\mu _t\)  corresponds:
\(\frac{1}{h^q}\lim _{\varepsilon \rightarrow 0}\int _0^T\!\!\!\int _{L^q(D)}\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x+hej) |u(x)-u(y)|q  dy  dx dt(u)  dt
= 1hq00TD







B(x+hej)U2 |1-2|q  dx,y2()  dy  dx  dt

Because of the equivalence theorem [2]} and assumption (REF ) combined with Corollary REF , we have
\(\lim _{\varepsilon \rightarrow 0}\frac{1}{h^q}\int _0^T\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x+hej)U2|1-2|q  dx,y2()  dy  dx  dt<,

uniformly in \(0<h<h_0\)  for a small enough \(h_0\) .
Since
\(&\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x+hej)U2|1-2|q  dx,y2()  dy  dx
      C D







B(x+hej)U2(|1|q+|2|q)  dx,y2()  dy  dx
consistency=2 C DU|1|q  dx1(1)  dx C<    for almost every  t,

we can apply the dominated convergence theorem and pass the limit in \(\varepsilon \)  inside:
\(\lim _{\varepsilon \rightarrow 0}\frac{1}{h^q}\int _0^T\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x+hej)U2|1-2|q  dx,y2()  dy  dx  dt
= 1hq0T0D







B(x+hej)U2|1-2|q  dx,y2()  dy  dx  dt
Now fix \(\varepsilon >0\) , a `good' \(t\in [0,T]\) ,
and note that for almost every \(x\in D\) ,
\(&\left|\left(\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\right.\right.{{\textstyle {\scriptstyle -}{\int } }\hbox{$\scriptstyle -$ }}\hspace{0.0pt}\) 





B(x+hej) U2 |1-2|q  dx,y2()  dy  dx)1/q - (D U2 |1-2|q  dx,x+hej2()  dx)1/q |
      consistency=|(D







B(x+hej) U3 |1-2|q  dx,y,x+hej3()  dy  dx)1/q
      consistency=|-(D







B(x+hej) U3 |1-3|q  dx,y,x+hej3() dy  dx)1/q|
      -ineq.(D







B(x+hej) U3 |3-2|q  dx,y,x+hej3()  dy  dx)1/q
      consistency=(D







B(x+hej) U2 |2-1|q  dy,x+hej2()  dy  dx)1/q

The last term goes to zero by the diagonal continuity property of \(\mathbf {\nu }\)  (c.f. [2]}). Hence, for almost every \(t\in [0,T]\) ,
\(\lim _{\varepsilon \rightarrow 0}\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x+hej) U2 |1-2|q  dx,y2()  dy  dx =D U2 |1-2|q  dx,x+hej2()  dx.

Thus, combining (REF ) and (REF ),
\(&\int _0^T \int _{L^q(D)}\int _{D}\frac{|u(x)-u(x+h\mathbf {e}_j)|^q}{h^q} \,dx \,d\mu _t(u) \,dt \\&\qquad =\lim _{\varepsilon \rightarrow 0}\frac{1}{h^q}\int _0^T\int _{L^q(D)}\int _{D}\mathchoice{{\displaystyle {\textstyle -}{\int } }\hbox{$\textstyle -$ }}{\hspace{0.0pt}}{-}{.}56\) 







B(x+hej)|u(x)-u(y)|q  dy  dx  dt(u)  dt
      =0 1hq 0TD







B(x+hej)U2|1-2|q  dx,y2()  dy  dx  dt
      = 1hq0TDU2 |1-2|q  dx,x+hej2() dx  dt.

Since this last identity holds for any \(h>0\) , and the limit is bounded by assumption (REF ), we can pass to the limit \(h\rightarrow 0\)  and obtain, using Lemma REF ,
\(\int _0^T\int _{L^q(D)}\int _{D}|\partial _j u|^q \,dx \,d\mu _t(u) \,dt & = \lim _{h\rightarrow 0} \int _0^T \int _{L^q(D)}\int _{D}\frac{|u(x)-u(x+h\mathbf {e}_j)|^q}{h^q} \,dx \,d\mu _t(u) \,dt\\& = \lim _{h\rightarrow 0}\frac{1}{h^q}\int _0^T\int _{D}\int _{U^{2}} |\xi _{1}-\xi _{2}|^q \,d\nu _{x,x+h\mathbf {e}_j}^2(\xi )\,dx \,dt.\) 

Summing over \(j=1,\dots , d\) , we obtain the claim.
Lemma 2.4 (Reverse direction of Lemma REF ) 
Let \(\mathbf {\nu }=(\mathbf {\nu }_t)_{0\leqslant t\leqslant T}=((\nu ^1_t,\nu _t^2\dots ))_{0\leqslant t\leqslant T}\)  be a correlation measure as in [2]} that satisfies
\(\lim _{h\rightarrow 0}\frac{1}{h^q}\int _0^T \int _{D}\sum _{j=1}^d\int _{U^2}|\xi _{1}-\xi _2|^q \,d\nu ^2_{x,x+h\mathbf {e}_j}(\xi ) \,dx \,dt\leqslant C<\infty .\) 
Let \(D\)  be the \(d\) -dimensional torus. Then \(\mathbf {\nu }\)  corresponds to a family of probability measures \(\mu =(\mu _t)_{0\leqslant t\leqslant T}\)  on \(L^q(D)\)  that satisfy
\(\int _0^T\int _{L^q(D)} \left\Vert  \nabla u \right\Vert _{L^q(D)}^q \,d\mu _t(u)\,dt = \lim _{h\rightarrow 0}\frac{1}{h^q}\int _0^T \int _{D}\sum _{j=1}^d\int _{U^2}|\xi _{1}-\xi _2|^q \,d\nu ^2_{x,x+h\mathbf {e}_j}(\xi ) \,dx \,dt\leqslant \infty .\) 

The proof follows from the proof of Lemma REF  by going in the reverse direction and checking that all the steps are equivalences.
Lemma 2.5 
Let \(\Phi :H^1(D)\rightarrow \mathbb {R}\)  be a function of the form
\(\Phi (u)= \varphi \left((u,\nabla g_1),\dots , (u,\nabla g_k)\right),\) 
where \(\varphi :\mathbb {R}^k\rightarrow \mathbb {R}\)  is bounded and continuous with \(\varphi (0)=0\)  and \(g_i\in H^1(D)\) . Then
\(\int _{L^2(D)}\Phi (u) \,d\mu _t(u) =0,\) 
for any \(\mu \)  supported on a bounded set of \(L^2(D)\)  corresponding to a family of correlation measures \(\mathbf {\nu }=(\nu _1,\nu _2,\dots )\)  satisfying (REF ), i.e.,
\(\int _{D^{k}}\int _{U^{k}} \xi _1\otimes \dots \otimes \xi _{\ell }\otimes \alpha _{\ell +1}(\xi _{\ell +1})\otimes \dots \otimes \alpha _k(\xi _k) \,d\nu _{t,x}^k(\xi )\cdot \nabla _{x_1,\dots , x_\ell }\varphi (x) \,dx = 0,\) 
where \(\nabla _{x_1,\dots ,x_\ell }=(\nabla _{x_{1}},\dots ,\nabla _{x_\ell })^\top \) , \(1\leqslant \ell \leqslant k\in \mathbb {N}\) , for all \(\varphi \in H^1(D^k;U^{k-\ell })\) , \(\alpha _j\in C(U;U)\) , with \(\alpha _j(v)\leqslant C(1+|v|^2)\)  for all \(j=1,\dots , k\) .

Since \(\varphi \)  is bounded and continuous, we can approximate it on every compact subset \(B\subset \mathbb {R}^k\)  by polynomials \(p_n\)  such that
\(\left\Vert  \varphi -p_n \right\Vert _{C(B)}\leqslant \frac{1}{n},\) 
\(n\in \mathbb {N}\) . Since \(\varphi (0)=0\) , we may assume that \(p_n(0)=0\)  and hence the constant term of the polynomial is zero. It is therefore of the form
\(p_n(\zeta ) = \sum _{|\underline{j}|=1}^{N_n} \beta _{\underline{j}} \zeta _1^{j_1}\dots \zeta _k^{j_k},\) 
where \(N_n\in \mathbb {N}\)  large enough, \(\zeta =(\zeta _1,\dots ,\zeta _k)\) , \(\underline{j}=(j_1,\dots ,j_k)\) , \(j_i\geqslant 0\) , and \(|\underline{j}| =j_1+\dots +j_k\) .
Hence
\(\int _{L^2(D)}p_n((u,\nabla g_1),\dots ,(u,\nabla g_k)) \,d\mu _t(u)= \sum _{|\underline{j}|=1}^{N_n} \beta _{\underline{j}} \int _{L^2(D)}(u,\nabla g_1)^{j_1}\cdots (u,\nabla g_k)^{j_k} \,d\mu _t(u).\) 
Consider one of the terms in the sum:
\(& \int _{L^2(D)}(u,\nabla g_1)^{j_1}\cdots (u,\nabla g_k)^{j_k} \,d\mu _t(u) \\& = \int _{L^2(D)}\int _D u(x_{1})\cdot \nabla _{x_1} g_1(x_1) \,dx_1\cdots \int _D u(x_{j_1})\cdot \nabla _{x_{j_1}} g_1(x_{j_1}) \,dx_{j_1} \\&\qquad \cdot \int _D u(x_{j_1+1})\cdot \nabla _{x_{j_1+1}} g_2(x_{j_1+1}) \,dx_{j_1+1} \cdots \int _D u(x_{|\underline{j}|})\cdot \nabla _{x_{|\underline{j}|}} g_k(x_{|\underline{j}|}) d x_{|\underline{j}|} \,d\mu _t(u)\\& = \int _{D^{|\underline{j}|}}\int _{(\mathbb {R}^{N})^{|\underline{j}|}}\xi _1\otimes \dots \otimes \xi _{|\underline{j}|} \,d\nu ^{|\underline{j}|}_x(\xi ) \nabla _{x_{1}}g_1(x_1)\dots \nabla _{x_{|\underline{j}|}} g_k(x_{|\underline{j}|}) \,dx =0,\) 

where the last equality follows from (REF ). Hence
\(\int _{L^2(D)}p_n((u,\nabla g_1),\dots ,(u,\nabla g_k)) \,d\mu _t(u)= 0.\) 
By the approximation property of the \(p_n\) , after passing \(n\rightarrow \infty \) , we obtain the result for arbitrary continuous \(\varphi \) .

An identity by T. Drivas for second order structure functions
We prove the following lemma discovered by Theo Drivas (Lemma REF ):
Lemma 3.1 (Weak anisotropy) 
Let \(\mu _t\)  be a statistical solution of the Navier–Stokes equation. Then \(\mu \)  satisfies
\(\begin{split}&3\int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)}_{\partial B_r(0)} (\delta _{rn} u\cdot n)^2dS(n)\,dx \,d\mu _t(u)\,dt\\&=\int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)}_{B_r(0)} |\delta _\ell u(x)|^2 d\ell \,dx \,d\mu _t(u)\,dt .\end{split}\) 

We let \(\omega \in C^\infty _c(\mathbb {R}^d)\) , \(\omega \geqslant 0\) , compactly supported in \(B_1(0)\)  and \(\omega _\eta (y)\eta ^{-d}\omega (y/\eta )\)  for \(\eta >0\)  its rescaled version. For some fixed \(\ell \in \mathbb {R}^3\)  and \(x\in D\) , we let \(\varphi (t,y) = \ell \omega _\eta (x-y)\) . We use this as a test function in (REF ) for \(k=1\) . Using the Einstein summation convention we then have
\(\begin{split}&\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} u(y)\cdot \ell \omega _\eta (x-y)\theta ^{\prime }(t) \,dy \,d\mu _t(u)\,dt\\&-\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} u^i(y)u^j(y)\ell ^j \partial _{y^i}\omega _\eta (x-y)\theta (t)\,dy \,d\mu _t(u)\,dt \\&+ \varepsilon \int _0^T\int _{L^2_{\operatorname{div}}(D;U)} \int _{D} u^i(y)\ell ^i\Delta _y\omega _\eta (x-y)\theta (t) \,dy\,d\mu _t(u) \,dt \\&+ \int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} u^i(y)u^j(y) \partial _{y^i}\partial _{y^j}\psi _{\omega }(x-y)\theta (t)\,dy \,d\mu _t(u)\,dt=0.\end{split}\) 
where the last term comes from the fact that the test function is not divergence free, c.f. Remark REF . We define \(u_\eta (x)\int _{D} u(y)\omega _{\eta }(x-y) \,dy\)  and rewrite the above in terms of \(u_\eta \) ,
\(\begin{split}&\int _0^T\int _{L^2_{\operatorname{div}}(D;U)} u_\eta (x)\cdot \ell \theta ^{\prime }(t) \,d\mu _t(u)\,dt+\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}u^i_\eta (x)\partial _{x^i} u^j_\eta (x)\ell ^j \,d\mu _t(u)\theta (t)\,dt \\&+ \varepsilon \int _0^T\int _{L^2_{\operatorname{div}}(D;U)} \Delta _x u^i_\eta (x)\ell ^i\theta (t) \,d\mu _t(u) \,dt \\= &-\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} u^i(y)u^j(y) \partial _{y^i}\partial _{y^j}\psi _{\omega }(x-y)\theta (t)\,dy \,d\mu _t(u)\,dt\\&+\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}[u^i_\eta (x)\partial _{x^i} u^j_\eta (x)- \partial _{x^i}((u^iu^j)_{\eta }(x))]\ell ^j \,d\mu _t(u)\theta (t)\,dt .\end{split}\) 
We subtract this identity from the same identity at \(x+\ell \)  (and use the notation \(\delta _\ell u=u(x+\ell )-u(x)\) ):
\(\begin{split}&\int _0^T\int _{L^2_{\operatorname{div}}(D;U)} \delta _\ell u_\eta (x)\cdot \ell \theta ^{\prime }(t) \,d\mu _t(u)\,dt\\&+\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}[u^i_\eta (x+\ell )\partial _{x^i} u^j_\eta (x+\ell )-u^i_\eta (x)\partial _{x^i} u^j_\eta (x)]\ell ^j \,d\mu _t(u)\theta (t)\,dt \\&+ \varepsilon \int _0^T\int _{L^2_{\operatorname{div}}(D;U)} \Delta _x\delta _\ell u^i_\eta (x)\ell ^i\theta (t) \,d\mu _t(u) \,dt \\= &-\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} u^i(y)u^j(y) \partial _{y^i}\partial _{y^j}\psi _{\omega }(x+\ell -y)\theta (t)\,dy \,d\mu _t(u)\,dt\\&+\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} u^i(y)u^j(y) \partial _{y^i}\partial _{y^j}\psi _{\omega }(x-y)\theta (t)\,dy \,d\mu _t(u)\,dt\\&+\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}[u^i_\eta (x+\ell )\partial _{x^i} u^j_\eta (x+\ell )- \partial _{x^i}((u^iu^j)_{\eta }(x+\ell ))]\ell ^j \,d\mu _t(u)\theta (t)\,dt \\&-\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}[u^i_\eta (x)\partial _{x^i} u^j_\eta (x)- \partial _{x^i}((u^iu^j)_{\eta }(x))]\ell ^j \,d\mu _t(u)\theta (t)\,dt .\end{split}\) 
Denote the above terms as \(E_1,E_2,\dots ,E_7\) , which are all functions of \(x\) . We now integrate (REF ) over \(x\in D\)  and consider each term separately. As we will see, all of the terms in (REF ) vanish. We obtain
\(\int _D E_1\,dx = \int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} \delta _\ell u_\eta (x)\cdot \ell \,dx\theta ^{\prime }(t) \,d\mu _t(u)\,dt = 0,\) 
after changing the integration from \(x\)  to \(x+\ell \)  in one of the integrals. Next,
\(\int _D E_3\,dx = \varepsilon \int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} \Delta _x u^i_\eta (x)\ell ^i \,dx\theta (t) \,d\mu _t(u) \,dt =0,\) 
again because of periodic boundary conditions. We also get
\(\begin{split}\int _DE_4\,dx &=\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\iint _{D^2} u^i(y)u^j(y) \partial _{y^i}\partial _{y^j}\psi _{\omega }(x+\ell -y)\theta (t)\,dy\,dx \,d\mu _t(u)\,dt\\&=\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\iint _{D^2} u^i(y)u^j(y) \partial _{x^i}\partial _{x^j}\psi _{\omega }(x+\ell -y)\theta (t)\,dy\,dx \,d\mu _t(u)\,dt\\&=\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D} u^i(y)u^j(y) \int _{D}\partial _{x^i}\partial _{x^j}\psi _{\omega }(x+\ell -y)\,dx\theta (t)\,dy \,d\mu _t(u)\,dt=0,\end{split}\) 
and similarly, \(\int _DE_5\,dx = 0\) . We also have
\(\begin{split}\int _DE_6\,dx &=\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D}[u^i_\eta (x+\ell )\partial _{x^i} u^j_\eta (x+\ell )- \partial _{x^i}(u^iu^j)_{\eta }(x+\ell )]\ell ^j\,dx \,d\mu _t(u)\theta (t)\,dt \\&=\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D}\partial _{x^i}\big [\big (u^i_\eta (x+\ell ) u^j_\eta (x+\ell )\big )- (u^iu^j)_{\eta }(x+\ell )\big ]\ell ^j\,dx \,d\mu _t(u)\theta (t)\,dt \\&= 0\end{split}\) 
from the fact that \(u\)  is divergence free and the identity \(\int _D \partial _{x^i} v\,dx = 0\)  for every \(v\in C^1(D)\) , which follows from the periodic boundary conditions. In the same manner we obtain \(\int _DE_7\,dx = 0\) . Since \(E_2\)  is the only remaining term in (REF ), we obtain also
\(\int _DE_2\,dx = 0.\) 
We next rewrite \(E_2\)  as
\(\begin{split}E_2&=\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\big [u^i_\eta (x+\ell )\partial _{x^i} u^j_\eta (x+\ell )-u^i_\eta (x)\partial _{x^i} u^j_\eta (x)\big ]\ell ^j \,d\mu _t(u)\theta (t)\,dt \\& =\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\big [u^i_\eta (x)\partial _{x^i} \delta _\ell u^j_\eta (x) + \delta _\ell u^i_\eta (x)\partial _{x^i} u^j_\eta (x+\ell )\big ]\ell ^j \,d\mu _t(u)\theta (t)\,dt \\& =\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\big [u^i_\eta (x)\partial _{x^i} \delta _\ell u^j_\eta (x) + \delta _\ell u^i_\eta (x)\partial _{\ell ^i} u^j_\eta (x+\ell )\big ]\ell ^j \,d\mu _t(u)\theta (t)\,dt \\& =\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\big [u^i_\eta (x)\partial _{x^i} \delta _\ell u^j_\eta (x) + \delta _\ell u^i_\eta (x)\partial _{\ell ^i} \delta _\ell u^j_\eta (x)\big ]\ell ^j \,d\mu _t(u)\theta (t)\,dt \\& =\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\big [u^i_\eta (x)\partial _{x^i} \delta _\ell u^j_\eta (x)\ell ^j + \delta _\ell u^i_\eta (x)\partial _{\ell ^i}\big (\delta _\ell u^j_\eta (x)\ell ^j\big )-\delta _\ell u^i_\eta (x) \delta _\ell u^j_\eta (x)\partial _{\ell ^i}\ell ^j\big ] \,d\mu _t(u)\theta (t)\,dt \\& =\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\big [u^i_\eta (x)\partial _{x^i} \delta _\ell u^j_\eta (x)\ell ^j + \delta _\ell u^i_\eta (x)\partial _{\ell ^i}\big (\delta _\ell u^j_\eta (x)\ell ^j\big )-\big |\delta _\ell u_\eta (x)\big |^2\big ] \,d\mu _t(u)\theta (t)\,dt\\&= E_{2,1}+E_{2,2}+E_{2,3}.\end{split}\) 
Integrating these terms, we have
\(\begin{split}\int _DE_{2,1}\,dx &=\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D}u^i_\eta (x)\partial _{x^i} \delta _\ell u^j_\eta (x)\ell ^j \,dx\,d\mu _t(u)\theta (t) \,dt\\& = \int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D}\partial _{x^i}\big (u^i_\eta (x) \delta _\ell u^j_\eta (x)\ell ^j\big ) \,dx\,d\mu _t(u)\theta (t) \,dt=0,\end{split}\) 
where the first equation follows because \(\mu \)  is supported on divergence free functions. Again using the divergence free property we can write
\(\int _DE_{2,2} \,dx = \int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)} \partial _{\ell ^i}\big (\delta _\ell u^i_\eta (x) \delta _\ell u^j_\eta (x)\ell ^j\big ) \,d\mu _t(u)\theta (t)\,dx \,dt.\) 
Now we integrate in \(\ell \)  over a ball of radius \(r\) , take \(\theta (t)\equiv 1\) , and use Gauss' divergence theorem:
\(\begin{split}_{B_r(0)} E_{2,2}\,d\ell &=\int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)}_{B_r(0)} \partial _{\ell ^i}\big (\delta _\ell u^i_\eta (x) \delta _\ell u^j_\eta (x)\ell ^j\big ) d\ell \,d\mu _t(u)\,dx\theta (t)\,dt\\&=3\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D}_{\partial B_r(0)} \delta _{rn} u^i_\eta (x)n^i \delta _{rn} u^j_\eta (x)n^j dS(n)\,dx \,d\mu _t(u)\,dt\\&=3\int _0^T\int _{L^2_{\operatorname{div}}(D;U)}\int _{D}_{\partial B_r(0)} \big (\delta _{rn} u_\eta \cdot n\big )^2 dS(n)\,dx \,d\mu _t(u)\,dt,\end{split}\) 
the factor 3 coming from the relation \(|B_r(0)| = 3^{-1}r|\partial B_r(0)|\) . Since \(\mu \)  has bounded support in \(L^2\) , we can let \(\eta \rightarrow 0\)  in (REF ) and obtain
\(\begin{split}&3\int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)}_{\partial B_r(0)} (\delta _{rn} u\cdot n)^2dS(n)\,dx \,d\mu _t(u)\,dt\\&=\int _0^T\int _{D}\int _{L^2_{\operatorname{div}}(D;U)}_{B_r(0)} |\delta _\ell u(x)|^2 d\ell \,dx \,d\mu _t(u)\,dt,\end{split}\) 
which is (REF ).
