Operator product expansion and multi-point correlations in turbulent energy cascades

Operator product expansion (OPE), also known as the 'fusion of correlations', is a powerful general tool for analyzing multi-point correlations of random fields. It claims, essentially, that for local quantities, whose correlations decay with distances, the correlations between values of the fields taken at a set of close points can be separated from their correlations with the values of the fields at distant points. That makes it possible to replace a product of closely located fields by a sum of composite local quantities (see e.g. [1]).

Originally, OPE was introduced in quantum field theory (QFT) [2, 3] and in (mathematically equivalent) equilibrium statistical mechanics of classical fields [4, 5]. In these contexts, OPE have proven to be a remarkably useful analytic tool, especially in conformal field theories, where it can even be promoted into the basic dynamical principle [5, 6]. In any case, OPE provides a natural way to separate local properties of the field statistics from global aspects such as the geometry of the space and the boundary conditions.

It is natural to ask if and to which degree OPE can be extended to more general cases of fluctuating fields, away from thermal equilibrium, such as the velocity field of developed turbulence. We are not aware of any first-principle derivation of the (equal-time) OPE in turbulence: while standard justifications of the equilibrium OPE are based on the Renormalization group analysis, no appropriate RG scheme in fluid turbulence is yet developed. The question of validity of OPE is important in particular in view of recent discovery of the traces of conformal invariance in the statistics of isolines in 2d inverse turbulent cascades in the whole family of fluid-mechanical models [7–9], including Euler hydrodynamics. In fact, conformal invariance of the statistics of fields in 2d turbulence was conjectured by Polyakov over 20 years ago [10]. We note that, regardless of the validity of this conjecture, even its general formulation requires the notion of OPE. More generally, in turbulence different symmetries (Galilean and scale invariance, time reversibility) are broken spontaneously i.e. at the level of correlation functions, see e.g. [11]. Those symmetries, however, may be present on a deeper level of OPE. For these reasons, we believe that establishing (or disproving) the validity of OPE in non-equilibrium statistics, in particular in turbulence, would be an important step in understanding its structure and symmetries [11–15].

There are several models of non-equilibrium field statistics for which OPE is known to hold. Most important is the Kraichnan model of passive scalar turbulence, which admits closed description in terms of equal-time correlation functions and analytic solution. The latter is fully consistent with OPE (see [12, 13, 16–18] and appendix). Another class of models is provided by many cases of weak wave turbulence. Here stationary distribution is essentially gaussian, and the associated correlation functions obey OPE trivially [19]. OPE for the velocity field exists in Burgers turbulence (see appendix). Yet another case amenable to straightforward perturbative and RG treatment is quasi-linear dynamics dominated at all scales by an external force with power-law spectrum [16, 20], where OPE can be derived as well.

In this situation, our strategy will be to assume validity of OPE for the velocity field and try to extract testable predictions. We believe that the natural setting to test this strategy is the stationary energy cascade in incompressible fluid in two or three spatial dimensions.

Mathematically, OPE may be formulated as follows. One assumes existence of a 'complete' set of local fluctuating fields $\{{{\mathcal{O}}_{a}}\}$, such that any product expands as

$$\begin{eqnarray}&&{{\mathcal{O}}_{a}}({{x}_{1}}){{\mathcal{O}}_{b}}({{x}_{2}})=\mathop{\sum }\limits_{c}C_{ab}^{c}({{x}_{12}}){{\mathcal{O}}_{c}}({{x}_{2}}).\end{eqnarray} \tag{ 1 }$$

Here $C_{ab}^{c}({{x}_{12}})$ are c-number (non-fluctuating) coefficients which depend on the spatial separation ${{x}_{12}}={{x}_{1}}-{{x}_{2}}$. Irrespective of the number of independent local degrees of freedom, the set $\{{{\mathcal{O}}_{a}}\}$ is always infinite—at the very least, it includes any field together with all its coordinate derivatives. It can be regarded as the extension of the notion of Taylor expansions to fields with fractal properties. Generally, the expansions (1) are understood as the relations between the correlation functions

$$\begin{eqnarray}&&\langle \;{{\mathcal{O}}_{a}}({{x}_{1}})\;{{\mathcal{O}}_{b}}({{x}_{2}})\ldots \rangle =\mathop{\sum }\limits_{c}C_{ab}^{c}({{x}_{12}})\;\langle \;{{\mathcal{O}}_{c}}({{x}_{2}})\ldots \rangle ,\end{eqnarray} \tag{ 2 }$$

where ⋯ stands for any product of local fields ${{\mathcal{O}}_{{{a}_{1}}}}({{y}_{1}})...{{\mathcal{O}}_{{{a}_{N}}}}({{y}_{N}})$. The sum (2) is assumed convergent, or at least asymptotic, to the correlation function in the lhs as long as the separation ${{x}_{12}}$ is sufficiently small.

Formally (i.e. ignoring the issues of the convergence and analyticity), successive applications of (1) allows one to reduce generic N-point correlation functions to one-point expectation values

$$\begin{eqnarray}&&\langle {{\mathcal{O}}_{{{a}_{1}}}}({{x}_{1}})\ldots {{\mathcal{O}}_{{{a}_{N}}}}({{x}_{N}})\rangle =\mathop{\sum }\limits_{b}C_{{{a}_{1}}\ldots {{a}_{N}}}^{b}({{x}_{1}},\ldots ,{{x}_{N}})\langle {{\mathcal{O}}_{b}}\rangle .\end{eqnarray} \tag{ 3 }$$

This relation clearly separates the 'local dynamics' of the theory, which is encoded in the coefficient functions, from global properties of the state, incorporated in the one-point expectation values.

Intuitively, OPE may be explained on the basis of some sort of RG procedure. Averaging over small-scale fluctuations with the wave numbers larger then $1/r$, with subsequent re-scaling $x\to x/r$, leads to r-dependent 'renormalization' of interactions between the remaining long-range degrees of freedom, while transforming original operators ${{\mathcal{O}}_{a}}(x)$ to certain linear combinations ${{\bar{\mathcal{O}}}_{a}}(x,r)=\mathop{\sum }\limits_{{{a}^{{\prime} }}}Z_{a}^{a^{\prime} }(r)\;{{\mathcal{O}}_{a^{\prime} }}(x)$, with c-number factors $Z(r)$ [37]. This allows one to relate correlations at different scales. Consider the lhs of (2), and apply the RG transformation with $r\sim |{{x}_{1}}-{{x}_{2}}|/\ell$, where $\ell$ is the microscopic scale. This relates it to the correlation function with 'microscopic' separation between the points 1 and 2, which obviously obeys OPE. Comparing with the RG transformed right-hand side (rhs) of (2), and assuming that the matrices $Z_{a}^{a^{\prime} }(r)$ are invertible, one concludes that (2) should hold at all scales, as long as $|{{x}_{12}}|$ is much smaller then all other separations $|{{x}_{2}}-{{y}_{i}}|$. While this argument (borrowed from similar consideration in equilibrium statistics, see e.g. [4, 21, 37], and used in [15]) does not constitute a proof of OPE in turbulence, it makes it tempting to assume OPE for the velocity field and other fluctuating fields in turbulent flows in 3d and 2d, and explore its possible form and consequences. This is the subject of the present note.

In this discussion we always have in mind stochastic fields in stationary, homogeneous, and isotropic regime. Thus, in the case of fluid turbulence that means pumped turbulence in the inertial interval. Spatial homogeneity was already taken into account in (1) when we claimed that the functions $C_{ab}^{c}({{x}_{12}})$ depend on the separations ${{x}_{12}}={{x}_{1}}-{{x}_{2}}$ only. Also, in the isotropic case the fields ${{\mathcal{O}}_{a}}$ can be classified according to representations of $O(d)$, the group of rotations of the d-dimensional space. Thus, the fields ${{\mathcal{O}}_{a}}$ may be scalars, vectors, or irreducible tensors of the higher ranks. In the cases of the vectors and tensors the coefficient functions $C_{ab}^{c}(x)$ will have tensor structures such that the whole OPE (1) is rotationally invariant.

In many physical situations (and those include fluid turbulence) the statistics of the random fields exhibits certain scaling properties. Assumption of the OPE allows one to incorporate those in a simple and unified way. We say that the OPE (1) is scale invariant if the fields ${{\mathcal{O}}_{a}}$ can be assigned scale dimensions d_a (which determine how the correlation functions depend on the large scale L), so that the OPE does not change under the space dilations accompanied by suitable rescalings of the fields, ${{\mathcal{O}}_{a}}(x)\to {{\tilde{O}}_{a}}(x)={{\lambda }^{{{d}_{a}}}}\;{{\mathcal{O}}_{a}}(\lambda x)$. It follows that the coefficient functions $C_{ab}^{c}(x)$ are homogeneous functions of the separation x,

$$\begin{eqnarray}&&C_{ab}^{c}(\lambda x)={{\lambda }^{{{d}_{c}}-{{d}_{a}}-{{d}_{b}}}}C_{ab}^{c}(x).\end{eqnarray} \tag{ 4 }$$

Thus, if the three fields ${{\mathcal{O}}_{a}},{{\mathcal{O}}_{b}},{{\mathcal{O}}_{c}}$ are all scalars, then $C_{ab}^{c}(x)$ is fixed up to a single overall constant, $C_{ab}^{c}(x)=C_{ab}^{c}\;|x{{|}^{{{d}_{c}}-{{d}_{a}}-{{d}_{b}}}}$, but when vectors and tensors are involved one may have several tensor structures. In all cases the scale invariance fixes $C_{ab}^{c}(x)$ up to a few numbers.

Physically, the scale invariant regime is limited to certain interval of scales,

$$\begin{eqnarray}&&\ell \ll x\ll L.\end{eqnarray} \tag{ 5 }$$

In equilibrium statistics, the 'UV' scale $\ell$ is usually the microscopic scale, such as the lattice size or interatomic distance, while the 'IR' scale is either the correlation length R_c, or the size of the system. In the case of turbulence, the origin of both UV and IR scales can be quite different. In the regime of the direct cascade (for instance, in incompressible 3d turbulence), the energy is injected by a large scale stirring and L is the pumping scale. The injected energy then 'propagates' through the inertial interval (5), to be dissipated at the viscous UV scale $\ell$. In 2d, the energy cascade is inverse: $\ell$ is the force scale, while L is either the box size or the scale where the energy is dissipated (say, by uniform friction). It is important to note that one cannot treat $\ell$ as the smallest scale even considering an inverse cascade, the viscous dissipative scale must be assumed the smallest scale; if one subject fluid to thermal driving (short-correlated random force and viscous friction related by the fluctuation–dissipation theorem), the resulting state is equilibrium rather than inverse turbulent cascade.

In scale-invariant theories the coefficients $C_{ab}^{c}(x)$ can be made independent on both the scales $\ell$ and L. Generally, the structure of the OPE (1) is thought to be determined by local properties of fluctuations, so L-independence is expected. In equilibrium statistical mechanics at criticality all properties of the scale invariant theory are determined by the associated fixed point of the RG flow [37], while all dependence of the correlation functions on the UV scale $\ell$ can be eliminated by (generally $\ell$-dependent) renormalization of the fields ${{\mathcal{O}}_{a}}$. We assume that these properties extend to the OPE in turbulence.

However, even if the OPE exhibit scale invariance according to (4), the correlation functions may not be scale invariant. This happens, for instance, if $\{{{\mathcal{O}}_{a}}\}$ includes scalar fields S of negative dimensions ${{d}_{S}}\lt 0$, having non-zero expectation values

$$\begin{eqnarray}&&\langle S\rangle \ \sim \ {{B}_{S}}\;{{L}^{-{{d}_{S}}}},\end{eqnarray} \tag{ 6 }$$

where the constant B_S depends not only on the field S but also on the details of the IR cutoff, that is on the form and the boundary conditions of the box containing the system, or, as is the case in the direct cascade, on the detailed statistics of the pumping force. In such cases the correlation functions are dominated by the negative dimension terms in (3), which diverge with the IR scale L. Formally, the situation is analogous to spontaneous symmetry breaking in field theory, and here we refer to this possibility as 'broken scale symmetry'. In all known cases of turbulence (energy cascades, vorticity cascade, turbulent convection etc), the velocity field itself has a negative dimension between 0 and −1, so that its derivatives, like vorticity, have positive dimensions. While velocity statistics necessarily depends on L, whether the correlation functions of the velocity derivatives are L-independent seems to be determined by the direction of the cascade: self-similarity is broken and the moments contain extra factors ${{(L/\ell )}^{\Delta }}$ with the whole spectrum of anomalous dimensions Δs in all direct cascades [12]. Particularly, the local energy flux $\epsilon (x)={{{\rm lim} }_{y\to 0,\ell \to 0}}{{\{[v(x+y)-v(x)]\cdot y\}}^{3}}/{{y}^{4}}$ is strongly fluctuating in 3d energy cascade and its pair correlation function decays as a power law in the inertial interval: $\langle \epsilon (0)\epsilon (x)\rangle \propto {{(L/x)}^{\mu }}$ with $\mu \approx 0.2$. On the contrary, there are claims that velocity derivatives and differences in the inverse cascade have self-similar probability distribution [11, 22], but the subject deserves more careful study.

Note that for a direct turbulent cascade of the passive scalar the OPE is scale-invariant and the correlation functions are not, while for an inverse cascade, the correlation functions are scale invariant, see [12, 13, 16–18].

In fluid turbulence, the main object of interest is the statistics of the velocity field ${{v}^{i}}(x)$, i.e. the collection of the correlation functions

$$\begin{eqnarray}&&\langle {{v}^{{{i}_{1}}}}({{x}_{1}})\cdots {{v}^{{{i}_{N}}}}({{x}_{N}})\rangle .\end{eqnarray} \tag{ 7 }$$

In this discussion we always have in mind incompressible flow, so that ${{\nabla }_{i}}{{v}^{i}}(x)=0$.

We are interested in the behavior of the correlation functions in the inertial interval (5), that is the domain where $\ell \ll |{{x}_{AB}}|\ll L$ for all the separations ${{x}_{AB}}={{x}_{A}}-{{x}_{B}}$ in (7). Known data (experimental and numerical) suggest that the velocity field in the inertial interval has the following properties:

(a) The correlation functions (7) with x_AB fixed admit finite limit $\ell \to 0$. Assuming this property, in what follows we discuss the correlation functions with the limit $\ell \to 0$ taken.

(b) The velocity field remains continuous in this limit, i.e. even after $\ell$ is sent to zero, the correlation functions (7) approach finite limits when two or more of the points x_A are brought together, ${{x}_{AB}}\to 0$. Note that this property implies existence of the composite field ${{u}^{ij}}(x)={{v}^{i}}(x){{v}^{j}}(x)$, whose dimension is twice the dimension d_v of the velocity field. Similarly, multiple powers of v can be defined, and their dimensions are multiples of d_v.

As is well known, the fluid dynamics is invariant with respect to Galilean transformations. In the stationary regime of forced turbulence the external forcing fixes a preferable frame, characterized by the condition

$$\begin{eqnarray}&&{{\bar{v}}^{i}}\equiv \langle {{v}^{i}}(x)\rangle =0.\end{eqnarray} \tag{ 8 }$$

The OPE is determined by the intrinsic dynamics of the turbulent flow, and therefore must be invariant with respect to the Galilean transformations. On the other hand, the correlation functions (7) can not be Galilean invariant, even in the inertial interval, where the forcing and dissipation manifest themselves only via the cascade. The symmetry is broken at the level of correlation functions, due to the non-invariance of the one-point expectation values in (3). The situation is formally similar to the spontaneous breaking of symmetry in phase transitions in statistical mechanics and in QFT. Possibility of breaking of scale symmetry in a similar manner was pointed out in the previous section.

The two-point correlation function has the following well known form

$$\begin{eqnarray}&&\langle {{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}})\rangle =\langle {{v}^{i}}{{v}^{j}}\rangle -{{D}^{ij}}({{x}_{12}})\;{{\bar{\epsilon }}_{2}},\end{eqnarray} \tag{ 9 }$$

where the first term involves the expectation value of the field ${{u}^{ij}}(x)={{v}^{i}}(x){{v}^{j}}(x)$, which is independent of ${{x}_{12}}$ but depends on the IR scale L. From now on, we consider energy cascades determined by the mean energy dissipation rate per unit mass $\bar{\epsilon }$. In this case it follows from isotropy and dimensional reasoning that

$$\begin{eqnarray}&&\langle \;{{v}^{i}}{{v}^{j}}\;\rangle =\frac{{{\delta }^{ij}}}{d}\langle \;{{v}^{2}}\;\rangle \sim {{\delta }^{ij}}\;{{(\bar{\epsilon }L)}^{2/3}},\end{eqnarray} \tag{ 10 }$$

i.e. the dimension d_v of the velocity field is exactly $-1/3$. The second term in (9) is as follows:

$$\begin{eqnarray}&&{{D}^{ij}}(x)=\left[ {{\delta }^{ij}}-\frac{{{\zeta }_{2}}}{{{\zeta }_{2}}+d-1}\;\frac{{{x}^{i}}{{x}^{j}}}{{{x}^{2}}} \right]\;|{{x}_{12}}{{|}^{{{\zeta }_{2}}}}.\end{eqnarray} \tag{ 11 }$$

The factor ${{\bar{\epsilon }}_{2}}$ in the second term in (9) generally is L-dependent, ${{\bar{\epsilon }}_{2}}\sim {{L}^{2/3-{{\zeta }_{2}}}}$, and in our treatment will be interpreted as the expectation value of a certain local field ${{\epsilon }_{2}}$, the exponent ${{\zeta }_{2}}$ being related to its dimension, ${{d}_{{{\epsilon }_{2}}}}={{\zeta }_{2}}-2/3$. Kolmogorov and Obukhov assumed that ${{\zeta }_{2}}=-2{{d}_{v}}=2/3$. Observed values are different in different d and in different regimes, but in all known cases ${{\zeta }_{2}}$ turns out to be close to this value. It can be related to the law of de-correlation of the velocities of two fluid particles or the scaling of the respective statistical Lagrangian integral of motion [12, 23].

The correlation function (9) is the special case n = 2 of the n momenta of the velocity difference

$$\begin{eqnarray}&&\langle {{(v(x)-v(0))}^{n}}\rangle \ \sim \ |x{{|}^{{{\zeta }_{n}}}}.\end{eqnarray} \tag{ 12 }$$

Here in general ${{(v(x)-v(0))}^{n}}$ stands for the tensor product of the velocity difference vectors, so that the rhs is a symmetric n-tensor built from the components of xⁱ. The exponents ${{\zeta }_{n}}$ may depend on d, and also may be different in the direct and inverse cascades in 2d, but ${{\zeta }_{3}}$ always equals 1, as follows from the flux constancy. With the tensor structure restored, this linear scaling of the third moment gives the exact correlation function

$$\begin{eqnarray}&&\langle {{v}^{i}}(x){{v}^{j}}(x){{v}^{k}}(0)\rangle =\frac{{\bar{\epsilon }}}{2}\left[ {{\delta }^{ik}}{{x}^{j}}+{{\delta }^{jk}}{{x}^{i}}-\frac{2}{d}{{\delta }^{ij}}{{x}^{k}} \right],\end{eqnarray} \tag{ 13 }$$

where $\bar{\epsilon }$ is a dimensional constant identified with the energy dissipation rate. In particular,

$$\begin{eqnarray}&&\langle {{v}^{i}}(x){{v}^{k}}(x){{v}^{k}}(0)\rangle =-\frac{1}{2}\langle {{v}^{i}}(x)|v(x)-v(0){{|}^{2}}\rangle =\bar{\epsilon }{{x}^{i}},\end{eqnarray} \tag{ 14 }$$

since $\langle {{v}^{2}}(x){{v}^{k}}(0)\rangle =0$. The (non-trivial) generalization of the flux relation for a curved space can be found in [24].

The form (14) is regarded as a manifestation of the energy cascade for two points: when the squared velocity difference between two points is larger than average, then the velocity difference is directed inwards for $\bar{\epsilon }\gt 0$ (the direct cascade) and outwards for $\bar{\epsilon }\lt 0$ (the inverse cascade). Further manifestations of the cascade behavior must exist at the level of multi-point correlation functions. In particular, one may ask how correlations of the velocities in three different points distinguish between direct and inverse cascade? Very little is known about the three point function from the first principles. Here we attempt to address this question by assuming the validity of OPE for the velocity field, and exploring its consequences.

Let us assume that the velocity field ${{v}^{i}}(x)$, its derivatives and infinitely many other fields, constitutes the set $\{{{\mathcal{O}}_{a}}\}$ closed with respect to the OPE (1). Consider the product ${{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}})$; by the assumption it admits the expansion

$$\begin{eqnarray}&&{{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}})=\mathop{\sum }\limits_{a}\;C_{a}^{ij}({{x}_{12}})\;{{\mathcal{O}}_{a}}({{x}_{2}}).\end{eqnarray} \tag{ 15 }$$

As was explained above, the fields ${{\mathcal{O}}_{a}}$ may have tensor indices; here we incorporate those into the index a. The coefficient functions $C_{a}^{ij}(x)$ must satisfy a number of conditions, some of them general properties of OPE, and some specific for the velocity field in turbulence. The most important, but also the most difficult, of the general conditions is the associativity of the algebra of OPE. With the exception of CFT, no efficient way of analysis of this condition is known, and we will not address this question here, limiting our attention to elementary properties.

Those include the symmetry with respect to the interchange ${{x}_{1}}\leftrightarrow {{x}_{2}},i\leftrightarrow j$, which the lhs of (15) obviously enjoys. In the rhs one can replace ${{\mathcal{O}}_{c}}({{x}_{2}})$ by ${{\mathcal{O}}_{c}}({{x}_{1}})$ at the price of adding terms with the derivatives of ${{\mathcal{O}}_{c}}$, according to the Taylor expansion, with the result consistent with the symmetry of the lhs. When ${{\mathcal{O}}_{\alpha }}$ is an independent field (i.e. is not a derivative of another local field), this simple condition demands

$$\begin{eqnarray}&&C_{\alpha }^{ij}(x)=C_{\alpha }^{ji}(-x),C_{\alpha k}^{ij}(x)=C_{\alpha k}^{ji}(-x)-C_{\alpha }^{ij}(x){{x}^{k}},\end{eqnarray} \tag{ 16 }$$

where $C_{\alpha k}^{ij}(x)$ is the coefficient at ${{\nabla }_{k}}{{\mathcal{O}}_{\alpha }}({{x}_{2}})$ in (15).

There are several conditions specific for velocity:

(i) Incompressibility, which demands

$$\begin{eqnarray}&&{{\nabla }_{i}}C_{a}^{ij}(x)=0.\end{eqnarray} \tag{ 17 }$$

Note that the divergence ${{\nabla }_{j}}\;C_{a}^{ij}(x)$ does not necessarily vanish for every field ${{\mathcal{O}}_{a}}$; however, it vanishes for independent (non-derivative) operators ${{\mathcal{O}}_{\alpha }}$.

(ii) Continuity of the velocity field. This condition demands that, when the limit $\ell \to 0$ is already taken, the limit ${{{\rm lim} }_{y\to 0}}{{v}^{i}}(x+y){{v}^{j}}(x)$ still exists, and defines the product field ${{u}^{ij}}(x)={{v}^{i}}(x){{v}^{j}}(x)$, which necessary has the dimension $2{{d}_{v}}=-2/3$. The field ${{u}^{ij}}({{x}_{2}})$ enters the OPE (15) as the dominating term, with the coefficient 1.

(iii) Galilean invariance: under the fixed-time Galilean transformation ${{v}^{i}}(x)\to {{v}^{i}}(x)+{{V}^{i}}$ the lhs of (15) changes as

$$\begin{eqnarray*}&&{{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}})\to {{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}})+{{V}^{i}}{{v}^{j}}({{x}_{2}})+{{V}^{j}}{{v}^{i}}({{x}_{1}}),\end{eqnarray*}$$

where we assumed Vⁱ to be infinitesimal. The invariance demands that the Galilean transformation of the rhs must generate the same terms. Similar but not identical terms are generated by the leading term ${{v}^{i}}({{x}_{2}}){{v}^{j}}({{x}_{2}})$ in the rhs of (15). The mismatch ${{V}^{j}}[{{v}^{i}}({{x}_{1}})-{{v}^{i}}({{x}_{2}})]$ can be written as the Taylor series in ${{x}_{12}}$, where the coefficients are the derivatives of ${{v}^{i}}(x)$ evaluated at the point x₂. Therefore the sum in (15) must include additional non-invariant fields which generate linear in V terms needed to fix the mismatch. Since the mismatch is a power series (PS below) in ${{x}_{12}}$, the non-invariant fields must come with the coefficients involving only integer powers of ${{x}_{12}}$. Emergence of such terms is not unexpected. Indeed, the OPE (15) must include, along with the leading term ${{u}^{ij}}({{x}_{2}})$, all the spatial derivatives of this field. However, the above mismatch can not always be closed with the help of these derivatives alone; in general one needs to introduce an additional non-derivative tensor field of the dimension $2{{d}_{v}}+1=\;+1/3$ (see below).

With these properties taken into account, the OPE (15) should have the following structure

$$\begin{eqnarray}&&{{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}})={\rm PS}+{\rm GIT}.\end{eqnarray} \tag{ 18 }$$

where PS and GIT stand for 'power series' (in ${{x}_{12}}$) and 'Galilean Invariant Terms', analyzed separately below.

Let us analyze power series part first. Its leading terms look as follows:

$$\begin{eqnarray}&&{\rm PS}={{u}^{ij}}({{x}_{2}})+\frac{x_{12}^{k}}{2}\left[ {{\nabla }_{k}}{{u}^{ij}}({{x}_{2}})+u_{k}^{ij}({{x}_{2}}) \right]+\cdots \;.\end{eqnarray} \tag{ 19 }$$

The first term in the bracket involves the derivative of the product field ${{u}^{ij}}(x)={{v}^{i}}(x){{v}^{j}}(x)$, and is symmetric in the indices $i,j$. The new field $u_{k}^{ij}(x)=-u_{k}^{ji}(x)$ in the second term is antisymmetric in $(i,j)$ and transforms as

$$\begin{eqnarray}&&u_{k}^{ij}(x)\to u_{k}^{ij}(x)+{{V}^{i}}{{\nabla }_{k}}{{v}^{j}}(x)-{{V}^{j}}{{\nabla }_{k}}{{v}^{i}}(x).\end{eqnarray} \tag{ 20 }$$

under the Galilean transformation. The equation (17) demands that

$$\begin{eqnarray}&&u_{i}^{ij}(x)=-{{\nabla }_{i}}\left( {{v}^{i}}(x){{v}^{j}}(x) \right).\end{eqnarray} \tag{ 21 }$$

This new field can be understood as follows. As the velocity field is not smooth (not continuously differentiable), generally the limit ${{{\rm lim} }_{{{x}_{12}}\to 0}}{{v}^{j}}({{x}_{2}}){{\nabla }_{k}}{{v}^{i}}({{x}_{1}})$ does not exist. That is clear, in particular, from the fact that the correlation function $\langle {{v}^{j}}({{x}_{2}}){{v}^{i}}({{x}_{1}})\rangle$ contains terms with ${{x}_{12}}$ in the power less than 1. It is possible, however, that a non-zero limit ${{x}_{12}}\to 0$ exists in the antisymmetric combination

$$\begin{eqnarray}&&{{v}^{j}}({{x}_{2}}){{\nabla }_{k}}{{v}^{i}}({{x}_{1}})-{{v}^{i}}({{x}_{2}}){{\nabla }_{k}}{{v}^{j}}({{x}_{1}}).\end{eqnarray} \tag{ 22 }$$

Our arguments here suggest that consistency of OPE with the Galilean invariance demands this to be the case, and the limit ${{x}_{12}}\to 0$ of (22) then defines the field $u_{k}^{ij}(x)$ in (19). By the transformation property (20) this field is identical to the microscopic field ${{v}^{i}}{{\nabla }_{k}}{{v}^{j}}-{{v}^{j}}{{\nabla }_{k}}{{v}^{i}}$. Generally, it is possible to show that the full Galilean invariance of the OPE (18) can be satisfied by the contributions from the fields ${{u}^{ij}}$, $u_{k}^{ij}$, and their derivatives. Moreover, it is not difficult to collect all these contributions explicitly; we do not present the full expression here as it is not used below.

The cases d = 1 and d = 2 are special. In d = 1 (Burgers equation) the field $u_{k}^{ij}$ vanishes identically, and the terms with the derivatives of ${{v}^{2}}(x)$ are sufficient to satisfy the Galilean invariance. This case is briefly described in the appendix. In d = 2 the filed $u_{ij}^{k}$ is not zero, but its components can be reduced to derivatives of ${{u}^{ij}}$ using the identity

$$\begin{eqnarray}&&u_{k}^{ij}={{\nabla }_{l}}\left[ \delta _{k}^{j}{{u}^{li}}(x)-\delta _{k}^{i}{{u}^{lj}}(x) \right],\end{eqnarray} \tag{ 23 }$$

verified component-wise.

Let us note that the separation between PS and GIT in (18) is ambiguous. Thus, the fields ${{\nabla }_{l}}u_{k}^{ij}-{{\nabla }_{k}}u_{l}^{ij}$ and ${{\nabla }_{i}}{{\nabla }_{k}}{{u}^{ik}}$ which may appear in the PS term (18) at the order $x_{12}^{2}$ are Galilean invariant, and their contributions can well be attributed to the GIT part. The coefficients in front of these and similar higher order regular terms are not fixed by Galilean symmetry.

The GIT generally involve fractional powers of $|{{x}_{12}}|$. The associated Galilean invariant fields may carry any representations of $O(d)$, but it is natural to expect that the leading contributions come from low spins. Including only scalar, vector and second-rank tensor fields, the OPE would look like this

$$\begin{eqnarray}\begin{array}{rcl} {\rm GIT} & = & {{D}^{ij}}({{x}_{12}})\;{{\epsilon }_{2}}({{x}_{2}})+C_{k}^{ij}({{x}_{12}})\;{{W}^{k}}({{x}_{2}}) \\ {} & {} & +G_{kl}^{ij}({{x}_{12}}){{U}^{kl}}({{x}_{2}})+H_{kl}^{ij}({{x}_{12}}){{A}^{kl}}({{x}_{2}})+...\;. \\ \end{array}\end{eqnarray} \tag{ 24 }$$

Here ${{\epsilon }_{2}}(x)$ denotes the scalar field of the lowest dimension ${{\zeta }_{2}}-2/3$, whose expectation value ${{\bar{\epsilon }}_{2}}=\langle {{\epsilon }_{2}}\rangle$ generates the second term in the correlation function (9). Note that ${{\bar{\epsilon }}_{2}}\sim {{\bar{\epsilon }}^{2/3}}{{L}^{{{\zeta }_{2}}-2/3}}$. Other terms explicitly written in (24) exemplify a priory admissible contributions of a vector field ${{W}^{k}}(x)$, and antisymmetric and irreducible symmetric tensors, ${{A}^{kl}}(x)=-{{A}^{lk}}(x)$ and ${{U}^{kl}}(x)={{U}^{lk}}(x),{{U}^{kk}}(x)=0$, respectively. Further terms may involve additional scalars, vectors and tensors, as well as independent fields with higher spins. Note that in homogeneous isotropic flow all derivatives and all fields with non-zero spin have zero expectation values, hence bring no contribution to the two-point function (9) (the situation changes when the isotropy is broken, as we discuss below). Contributions from additional scalar fields S of higher dimensions are possible, but assuming their dimensions d_S are positive, their expectation values (6) bring negligible contribution in the inertial interval ${{x}_{12}}\ll L$.

The coefficient functions $C_{k}^{ij}(x)$, $G_{kl}^{ij}(x)$, $H_{kl}^{ij}(x)$ are severely restricted by the conditions (16) and (17). Let us first show that the non-vanishing contributions from W^k and ${{A}^{ij}}$ are inconsistent with these conditions. Indeed, by isotropy and (16) the coefficient $C_{k}^{ij}(x)$ must have the form $\left[ {{\delta }^{ik}}{{x}^{j}}-{{\delta }^{jk}}{{x}^{i}} \right]\;C({{x}^{2}})$ with scalar $C({{x}^{2}})$, and then the incompressibility condition (17) leads to $C({{x}^{2}})=0$. Next, for the coefficient $H_{kl}^{ij}(x)$ the antisymmetry in $(k,l)$ and the condition $H_{kl}^{ij}(x)=H_{kl}^{ji}(-x)$, which follows from (16), allow for only one structure $H_{kl}^{ij}(x)\sim \left[ {{\delta }^{ik}}{{x}^{j}}{{x}^{l}}+{{\delta }^{jk}}{{x}^{i}}{{x}^{l}}-(k\leftrightarrow l) \right]\;{{({{x}^{2}})}^{b}}$, and then (17) fixes b to be $-(d+2)/2$, which is not acceptable because of singularity at coinciding points. We conclude that in the OPE (24) contributions from vector and antisymmetric tensor fields, which are not derivatives of other local fields, are forbidden. Note that the derivatives like ${{\nabla }_{l}}{{W}^{k}},{{\nabla }_{m}}{{A}^{kl}}$ may well be present. Similar analysis of the term with the symmetric tensor ${{U}^{kl}}$ allows one to fix the coefficient function $G_{kl}^{ij}(x)$ up to two constants,

$$\begin{eqnarray}\begin{array}{rcl} G_{kl}^{ij}(x) & = & B{{({{x}^{2}})}^{a-2}}\left\{ [2{{\delta }^{ij}}{{\delta }^{kl}}-{{\delta }^{jk}}{{\delta }^{il}}-{{\delta }^{ik}}{{\delta }^{jl}}]{{x}^{4}} \right. \\ {} & {} & +2a\;[(d+1+2a){{\delta }^{ij}}{{x}^{k}}{{x}^{l}}-{{\delta }^{jk}}{{x}^{i}}{{x}^{l}}-{{\delta }^{jl}}{{x}^{i}}{{x}^{k}} \\ {} & {} & -\left. {{\delta }^{ik}}{{x}^{j}}{{x}^{l}}-{{\delta }^{jl}}{{x}^{i}}{{x}^{k}}]{{x}^{2}}-4a(a-1)\;{{x}^{i}}{{x}^{j}}{{x}^{k}}{{x}^{l}} \right\}, \\ \end{array}\end{eqnarray} \tag{ 25 }$$

where a relates to the dimension of the field U, $2a={{d}_{U}}-2{{d}_{v}}$, and the numerical factor depends on the normalization of the field ${{U}^{kl}}$.

The presence and nature of one such tensor field can be inferred from what is known about anisotropic turbulence. Indeed, the term with a symmetric tensor should describe, in particular, the response of the pair correlation function to the strain of a mean flow. It is reasonable to assume and is known empirically that the correction to the spectrum is proportional to the gradient of the mean-flow velocity, which suggests that the OPE (18) must involve terms with the first derivatives of the velocity field, and allows to identify

$$\begin{eqnarray}&&{{U}^{kl}}={{\nabla }_{k}}{{v}_{l}}+{{\nabla }_{l}}{{v}_{k}}.\end{eqnarray} \tag{ 26 }$$

Then $G_{kl}^{ij}(x)$ must scale as a typical turn-over time at the scale x, that is in Kolmogorov phenomenology $a=2/3$, consistent with ${{d}_{U}}=2/3$. That scaling was suggested for the anisotropic part of the spectrum, which carries the momentum flux $\langle {{v}_{x}}{{v}_{z}}\rangle$ in [25–27] through a dimensional argument. Experimental evidence for it was found in [28] and [29], among others. Similar identification for the case of power-law force was made in [20]. With this in mind, here we assume that the tensor field of the lowest dimension contributing to (24) is the strain field (26). It is worth stressing, that the Kolmogorov scaling for anisotropic correction to the second moment within OPE follows from the assumption that it is linearly proportional to the external gradient i.e. from the identification (26), and it is independent of the scaling exponent ${{\zeta }_{2}}$ (whether Kolmogorov or not) of the isotropic part of the second moment. To avoid misunderstanding, note that all our results are valid for isotropic homogeneous turbulence only; we briefly invoked the data on turbulence anisotropy only to justify (26).

The overall factor B in (25) depends on $\bar{\epsilon }$. The character of this dependence can be inferred from the fact that the single-point value of the cross-correlation function $\langle {{v}_{i}}{{v}_{j}}\rangle$ can be related to the gradient of the mean flow $\partial \langle {{v}_{i}}\rangle /\partial {{x}_{j}}$ and the energy flux $\bar{\epsilon }$. Indeed, the energy exchange rate between turbulence and the mean flow is $\langle {{v}_{i}}{{v}_{j}}\rangle \partial \langle {{v}_{i}}\rangle /\partial {{x}_{j}}=\bar{\epsilon }$, which is positive in 3d [30] and negative in 2d [31]. If we now assume that the different-point correlation function $\langle {{v}_{i}}(0){{v}_{j}}(x)\rangle$ decreases with the distance x, we must require B and $\bar{\epsilon }$ to be of the same sign, i.e. up to a positive order-unity factor $B\simeq \bar{\epsilon }/|\bar{\epsilon }{{|}^{2/3}}$.

With all the non-vanishing contributions found above, the OPE (15) has the form

$$\begin{eqnarray}\begin{array}{rcl} {{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}}) & = & {{v}^{i}}{{v}^{j}}+{{D}^{ij}}({{x}_{12}}){{\epsilon }_{2}}+\frac{x_{12}^{k}}{2}\left[ {{\nabla }_{k}}({{v}^{i}}{{v}^{j}}) \right. \\ {} & {} & +\left. ({{v}^{i}}{{\nabla }_{k}}{{v}^{j}}-{{v}^{j}}{{\nabla }_{k}}{{v}^{i}}) \right]+G_{kl}^{ij}({{x}_{12}})({{\nabla }_{k}}{{v}^{l}}+{{\nabla }_{l}}{{v}^{k}}). \\ \end{array}\end{eqnarray} \tag{ 27 }$$

where all fields v and ${{\epsilon }_{2}}$ in the rhs are taken at the point x₂, and we omitted terms which involve ${{x}_{12}}$ in powers higher then 4/3. According to (11) and (25), the coefficients ${{D}^{ij}}({{x}_{12}})\sim x_{12}^{{{\zeta }_{2}}}$ and $G_{kl}^{ij}({{x}_{12}})\sim x_{12}^{4/3}$ involve fractional powers of ${{x}_{12}}$. One can say that while the first and the third terms in the rhs represent smooth component of the field $v(x)$, the second and the fourth terms there should be regarded as the manifestation of its fractal nature.

The OPE (27) suggests certain behavior of the correlation functions (7) in the situation when the separation between two points, say x₁ and x₂, is much smaller then all other separations. Here we limit attention to the three point function

$$\begin{eqnarray}&&\langle {{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}}){{v}^{k}}({{x}_{3}})\rangle .\end{eqnarray} \tag{ 28 }$$

At zero separation ${{x}_{12}}=0$ we have the exact form (13). Conventional derivation of this form implies that the limit ${{x}_{12}}\to 0$ in (28) is taken before the limit $\ell \to 0$. However, the OPE (15) suggests that this form is retained if the limits are taken in the opposite order, i.e. $\ell$ is sent to zero first, while all separations in (28) are kept within the inertial interval, see also [32]. The form (13) represents the contribution of the leading term ${{u}^{ij}}({{x}_{2}})$ in (27), which is the only survivor in the limit ${{x}_{12}}\to 0$. The behavior of (28) at ${{x}_{12}}$ finite but small as compared to ${{x}_{23}}$ is determined by further terms in (27). There is no contribution from the term with the scalar field ${{\epsilon }_{2}}$ because the two-point function $\langle {{\varepsilon }_{2}}({{x}_{2}}){{v}^{k}}({{x}_{3}})\rangle$ vanishes. Indeed, by isotropy and incompressibility (17) the two-point function must have the form Const $\times {{x}^{k}}/|x{{|}^{d}}$ at all scales, including the domain $|x|\lt \ell$. Then absence of singularity at x = 0 demands Const = 0. Therefore the leading corrections are determined by the term $\sim {{x}^{k}}$ and the contribution of the symmetric tensor ${{U}^{kl}}$ in (27). To find these contributions we need the two-point correlations of the corresponding composite operators at x₂ with ${{v}^{k}}({{x}_{3}})$. From the flux law (13) we have

$$\begin{eqnarray}&&\langle {{\nabla }_{l}}{{u}^{ij}}({{x}_{2}}){{v}^{k}}({{x}_{3}})\rangle =\frac{{\bar{\epsilon }}}{2}\left( {{\delta }_{ik}}{{\delta }_{jl}}+{{\delta }_{jk}}{{\delta }_{il}}-\frac{2}{d}\;{{\delta }_{kl}}{{\delta }_{ij}} \right).\end{eqnarray} \tag{ 29 }$$

The correlation function $\langle {{v}^{k}}(x)u_{l}^{ij}(0)\rangle$ is also dimensionless, and the conditions (17) and (21) fix its form up to single constant A,

$$\begin{eqnarray}\begin{array}{rcl} \langle u_{l}^{ij}(x){{v}^{k}}(0)\rangle & = & {{A}_{0}}\left( {{\delta }^{il}}{{\delta }^{jk}}-{{\delta }^{jl}}{{\delta }^{ik}} \right)+A[{{\delta }^{il}}{{\delta }^{jk}}-{{\delta }^{jl}}{{\delta }^{ik}}+{{e}^{i}}{{e}^{k}}{{\delta }^{jl}}-{{e}^{j}}{{e}^{k}}{{\delta }^{il}} \\ {} & {} & \left. -(d-1)\;({{e}^{i}}{{e}^{l}}{{\delta }^{jk}}-{{e}^{j}}{{e}^{l}}{{\delta }^{ik}}) \right], \\ \end{array}\end{eqnarray} \tag{ 30 }$$

where ${{A}_{0}}=(\bar{\epsilon }/2)(d+1-2/d)/(d-1)$, and ${{e}^{i}}={{x}^{i}}/|x|$ is the unit vector in the direction of xⁱ.

Remarkably, in 2d the square brackets in (30) vanishes identically and this correlation function is fixed completely. This relates to the identity (23), but can be verified directly: note that at d = 2 the form is antisymmetric in both ij and kl, thus it has only one independent component. Taking, say $(ij)=(12),(kl)=(12)$ one finds that the three terms in (30) combine to zero.

Finally, with the identification (26), the correlation function $\langle {{U}^{il}}(x){{v}^{k}}(0)\rangle$ is equal to

$$\begin{eqnarray}\begin{array}{ll} {} & {} & {{{\bar{\epsilon }}}_{2}}\left[ {{\nabla }^{i}}{{D}^{lk}}(x)+{{\nabla }^{l}}{{D}^{ik}}(x) \right]=\frac{{{{\bar{\epsilon }}}_{2}}{{\zeta }_{2}}|x{{|}^{{{\zeta }_{2}}-1}}}{{{\zeta }_{2}}+d-1}\left[ -2{{\delta }^{il}}{{x}^{k}} \right. \\ {} & {} & +\left. ({{\zeta }_{2}}+d-2)({{\delta }^{ik}}{{x}^{l}}+{{\delta }^{lk}}{{x}^{i}})+4(1-{{\zeta }_{2}})\frac{{{x}^{i}}{{x}^{l}}{{x}^{k}}}{{{x}^{2}}} \right]. \\ \end{array}\end{eqnarray} \tag{ 31 }$$

Here ${{D}^{lk}}$ is given by (11). To write the leading terms of the expansion of (28) at $|{{x}_{12}}|\ll |{{x}_{23}}|$ in a compact form, it is useful to split it into parts symmetric and antisymmetric in $i,j$. The contribution of ${{U}^{kl}}$ is symmetric in $i,j$, and thus drops out from the antisymmetric part. As a result, we have

$$\begin{eqnarray}\begin{array}{lll} {} & {} & \langle \left[ {{v}^{i}}({{x}_{1}}){{v}^{j}}({{x}_{2}})-{{v}^{j}}({{x}_{1}}){{v}^{i}}({{x}_{2}}) \right]{{v}^{k}}({{x}_{3}})\rangle ={{A}_{0}}\left( {{\delta }^{ik}}x_{12}^{j}-{{\delta }^{jk}}x_{12}^{i} \right) \\ {} & {} & +2A\left[ ({{\delta }^{ik}}x_{12}^{j}-{{\delta }^{jk}}x_{12}^{i})+{{n}^{i}}{{n}^{k}}x_{12}^{j}-{{n}^{j}}{{n}^{k}}x_{12}^{i} \right. \\ {} & {} & +\left. (d-1)({{\delta }^{ik}}{{n}^{j}}-{{\delta }^{jk}}{{n}^{i}})({{x}_{12}}\cdot n) \right]+..., \\ \end{array}\end{eqnarray} \tag{ 32 }$$

where ${\bf n}={{{\bf x}}_{23}}/|{{x}_{23}}|$, and dots are $o(|{{x}_{12}}{{|}^{4/3}})$. It is tempting to assume that $A=0$ at all d, but it needs to be checked in numerical simulations.

To treat the symmetric part, consider the velocity difference $\delta v_{12}^{i}={{v}^{i}}({{x}_{1}})-{{v}^{i}}({{x}_{2}})$ and note that the product $\delta v_{12}^{i}\delta v_{12}^{j}$, being evidently Galilean invariant, receives contributions from (24) only. The leading contribution comes from the term with ${{U}^{kl}}$,

$$\begin{eqnarray}&&\langle \delta v_{12}^{i}\delta v_{12}^{j}\;{{v}^{k}}({{x}_{3}})\rangle =2G_{kl}^{ij}({{x}_{12}})\langle {{U}^{ml}}({{x}_{23}}){{v}^{m}}(0)\rangle +...,\end{eqnarray} \tag{ 33 }$$

where $\langle {{U}^{ml}}({{x}_{23}}){{v}^{m}}(0)\rangle$ is given by (31), and the dots represent further terms which involve powers of ${{x}_{12}}$ higher then $4/3$. In particular

$$\begin{eqnarray}\begin{array}{ll} {} & {} & \langle \delta v_{12}^{2}{{v}^{k}}({{x}_{3}})\rangle =-2\langle {{v}^{i}}({{x}_{1}}){{v}^{i}}({{x}_{2}}){{v}^{k}}({{x}_{3}})\rangle \simeq {{\left( \frac{{{x}_{12}}}{{{x}_{13}}} \right)}^{1/3}}\bar{\epsilon }\left[ x_{13}^{k}(4{{{\rm cos} }^{2}}\theta -3)\frac{{{x}_{12}}}{{{x}_{13}}} \right. \\ {} & {} & +\left. x_{12}^{k}(3d-4){\rm cos} \theta \right]. \\ \end{array}\end{eqnarray} \tag{ 34 }$$

Here θ is the angle between the vectors ${{{\bf x}}_{12}}$ and ${{{\bf x}}_{23}}$. Note that the terms (25) did not contribute. It is interesting that if we take the longitudinal velocity difference ${{u}_{l}}=({\boldsymbol{ \delta }} {{{\bf v}}_{12}}\cdot {{{\bf x}}_{12}})/{{x}_{12}}$ then the moment has an opposite sign in 2d:

$$\begin{eqnarray}&&\frac{\langle |{{u}_{l}}{{|}^{2}}{{v}^{k}}({{x}_{3}})\rangle }{\langle \delta v_{12}^{2}{{v}^{k}}({{x}_{3}})\rangle }=\frac{3(2d-5)}{2(3{{d}^{2}}+7d-10)}.\end{eqnarray} \tag{ 35 }$$

Analytical closures predict that fusing points 1 and 2 inside the triple velocity moment gives the second moment times interaction time i.e. $x_{12}^{4/3}$ [27, 33, 34], which is supported by the numerical data [35, 36] and agrees with the leading term of (34).

Two-point manifestation of the energy cascade, $\nabla \langle \delta {{{\bf v}}_{12}}\delta v_{12}^{2}\rangle =-4\epsilon$, has positive/negative for direct/inverse cascade. The physical meaning is quite simple: u_l is negative for larger $\delta v$ in 3d, i.e. when there is more energy at a given scale the points tend to approach so that the energy goes to shorter distance. It is opposite in 2d: two fluid particles having large velocity difference tend to separate. Three-point correlation function (34) tells us more intimate information about the energy cascade. The sign of the correlation function (34) depends not only on the sign of the energy flux but on the geometry as well. Consider first the simplest situation when all three points are on a line (i.e. $\theta =0$). Apparently, the velocity at the distant point is more correlated with the velocity at the point which is closer. That means that the velocity at a distant point has the same sign as u_l. When there is a strong vortex in some region (large $|v({{x}_{1}})-v({{x}_{2}}){{|}^{2}}$), the velocity far away is predominantly directed towards it in the direct cascade and away from it in the inverse cascade. For arbitrary angle, the projection of this correlation function on the long direction is proportional to $2{{{\rm cos} }^{2}}\theta -1$ in 2d, i.e. changes sign at $\theta =\pi /4$ and has zero average over the angle due to incompressibility (generally $d{{{\rm cos} }^{2}}\theta -1$). In other words, when $\theta =\pi /2$ the velocity is opposite to that at $\theta =0$. The projection on the short direction is proportional to ${\rm cos} \theta (4{{{\rm cos} }^{2}}\theta -1)$ and turns into zero for $\theta =\pi /6,\pi /2$ (generally $(d-4/3){\rm sin} 2\theta$).

The general expansion equation (27) determines behavior of any correlation function (7) when two points are brought close together, much closer then the distances to all remaining points. This behavior was derived on the basis of two assumptions: the general assumption of the validity of equal-time OPE in fluid turbulence, and the identification of the leading tensor contribution with the strain (26). We believe that elucidating validity of these assumptions, especially the first one, would be an important step in understanding statistics of velocity field in turbulence. The expansion leads to the asymptotic forms (32)–(34) of the three-point correlation function, providing testable predictions of the above assumptions.

The work of GF was supported by the grants of the BSF, ISF and Minerva foundations and by the RScF grant 14-22-00259 (analytic theory). GF acknowledges kind hospitality of the Rutgers University and AZ acknowledges kind hospitality of the Weizmann Institute of Science, where most of this work was done. Research of AZ was supported in part by DOE grant DE-FG02-96ER40959, and by the BSF grant.

A.1. OPE in Burgers turbulence

For one-dimensional Burgers model in the limit of large Reynolds numbers, it is known that the statistics of velocity differences and derivatives is largely determined by dilute shocks. In particular, one can write the Burgers OPE for $v(0)v(x)$ regarding $v(x)$ as smooth field distorted by shocks. Contributions non-analytic in the inviscid limit are determined by a single shock $v(x)=u{\rm tanh} (ux/2\nu )$ in the main order in x. Typical size of the shock $\ell =\nu /\bar{v}$ (which is determined by viscosity ν and the mean velocity $\bar{v}$) plays the role of the UV scale in (5). In the inertial interval $\ell \ll {{x}_{12}}$ we have

$$\begin{eqnarray*}&&v(0)v(x)={{v}^{2}}(x)-2|x|v_{s}^{2}(x)+\frac{x}{2}{{\partial }_{x}}({{v}^{2}})(x)+\ldots ,\end{eqnarray*}$$

where $v_{s}^{2}(x)$ is related to 'shock operator' [38, 39], namely $v_{s}^{2}(x)$ returns the squared velocity of the shock (${{({{v}_{+}}+{{v}_{-}})}^{2}}/4$ where ${{v}_{\pm }}$ are the velocity immediately to the right and to the left of the shock) if a shock is present in the interval $(x,x+l)$, and zero otherwise. Its expectation value is $\langle v_{s}^{2}\rangle ={{\bar{v}}^{2}}/L\sim {{\bar{\epsilon }}^{2/3}}/{{L}^{1/3}}$, where L is the mean distance between the shocks, which plays the role of the upper bound for the inertial interval (5).

One can develop further terms in the above OPE. In the order x² we have the terms ${{x}^{2}}[v(x){{v}_{xx}}(x)/2+({{v}_{x}})_{s}^{2}(x)]$, where the operator $({{v}_{x}})_{s}^{2}(x)$ measures the square of the velocity gradient immediately outside the shock if the shock is present in $(x,x+l)$.

One can write OPE of another type, considering ${{v}_{j}}({{x}_{2}}){{\partial }_{t}}{{v}_{i}}({{x}_{1}})$. Namely, we can expect relation like $({{v}_{1}}\cdot {{\nabla }_{1}})v_{1}^{i}v_{2}^{j}=a\epsilon {{\delta }_{ij}}$ + time derivative + space derivative and the coefficient a can be determined exactly. There are infinity of such relations [40, 41].

A.2. Passive scalar model. 

The scalar field $\theta ({\bf r},t)$ satisfies the equation

$$\begin{eqnarray}&&{{\partial }_{t}}\theta +({\bf v}\nabla )\theta =f+\kappa \Delta \theta .\end{eqnarray} \tag{ 36 }$$

Here the random force is assumed to be Gaussian with the variance $\langle f(0,0)f({\bf r},t)\rangle =\delta (t)\chi (r)$ where $\chi (r)$ decays on the scale L. In the Kraichnan model, we also consider the velocity field incompressible white Gaussian with $\langle {{v}_{i}}(0,0){{v}_{j}}({\bf r},t)\rangle =\delta (t)[{{v}^{2}}{{\delta }_{ij}}-{{K}_{ij}}(r)]$, where ${{K}_{ij}}(r)={{r}^{2-\gamma }}[(d+1-\gamma ){{\delta }_{ij}}-(2-\gamma ){{r}_{i}}{{r}_{j}}{{r}^{-2}}]$. Here γ is the degree of non-smoothness, Kolmogorov cascade would correspond to $\gamma =2/3$.

In the model, the passive scalar correlation functions ${{F}_{n}}=\langle \theta ({{x}_{1}})\cdots \theta ({{x}_{n}})\rangle$ satisfy the equation

$$\begin{eqnarray}&&\mathcal{L}{{F}_{n}}=\mathop{\sum }\limits_{AB}\chi ({{x}_{AB}}){{F}_{n-2}},\end{eqnarray} \tag{ 37 }$$

where $\mathcal{L}$ is the n-point differential operator [12]

$$\begin{eqnarray}&&\mathcal{L}({{x}_{1}},...,{{x}_{n}})=\mathop{\sum }\limits_{A,B=1}^{n}{{K}_{ij}}({{x}_{AB}})\nabla _{A}^{i}\nabla _{B}^{j}.\end{eqnarray} \tag{ 38 }$$

Building OPE is straightforward in this case, see also [12, 13, 18, 42]. The OPE must be invariant with respect to the following two symmetries: (a) shift $\theta (x)\to \theta (x)$ + constant, and (b) renormalization symmetry $\theta (x)\to \Lambda \;\theta (x)$. By the last symmetry, generic operator is characterized by its 'degree' n related to the transformation property $\mathcal{O}(x)\to {{\Lambda }^{n}}\;\mathcal{O}(x)$, and the degrees of the operators in the rhs of the (1) must agree with the degree of the lhs. The OPE also have the scale symmetry (4).

Operators can be described explicitly in the case $\gamma =2$. In this case we just have all polynomials of θ and derivatives. The dimensions are $d=K-n$, where n is the degree and K is the total number of derivatives. Thus, in this case the spectrum of dimensions is highly degenerate. At $\gamma \lt 2$ the space of operators does not change, but the degeneracy of dimensions is lifted. The spectrum of the dimensions changes continuously with γ, therefore it is natural to assume that for fixed n the spectrum of dimensions remains bounded from below, at least in some range of γ. From the differential equations (37) certain general properties can be deduced [12]:

• (1)  
  The dimension of θ is exactly $-\gamma /2$.
• (2)  
  The zero-dimension operator can be defined as the flux $\epsilon ={{{\rm lim} }_{y\to 0}}{{(\theta (x+y)-\theta (x))}^{2}}/{{y}^{\gamma }}$, where it is assumed that the limit $\ell \equiv {{\kappa }^{1/\gamma }}$ is taken first, so that y is always kept within the inertial interval $\ell \ll |y|$. Alternatively, starting with finite κ one can define the flux operator as suitably renormalized square of the gradient of θ, i.e. $\epsilon ={{{\rm lim} }_{\kappa \to 0}}\kappa \;{{(\nabla \theta )}^{2}}$.
• (3)  
  For any operator $\mathcal{O}$ of degree n and dimension d there is at least one operator of the same degree and spin, and dimension $d+\gamma$. We call such operators 'descendants'. Thus, for each degree n there are 'primary' operators (denoted below as ${{\Theta }_{n,\alpha }}$, where α distinguishes between different primaries of the same degree), each accompanied by a string of its descendants, with the dimensions ${{d}_{n,\alpha }}+m\gamma$. For instance, is the first descendant of ${{\theta }^{2}}$.

There is the correspondence between the primary operators of degree n and zero modes ${{Z}_{n,\alpha }}({{x}_{1}},\ldots ,{{x}_{n}})$ of the operator $\mathcal{L}$. The zero modes are homogeneous functions with the degree of homogeneity ${{d}_{n,\alpha }}+n\gamma$. The n-fold product of θ's can be written as

$$\begin{eqnarray}&&\theta ({{x}_{1}})\ldots \theta ({{x}_{n}})=\mathop{\sum }\limits_{\alpha }{{Z}_{n,\alpha }}({{x}_{1}},\ldots ,{{x}_{n}}){{\Theta }_{n,\alpha }}({{x}_{n}})+\ldots ,\end{eqnarray} \tag{ 39 }$$

where dots stand for contribution of the descendants of ${{\Theta }_{n,\alpha }}$ and their derivatives.

For n = 2 there is only one scalar primary field ${{\theta }^{2}}$ of the dimension $-\gamma$, and all other scalar operators are its descendants. The mth descendant can be understood as the regularized composite fields ${{({{\nabla }^{m}}\theta )}^{2}}$ (it reduces to this product at $\gamma =2$). In addition, there are primary fields of any spin

For n = 4 we have infinitely many scalar primaries (and of course primaries with any spin), but only three have negative dimensions: the operator ${{\theta }^{4}}$ has the dimension $-2\gamma$, while the dimensions of the other two—the regularized products ${{\epsilon }^{2}}$ and $\epsilon \theta \Delta \theta -{\rm const}{{\epsilon }^{2}}$—are related to the exponents ${{\Delta }_{4,4}}$ and ${{\Delta }_{3,4}}$ of [12, 18]. For general n, the operator with the lowest negative dimension $-{{\Delta }_{n,n}}$ is shift invariant, and can be understood as the regularized ${{\epsilon }^{n}}$. It corresponds to the dominant zero mode of the n-point operator $\mathcal{L}$ [12].

Computation of the two-point OPE can be done, for instance, by using the decomposition of the operator (38) from [18], $\mathcal{L}={{\mathcal{L}}_{0}}(\rho )+{{\mathcal{L}}_{1}}(\rho ,R)$, where

$$\begin{eqnarray*}&&{{\mathcal{L}}_{0}}={{K}_{ij}}({{x}_{12}})\nabla _{1}^{i}\nabla _{2}^{j}={{\rho }^{1-d}}{{\partial }_{\rho }}{{\rho }^{d+1-\gamma }}{{\partial }_{\rho }},\end{eqnarray*}$$

and R is the short-hand for ${{x}_{2}},\ldots ,{{x}_{n}}$. Assume that $\rho =|{{x}_{12}}|$ is much smaller then all other separations in the correlation function F_n. Then equation (37), written as ${{\mathcal{L}}_{0}}F_{n}^{(0)}(\rho ,R)=\chi (\rho ){{F}_{n-2}}(R)-{{\mathcal{L}}_{1}}{{F}_{n}}(\rho ,R)$, can be solved iteratively, starting with the zeroth-order approximation $F_{n}^{(0)}(\rho ,R)={{F}_{2}}(\rho ){{F}_{n-2}}(R)$. The iterations generate series which is interpreted in terms of the OPE. To put it simply, for passive scalar in the Kraichnan model, the OPE follows from the two facts: (i) the main ρ-dependent part of the operator $\mathcal{L}$ separates and (ii) the rhs of (37) is a sum of reducible contributions.

Realization of that scheme requires specifying whether distances between points are larger or smaller than the correlation radius L of the pumping correlation function χ. For an incompressible velocity, which we consider here, the cascade of the scalar is direct i.e. takes place at $|{{x}_{AB}}|\lt L$. For such non-equilibrium (flux) state, we have

$$\begin{eqnarray}&&\theta (x)\theta (0)=[{{\theta }^{2}}(x)+{{\theta }^{2}}(0)]/2+|x{{|}^{\gamma }}\;\epsilon (0)+\ldots ,\end{eqnarray} \tag{ 40 }$$

where further terms involve the descendants of as well as operators of higher spin. Non-zero-dimensions and non-zero expectation values of the powers of mean that the scale invariance is spontaneously broken at the level of correlation functions.

One can also build an OPE for the fluxless state at $\rho \gt L$ (it tends to equilibrium as $r\to \infty$). This OPE is simpler. The main contributions into correlation functions have normal scaling. There are anomalous dimensions in subleading terms. The leading part $F_{n}^{(0)}(\rho ,R)={{F}_{2}}(\rho ){{F}_{n-2}}(R)$ gives the leading contribution due to the identity operator I, normalized by $\bar{\chi }=\int \chi (r){{r}^{d-1}}{\rm d}r$; the next iteration yields $F_{n}^{(1)}(\rho ,R)$, which gives the first subleading contribution: $\theta (x)\theta (0)=|x{{|}^{\gamma -d}}I+|x{{|}^{2\gamma -2-d}}{{x}^{i}}{{x}^{j}}{{\nabla }_{i}}{{\nabla }_{j}}{{O}_{1}}(0)+\ldots$, where $\langle {{O}_{1}}\rangle =0$ and ${{d}_{{{O}_{1}}}}=2-\gamma$, see [43] for the details.