Emergence of skewed non-Gaussian distributions of velocity increments in isotropic turbulence

TL;DR

This paper introduces an exactly solvable stochastic model that explains the emergence of skewed, heavy-tailed velocity increment distributions in isotropic turbulence, aligning well with numerical data across scales and Reynolds numbers.

Contribution

The authors develop a novel multiscale stochastic model that captures the skewness and heavy tails of velocity increments in turbulence, extending understanding of non-Gaussian behaviors.

Findings

Model accurately fits numerical turbulence data across scales.

Single scale limit fails to reproduce background density.

Multiscale mechanism is crucial for distribution shape.

Abstract

Skewness and non-Gaussian behavior are essential features of the distribution of short-scale velocity increments in isotropic turbulent flows. Yet, although the skewness has been generally linked to time-reversal symmetry breaking and vortex stretching, the form of the asymmetric heavy tails remain elusive. Here we describe the emergence of both properties through an exactly solvable stochastic model with a scale hierarchy of energy transfer rates. From a statistical superposition of a local equilibrium distribution weighted by a background density, the increments distribution is given by a novel class of skewed heavy-tailed distributions, written as a generalization of the Meijer $G$-functions. Excellent agreement in the multiscale scenario is found with numerical data of systems with different sizes and Reynolds numbers. Remarkably, the single scale limit provides poor fits to the background density, highlighting the central role of the multiscale mechanism. Our framework can be also applied to describe the challenging emergence of skewed distributions in complex systems.

Figures

Equations

$$P(\delta v_r)={\int }_0^{\infty }P(\delta v_r|{\epsilon }_{\ell })f({\epsilon }_{\ell })d{\epsilon }_{\ell },$$

$${\epsilon }_{\ell }(x)=15\nu \left[\frac{1}{\ell }{\int }_x^{x+\ell }{\left(\frac{\partial v}{\partial x^{\prime }}\right)}^{2}d{x}^{\prime }-{\left(\frac{\delta v_{\ell }}{\ell }\right)}^2\right],$$

$$P(\delta v_r)={\int }_0^{\infty }P(\delta v_r|{ϵ}_{\ell })f({ϵ}_{\ell })d{ϵ}_{\ell }={\int }_0^{\infty }\frac{1}{\sqrt{2\pi {ϵ}_{\ell }}}\exp \left[-\frac{{\left(\delta v_r-⟨\delta v_r|{ϵ}_{\ell }⟩\right)}^2}{2{ϵ}_{\ell }}\right]f({ϵ}_{\ell })d{ϵ}_{\ell },$$

$$⟨\delta v_r|{ϵ}_{\ell }⟩=\mu \left({ϵ}_{\ell }-⟨{ϵ}_{\ell }⟩\right),$$

$$d{ϵ}_i=-{\gamma }_i\left({ϵ}_i-{ϵ}_{i-1}\right)\left(1+{\alpha }^2\frac{{ϵ}_{i-1}}{{ϵ}_i}\right)dt+{\kappa }_i\sqrt{{ϵ}_i{ϵ}_{i-1}}d{W}_i,$$

$$f({ϵ}_i|{ϵ}_{i-1})=\frac{({ϵ}_i/{ϵ}_{i-1}{)}^{p-1}}{2{ϵ}_{i-1}{\alpha }^p{K}_p(\omega )}\exp \left(-\frac{\beta {ϵ}_i}{{ϵ}_{i-1}}-\frac{\beta {\alpha }^2{ϵ}_{i-1}}{{ϵ}_i}\right),$$

$$f({ϵ}_N)={\int }_0^{\infty }...{\int }_0^{\infty }f({ϵ}_N|{ϵ}_{N-1})\prod _{i=1}^{N-1}[f({ϵ}_i|{ϵ}_{i-1})d{ϵ}_i].$$

$$f(\epsilon_{N})=\frac{1}{\epsilon_{0}\left[\alpha K_{p}(\omega)\right]^{N}}R_{0,N}^{N,0}\left(\begin{array}[]{c}{-}\\ {(\boldsymbol{p}-\boldsymbol{1},\boldsymbol{\omega}/2)}\end{array}\bigg{|}{\beta}^{N}\frac{\epsilon_{N}}{\epsilon_{0}}\right),$$

$$P_{N}(\delta v_{r})=ce^{\mu y}R_{0,N+1}^{N+1,0}\left(\begin{array}[]{c}{-}\\ {[(0,\boldsymbol{p}-\frac{\boldsymbol{1}}{\boldsymbol{2}}),[(\frac{|\mu y|}{2},\frac{\boldsymbol{\omega}}{2})]}\end{array}\bigg{|}\frac{{\beta}^{N}y^{2}}{2\epsilon_{0}}\right),$$

$$P_N(\delta v_r)\sim |y{|}^{\theta }\exp \left[-\beta N{\left(\frac{y}{{ϵ}_0|\mu |}\right)}^{1/N}\right]g(\delta v_r),$$

$$\frac{⟨{ϵ}_{\ell }^q⟩}{{ϵ}_0^q}=\prod _{j=1}^{q-1}{\left(1+\frac{j}{\beta }\right)}^N.$$

$$S_{2q}=\frac{⟨({\partial }_{x}v{)}^{2q}⟩}{⟨{\partial }_x{v}^2{⟩}^q}\sim \frac{⟨{ϵ}_{\ell }^q⟩}{{ϵ}_0^q}.$$

$$f({ϵ}_N)=\underset{0}{\overset{\infty }{\int }}...\underset{0}{\overset{\infty }{\int }}f({ϵ}_N|{ϵ}_{N-1})\prod _{i=1}^{N-1}[f({ϵ}_i|{ϵ}_{i-1})d{ϵ}_i],$$

$$f({ϵ}_i|{ϵ}_{i-1})=\frac{({ϵ}_i/{ϵ}_{i-1}{)}^{p-1}}{2{ϵ}_{i-1}{\alpha }^p{K}_p(\omega )}\exp \left(-\frac{\beta {ϵ}_i}{{ϵ}_{i-1}}-\frac{\beta {\alpha }^2{ϵ}_{i-1}}{{ϵ}_i}\right),$$

$$f({ϵ}_i|{ϵ}_{i-1})d{ϵ}_i=g_i(x_i)d{x}_i,$$

$$g_i(x_i)=\frac{x_i^{p-1}}{2{\alpha }^p{K}_p(\omega )}\exp \left(-\beta x_i-\frac{\beta {\alpha }^2}{x_i}\right).$$

$${ϵ}_N$$

$$={ϵ}_0\prod _{j=1}^N{x}_j.$$

$$\tilde{f}(s)$$

$$=\underset{0}{\overset{\infty }{\int }}...\underset{0}{\overset{\infty }{\int }}{ϵ}_N^{s-1}\prod _{i=1}^N[f({ϵ}_i|{ϵ}_{i-1})]d{ϵ}_N\dots d{ϵ}_1$$

$$={ϵ}_0^{s-1}\prod _{i=1}^N\left[\underset{0}{\overset{\infty }{\int }}{x}_i^{s-1}{g}_i(x_i)d{x}_i\right]$$

$$={ϵ}_0^{s-1}\prod _{i=1}^N{\tilde{g}}_i(s).$$

$${\tilde{g}}_i(s)={\alpha }^{s-1}\frac{K_{s+p-1}(\omega )}{K_p(\omega )}.$$

$$\tilde{f}(s)$$

$$={\left({ϵ}_0{\alpha }^N\right)}^{s-1}{\left[\frac{K_{s+p-1}(\omega )}{K_p(\omega )}\right]}^N.$$

$$f({ϵ}_N)=\frac{1}{{ϵ}_0{\left[\alpha K_p(\omega )\right]}^N}\frac{1}{2\pi i}\underset{\Gamma }{\int }{\left(\frac{{ϵ}_N}{{ϵ}_0{\alpha }^N}\right)}^{-s}{\left[K_{s+p-1}(\omega )\right]}^{N}ds.$$

$$R_{p,q}^{m,n}\left(\begin{matrix}\boldsymbol{a},\boldsymbol{A}\\ \boldsymbol{b},\boldsymbol{B}\end{matrix}\bigg{|}x\right)=\frac{1}{2\pi i}\displaystyle\int\limits_{\Gamma}x^{-s}\widetilde{R}_{p,q}^{m,n}\left(\begin{matrix}\boldsymbol{a},\boldsymbol{A}\\ \boldsymbol{b},\boldsymbol{B}\end{matrix}\bigg{|}s\right)ds,$$

$$\widetilde{R}_{p,q}^{m,n}\left(\begin{matrix}\boldsymbol{a},\boldsymbol{A}\\ \boldsymbol{b},\boldsymbol{B}\end{matrix}\bigg{|}s\right)=\frac{\displaystyle\prod_{j=1}^{m}B^{s}_{j}K_{{b_{j}}+s}(2B_{j})\displaystyle\prod_{k=1}^{n}A^{-s}_{k}K_{1-a_{k}-s}(2A_{k})}{\displaystyle\prod_{k=n+1}^{p}A^{s}_{k}K_{a_{k}+s}(2A_{k})\displaystyle\prod_{j=m+1}^{q}B^{-s}_{j}K_{1-b_{j}-s}(2B_{j})},$$

$$R_{0,N}^{N,0}\left(\begin{matrix}-\\ \boldsymbol{b},\boldsymbol{B}\end{matrix}\bigg{|}x\right)=\frac{1}{2\pi i}\displaystyle\int\limits_{\Gamma}x^{-s}\prod_{j=1}^{N}B^{s}_{j}K_{{b_{j}}+s}(2B_{j})ds,$$

$$f(\epsilon_{N})=\frac{1}{\epsilon_{0}\left[\alpha K_{p}(\omega)\right]^{N}}R_{0,N}^{N,0}\left(\begin{array}[]{c}{-}\\ {(\boldsymbol{p}-\boldsymbol{1},\boldsymbol{\omega}/2)}\end{array}\bigg{|}\beta^{N}\frac{\epsilon_{N}}{\epsilon_{0}}\right),$$

$$P(\delta v_{r}|\epsilon_{\ell})=\frac{e^{\mu\,y}}{\sqrt{\pi}}\left(\frac{2\,\mu}{y}\right)^{\frac{1}{2}}R_{1,0}^{0,1}\left(\begin{array}[]{c}{\frac{1}{2},\frac{|\mu\,y|}{2}}\\ {-}\end{array}\bigg{|}\frac{2\,\epsilon_{\ell}}{y^{2}}\right),$$

$$P_{N}(\delta v_{r})=ce^{\mu y}R_{0,N+1}^{N+1,0}\left(\begin{array}[]{c}{-}\\ {[(0,\boldsymbol{p}-\frac{\boldsymbol{1}}{\boldsymbol{2}}),[(\frac{|\mu y|}{2},\frac{\boldsymbol{\omega}}{2})]}\end{array}\bigg{|}\frac{{\beta}^{N}y^{2}}{2\epsilon_{0}}\right),$$

Full text

I Introduction

The phenomenon of turbulence is plentiful of challenging features that still remain elusive after decades of efforts frisch ; review . In particular, the negative skewness and non-Gaussian behavior of the distribution of velocity increments between close points in a homogeneous and isotropic turbulent flow have long figured among the most intriguing ones. Though the negative asymmetry can be derived from the Navier-Stokes equations and has been connected to the time-reversal symmetry breaking xuetal14 , elucidating its physical origins and determining the form of the heavy tails persist as long-standing open questions.

Indeed, understanding the statistical properties of velocity fluctuations has always been, and remains, an essential issue in turbulence. A significant step in this direction was Kolmogorov’s theory of turbulence frisch . One of its few exact results is the so-called 4/5-law: $\langle(\delta v_{r})^{3}\rangle=-\frac{4}{5}\langle\varepsilon\rangle r$, where $\delta v_{r}=v(x+r)-v(x)$ represents the longitudinal velocity increment and $\langle\varepsilon\rangle$ is the mean energy dissipation rate. For homogeneous and isotropic turbulent flows, in which $\langle\delta v_{r}\rangle=0$, Kolmogorov’s 4/5-law implies negative skewness and non-Gaussian statistics of velocity increments. Considerable effort has been also devoted to investigate the scaling properties of higher-order structure functions, $\langle(\delta v_{r})^{n}\rangle\sim r^{\zeta_{n}}$, $n>3$, for which no exact results are known frisch ; PRF2019 . Moreover, a renewed interest has arisen as well in the study of the increments distribution itself, rather than its set of moments andrews_1989 ; Kailetal92 ; castaing_PhysD90 ; naert-1998 ; chevetal03 ; chevetal05 ; chevetal06 . In particular, it has long been known that velocity increments for large separations tend to be Gaussian distributed, whereas non-Gaussian behavior is observed at short scales frisch . In this context, a more recent work pnas2014 found that short-scale non-Gaussian effects appear at Reynolds numbers much smaller than initially thought.

Here we report on a statistical approach to the distribution of short-scale velocity increments in isotropic turbulent flows that describes the emergence of both the negatively skewed asymmetry and non-Gaussian heavy tails, with very nice agreement with numerical turbulence data of systems featuring distinct sizes and Reynolds numbers. Our work is based on two central tenets of turbulence theory frisch ; review , namely the intermittency phenomenon and the concept of energy cascade, whereby energy is transferred from large to small eddies until dissipation by viscous forces at the shortest (Kolmogorov) scale.

Our intermittency model is built upon a hierarchy of multiple coupled scales of energy transfer rates pre2010 ; pre2017 ; pre2018 . The marginal distribution of short-scale velocity increments $P(\delta v_{r})$ is related to the energy transfer rate $\varepsilon_{\ell}$ at a larger scale $\ell$ [see Eq. (3) for a formal definition of $\varepsilon_{\ell}$] through a statistical superposition of the conditional distribution $P(\delta v_{r}|\varepsilon_{\ell})$, weighted by a background distribution $f(\varepsilon_{\ell})$ obtained in exact closed form from our model. By considering $P(\delta v_{r}|\varepsilon_{\ell})$ as a Gaussian with nonzero mean characterized by an asymmetry parameter $\mu$, we obtain an exact $P(\delta v_{r})$ in the form of a novel class of skewed functions with stretched exponential heavy tails. These newly defined functions constitute a generalization of the Meijer-$G$ functions and, to our knowledge, have never been considered in the literature.

The theoretical predictions emerging from this multiscale scenario are found to be in excellent agreement with turbulence data from two extensive and independent numerical simulations of the Navier-Stokes equations. Remarkably, a poor agreement is found if only a single scale is considered. Also, the origin of the stretched exponential heavy tails is shown to be directly related to the multiscale behavior, since a simple exponential decay would result if only a single scale were present. Therefore, our results highlight the crucial role of the interplay of multiple coupled scales of energy transfer rates, advancing on the multiscale modeling of turbulent systems in an alternative way to other approaches, such as multiplicative cascades frisch , shell bif03 and Lagrangian bifetal07 ; JohMen17 models. Moreover, our framework can be also applied to investigate the emergence of skewed distributions in other complex systems, such as financial markets complex1 and biological systems complex2 .

II Theoretical Background

We work under the formalism of a unified hierarchical approach to describe the statistics of fluctuations in multiscale complex systems pre2010 ; pre2017 ; pre2018 . This framework, called H-theory, is an extension to multiscale systems of the compounding andrews_1989 ; castaing_PhysD90 or superstatistics beck approaches to describe complex fluctuating phenomena. In this formalism, the probability distribution of the relevant signal—say, the velocity increments—at short scales is given by a statistical superposition of a large-scale conditional distribution weighted by the distribution of certain internal degrees of freedom related to the slowly fluctuating environment,

[TABLE]

where the variable $\varepsilon_{\ell}$ characterizes the local equilibrium at scale $r$. The large-scale conditional distribution is assumed to be known, so that the complex statistical properties of the turbulent state are entirely captured by the weighting density $f(\varepsilon_{\ell})$, which incorporates the effect of the fluctuating energy flux (intermittency). In turbulence modelling, the conditional distribution $P(\delta v_{r}|\varepsilon_{\ell})$ in Eq. (1) is often chosen to be a Gaussian with zero mean, while several different weighting distributions have been used, such as the gamma andrews_1989 , lognormal castaing_PhysD90 ; gagneetal94 ; chabaud_PRL94 ; naert-1998 ; yakhot2006 and inverse-gamma beck distributions. A distinctive feature of our formalism, however, is that the distribution $f(\varepsilon_{\ell})$ in (1) is not prescribed a priori—as in these previous works—, but rather is calculated from a hierarchical intermittency model; see below.

One important physical assumption built into Eq. (1) is the separation of time and length scales: the background variable $\varepsilon_{\ell}$ is supposed to vary more slowly (in time and space) than the signal $\delta v_{r}$ beck , thus allowing it to reach a quasi-equilibrium distribution $P(\delta v_{r}|\varepsilon_{\ell})$. In the statistical mechanics language, structures of size $\ell$ act as a ‘heat bath’ for the fast fluctuating quantity $\delta v_{r}$ pre2012 . In the turbulence context, $\varepsilon_{\ell}$ can be associated with the energy transfer rate from scale $\ell$ towards smaller scales, where $\ell\gg r$ in view of the assumed scale separation.

Following gagneetal94 , we consider the energy transfer rate $\varepsilon_{\ell}$ at scale $\ell$ as defined by

[TABLE]

where $\nu$ is the viscosity. The first term in the right-hand side of (2) is the space average of the dissipation rate over a volume of size $\ell$, which is Obukov’s proposal for estimating the rate of energy transfer frisch , whereas the term $15\nu\left(\frac{\delta v_{\ell}}{\ell}\right)^{2}$ takes into account the energy dissipation at the scale $\ell$ itself gagneetal94 . For large $\ell$ (say, in the inertial range), the second term is negligible and so $\varepsilon_{\ell}$ agrees with Obukov’s prescription for the energy transfer rate. In Ref. gagneetal94 it is argued that the energy transfer rate $\varepsilon_{\ell}$ defined in (2) can be approximated by $\varepsilon_{\ell}\approx 15\nu\epsilon_{\ell}/r^{2}$, where $\epsilon_{\ell}=\langle(\delta v_{r})^{2}\rangle-\langle\delta v_{r}\rangle^{2}$ is the variance of $\delta v_{r}$ at the scale $\ell$, meaning that the averages $\langle(\cdots)\rangle$ are performed over windows of size $\ell$. Here we shall make a similar assumption and take the variance, $\epsilon_{\ell}$, of $\delta v_{r}$ over a region of size $\ell$ as a proxy measure for the energy transfer rate $\varepsilon_{\ell}$. We note, however, that in our approach the scale $\ell$ is not initially known and must be determined from the velocity data, as explained in Section III.

Experimental and theoretical studies on homogeneous and isotropic turbulent flows indicate andrews_1989 ; castaing_PhysD90 ; naert-1998 ; chevetal05 ; gagneetal94 ; chabaud_PRL94 ; SKS that the conditional distribution $P(\delta v_{r}|\epsilon_{\ell})$ is given by a Gaussian with variance $\epsilon_{\ell}$. For the sake of simplicity, a Gaussian with zero mean is often considered in theoretical turbulence models andrews_1989 ; pre2010 ; pre2017 ; beck ; pre2012 , leading to symmetric (i.e., non-skewed) distributions $P(\delta v_{r})$.

Here we introduce a model for $P(\delta v_{r}|\epsilon_{\ell})$ that yields an asymmetric (skewed) distribution $P(\delta v_{r})$ which can be written in exact closed form in terms of certain special functions, see below. More specifically, we consider

[TABLE]

where the conditional mean of velocity increments $\langle\delta v_{r}|\epsilon_{\ell}\rangle$ is a function of $\epsilon_{\ell}$ with the constraint of null global average, i.e., $\langle\delta v_{r}\rangle=0$, as required for homogeneous and isotropic turbulence. We thus make the choice

[TABLE]

where $\mu$ is a flow-related asymmetry parameter so that $\langle\delta v_{r}\rangle=0$ is ensured for any $\mu$, with the advantage that it renders possible a closed analytical form for $P(\delta v_{r})$. We shall see below that the parameter $\mu$ controls the overall asymmetry of the resulting distribution $P(\delta v_{r})$. (We can also introduce a dimensionless parameter $b=|\mu|\sqrt{\langle\epsilon_{\ell}\rangle}$, but for our purposes here it is more convenient to work with $\mu$ itself; see below.) In Fig. 1 we show qualitatively how a weighted mixture of Gaussians with nonzero mean (lower curves) can yield an asymmetric, heavy-tailed distribution (uppermost curve).

The possibility of producing asymmetric distributions by compounding Gaussian distributions with nonzero mean as indicated in (3) has been generally discussed, e.g., in Refs. naert-1998 ; gagneetal94 ; Dubrulle2000 ; valvoetal15 , but with no specific models for the resulting distribution. In Ref. castaing_PhysD90 a particular non-Gaussian model was also proposed, although it did not lead to a closed form solution and had the drawback of producing a marginal distribution with nonzero mean. A model for non-Gaussian statistics and intermittency based on an ensemble of Gaussian fields—albeit with zero mean—has been also considered in the literature wilczek16 . To the best of our knowledge, the prescription given in Eq. (4) for the conditional mean velocity $\langle\delta v_{r}|\epsilon_{\ell}\rangle$ has not been used before. This is a crucial ingredient that allows us to obtain an analytic solution for the skewed marginal distribution $P(\delta v_{r})$.

We now turn to the calculation of the background distribution $f(\epsilon_{\ell})$ in Eq. (3). The scale $\ell$ is assigned to the $N$-th level of the turbulence hierarchy $(\epsilon_{\ell}\leftrightarrow\epsilon_{N})$, that is, $\ell=L/2^{N}$, where $L$ is the integral scale and $N$ is the number of levels in the cascade down from $L$ to $\ell$. Our hierarchical intermittency model is defined by the following set of $N$ stochastic differential equations:

[TABLE]

for $i=1,\dots,N$, where $\epsilon_{i}\geq 0$ represents the energy transfer rate from the hierarchy level $i$ to smaller scales, $\gamma_{i}>0$ is a relaxation rate, $\kappa_{i}>0$ characterizes the strength of the multiplicative noise (and hence of the intermittency) in the hierarchical level $i$, and $W_{i}$ denotes a Wiener process. The intermittency model (5) with $\alpha=0$ has been introduced in pre2017 . The generalization above (with $\alpha\neq 0$) is important to consider because the parameter $\alpha>0$ can be associated with a residual dissipation in the inertial range (see below), which is usually neglected in phenomenological cascade models.

Physically, the deterministic term in Eq. (5) represents the coupling between adjacent scales, whereas the stochastic term emerges from the complex interactions among all scales and is necessary for intermittency pre2017 . We further observe that a rescaling of variables $\epsilon_{i}\to\zeta\epsilon_{i}$ properly leaves the model dynamics unchanged, which is a required property for a multiplicative cascade model jimenez2000 in the sense that it implies $f(\epsilon_{i}|\epsilon_{i-1})d\epsilon_{i}=g(x)dx$, for $x=\epsilon_{i}/\epsilon_{i-1}$, where $f(\epsilon_{i}|\epsilon_{i-1})$ is the conditional distribution for $\epsilon_{i}$ with $\epsilon_{i-1}$ fixed and $g(x)$ is some function of $x$. Moreover, one can verify that if $\alpha=0$ then $\langle\epsilon_{i}\rangle=\epsilon_{0}$ for $t\to\infty$, whereas for $\alpha\neq 0$ it can be shown [see Appendix A, Eq. (31)] that $\langle\epsilon_{i}\rangle/\langle\epsilon_{i-1}\rangle=1-\alpha^{2}$, as $\alpha\to 0$, thus showing that the energy flux leaving the scale $i$ is actually smaller than that entering it. In this sense, it is thus expected that $\alpha$ becomes negligible for very large Reynolds number. The model above is perhaps the simplest stochastic dynamical model of intermittency that allows for an analytic solution (see below) and incorporates a small degree of dissipation in the cascade, so that it can describe intermittency even at not so high Reynolds numbers where residual dissipation might be relevant. It is interesting to notice that the nonlinear relaxation term in (5) is similar to the anomalous drift coefficient discussed in Nature_Lutz to model friction in the context of the unusual transport of cold atoms in dissipative optical lattices. (Higher-order terms could in principle be added in (5) but they should not affect our findings significantly and, besides, destroy the exact solvability of the model. Other non-exactly-solvable stochastic models of intermittency were considered, e.g., in eggers_1992 ; japanese2003 .)

We assume that the time scales within the cascade are largely separated, with faster dynamics at smaller scales, i.e., $\gamma_{N}\gg\gamma_{N-1}\gg\cdots\gg\gamma_{1}$. We consider furthermore that $\kappa_{N}\gg\kappa_{N-1}\gg\cdots\gg\kappa_{1}$, which is reasonable since one expects stronger intermittency at smaller scales, in such a way that the dimensionless ratio $\beta\equiv 2\gamma_{i}/\kappa_{i}^{2}$ remains invariant across scales. Under these assumptions, the stationary solution of the Fokker-Planck equation associated with (5) under Itô prescription is given by

[TABLE]

where $p=\beta(1-\alpha^{2})$, $\omega=2\alpha\beta$ and $K_{p}(x)$ is the modified Bessel function of second kind. We notice that the density function (6) has the form of a generalized inverse Gaussian (GIG) distribution, which has been applied to describe diverse fluctuation phenomena gig1 .

By denoting $f(\epsilon_{N})\equiv f(\epsilon_{\ell})$ in Eq. (1), we write

[TABLE]

Notably, these integrals can be performed exactly to give

[TABLE]

where ${\boldsymbol{p}}\equiv(p,...,p)$, ${\boldsymbol{\omega}}\equiv(\omega,...,\omega)$, and $R_{p,q}^{m,n}$ is a new special function defined in Appendix A. The function $R_{p,q}^{m,n}$ can be viewed as a generalization of the Meijer $G$-function $G_{p,q}^{m,n}$, in which the gamma functions $\Gamma(\nu)$ are essentially replaced by the Bessel functions $K_{\nu}(x)$ in the Mellin transform mathai .

Finally, substituting Eq. (8) into Eq. (3) and using some properties of the $R$-functions (see Appendices A and B), we obtain

[TABLE]

with $y=\delta v_{r}+\mu\langle\epsilon_{N}\rangle$ and $c=(2/\pi\epsilon_{0}\alpha^{N})^{1/2}/[K_{p}(\omega)]^{N}$. For a given $N$, the distribution above has four parameters, namely: $\alpha$, $\beta$, $\epsilon_{0}$, and $\mu$. As discussed above, the parameter $\alpha$ is physically related to a residual energy dissipation in the inertial range. On the other hand, the dimensionless constant $\beta$ together with $\epsilon_{0}$ define a typical scale $(2\epsilon_{0}/\beta^{N})^{1/2}$ for the fluctuations of the velocity increments $\delta v_{r}$, so that a larger relative noise (intermittency) strength and/or a lower relaxation rate consistently yields a broader distribution $P_{N}(\delta v_{r})$. Lastly, the parameter $\mu$ controls the asymmetry of the distribution, as already mentioned.

At this point, we emphasize that, although our model has four free parameters, they are determined in pairs—first $\epsilon_{0}$ and $\mu$, then $\alpha$ and $\beta$—in a two-step procedure involving the background distribution $f(\epsilon_{N})$, which is a more stringent constraint than a direct fit of $P_{N}(\delta v_{r})$; see Section III. Indeed, the fact that the background distribution $f(\epsilon_{N})$ is available to fit the empirical data in an unambiguous way, as seen below, actually proves to be an important feature of our method, since it is known that the distribution of velocity increments $P_{N}(\delta v_{r})$ can be almost equally well fitted by different theoretical expressions, thus making it difficult to select between competing models pre2017 .

We note that the single-scale case, i.e., $N=1$, in Eq. (9) corresponds to the generalized hyperbolic distribution, as the distribution $P_{1}(\delta v_{r})$ in this case reduces to a Gaussian variance-mean mixture where the mixing distribution is the GIG distribution; see Eqs. (1)-(4) and (6). The generalized hyperbolic distribution has found many applications, including in the analysis of turbulent velocity increments NIG . It appears, however, that the $N>1$ multiscale scenario and the corresponding $R$-distribution defined in (9) have not been considered before in the literature. We anticipate here that the multiscale behavior ($N>1$) is crucial to generate heavy tails, as the case $N=1$ yields only semi-heavy tails; see below. We also highlight that $P_{N}(\delta v_{r})$ given by Eq. (9) is negatively (positively) skewed for $\mu<0$ ($\mu>0$), whereas for $\mu=0$ a symmetric (non-skewed) distribution arises.

The large-$|\delta v_{r}|$ behavior of $P_{N}(\delta v_{r})$ evidences the presence of non-Gaussian tails. Indeed, for $N>1$ and negative asymmetry, $\mu<0$, we obtain

[TABLE]

where $\theta=p+1/(2N)-3/2$ and $g(\delta v_{r})=1$ for $\delta v_{r}\to-\infty$ and $g(\delta v_{r})=e^{-2|\mu|y}$ for $\delta v_{r}\to+\infty$. The negatively skewed marginal distribution displays an asymptotic behavior to the right $(\delta v_{r}\to+\infty)$ with exponential decay, while the left tail is heavier, in the form of a modified stretched exponential. In contrast, for $N=1$ modified exponential tails emerge on both sides: $P_{N=1}(\delta v_{r})\sim z^{p-1}e^{\mu y-\kappa z}$, where $\kappa=\sqrt{\mu^{2}+2\beta/\epsilon_{0}}$ and $z=\sqrt{y^{2}+2\alpha^{2}\beta\epsilon_{0}}$ for $\delta v_{r}\to\pm\infty$. Stretched exponentials have for long been used to fit turbulence data Kailetal92 despite the lack of a theoretical basis for this. Our model thus provides a reasonable physical framework for the emergence of such heavy tailed distributions.

III Data analysis

We now describe how to apply the above formalism to the data analysis of turbulent flows.

Consider a large dataset $\{\delta v_{r}(j)\}$ of longitudinal velocity increments, with $j=1,...,N_{v}$. As a first step, we need to determine the optimal window size $M$ over which the variance of $\delta v_{r}(j)$ is supposed to remain approximately constant. By dividing the original series into overlapping intervals of size $M$, we define pre2017 ; nature2017 an estimator of the local variance for each interval as $\epsilon(k)=\sum_{j=1}^{M}[\delta v(k-j)-\overline{\delta v}(k)]^{2}/M$, where $\overline{\delta v}(k)=\sum_{j=1}^{M}\delta v(k-j)/M$, with $k=M,...,N_{v}$. As discussed in Section II, we take the variance of $\delta v_{r}$ over a region of size $\ell=Mr$ as a proxy measure for the energy transfer rate from scale $\ell$ to smaller scales naert-1998 ; gagneetal94 . For various choices of $M$ and varying the asymmetry parameter $\mu$ for each $M$, we numerically compound the empirical distribution of the variance series $\{\epsilon(k)\}$ with the Gaussian as given in (3), and select the optimal parameters $M$ and $\mu$ for which the compounding integral (1) best fits the distribution of velocity increments computed from the original data. (See, e.g., guhr-beck-swinney ; guhr-beck-swinney2 ; guhr-beck-swinney3 for other methods to estimate the optimal window size for the variance series in the case of superposition of Gaussians with zero mean, therefore not corresponding to our context.) The knowledge of $M$ then allows to express $\epsilon_{0}$ in terms of the mean $\langle\epsilon_{N}\rangle$ of the variance series and the parameters $\alpha$ and $\beta$ (see Appendix C), thus leaving only $\alpha$ and $\beta$ to be determined.

Once $M$ is set, we estimate the number $N$ of scales in the cascade by $N=\log_{2}(L/\ell)=\log_{2}(L/Mr)$, where $L$ is the integral scale; see discussion preceding Eq. (5). After obtaining the variance series $\{\epsilon_{N}(k)\}$, we fit the empirical distribution $f(\epsilon_{N})$ to Eq. (8) to determine $\alpha$ and $\beta$. Finally, the theoretical distribution of velocity increments $P_{N}(\delta v_{r})$ is computed by inserting the parameters in (9). Therefore, we remark that the setting of parameters is completed prior to the calculation of $P_{N}(\delta v_{r})$.

Let us now apply this procedure to the analysis of isotropic turbulence data jhu08 generated by the extensive direct numerical simulation (DNS) of the Navier-Stokes equations for a system with $1024^{3}$ lattice points in a periodic cube and Taylor-based Reynolds number $\mbox{Re}_{\lambda}\approx 433$. The dataset was obtained from the Johns Hopkins University turbulence research group’s database jhu08 . The simulation spans five large eddy turnover times, from which we considered $\approx 3\times 10^{8}$ points for our statistics. To test our intermittency model and show that it applies well to turbulence data, we shall analyze here the velocity increments $\delta v_{r}$ computed at the smallest resolved scale $r$, which lies in the near dissipation range as $r\approx 2.14\eta$ jhu08 , where $\eta$ denotes the Kolmogorov scale. A more complete analysis including other scales $r$ will be left for future studies.

We begin by analyzing the conditional distribution $P(\delta v_{r}|\epsilon_{\ell})$, which requires computing first the joint distribution $P(\delta v_{r},\epsilon_{\ell})$ castaing_PhysD90 ; gagneetal94 ; naert-1998 ; homannetal11 . To this end, we adopt here the following ad-hoc prescription: for each window of size $M$ of the dataset $\{\delta v_{r}(j)\}$ we compute the corresponding variance $\epsilon_{N}(k)$ and associate it with the velocity increment $\delta v_{r}$ at the center of the respective window. The variance series $\{\epsilon_{N}(k)\}$ thus generated is then ‘binarized’ and for each bin we compute the respective histogram $P(\delta v_{r}|\epsilon_{N})$ of velocity increments. In Fig. 2 we show that the empirical conditional distributions $P(\delta v_{r}|\epsilon_{N})$ obtained for $M=19$ are indeed well described by Gaussians, thus validating the assumption (3), with the upper (lower) inset displaying the conditional mean (variance). Although the observed behavior of $\langle\delta v_{r}|\epsilon_{N}\rangle$ vs. $\epsilon_{N}$ is only approximately linear, the important point to note is that $\langle\delta v_{r}|\epsilon_{N}\rangle$ decreases from a positive value to a negative one as $\epsilon_{N}$ increases, thus implying $\mu<0$. (Models with a nonlinear mean-variance relationship could in principle be introduced but the distributions may not be given in analytical form.) A similar trend as that seen in the upper inset of Fig. 2 has been observed before, e.g., in Refs. hosokawa1994 ; gagneetal94 ; naert-1998 , although there the velocity increments $\delta v_{r}$ and the variance $\epsilon_{r}$ are computed over the same scale $r$, while in our case $\epsilon_{\ell}$ is defined over a larger scale $\ell=Mr$; see the discussion after Eq. (1).

Our analysis goes further, however, in that it shows mathematically that such ‘local behavior’ of the average velocity increment is linked to both the ‘global asymmetry’ and the non-Gaussian tails of the marginal distribution of the velocity increments; see the role of $\mu$ in Eqs. (9) and (10). Physically, the change in $\langle\delta v_{r}|\epsilon\rangle$ from positive to negative values for increasing energy dissipation rate $\epsilon$, as inferred from Fig. 2, is a clear evidence that the emergence of skewness is directly related to intermittency: in regions of small (large) $\epsilon$ the fluid particle is more likely to accelerate (decelerate) from one point to the next, resulting in a positive (negative) local average $\langle\delta v_{r}|\epsilon\rangle$, so that the long-time statistics of $\delta v_{r}$ has zero mean but negative skewness. We remark that a link between intermittency and skewness governed by a single parameter was also found in a recent stochastic model for the turbulent velocity field JFM16 , but no explicit distribution is obtained there.

We now proceed to further test the model. By applying the optimization procedure described above to select $M$ and $\mu$ we obtain $\mu=-1.82$ and $M=19$ (yielding $\langle\epsilon_{N}\rangle=1.09\times 10^{-3}$), which leads to the nice agreement in the inset of Fig. 3(b) between the numerical compounding (solid line) of the Gaussian with the empirical $f(\epsilon_{N})$, see Eq. (3), and the velocity increments distribution from the DNS data (circles). We note furthermore that the scale $\ell=Mr$ belongs to the inertial range, since $\ell\approx 40.7\eta$ nearly coincides with the Taylor scale $\lambda\approx 41.1\eta$ jhu08 , thus confirming the separation of scales anticipated in the discussion of Eq. (1). Using that the integral scale in this case jhu08 is $L=104.7\eta=224r$, we estimate the number of scales in the model hierarchy: $N=\log_{2}(L/\ell)=\log_{2}(224r/19r)\approx 4$.

Figure 3(a) and the main panel of Fig. 3(b) display, respectively, the marginal distribution $P_{N}(\delta v_{r})$ and background density $f(\epsilon_{N})$ for $N=4$. The theoretical results are shown in solid lines and the empirical data are depicted in circles, with excellent agreement observed in both $P_{N}(\delta v_{r})$ and $f(\epsilon_{N})$. The best fit parameters are $\alpha=0.17$ and $\beta=2.72$. For comparison, we also plot the best fit using $N=1$ (single scale), which clearly does not perform so well as the one with $N=4$. This result evidences that this DNS dataset cannot be properly described with only a single scale. We also confirmed that the $N=4$ case indeed produces a better fit than $N=2,3,5$. This is depicted in Fig. 4, in which we show the relative squared error (solid circles) of the fitted background distribution $f(\epsilon_{N})$ for different hierarchy levels $N$.

We stress that no curve fitting was performed in Fig. 3(a); the fit was done only in the main panel of Fig. 3(b) to obtain the parameters $\alpha$ and $\beta$ entering the density $f(\epsilon_{N})$. Once these parameters were known, we simply plotted $P_{N}(\delta v_{r})$ using Eq. (9) and superimposed it with the empirical histogram for the velocity increments $\delta v_{r}$. Thus, the nice agreement exhibited for $N=4$ in Fig. 3(a) using the parameters determined from Fig. 3(b) attests to the method’s self-consistency.

We now turn to analyze more recent turbulence data 4096 from the DNS of the Navier-Stokes equations for a larger system with 40963 points and higher $\mbox{Re}_{\lambda}\approx 600$. The dataset consists of only one snapshot in time from which we took $\approx 3\times 10^{8}$ points, with the smallest resolved scale being $r\approx 1.11\eta$ and the integral scale $L=907r$ 4096 . (Here again, the analysis of other scales lies out of the scope of this work.)

The theoretical results (solid lines) and DNS data (circles) for $P_{N}(\delta v_{r})$ (main panel) and $f(\epsilon_{N})$ (inset) are shown in Fig. 5. Here we find $\mu=-1.50$ and $M=27$ (yielding $\langle\epsilon_{N}\rangle=9.06\times 10^{-4}$), whereas $\alpha\approx 0$ and $\beta=2.55$. From the data in 4096 we obtain a larger number of scales $N=\log_{2}(907r/27r)\approx 5$. Indeed, for $N=5$ a remarkable agreement with the empirical data is observed for both $P_{N}(\delta v_{r})$ and $f(\epsilon_{N})$, as seen in Fig. 5(a) and the the main panel of Fig. 5(b), respectively (red curves). As in the previous analysis, the fit (green curve) using only a single scale $(N=1)$ is not as accurate as that with $N=5$. Accordingly, we found that the cases $N=2,3,4,6$ also led to poorer fits when compared to $N=5$ (see Fig. 4).

Let us now briefly examine the behavior of the model parameters with Reynolds number. First, note that the larger $N$ obtained for the second dataset, which has a higher $\mbox{Re}_{\lambda}$, is consistent with the fact that $L/\eta$ increases with $\mbox{Re}_{\lambda}$, and so we expect more steps in the cascade (hence a larger $N$) as $\mbox{Re}_{\lambda}$ enhances. Furthermore, the fact that $\alpha\approx 0$ for $N=5$ also agrees with the suggested interpretation that the $\alpha$-term in Eq. (5) represents a residual dissipation in the inertial range, which is expected to become negligible for very large $\mbox{Re}_{\lambda}$, as commented above. Note also that the asymmetry parameter $\mu$ is smaller in magnitude for the second dataset, as expected, since this case corresponds to higher $\mbox{Re}_{\lambda}$ and smaller $r$. (Recall that in the second dataset $r/L$ decreases by a factor of four and $r/\eta$, by a factor of 2.)

Moreover, the parameter $\beta$, which controls the shape of the background distribution $f(\epsilon_{\ell})$, was found to decrease slightly in the second dataset, implying that $f(\epsilon_{\ell})$ is broader in this case; compare the main plots of Figs. 3(b) and 5(b). This behavior is consistent with the expected “amplification of intermittency” chevetal05 as $r$ decreases. If this trend persists for higher $\mbox{Re}_{\lambda}$ and smaller $r$, then the normalized moments of the velocity derivative distribution should diverge for $\mbox{Re}_{\lambda}\to\infty$ (see below). At present, however, one cannot rule out the possibility that $\beta$ eventually becomes an increasing function of $N$ as $r$ gets smaller and $\mbox{Re}_{\lambda}$ larger, which would lead to constant normalized moments. It thus follows from this discussion that how $\beta$ varies with $N$ is crucial to determine the statistical properties of velocity increments at small scales. To see this more explicitly, we recall that for $\alpha\to 0$ the intermittency model given in Eq. (2) recovers that described in Ref. pre2017 for which the normalized moments of the background variable $\epsilon_{\ell}$ are

[TABLE]

If we assume that our hierarchical model remains valid at very small scales, it follows from (3) [in the limit that $\mu\to 0$] that

[TABLE]

One then sees from Eqs. (11) and (12) that as $N\to\infty$ two quite distinct scenarios arise: if $N/\beta\to\infty$ then the normalized derivative moments $S_{2q}$ diverge, as in Kolmogorov’s 1962 theory (K62) Kol62 , whereas if $N/\beta\to 0$ then $S_{2q}=\mbox{constant}$, as predicted by Kolmogorov’s original 1941 theory (K41) K41a ; K41b .

In this context, it is interesting to point out that it has recently been suggested PRF2019 that one should observe an approach towards the predictions of K41 (rather than K62) as $\mbox{Re}_{\lambda}$ continues to increase. Our hierarchical theory described above thus suggests an alternative way to assess this claim through a careful analysis of the behavior of $\beta$ for increasing $\mbox{Re}_{\lambda}$. If K41 is indeed to be obtained in such limit one should observe a faster growth of $\beta$ in comparison to $N$, whereas if $\beta$ continues to decrease (or eventually increases but slower than $N$) then $K62$ is favored. This interesting possibility certainly deserves further investigation.

In summary, we have seen from the preceding discussion that our model is rather versatile in that the changing behavior of the distribution of velocity increments with varying $r$ can be well accommodated in the background distribution $f(\epsilon_{N})$. Further investigation of the model at more scales, as well as of the dependence of the model parameters with the Reynolds number, will be left for future studies.

IV Conclusions

We have developed a hierarchical model to investigate the emergence of the negative skewness and non-Gaussian behavior of the distribution of short-scale velocity increments in isotropic turbulence. The fine agreement between the theoretical distributions, given in terms of a newly-defined transcendental function, and the empirical histograms from two independent numerical datasets highlights the crucial role of the multiple scales of the intermittent energy cascade.

The general character and plasticity of our formalism make it readily adaptable to investigate the emergence of skewed (non-Gaussian) statistics in other complex systems. For example, the symmetric version of our theory has been successfully applied to explain the emergence of turbulence in a photonic random laser, as recently reported in nature2017 ; and so we expect that the asymmetric model introduced here should also have great applicability. Indeed, the compounding of an asymmetric conditional Gaussian distribution with a background density built from a hierarchical stochastic model might be a common feature in contexts as diverse as financial markets complex1 and biological systems complex2 .

In conclusion, we emphasize that the formalism presented here not only advances on the statistical description of turbulent phenomena but can also be applied to investigate the origin of skewed non-Gaussian distributions in other complex systems.

Acknowledgements.

This work was supported in part by the following Brazilian agencies: Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), under Grants No. 303772/2017-4, No. 305062/2017-4 and No. 311497/2015-2, Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), and Fundação de Amparo à Ciência e Tecnologia do Estado de Pernambuco (FACEPE), under Grants No. APQ-0073-1.05/15 and No. APQ-0602-1.05/14.

Appendix A Derivation of the background distribution

We start by providing additional details on the exact calculation of the background probability density, $f(\epsilon_{N})\equiv f(\epsilon_{\ell})$, which incorporates the crucial effect of the fluctuating energy flux (intermittency) on the turbulence properties. As indicated in Eq. (1), in this scenario the marginal distribution of short-scale velocity increments, $P_{N}(\delta v_{r})\equiv P(\delta v_{r})$, is obtained by compounding $f(\epsilon_{\ell})$ with the Gaussian conditional distribution $P(\delta v_{r}|\epsilon_{\ell})$.

Our starting point is the multiple integral representation of the background density, Eq. (7),

[TABLE]

in which the generalized inverse Gaussian (GIG) distribution,

[TABLE]

arises as the solution of the system of stochastic differential equations, Eq. (5), with $\kappa_{i}=\sqrt{{2\gamma_{i}}/{\beta}}$, $p=\beta(1-\alpha^{2})$, $\omega=2\alpha\beta$, and $K_{\nu}(x)$ as the modified Bessel function of second kind. Introducing the new variable $x_{i}=\epsilon_{i}/\epsilon_{i-1}$, we write

[TABLE]

where

[TABLE]

We proceed by observing that

[TABLE]

Next we recall that the Mellin transform mellin of a function $f(x)$ is defined by

[TABLE]

which implies the following relation between the Mellin transforms of $f(\epsilon_{N})$ and $g_{i}(x)$:

[TABLE]

The Mellin transform of Eq. (16) is mellin

[TABLE]

Inserting Eq. (20) into Eq. (19), we see that the Mellin transform of $f(\epsilon_{N})$ is

[TABLE]

Now, using Eq. (21) and the formula of the inverse Mellin transform, we can write $f(\epsilon_{N})$ as the contour integral

[TABLE]

Further progress can be made by introducing a generalization of the Meijer-$G$ function (mathai, ), which we shall refer to as the $R$-function, in terms of the following Mellin-Barnes integral,

[TABLE]

where

[TABLE]

and $\boldsymbol{a}=(a_{1},\allowbreak\dots,\allowbreak a_{n},\allowbreak a_{n+1},\allowbreak\dots,\allowbreak a_{p})$, $\boldsymbol{A}=(A_{1},\allowbreak\dots,\allowbreak A_{n},\allowbreak A_{n+1},\allowbreak\dots,\allowbreak A_{p})$, $\boldsymbol{b}=(b_{1},\allowbreak\dots,\allowbreak b_{m},\allowbreak b_{m+1},\allowbreak\dots,\allowbreak b_{q})$, $\boldsymbol{B}=(B_{1},\allowbreak\dots,\allowbreak B_{m},\allowbreak B_{m+1},\allowbreak\dots,\allowbreak B_{q})$. The contour path $\Gamma$ is chosen so that the conditions for the existence of the inverse Mellin transform are satisfied mellin , since $R$ and $\widetilde{R}$ are Mellin pairs (see also below). From Eq. (24) we see that

[TABLE]

where $\boldsymbol{b}=(b_{1},\dots,b_{N})$ and $\boldsymbol{B}=(B_{1},\dots,B_{N})$.

The newly defined special function $R_{p,q}^{m,n}$ can be viewed as a generalization of the Meijer-$G$ function $G_{p,q}^{m,n}$, in which the gamma functions $\Gamma(\nu)$ are essentially replaced by $K_{\nu}(x)$ Bessel functions in the Mellin-Barnes integral above, Eqs. (24) and (25). Indeed, by using the limit form $K_{\nu}(x)\to\Gamma(\nu)2^{\nu-1}x^{-\nu}$, $\nu>0$, $x\to 0$, we observe that $R_{p,q}^{m,n}\to a\hskip 0.85355ptG_{p,q}^{m,n}$, where $a$ is a constant, when the argument of the Bessel functions tends to zero.

Finally, by comparing Eq. (22) with Eq. (25) we obtain the expression for the background density, Eq. (8):

[TABLE]

where ${\boldsymbol{p}}\equiv(p,...,p)$ and ${\boldsymbol{\omega}}\equiv(\omega,...,\omega)$.

The last step consists in compounding Eq. (26) with the Gaussian conditional distribution via Eq. (3) to obtain exactly the marginal distribution of short-scale velocity increments, $P_{N}(\delta v_{r})$, Eq. (9), which is also given in terms of an $R$-function. To see this, note that the Gaussian distribution in Eq. (3) can be written as

[TABLE]

with $y=\delta v_{r}+\mu\langle\epsilon_{\ell}\rangle$. Thus, the statistical composition of Eqs. (26) and (27) is performed using the integral involving the product of two $R$-functions (see property (36) below). We thus find

[TABLE]

with $c=(2/\pi\epsilon_{0}\alpha^{N})^{1/2}/[K_{p}(\omega)]^{N}$.

It follows from Eq. (19) that the mean of $f_{N}(\epsilon_{N})$ is obtained by setting $s=2$ in (21):

[TABLE]

which implies

[TABLE]

Now, using $K_{\nu}(z)\approx\Gamma(\nu)2^{\nu-1}z^{-\nu}$, for $z\to 0$, $\nu>0$, it then leads to

[TABLE]

Recursive application of this relation yields

[TABLE]

We lastly remark that the novel transcendent $R$-function, which emerges from our $N$-scale intermittency model, seems to have never been previously considered in the literature.

Appendix B Properties of the $R$-function

The general usefulness of the $R$-function representation arises from a number of identities that can be derived from extensions of related identities of the Meijer-$G$ function. Therefore, we give below a short list of some general properties of the $R$-function.

• •

Mellin transform.

[TABLE]

• •

Argument inversion.

[TABLE]

• •

Power absorption.

[TABLE]

• •

Integral involving the product of two R-functions.

[TABLE]

Appendix C Numerical Procedure

We now provide further details on the numerical procedure to apply our theoretical formalism to the analysis of general (i.e., either numerical or experimental) turbulence data.

The first step is to determine the optimal window size $M$ to compute the background series of variance estimators $\left\{\epsilon(k)\right\}$ built from the dataset as described in Section III. This is done simultaneously to the fitting of the asymmetry parameter $\mu$.

The general idea is to search for the optimal pair $(M,\mu)$ that yields the best agreement between the distribution computed numerically from Eq. (3), using the empirical density $f(\epsilon)$, and the empirical distribution of velocity increments. In practice, we compute the integral in Eq. (3) as a Monte Carlo sum,

[TABLE]

where $\langle\epsilon\rangle=\sum_{i}\epsilon_{i}/N_{M}$ and $N_{M}=N_{v}-M$ is the number of windows of size $M$. If this step is successful then one guarantees that a proper modeling of the background density will lead to a good theoretical description of the increments distribution, as described below.

We therefore note that the window size $M$ is not a free parameter in the usual sense, but it rather represents an internal length scale that needs to be obtained from the data. Other methods to estimate $M$ for Gaussians with zero mean have been proposed, e.g., in Refs. [30-32], but they do not apply to our case since our conditional Gaussians have nonzero mean, and so it was necessary to find both $M$ and $\mu$ simultaneously.

The next step is to compute the background distribution of the variance series for the optimal value of $M$ and proceed to the fitting of the theoretical prediction, Eq. (8). Through the Mellin transform formula (19) with $s=2$, yielding Eq. (29), we can relate the $\varepsilon_{0}$ parameter to the first statistical moment of the distribution (26), which is measured from the variance series, and the parameters $\alpha$ and $\beta$. This means that $\varepsilon_{0}$ is not a free parameter, so that the only two free parameters in (26) are $\alpha$ and $\beta$. These two parameters are then fitted using the value of $N$ estimated according to the description in Section III. (For comparison, we also analyze fits for other values of $N$; see main text.)

To perform the fit to Eq. (26), we must calculate the $R$-function. We note that for $N$ from 1 up to 6 the multiple integral (13) may be the most efficient way. As mentioned in Section II, the $N=1$ case is a generalized hyperbolic distribution. Interestingly, the case $N=2$ also allows for an exact integration, and, in fact, for every two new hierarchy levels—and hence two additional integrals in (13)—one integral can be executed exactly, reducing at least by half the number of integrals to be computed numerically.

On the other hand, it is also possible to compute numerically the complex integral (25). In this sense, a striking fact is that the aforementioned generalization of the Meijer-$G$ function through the substitution of the gamma functions $\Gamma(\nu)$ by the Bessel functions $K_{\nu}(x)$ in the Mellin-Barnes integral, Eqs. (23) and (24), greatly simplifies the structure of poles of the integrand. Regarding the index $\nu$, the Bessel function for a fixed $x>0$ has a pole only at infinity, and decays to zero for $\nu=c\pm i\infty$. Thus, any vertical contour in the complex plane satisfies the conditions of the Mellin inversion theorem and is suitable for the computation. The function grows very rapidly away from $\nu=0$, developing strong oscillations in the real and imaginary parts, which led us to choose a contour that passes through $\nu=0$ in the real line to attain fast numerical convergence. For a purely imaginary $\nu$ the function $K_{\nu}(x)$ is real for $x>0$, so that for a single $K$-function the integral in Eq. (23) is real. For a product of $K$-functions with different indexes, which happens for any $N>1$, the contour should pass as close as possible to the zeros of these indexes to provide convergence and stability.

Lastly, with all parameters in hand, we plot the model prediction for the distribution of velocity increments, which depends on another $R$-function, as given by Eq. (9), and compare with the one from the original empirical turbulence data.