Explicit construction of joint multipoint statistics in complex systems

Complex systems often involve random fluctuations for which self-similar properties in space and time play an important role. Fractional Brownian motions, characterized by a single scaling exponent, the Hurst exponent $H$, provide a convenient tool to construct synthetic signals that capture the statistical properties of many processes in the physical sciences and beyond. However, in certain strongly interacting systems, e.g., turbulent flows, stock market indices, or cardiac interbeats, multiscale interactions lead to significant deviations from self-similarity and may therefore require a more elaborate description. In the context of turbulence, the Kolmogorov-Oboukhov model (K62) describes anomalous scaling, albeit explicit constructions of a turbulent signal by this model are not available yet. Here, we derive an explicit formula for the joint multipoint probability density function of a multifractal field. To this end, we consider a scale mixture of fractional Ornstein-Uhlenbeck processes and introduce a fluctuating length scale in the corresponding covariance function. In deriving the complete statistical properties of the field, we are able to systematically model synthetic multifractal phenomena. We conclude by giving a brief outlook on potential applications which range from specific tailoring or stochastic interpolation of wind fields to the modeling of financial data or non-Gaussian features in geophysical or geospatial settings.


I. INTRODUCTION
The theory of fractals has provided a unifying view on processes that occur in complex systems [1]. In particular, the concept of self-similarity has numerous applications, ranging from the solutions of ordinary and partial differential equations to the description of complex processes in nature [2][3][4]. Typical examples include subsurface hydrology [5,6], turbulence measurements in the solar wind [7], random amoeboid motion [8], and heart interbeat fluctuations [9]. In this context, it is convenient to describe fluctuating time series or spatial fields as fractional Brownian motion (fBm) where the roughness of the signal is determined by a single scaling exponent, the Hurst exponent H [10][11][12]. Nonetheless, many strongly interacting systems exhibit deviations from statistical self-similarity. Perhaps one of the most emblematic examples for such multifractal properties is the phenomenon of intermittency in turbulence [13]. The latter manifests itself by the occurence of strong velocity fluctuations at small scales which follow a non-Gaussian distribution. Similar phenomena can also be encountered in other systems, e.g., stock market fluctuations [14], hydraulic conductivity measurements in geophysics [15], or urban rent price fields [16,17]. Hence, in contrast to fractional Brownian motion, which is a Gaussian process and thus fully determined by its covariance function, multifractal features are inherently non-Gaussian and, therefore, fundamentally more difficult to describe. Accordingly, numerous scholars have been focusing on random multifractal models, which include multiplicative or hierarchical cascade processes [18][19][20][21] and random multifractal walks [22] (see also recent works in Refs. [23][24][25][26]). The corresponding multifractal models are oftentimes formulated as non-Gaussian generalizations of fBm but a statistical characterization of their moments or probability distribution functions is rather difficult to assess.
Furthermore, the multifractal description has been devised in order to infer the scaling behavior of moments of certain quantities, such as velocity or temperature increments and is -with a few exceptions [27,28] -restricted to a single-scale separation or two-point analysis [13]. Nonetheless, to explore the spatial and temporal structures of the fluctuating field in more detail, it is useful to go beyond considerations of mere scaling aspects. This can be done by determining the probability of field fluctuations at several points x i [29]. Hence, the most general quantity in such a statistical description of a random field u(x) (e.g., a turbulent velocity field) is the joint n-point probability density function (PDF) where brackets denote ensemble averages. The knowledge of the joint n-point PDF is particularly useful to unravel the complex behavior resulting from the interplay between spatial scales and is needed in many applications such as time series reconstructions [30][31][32][33][34], for spherical n-point correlations in 2D string theory [35], many body problems in strongly correlated quantum field theory [36] or solid state physics [37][38][39], and statistical turbulence theory [33,[40][41][42]. Nevertheless, explicit determination of the multipoint PDF (1) is rather difficult if the field u(x) exhibits non-Gaussian properties and one typically has to resort to approximations, e.g., the assumption of a Markov process in scale which reduces Eq. (1) to a product of three-point quantities [31,32,43].
We introduce here an explicit construction of the joint multipoint PDF (1) with multifractal scaling that belongs to the class of the Kolmogorov-Oboukhov (K62) model of turbulence [44,45]. The method is based on an ensemble of fractional Ornstein-Uhlenbeck processes which have been modified by the introduction of fluctuating length scales. It can thus be considered as a multipoint generalization of the two-point statistics proposed by Kolmogorov and Oboukhov [44,45]. In the context of the present work, we notice that the K62 model belongs to the class of "superstatistics" put forth by Beck [46], Castaing [47][48][49] and Yakhot [50]. The concept of superstatistics has also been invoked to describe non-Gaussian behavior, e.g., the anomalous diffusion of nanospheres in an F-actin network [51], the spreading process of malignant cancer cells [52], or the modeling of financial data [53] (see also [54] for further reviews on superstatistics). The main novelty of our method is that it allows for a full statistical characterization (i.e., the determination of the n-point PDF (1)) of the multifractal field. Moreover, individual field realizations can be drawn from a so-called Gaussian scale mixture which is thus particularly suitable for modeling purposes, e.g., in the context of mesoscale atmospheric modeling [21] or synthetic wind field modeling [55].
The paper is organized as follows: Sec. II discusses the ensemble of fractional Ornstein-Uhlenbeck processes and derives their respective correlation functions. Furthermore, it is shown through the example of a turbulent velocity field that a lognormal distribution of the fluctuating scale parameters reproduces the Komogorov-Oboukhov scaling of turbulence. In Sec. III we introduce the multipoint statistics of the fractional Ornstein-Uhlenbeck scale mixture and explicitly consider one-, two-, and three-point statistical quantities. These quantities are also compared to empirical findings such as the PDF of velocity multipliers [56]. Further extensions of our method, such as the generalization to a threedimensional field, are outlined in Sec. IV. Sec. V briefly discusses potential applications of the proposed multipoint statistics in the fields of wind energy, geophysics, and also as a potential ansatz for a multi-point hierarchy in turbulence [40].

II. THE KOLMOGOROV-OBUKHOV MODEL OF TURBULENCE AS A MIXTURE OF FRACTIONAL ORNSTEIN-UHLENBECK PROCESSES
The methodology that is outlined in this section borrows from the superstatistical approach, which describes non-equilibrium systems in terms of a superposition of equilibrium statistics [57]. Hence, we propose a statistical description of a turbulent velocity field u(x) as an ensemble of Gaussian distributed velocity fields w ξ (x). To this end, we use a functional formulation of the velocity field statistics, which was first introduced by Hopf [58]. In the following, we restrict ourselves to the case of a one-dimensional field u(x) whereas a generalization to multiple dimensions is discussed in Sec. IV. We are thus led to consider the characteristic functional in terms of a superposition of characteristic functionals ϕ ξ [α] belonging to Gaussian sub-ensembles Here, g(ξ) denotes the PDF of a parameter ξ that is arbitrary at this stage. Due to the Gaussianity of the subensemble, the characteristic functional can be calculated explicitly (see for instance [29,59,60]) as where we assume that the sub-ensemble possesses no mean and is thus fully characterized by the correlation function C ξ (r) = w ξ (x + r)w ξ (x) of a Gaussian random field w ξ (x). Furthermore, we make use of the assumption of homogeneity which entails that the two-point correlation function is a function of the spatial separation r only.
In the next step, we determine the correlation function C ξ (r). To this end, we consider a fractional Ornstein-Uhlenbeck process w(x) which obeys the following ordinary stochastic differential equation where dB H (x) denotes the increment of fractional Gaussian noise with covariance Here, H is the Hurst exponent which lies in between 0 and 1. The explicit calculation of the correlation function C(r) = w(x + r)w(x) of the stationary fractional Ornstein-Uhlenbeck process can be found in Mardoukhi et al. [61] and is reproduced in Appendix A. In order to introduce fluctuations at multiple points of the scale mixture (3), we replace the scale r in the correlation function C(r) by a parametrized scale ρ ξ (r). Hence, we consider a family of Gaussian processes w ξ (x) with zero mean, which are characterized by the correlation function where we imply that r ≥ 0. Furthermore, we introduce Kummer's confluent hypergeometric function as K H (r) = 1 F 1 (2H + 1, 2H + 2, r).
We thus parametrize the correlation function of the fractional Ornstein-Uhlenbeck process by allowing for a more general function of the scale separation r, namely ρ ξ (r). In the following, we consider an explicit form of this generalized scale where µ denotes the intermittency coefficient, x ≺ a smallscale cut-off, and x a large-scale quantity which does not necessarily have to coincide with the integral length scale L. Hence, by this re-parametrization of the scale r to ρ ξ (r), we have introduced a family of Gaussian processes with two-point correlations that are a varying function of the parameter ξ and thus deviate from the simple inertial range scaling C(r) ∼ r 2H of the ordinary fractional Ornstein-Uhlenbeck process w(r). Moreover, the logarithm in Eq. (9) is restricted to positive values only It should be noted that for µ = 0 and A = 0, ρ ξ (x) = r, and the correlation function (7) reduces to the (stationary) correlation function of the fractional Ornstein-Uhlenbeck process [61]. Furthermore, in the limit r → ∞, the asymptotic behavior of Kummer's confluent hypergeometric function and we obtain lim r→∞ C ξ (r) = 0. In the final step, we demand that the parameter ξ follows a lognormal distribution which results in a scale mixture (3) that exhibits strong correlations between multiple points.
We will now verify that the increment statistics derived from (3) belongs to the class of the K62 model of turbulence. To this end, we calculate the velocity increment Here, we can introduce the second-order structure function which yields .
and for large L, we can approximate S 2,ξ (r) ≈ σ 2 ρ ξ (r) 2H and obtain which, for H = 1/3 and in the limit of r x ≺ , reduces to the original K62 scaling with scaling exponents ζ n = n 3 − µ 18 n(n − 3). Moreover, for r L, it is possible to transform the one-increment PDF (14) into a form similar to the PDF proposed by Yakhot [50] which can be obtained from a Mellin transform of the structure functions (16). To this end, we substitute ξ = v σρ ξ (r) H for ξ and obtain Similar formulas for the increment PDF have also been proposed by Beck [46] as well as by Castaing [47] who allowed for a more general form of Eq. (17) which does not necessarily entail scaling solutions (16).
We want to end this section by a recapitulation of the three steps that leed to the reconstruction of the characteristic functional (3) which contains the K62 increment statistics in Eqs. (16) and (17): i.) we consider a Gaussian random process in form of a fractional Ornstein-Uhlenbeck process w(r) whose monofractal behavior is characterized by the Hurst exponent H and whose twopoint correlation function is given by C(r), ii.) in this correlation function, we introduce a re-parameterization of the scale parameter r → ρ ξ (r). We thus consider a family of Gaussian processes with fluctuating two-point correlation functions (7). We propose an explicit form for the re-parametrization in Eq. (9), and iii.) we sum over all fluctuations of the parameter ξ by weighing it according to a lognormal distribution (11). Hence, by virtue of these three step process, we introduce strong correlations which culminate in the non-Gaussian statistics of field increments (17). This will be further elaborated upon in the following section, where we derive the main result of this paper, i.e., the multipoint statistics of the Kolmogorov-Oboukhov-type velocity field.

III. MULTIPOINT STATISTICS OF THE GAUSSIAN SCALE MIXTURE
As it has been shown in the previous section, the Gaussian scale mixture (3) with correlation function (7) and the parameter distribution leads to the K62-scaling of velocity increments (16). In this section, we derive the joint n-point probability density function of the ensemble as a superposition of multivariate normal distributions where the parameter ξ of the covariance matrices σ ij,ξ = C ξ (x i , x j ) in Eq. (7) is distributed according to a lognormal distribution (11). Hence, the n-point PDF (19) belongs to the class of compound probability distributions [62]. Moreover, as the parameter ξ in Eq. (9) represents a scale parameter, the multipoint model can also be apprehended as a multivariate Gaussian scale mixture which, in this particular case, leads to heavy-tailed behavior (see Eq. (17)). In the following, we want to further examine the multipoint statistics governed by Eq. (19). In the general case (i.e., for n > 1), Eq. (19) has no closed form and one has to resort to numerical treatments such as Markov Chain Monte Carlo or collapsed Gibbs sampling. In the following, we draw samples from the n-point PDF (19) by means of a multivariate Gaussian mixture model using TensorFlow [63]. A typical realization is depicted in Fig. 1 with model parameters L = l = 1.01, η = 0.0005, H = 1/3, µ = 0.227, and a spatial resolution of N = 4096. Furthermore, in order to generate the Gaussian scale mixture, we consider 200 realizations of the lognormally distributed parameter ξ. In the following sections, we compare the statistics of the simulated velocity field realizations to explicitly derived one-, two-, and three-point quantities of the Gaussian scale mixture.

A. Single-point statistical quantities
Due to the stationarity of the underlying fractional Ornstein-Uhlenbeck process, the single-point PDF f 1 (u 1 , x 1 ) becomes independent of x 1 and reduces to From Eq. (7) we obtain C ξ (0) = σ 2 L 2H 2 Γ(2H + 1) which is independent of the parameter ξ and yields a Gaussian   single-point PDF Hence, the correlation function (7) of the fractional Ornstein-Uhlenbeck process ensures that the multipoint statistics is homogeneous. Fig. 2 depicts the single-point PDF (blue) obtained from 5000 realizations similar to the one in Fig. 1. The dashed line corresponds to Eq. (22) and agrees well with the numerical simulations. Thhe current framework thus entails a Gaussian single-point velocity PDF and non-Gaussian features emerge in higher-order statistics such as the two-point quantities discussed in the next section.

B. Two-point statistical quantities
The correlation function C(r) = u(x + r)u(x) can readily be derived from Eq. (19) and yields which implies that C(0) = σ 2 L 2H Γ(2H+1) 2 and thus coincides with the variance of the individual fractional Ornstein-Uhlenbeck process. By the same token, we can calculate the structure functions of the velocity incre- which obviously agree with the moments of the increment PDF derived in Eq. (15). The structure functions of the simulated Gaussian scale mixture model for n = 2, 4, 6 are depicted in Fig. 3 and agree well with the theoretical predictions of the K62 model (dashed lines) at small scales. At large scales the structure functions saturate which is a consequence of the underlying fractional Ornstein-Uhlenbeck processes. Indeed, due to the asymptotic behavior of Kummer's confluent hypergeometric function we obtain lim r→∞ S n (r) = (n − 1)!!2 n/2 u n rms , where u 2 rms = σ 2 L 2 Γ(2H+1) 2 . Furthermore, it has to be stressed that all odd-order structure functions are strictly zero which also entails that the multipoint statistics possesses no third-order structure function. This result is inconsistent with the so-called 4/5-law for the third-order structure function which can be derived from the Navier-Stokes equation and suggests S 3 (r) = − 4 5 ε r, where r lies within the inertial range of scales [13]. Here, ε is defined as the average local energy dissipation rate, which for the present one-dimensional velocity field is defined according to The statistics of velocity gradients, however, cannot be predicted by the present framework due to the nondifferentiability of the velocity field. This is also highlighted by the fact that structure functions in Fig. 3 exhibit no dissipative range which would imply that S n (r) ∼ r n for small r. In Sec. V, we will briefly discuss how skewness and differentiability can be incorporated in the current framework.
The velocity increment PDFs f (v, r) are depicted in Fig. 4 and agree quantitatively with the theoretical predictions which have been obtained from a numerical evaluation of Eq. (14). In particular, the increment PDF exhibits pronounced tails at small scales r.

C. Three-point statistical quantities
As it has been discussed in the previous section, the Gaussian scale mixture was devised in a way to ensure that two-point statistical quantities obey the K62 phenomenology in the inertial range. Therefore, structure functions as well as velocity increment PDFs in Figs. 3 and 4 agree with phenomenological predictions. Hence, in this section, we want to investigate the behavior of three-point statistical quantities. Such quantities are quite important in the context of multiscale models of turbulence, e.g., the operator product expansion/fusion rules [64,65] or a stochastic interpretation of the turbulent energy cascade as a Markov process of velocity increments in scale [31,32]. In the following, we want to consider the PDF of the ratio of velocity increments for r < r , which was first discussed by Kolmogorov [44] and is also referred to as the PDF of velocity increment multipliers [56]. On the basis of experimental data, Chen et al. [56] argued that the multiplier PDF is reproduced by a Cauchy distribution derived for increments of fractional Brownian motion with Hurst parameter H, namely with scale and location parameter given by In the following, we want to determine the multiplier PDF of the Gaussian scale mixture and estimate intermittency corrections. It will be shown that the multiplier PDF agrees fairly well with the prediction from fBm (27) and that intermittency corrections are largely dominated by self-similar features (the ratio of two Gaussian distributed random variables follows a Cauchy distribution [66]) and thus might not be observable in turbulence experiments or numerical simulations. Similar arguments have already been put forth by Siefert and Peinke [67] who reproduced the Cauchy distribution by an ordinary Ornstein-Uhlenbeck process. The multiplier PDF (26) can be obtained from the twoincrement PDF (30) according to We thus derive the two-increment PDF from the threepoint PDF (Eq. (19) for n = 3) according to which results in a compound bivariate normal distribution where the correlation ρ(r, r ) between two increments is defined as ρ(r, r ) = S 2,ξ (r) + S 2,ξ (r ) − S 2,ξ (r − r) 2 S 2,ξ (r)S 2,ξ (r ) .
Inserting the two-increment PDF into Eq. (31) yields p(w, r, r ) = dξg(ξ) 1 π γ ξ (r, r ) Hence, the multiplier PDF is a compound Cauchy distribution with location and scale parameter defined as and w 0,ξ (r, r ) = S 2,ξ (r) respectively. We have numerically evaluated the integral in Eq. (35) for fixed r = 0.25L and for different values of small-scale r. As shown in Fig. 5(a), the multiplier PDF of the Gaussian scale mixture is in good agreement with the prediction of the fBm model (dashed lines) given by Eq. (27). At larger scale separations, however, slight differences in the tails can be perceived. These differences manifest themselves in the cumulative distribution function P (w, r, r ) in Fig. 5(b) as well. In order to further quantify the intermittency corrections, we have computed the Kullback-Leibler divergence [68] K(r, r ) = dwp(w, r, r ) log p(w, r, r ) p f Bm (w, r, r ) , (38) between the compound Cauchy distribution and the prediction by the fBm model (27) for different intermittency coefficients µ and r (with fixed large-scale r = 0.25L).
Although moments of the Cauchy distributions are not defined, it is possible to compute the Kullback-Leibler divergence of two Cauchy distributions [69], which should be sufficient to justify using it for our purpose here. As expected, the Kullback-Leibler divergence depicted in Fig. 5(c) suggests deviations for large scale separations r − r, but decreases for smaller intermittency coefficients µ. For intermediate scale separations, the fBm approximation holds regardless of the magnitude of µ, whereas deviations reappear for r ≈ r , as the Cauchy distribution approaches the Dirac delta-function. Hence, although the multiplier PDF appears to be rather insensitive to intermittency corrections, the multiplier PDF is not fully monofractal as suggested by experimental studies [56,67]. Due to the Gaussian self-similar features contained in the velocity increments, the multiplier PDF is already so heavy-tailed that intermittency corrections are barely observable. For a proper multiscale description of the turbulent energy cascade, it might thus be appropriate to consider other observables such as transition PDFs [31] or multipoint moments [64,65,70].

IV. EXTENSIONS AND GENERALIZATIONS OF THE GAUSSIAN SCALE MIXTURE MULTIPOINT STATISTICS
In the previous sections we have presented a onedimensional multipoint statistics of a Gaussian scale mixture which exhibits non-Gaussian features. In the following, we focus on multi-dimensional generalizations as well as further potential improvements. The n-point PDF (19) can readily be generalized to the multipoint statistics of a three-dimensional velocity field according to where greek indices indicate spatial components α = 1, 2, 3, and where we introduce the covariance matrix σ ξ,i,α;j,β = C ξ,αβ (x i − x j ) with the correlation function C ξ,αβ (r) = (C ξ,rr (r) − C ξ,tt (r)) r α r β r 2 + C ξ,tt (r)δ αβ . (40) Here, we assume that either the longitudinal C ξ,rr (r) or the transverse correlation function C ξ,tt (r) obey Eq. (7). A relation between both functions can then be established by the incompressibility condition ∇ x · u(x) = 0 for the velocity field. This leads to the following relation: which can be considered as a generalization of the von-Kármán-Howarth relation [71]. In the context of turbulence, it remains unclear whether longitudinal and transverse statistical quantities exhibit different scaling properties [72][73][74][75][76][77][78]. Establishing connections between longitudinal and transverse statistics on the basis of the integral constraint (41) will be a task for future work. Nevertheless, other potential requirements for the statistics such as anisotropies due to the existence of mean fields can be incorporated in the framework by an alternative form of the defining correlation tensor in Eq. (40) [79,80]. For the case of other turbulent systems such as passive scalars, Rayleigh-Bénard convection or magnetohydrodynamic (MHD) turbulence, multipoint statistics of different fields can be generated by using cross correlation functions in the Gaussian ansatz, e.g., in the case of MHD, the cross correlation between velocity and magnetic field [81]. Moreover, spatiotemporal statistics of the velocity field u(x, t) could be introduced by generalizing the ansatz for the correlation functions to C ξ,ij (r, τ ). Depending on how spatial and temporal fluctuations are distributed, one could thus modify the correlation function in Eq. (7) accordingly. Furthermore, due to the underlying fractional Ornstein-Uhlenbeck processes, which dictate the monofractal behavior, the framework should also be general enough to include other phenomenological models of turbulence. The She-Leveque model [82], for instance, could be reproduced by choosing a log-Poisson distribution instead of the lognormal distribution (11), possibly in combination with another scale parametrization for ρ ξ (r) in Eq. (9).

V. OUTLOOK
We have proposed an explicit form of the joint multipoint PDF of a multifractal field. The proposed methodology relies on three different elements: i.) the stationary solution of a fractional Ornstein-Uhlenbeck process with Hurst exponent H, ii.) a re-parametrization r → ρ ξ (r) of the scale separation in the corresponding correlation function C ξ (r), and, iii.) introduction of fluctuations of the parameter ξ. As highlighted by the joint multi-point PDF (19), these three elements introduce fluctuating correlations which are determined by the parameter or mixing distribution g(ξ). The proposed Gaussian scale mixture (19) has been used to generate synthetic multifractal fields. Statistical quantities such as the one-point PDF (which is Gaussian due to the underlying fractional Ornstein-Uhlenbeck process), the increment PDF (two-point quantity, which exhibits non-Gaussian features at small scales), the PDF of increment multipliers (three-point quantity), and even the statistics of longitudinal and transverse increments have been reproduced and agree well with empirical findings from turbulence research. Hence, the joint multipoint PDF (39) defines an appropriate statistical model of a turbulent flow [29].
Future work should address the inclusion of skewness of the velocity increment distribution into the framework which would thus make it more truthful to the original turbulence problem (i.e., the 4/5-law derived from the Navier-Stokes equation). This can be done by the introduction of a mean in the Gaussian ansatz (4), which has already been discussed in the context of superstatistics [49,83]. Moreover, dissipative effects might by incorporated by embedding the Ornstein-Uhlenbeck process which effectively leads to the differentiability of the velocity field [84,85].