A Bregazzi

P F Griffin

A S Arnold

D P Burt

G Martinez

R Boudot

J Kitching

E Riis

J P Mcgilligan
email: james.mcgilligan@strath.ac.uk

A simple imaging solution for chip-scale laser cooling

We demonstrate a simple stacked scheme that enables absorption imaging through a hole in the surface of a grating magneto-optical trap (GMOT) chip, placed immediately below a micro-fabricated vacuum cell. The imaging scheme is capable of overcoming the reduced optical access and surface scatter that is associated with this chip-scale platform, while further permitting both trapping and imaging of the atoms from a single incident laser beam. The through-hole imaging is used to characterise the impact of the reduced optical overlap volume of the GMOT in the chip-scale cell, with an outlook to an optimised atom number in low volume systems.

Micro-fabricated physics packages based on the measurement of thermal atomic ensembles have been integrated in metrological instruments, ranging from magnetometers to interferometers and clocks [START_REF] Knappe | A microfabricated atomic clock[END_REF][START_REF] Kitching | Chip-scale atomic devices[END_REF][START_REF] Carlé | Exploring the use of Ramsey CPT spectroscopy for a microcell-based atomic clock[END_REF] . Most such instruments rely on spin transitions in atoms, and require the presence of buffer gas mixtures or cell wall coatings to reduce the relaxation rate due to wall collisions. For clocks, the presence of these buffer gases causes temperature dependent frequency offsets and require careful temperature stabilisation to achieve good medium-and long-term stability [START_REF] Kozlova | Temperature and pressure shift of the cs clock transition in the presence of buffer gases: Ne, N 2 , Ar[END_REF][START_REF] Hasegawa | Microfabrication of cesium vapor cells with buffer gas for MEMS atomic clocks[END_REF][START_REF] Seltzer | High-temperature alkali vapor cells with antirelaxation surface coatings[END_REF] . The use of cold atom ensembles avoids these difficulties and can result in an increased interaction time and improved absolute accuracy compared to their thermal atom counterparts.

While the transition to cold atoms offers clear advantages over thermal atom packages [START_REF] Eckel | Challenges to miniaturizing cold atom technology for deployable vacuum metrology[END_REF] , the additional experimental size and complexity associated with laser cooling has limited the deployability and application range of cold atom devices [START_REF] Rushton | Contributed review: The feasibility of a fully miniaturized magneto-optical trap for portable ultracold quantum technology[END_REF][START_REF] Marlow | A review of commercial and emerging atomic frequency standards[END_REF][START_REF] Muquans | Muclock data sheet[END_REF] . Significant efforts have been made on the miniaturisation of cold atom components to facilitate portability for next-generation atomic sensors [START_REF] Burrow | Centilitre-scale vacuum chamber for compact ultracold quantum technologies[END_REF][START_REF] Little | A passively pumped vacuum package sustaining cold atoms for more than 200 days[END_REF][START_REF] Nichols | Magneto-optical trap with millimeter ball lenses[END_REF][START_REF] Nshii | A surface-patterned chip as a strong source of ultracold atoms for quantum technologies[END_REF][START_REF] Ravenhall | High-flux, adjustable, compact coldatom source[END_REF][START_REF] Kang | Magneto-optic trap using a reversible, solid-state alkali-metal source[END_REF][START_REF] Mcgehee | Magneto-optical trapping using planar optics[END_REF][START_REF] Sitaram | Confinement of an alkaline-earth element in a grating magnetooptical trap[END_REF] . However, the fabrication complexity of many of these components remains unfavourable for mass production, preventing their adoption in commercial applications.

Recent work has demonstrated a dramatic reduction of the cooling platform to the chip-scale by combining an anodically bonded glass-silicon-glass vapour cell with a diffractive optical element that redirects a single incident beam into the components required for laser cooling [START_REF] Mcgilligan | Laser cooling in a chip-scale platform[END_REF] . However, such systems have demonstrated that the detection of cold atoms is made difficult by the reduced optical access of the cell and light scattering from the grating and cell surfaces. As such, the authors required adopting a non-trivial two-photon spectroscopy scheme for improved detection, greatly increasing both the size and complexity of the optical system [START_REF] Mcgilligan | Laser cooling in a chip-scale platform[END_REF][START_REF] Sheludko | State-selective imaging of cold atoms[END_REF] .

In this letter we demonstrate a simple imaging solution for chip-scale laser cooling platforms. The stackable structure of the apparatus provides simplicity in alignment, as well as enabling future scalability for device mass production. A hole is laser cut in the centre of the grating chip to enable on-axis ab-sorption imaging from the cooling beam without significantly degrading atom trapping. Although grating chips with a central hole have been used to reduce surface reflections [START_REF] Imhof | Two-dimensional grating magneto-optical trap[END_REF] , as a source for cold electrons [START_REF] Franssen | Compact ultracold electron source based on a grating magnetooptical trap[END_REF] , and for Zeeman slowing with an alkali [START_REF] Barker | Single-beam zeeman slower and magneto-optical trap using a nanofabricated grating[END_REF] and alkaline-earth metal source [START_REF] Sitaram | Confinement of an alkaline-earth element in a grating magneto-optical trap[END_REF] , the impact of the reduced grating surface area on the MOT number has not yet been quantified or used as an imaging axis. The ability to measure the trapped atom number in such a chip-scale system is used to characterise the impact of the 3 mm internal height constraint on the optical overlap volume, with an outlook to an optimised atom number.

A simplified schematic for the cooling and imaging of 87 Rb is shown in Fig 1 . The incident light is derived from a single volume-Bragg-grating laser (VBG) frequency stabilised using saturated absorption spectroscopy. A double-pass acoustooptic modulator (AOM) shifts the light frequency to be approximately 8 MHz red detuned from the 780 nm D 2 F = 2 → F = 3 cycling transition, while also enabling frequency and intensity control for the imaging process. Re-pumping is achieved by modulating a free-space electro-optic modulator (EOM) at 6.5 GHz to generate 5 % sidebands on the carrier frequency. The light is then fibre coupled into a single-mode, polarization maintaining fibre and passed to the cooling platform.

From the fibre, we expand an ≈30 mW beam to a 1/e 2 radius of ≈ 1.6 cm to flatten the intensity distribution at the grating surface. The trap beam is then circularly polarised with a quarter wave-plate and aligned onto the grating chip, mounted externally to the actively pumped chip-scale cell. A pair of anti-Helmholtz coils are used to produce a trapping field within the cell volume, with an axial gradient of ≈1.5 mT/cm (15 G/cm).

The cell, shown in Fig. 1, is composed of a 3 mm thick silicon wafer, sandwiched between two anodically bonded aluminosilicate glass wafers of 0.7 mm thickness [START_REF] Dellis | Low helium permeation cells for atomic microsystems technology[END_REF] . The cell was fabricated with a square central region of dimensions 2.5 cm×2.5 cm, to enable cooling with a 2 cm wide grating chip. The upper window of the cell is mechanically drilled with a 5 mm hole to enable active pumping from an ion pump through a copper pinch-off tube, adhered to the upper cell surface. Square cavities are cut in the silicon walls to house nonevaporable getters (NEGs) for future passive pumping measurements.

The grating chip, shown in Fig. 1 (a), has a 1100 nm period in a tri-segment geometry. The central region of the chip is laser cut to provide a through-hole axis for absorption imaging. This central region of the grating plays a limited role in forming the trap overlap volume while also providing a unfavourable zeroth order reflection, such that its removal has a minimal impact on atom number. A two-lens imaging telescope and camera are placed behind the grating hole to detect atomic absorption from the GMOT.

While extracting atom number with the through-hole system, fluorescence measurements from the MOT loading curves were simultaneously measured using a separate detection telescope with a spatially selective focal plane (described in Ref. [START_REF] Boudot | Enhanced observation time of magneto-optical traps using micro-machined non-evaporable getter pumps[END_REF] ) at an angle of around 30 degrees to the grating surface (not shown in Fig. 1). The loading curves were used for the extraction of background pressure at the MOT location using the relation between MOT lifetime and background rubidium pressure [START_REF] Burrow | Centilitre-scale vacuum chamber for compact ultracold quantum technologies[END_REF][START_REF] Arpornthip | Vacuum-pressure measurement using a magneto-optical trap[END_REF][START_REF] Moore | Measurement of vacuum pressure with a magnetooptical trap: A pressure-rise method[END_REF] , while also checking the validity of the through-hole measured atom number. We note that the atom number measured from fluorescence imaging at this angle well matched the atom number extracted from conventional orthogonal fluorescence imaging. Additionally, an average ratio of 1.6±0.4 was observed between the calculated atom number from the through-hole absorption and orthogonal fluorescence during our data sets. However, the alignment of the angled fluorescence imaging axis was complicated due to the surface scatter and diffracted orders further restricting the available imaging angles and contributing to an increased noise level on the atom number extraction. The complexity in aligning the camera at a position and angle that minimises the impact of surface scatter does not meet the needs of a simple, mass-producible device.

The experimental procedure is initiated by turning the trap coils on, followed by a 500 ms loading time at 8 mW/cm 2 peak intensity. The 500 ms load time enables resolving the vacuum pressure down to 1.4 ×10 -6 Pa (1.4×10 -8 mbar) [START_REF] Mcgilligan | Grating chips for quantum technologies[END_REF] . The trap coils are then turned off while concurrently decreasing the incident beam below saturation to around 115 µW/cm 2 and bringing the frequency on resonance. This serves the dual purpose of reducing image distortion due to diffraction effects within the atom cloud while also maximising the signal contrast. We note that the measured atom number from this method did not differ with the addition of a static magnetic field, provided from a Helmholtz pair along the imaging axis, to aid the optical pumping of the atoms. An initial image I 1 is taken with the MOT present for an exposure time of 25 µs. We then wait 100 ms such that the trapped atoms are no longer present before taking a second image I 2 . Finally the trap beam is turned off with an extinction ratio of 62 dB and a dark background image I 3 is acquired. These three images are then processed in terms of optical density (OD) using the equation OD = ln I 1 -I 3 I 2 -I 3 . An objective lens of focal length f 1 = 15 mm and an imaging lens f 2 = 30 mm, arranged in 2 f 1 : 2 f 2 configuration, is used to image the atoms, with the objective placed directly behind the hole laser cut in the centre of the grating. This arrangement of lenses provides an improved signal-to-noise ratio for data extraction [START_REF] Smith | Absorption imaging of ultracold atoms on atom chips[END_REF] .

In the present work, the stretched state saturation intensity I sat =1.67 mW/cm 2 was used for calculating the trapped atom number from the fluorescence and absorption methods, as is consistent with previous publications of our group [START_REF] Nshii | A surface-patterned chip as a strong source of ultracold atoms for quantum technologies[END_REF] . Additionally, we found good agreement between the calculated atom numbers from the two imaging methods with the stretched state saturation intensity, compared to the average of all polarizations and magnetic sub-levels, I sat =3.57 mW/cm 2 .

A saturation parameter, S = I I sat , sums over the intensities of the single input trap beam and the three diffracted beams from the grating to account for the total intensity the atoms experience. This inclusion modifies S such that S = (1 + nη 1 sec θ ) I im I sat ≈ 2.4I im I sat , where I im is the imaging light intensity, n is the number of diffracted first orders interacting with the atoms, η 1 is the efficiency of the first diffracted order, and θ is the angle of diffraction. An example of the obtained MOT image is shown in Fig. 1 (b), where 2×10 5 atoms are trapped in the chip-scale cell.

Our investigation of the impact of removing the central region of the grating chip was initially carried out in a conventional glass vacuum chamber (2.5 cm×2.5 cm×10 cm with 1 mm thick glass walls) using fluorescence imaging, with the results shown in Fig. 2. We fabricated 5 identical grating chips, and laser cut the central regions for a hole size ranging from 1 mm to 5 mm. For each subsequent measurement of the atom number, the grating chip was carefully implemented immediately below the same vacuum cell, under the same conditions as the previous grating chip. We found that a hole diameter up to 3 mm did not significantly degrade the measured atom number, shown in Fig. 2 (a), due to the minimal impact that such a hole size has on the total optical overlap volume, shown in Fig. 2 (b). The theoretical optical overlap volume is numerically integrated for a 2 cm beam diameter, incident upon a 1100 nm period tri-segmented grating chip with a central hole of varying diameter. As the grating hole was increased to 4 mm, the atom number decreases, reducing to the point that no MOT was detected for a 5 mm hole diameter, even with additional optimisation of the MOT parameters, such as detuning and intensity. The atom number drops faster than would be expected from the decrease in the overlap volume, likely due to the additional sensitivity of the radiation pressure balance from the grating chip. This is emphasized by the red curve shown in Fig. 2 (a), where the theoretically calculated overlap volume is scaled to the atom number through N ∝ V 1.2 and normalised to the measured atom number for the 1 mm hole. For a grating hole diameter of 4 mm, a factor of ∼2 difference between the expected and measured atom number is observed. However, the increased hole diameter aids an improved numerical aperture, also plotted in Fig. 2 (b) for a MOT position of 2 mm above the grating surface, such that a larger hole diameter is favourable for imaging with the small atom numbers that have been observed in passively pumped vacuum cells [START_REF] Boudot | Enhanced observation time of magneto-optical traps using micro-machined non-evaporable getter pumps[END_REF] . Taking this into account, a 3 mm hole diameter was selected for through-hole imaging.

Once the through-hole imaging was well aligned and in focus, the glass cell was replaced by the actively pumped chipscale vacuum cell. The chip-scale cell was pumped down to an initial pressure of 10 -6 Pa (10 -8 mbar), measured using an ion pump. The validity of the ion pump reading was later verified with the pressure calculated from MOT loading curves. Following pump down, a resistively heated dispenser within the larger vacuum chamber was used to provide a moderate rubidium density. With the grating overlap volume aligned within the chip-scale cell, an image of the cold atoms was extracted from the through-hole scheme, highlighted in Fig. 1. The initial atom number measurement yielded a total of 2×10 5 atoms, which was an order of magnitude lower than what had been observed for the glass cell.

To investigate the cause of the reduced atom number, the diameter of the incident cooling light onto the grating surface, d, was reduced while measuring the trapped atom number for both the glass and chip-scale cells. This measurement provides an insight to the impact of the reduced vacuum volume of the chip-scale cell on the grating chip's optical overlap volume, and therefore trapped atom number. During this process a comparable background vapour density of ≈ 3.1×10 6 cm -3 and total pressure of ≈ 3 × 10 -5 Pa (≈ 3 × 10 -7 mbar) was maintained in both vacuum systems. The resulting data set is shown in Fig. 3 (a) in blue and red for the glass and chip-scale cells, respectively.

At the largest values of d an order of magnitude difference in the atom number is observed between the glass and chipscale cell, with the glass cell data reaching a maximum around 10 6 atoms, in line with previous atom numbers measured with a similar incident intensity [START_REF] Mcgilligan | Phase-space properties of magneto-optical traps utilising micro-fabricated gratings[END_REF] . Concurrently, the red data set extracted from the chip-scale cell reaches a maximum around 10 5 atoms. When d is reduced from this maximum atom number, there is initially no impact on the optical overlap volume, since the majority of the overlap volume is outside the cell's internal thickness and therefore cannot play a role in the cooling process. As d is decreased below ≈ 17 mm, the overlap volume is reduced below the internal height of the chip-scale cell, contributing to a degradation of the atom number.

The impact of the reduced optical overlap volume is emphasised by the theoretical atom number expected from the grating MOT as a function of the incident beam diameter, shown with solid curves in Fig. 3 (a). This atom number, N, is derived from a numerical integration of the optical over-lap volume, V , for a 1100 nm period linear grating, where N ∝ V 1.2 ∝ d 3.6 14 . The red curve represents the overlap volume from the grating chip with a 3 mm diameter hole at its centre. A vertical restriction is placed on the integration to account for the 3 mm thickness of the silicon frame. The overlap volume is then converted to a normalised atom number with a maximum corresponding to the largest measured atom number from the glass cell. A total potential overlap volume of 0.23 cm 3 is possible for a 2 cm diameter incident upon a grating chip with a 3 mm diameter hole, held 1 mm from the cell window. However, when the grating is mounted similarly for the chip-scale cell, only 21% of the relative overlap volume coincides within the cell vacuum volume. This process is illustrated in Fig. 3 The dashed black line in Fig. 3 shows a d 3.6 scaling of the atom number, expected for a constant intensity of the incident beam [START_REF] Lindquist | Experimental and theoretical study of the vapor-cell zeeman optical trap[END_REF][START_REF] Gibble | Improved magneto-optic trapping in a vapor cell[END_REF] that expands in three-dimensions. Previous studies have demonstrated that significantly smaller overlap volumes for 6-beam MOTs [START_REF] Hoth | Atom number in magneto-optic traps with millimeter scale laser beams[END_REF] and pyramid MOTs [START_REF] Pollock | Characteristics of integrated magneto-optical traps for atom chips[END_REF] are reduced to a d 6 scaling, as the atom number is then limited by the stopping distance of the overlap volume rather than the cooling light intensity that dominates the process for larger beam diameters. We note that the measured atom number in the chip-scale cell is lower than the theoretical V 1.2 , shown as a solid red curve. The cause of this is currently unknown and will be investigated in future studies.

When limited by atom shot noise, the stability of cold atom sensors scales inversely with the square-root of the atom number. In this case, an order of magnitude improvement in the chip-scale cell would be favourable. Expanding the overlap volume simulation to a 5 mm thick silicon frame enables a 71% overlap between the optical overlap volume and the vacuum volume of the chip-scale cell when the grating is held 1 mm from the cell lower window. Using identical scaling parameters to the 3 mm chip-scale cell curve shown in Fig. 3 (a) results in an achievable atom number of 10 6 atoms for a 2 cm incident beam diameter in this thicker cell geometry. Thus, we plan to fabricate thicker silicon cells in the near-future to enable a wider range of applications for this chip-scale laser cooling platform. Additionally, a steeper angle of diffraction with a two-dimensional grating geometry would provide an overlap volume that is largest at the chip surface for an improved atom number within the reduced vacuum volume compared to the tri-segmented chip used here. The improved atom number will enable studies into passive pumping with NEGs in an aluminosilicate based chip-scale cell [START_REF] Boudot | Enhanced observation time of magneto-optical traps using micro-machined non-evaporable getter pumps[END_REF] .

In conclusion, we have demonstrated a simple imaging solution for atom number extraction in a chip-scale cooling platform. The laser cut hole in the grating chip has shown a negligible atom number degradation, while enabling a sufficient numerical aperture for imaging small atom numbers in future experiments. We have shown that the 3 mm height of the micro-fabricated cell reduces the achievable overlap volume, and hence atom number, when coupled with a micro- fabricated grating chip. However, moving to the machining of thicker silicon wafers will enable a work around to bring this platform to a competitive stance for quantum technologies. In addition it may also be advantageous to apply an anti-reflection coating to the cell windows during fabrication due to the increased etaloning observed in the through-hole images. While the effect of these fringes on the images was minimal, especially for larger MOTs, and can be mitigated through alignment there are also several fringe removal algorithms that can be applied in post-processing to further reduce their impact [START_REF] Niu | Optimized fringe removal algorithm for absorption images[END_REF][START_REF] Xiong | Enhanced principle component method for fringe removal in cold atom images[END_REF] . Glass cell Chip-scale cell d 3.6 (a)

FIG. 1 .

 1 FIG. 1. Experimental set-up. CCD: charge coupled device, AHC: anti-Helmholtz coils, λ /4: quarter-wave plate. (a): Image of the grating chip with a central hole removed. (b): An on-axis, through-hole absorption image of ≈ 2×10 5 atoms within the chip-scale platform. This image was taken with a peak beam intensity of 115 µW/cm 2 at zero detuning.

FIG. 2 .

 2 FIG. 2. Critical parameters as a function of the grating centre hole diameter for an incident 2 cm beam diameter. (a): The measured atom number from orthogonal fluorescence imaging in the conventional glass cell. The measured error bars are smaller than the shown data points. The red curve shows the normalized N ∝ V 1.2 for the overlap volume shown in (b). (b): Theoretically calculated optical overlap volume and numerical aperture from the grating chip assuming that the grating hole is the aperture stop of the imaging system.

  (b)-(d), where the numerically calculated geometry of the optical overlap volume is shown for different values of d. The region of the overlap volume occupied by the 3 mm tall chip-scale cell is highlighted to show the clipping of the overlap volume as a function of d.

FIG. 3 .

 3 FIG. 3. (a): Atom number measured from absorption imaging as a function of the incident beam diameter, d. Blue/red data points represent measurements made in a glass/chip-scale cell. Error bars were calculated from 5 subsequent atom number measurements at each beam diameter. Where error bars are not visible they are smaller than the data point. The error bars at lower values of d are noticeably larger due to the atom number approaching the detection noise floor. Solid curves are theoretical expectations of the atom number, determined from numerical integration of the optical overlap volume from the grating chip as a function of d. The normalised theoretical atom number is pinned to the maximum measured atom number for the glass cell. The red line shows an overlap volume that is restricted to a 3 mm height due to the vacuum volume of the chip-scale cell. The dashed line represents a d 3.6 fit to the atom number in the glass cell. (b)-(d) illustrate the optical overlap volume of the grating chip as a function of d with a highlighted 3 mm height of the chip-scale cell. The grey line highlights the grating chip position.

ACKNOWLEDGMENTS

The authors acknowledge funding from Defence Security and Technology Laboratory, Engineering and Physical Sciences Research Council (EP/T001046/1), and Defence and Security Accelerator. A. B. was supported by a PhD studentship from the Defence Science and Technology Labo-

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

ratory (Dstl). The authors gratefully acknowledge W. R. McGehee and D. S. Barker (NIST) for the careful reading of the manuscript before submission. Additionally, the authors would like to thank B. Lewis and M. Himsworth for useful conversations.

A. Conflict of interest

The authors declare they have no competing interest. Approved for Public Release, Distribution Unlimited.


Qing Chu
email: qing.chu@univ-mlv.fr

CONVERGENCE OF WEIGHTED POLYNOMIAL MULTIPLE ERGODIC AVERAGES

Keywords: Mathematics Subject Classification. 37A05, 37A30 Weighted ergodic averages, Universally good sequences, Wiener-Wintner Ergodic Theorem, Nilsequences

In this article we study weighted polynomial multiple ergodic averages. A sequence of weights is called universally good if any polynomial multiple ergodic average with this sequence of weights converges in L 2 . We find a necessary condition and show that for any bounded measurable function φ on an ergodic system, the sequence φ(T n x) is universally good for almost every x. The linear case was covered by Host and Kra.

introduction

In his innovative proof of Szemerédi's Theorem via ergodic theory, Furstenberg introduced certain multiple ergodic averages. There have been many results on these and other nonconventional ergodic averages, including the multiple ergodic theorems of Host and Kra [START_REF] Host | Nonconventional ergodic averages and nilmanifolds[END_REF], [START_REF] Host | Convergence of polynomial ergodic averages[END_REF], Leibman [START_REF] Leibman | Convergence of multiple ergodic averages along polynomials of several variables[END_REF], Ziegler [START_REF] Ziegler | A non-conventional ergodic theorem for a nilsystem[END_REF]. Recently Host and Kra studied weighted ergodic theorems for multiple averages along arithmetic progressions. We give a generalization of this result for polynomial averages, showing: Theorem 1.1. Let (Y, ν, S) be an ergodic system and φ ∈ L ∞ (ν). Then there exists Y 0 ⊂ Y with ν(Y 0 ) = 1 such that, for every y 0 ∈ Y 0 , every system (X, µ, T ), every r ≥ 1, all integer polynomials p 1 , . . . , p r and all functions f 1 , . . . , f r ∈ L ∞ (µ), the averages

1 N N -1 n=0 φ(S n y 0 )T p1(n) f 1 • • • T pr(n) f r converge in L 2 (µ) .
Throughout this article, by integer polynomial we mean a polynomial all of whose coefficients are integers.

The case of p i (n) = in was proved by Host and Kra [START_REF] Host | Uniformity seminorms on l ∞ and applications[END_REF].

Remark 1.2. One may wonder why in the theorem S n is not replaced by S p(n) for some integer polynomial p(n). In fact, the latter would be much harder to prove, even in the simplest case that r = 1 and p 1 (n) = n. This problem is equivalent to showing the convergence of the averages of φ(S p(n) y 0 )e 2πint for all t ∈ T and all y 0 ∈ Y 0 , where Y 0 does not depend on t. But this problem reduces to a question of almost everywhere convergence of multiple ergodic averages along polynomials, and this question is out of reach for the moment.

Note that the set Y 0 does not depend on X or on f i , i = 1, . . . , r. We say that for every y 0 ∈ Y 0 , the sequence φ(S n y 0 ) is universally good for the convergence in the mean of polynomial multiple ergodic averages.

For r = 1 and p(n) = n, the result follows immediately from the classical Wiener-Winter ergodic Theorem [START_REF] Wiener | Harmonic analysis and ergodic theory[END_REF] and a corollary of the Spectral Theorem. We follow a similar strategy, generalizing the proof in [START_REF] Host | Uniformity seminorms on l ∞ and applications[END_REF] along arithmetic progressions to polynomial progressions, but need to address some deeper technical issues.

We first recall some definitions, see [START_REF] Bergelson | with an Appendix by I. Ruzsa. Multiple recurrence and nilsequences[END_REF] and [START_REF] Host | Uniformity seminorms on l ∞ and applications[END_REF] for details. Let G be a k-step nilpotent Lie group and Γ ⊂ G be a discrete, cocompact subgroup of G. The compact manifold X = G/Γ is called a k-step nilmanifold. The Haar measure µ of X is the unique probability measure invariant under the left translations x → gx of G on X. Letting T denote left multiplication by a fixed element α ∈ G, we call (X, µ, T ) a k-step nilsystem. Let f : X → C be a continuous function, x 0 ∈ X, then the sequence (f (α n x 0 ) : n ∈ Z) is called a basic k-step nilsequence. The family of basic k-step nilsequences forms a subalgebra of l ∞ . Under the uniform norm || • || ∞ of l ∞ , we call a uniform limit of basic k-step nilsequences a k-step nilsequence.

The proof of Theorem 1.1 is broken down into two pieces. First we give a convergence criterion for weighted polynomial multiple ergodic averages.

Theorem 1.3 (Convergence Criterion for Weighted Averages). For any r, b ∈ N, there exists an integer K ≥ 1 with the following property: for any bounded sequence c = (c n : n ∈ Z), if the averages

1 N N -1 n=0
c n d n converge as N → ∞ for every K-step nilsequence d = (d n : n ∈ Z). Then for every system (X, µ, T ), all f 1 , . . . , f r ∈ L ∞ (X), and all integer polynomials p 1 , . . . , p r of degree ≤ b, the averages

1 N N -1 n=0 c n T p1(n) f 1 • T p2(n) f 2 • • • T pr(n) f r (1) converge in L 2 (X).
The bulk of this paper is devoted to the proof of this theorem. Then our main result follows from the following Generalized Wiener-Wintner Theorem proved by Host and Kra in [START_REF] Host | Uniformity seminorms on l ∞ and applications[END_REF]. The case of a polynomial version of the Wiener-Wintner theorem was proved by Lesigne ( [7], [START_REF] Lesigne | Spectre quasi-discret et théorème ergodique de Wiener-Wintner pour les polynômes[END_REF]).

Theorem 1.4 (Generalized Wiener-Wintner Theorem [START_REF] Host | Uniformity seminorms on l ∞ and applications[END_REF]). Let (X, µ, T ) be an ergodic system and φ be a bounded measurable function on X. Then there exists X 0 ⊂ X with µ(X 0 ) = 1 such that for every x ∈ X 0 , the averages

1 N N -1 n=0 φ(T n x)b n converge as N → ∞ for every x ∈ X 0 and every nilsequence b = (b n : n ∈ Z).
While nilsequences do not appear in the statement of Theorem 1.1, they are used as tools in its proof. Both Theorems 1.3 and 1.4 are of interest on their own, as results on nilsequences.

proof of theorem 1.3

Using a standard ergodic decomposition argument, it suffices to prove Theorem 1.3 for ergodic systems.

2.1. We first remind the reader of a definition from Leibman's paper [START_REF] Leibman | Pointwise convergence of ergodic averages for polynomial sequence of translation on a nilmanifold[END_REF]. We call a sequence {g(n

)} n∈Z with values in a nilpotent group G a polynomial sequence, if g(n) is of the form g(n) = a p1(n) 1 . . . a pm(n) m
, where a 1 , . . . , a m ∈ G and p 1 , . . . , p m are polynomials taking integer values on the integers.

Before stating the next proposition, we explain briefly its meaning: we can view the sequence of values of a continuous function along a polynomial sequence on a nilmanifold as the sequence of values of some other continuous function along an ordinary orbit of some other nilsystem.

Proposition 2.1. Let (X = G/Γ, T ) be a nilsystem, x 0 ∈ X, p be an integer polynomial, and f ∈ C(X). Then there exists a nilsystem (Y, S), y 0 ∈ Y , h ∈ C(Y ), such that f (T p(n) x 0 ) = h(S n y 0 ) for every n.

Proof. Let (X = G/Γ, T ) be a nilsystem. Suppose T x := αx, for some α ∈ G.

Then T p(n) x = α p(n) x. Let g(n) := α p(n) , then g is a polynomial sequence in G.
Let π : G → X be the factorization mapping. We will assume that

x 0 = π(1 G ); otherwise if x 0 = π(γ), γ ∈ G, we write g(n)x 0 = g(n)γπ(1 G ), and replace g(n) by g(n)γ.
Now we have a nilpotent Lie group G, a discrete cocompact subgroup Γ and a polynomial sequence g in G. By Proposition 3.14 in Leibman's paper [START_REF] Leibman | Pointwise convergence of ergodic averages for polynomial sequence of translation on a nilmanifold[END_REF], there exist a nilpotent Lie group G, a discrete cocompact subgroup Γ, an epimorphism η : G → G with η( Γ) ⊆ Γ, a unipotent automorphism τ of G with τ ( Γ) = Γ, and an element c ∈ G such that

g(n) = η( τ n (c)), n ∈ Z.
Let X = G/ Γ and let π : G → X be the factorization mapping. The epimorphism η : G → G factors to a map X → X, we also denote it by η, which is onto and satisfies, π

• η = η • π.
The map τ induces a homomorphism X → X, which we also denote it by τ . It satisfies,

τ • π = π • τ . Let x 0 = π(1 G ), then η( τ n (c x 0 )) = g(n)x 0 , n ∈ Z.
Let G be the extension of G by τ , then G is a nilpotent Lie group (see Proposition 3.9 in [START_REF] Leibman | Pointwise convergence of ergodic averages for polynomial sequence of translation on a nilmanifold[END_REF]). Let τ be the element in G representing τ , so that τ ( α) = τ α τ -1 for any α ∈ G, the multiplication of G is given by this formula. We have

G = { g τ n : g ∈ G, n ∈ Z} and G is open in G.
Let Γ be the subgroup of G spanned by Γ and τ . As τ ( Γ) = Γ, we have Γ = { γ τ n : γ ∈ Γ, n ∈ Z} and Γ ∩ G = Γ. By the definition of the relative topology, Γ is a discrete subgroup of G. Moreover, G/ Γ = ( G Γ/ Γ) can be identified with G/( Γ ∩ G) = G/ Γ = X. We write π : G → X for the quotient map.

Let x ∈ X and g ∈ G with π(g) = x. We have

τ ( x) = π( τ (g)) = π( τ g τ -1 ) = π( τ g) = τ x, because τ -1 ∈ Γ. So for every n, g(n)x 0 = η( τ n (c x 0 )) = η( τ n (c x 0 )).
Let Y = ( G/ Γ, S) = ( X, S), S x = τ x, and let h = f • η, and y 0 = c x 0 . This system and this function satisfy the announced properties.

2.2. We recall a few properties of the seminorms and the factors introduced in [START_REF] Host | Nonconventional ergodic averages and nilmanifolds[END_REF]. Let (X, µ, T ) be an ergodic system. For an integer k ≥ 0, we write X k] for the map T × T × . . . × T , taken 2 k times. We define by induction a probability measure µ [k] on X [k] that is invariant under T

[k] = X 2 k and T [k] : X [k] → X [
[k] . Set µ [0] = µ. For k ≥ 0, let I [k] be the σ-algebra of T [k] -invariant subsets of X [k] . Then µ [k+1] is the relatively independent product of µ [k] over I k , which means for F, F ′ ∈ L ∞ (X [k] ), X [k+1] F ⊗ F ′ dµ [k+1] = X [k] E(F |I [k] )E(F ′ |I [k] ) dµ [k] .
For a bounded measurable function f , we define

|||f ||| k =   X [k] ε∈{0,1} k C |ε| f (x ε ) dµ [k] (x)   1/2 k , (2)
where C :

C → C is the conjugacy map z → z, ε = ε 1 ε 2 . . . ε k with ε i ∈ {0, 1} and |ε| = ε 1 + ε 2 + • • • + ε k . It is shown in [2] that for every k ≥ 1, ||| • ||| k is a seminorm on L ∞ (µ).
Moreover, for every k ≥ 2, X admits a factor Z k-1 such that, for every f ∈ L ∞ (µ), |||f ||| k = 0 if and only if E(f |Z k-1 ) = 0. One of the main results of [START_REF] Host | Nonconventional ergodic averages and nilmanifolds[END_REF] is that, for every k, Z k is an inverse limit of k-step nilsystems. We call this result the Structure Theorem.

Let (Z = Z 1 (X), m, T ) be the Kronecker factor of (X, µ, T ). For s ∈ Z, we define a measure µ s on X × X by

X×X f (x)f ′ (x ′ ) dµ s (x, x ′ ) = Z E(f |Z)(z) • E(f ′ |Z)(sz) dm(z).
For every s ∈ Z the measure µ s is invariant under T × T and is ergodic for m-almost every s. The ergodic decomposition of µ × µ under T × T is

µ × µ = Z µ s dm(s).
For each s ∈ Z such that (X × X, µ s , T × T ) is ergodic, and for each integer k ≥ 1, a measure (µ s ) [k] on (X × X) [k] can be defined in the same way as µ [k] . Furthermore, a seminorm ||| • ||| k,s on L ∞ (µ s ) can be associated to this measure in the same way as the seminorm ||| • ||| k is associated to µ [k] . It follows from the definition (2) that for every f ∈ L ∞ (µ),

|||f ||| 2 k+1 k+1 = Z |||f ⊗ f ||| 2 k k,s dm(s). (3)
2.3. We return to the proof of theorem 1.3. We may assume that the polynomials p 1 , . . . , p r are nonconstant and essentially distinct, that is p i -p j = constant for i = j.

The following theorem will be proved in the next section.

Theorem 2.2. For any r, b ∈ N, there exists k ∈ N, such that for any nonconstant essentially distinct polynomials p 1 , . . . , p r : Z → Z of degree ≤ b, for every ergodic system (X, µ, T ), every f 1 , . . . , f r ∈ L ∞ (X) with |||f 1 ||| k = 0, and any bounded sequence c = (c n : n ∈ Z), one has

lim N →∞ 1 N N -1 n=0 c n T p1(n) f 1 • • • T pr(n) f r L 2 (X) = 0. (4)
2.4. Now we give the proof of Theorem 1.3.

Proof of Theorem 1.3 from Proposition 2.1 and Theorem 2.2. The proof is exactly the same as the proof of Theorem 2.24 in [START_REF] Host | Uniformity seminorms on l ∞ and applications[END_REF]. For any r, b ∈ N, let k ∈ N be the integer in Theorem 2.2, let Z k-1 be the (k -1)-th factor of (X, µ, T ), as given by the Structure Theorem. By definition, if E(f 1 |Z k-1 ) = 0, then |||f 1 ||| k = 0, and by Theorem 2.2, the averages (4) converge to zero in L 2 (X). We say that the factor Z k-1 is the characteristic for the convergence of these averages. Therefore, it suffices to prove the result when the functions are measurable with respect to the factor Z k-1 .

Since Z k-1 is an inverse limit of (k -1)-step nilsystems by density, we can assume that (X, µ, T ) is a (k -1)-step nilsystem and that the functions f 1 , . . . , f r are continuous.

But in this case, by Proposition 2.1, for every x ∈ X, and all polynomials p 1 , . . . , p r , there exist nilsystems (Y 1 , S 1 ), . . . , (Y r , S r ), y i ∈ Y i , and g i ∈ C(Y i ), such that f i (T pi(n) x) = g i (S n i y i ), i = 1, . . . , r. Let K be the maximal order of the nilsystems (Y i , S i ), i = 1, . . . , r.

Then the system (Y = Y 1 × • • • × Y r , S = S 1 × • • • × S r ) is a K-step nilsystem. Let g : Y 1 × • • • × Y r →
R be given by g(y) = g(y 1 , . . . , y r ) = g 1 (y 1 ) • . . . • g r (y r ). So the sequence

{f 1 (T p1(n) x) • f 2 (T p2(n) x) • . . . • f r (T pr(n) x)} n∈Z = {g 1 (S n 1 y 1 ) • g 2 (S n 2 y 2 ) • . . . • g r (S n r y r )} n∈Z = {g(S n y)} n∈Z
is a K-step nilsequence and by hypothesis, the averages (1) converge for every x ∈ X.

proof of theorem 2.2

Note that our goal is very similar to the main result in Leibman's paper [START_REF] Leibman | Convergence of multiple ergodic averages along polynomials of several variables[END_REF], the only difference being that in our case we are dealing with the weighted averages. In fact, we can deduce Theorem 2.2 by following very closely the arguments of Leibman in his paper to cover our weighted case with some modifications. But here we adopt another way that allows us to deduce it directly from Leibman's result.

Proof of Theorem 2.2. We prove the result by several steps.

(i) The following assertion follows immediately from Theorem 3 in Leibman's paper [START_REF] Leibman | Convergence of multiple ergodic averages along polynomials of several variables[END_REF]:

Let r, b ∈ N be fixed. There exists an integer k = k(r, b) such that for every family of nonconstant essentially distinct polynomials p 1 , . . . , p r : Z 2 → Z of degree ≤ b, we have: for every ergodic system (X, B, µ, T ), and every f

1 , . . . , f r ∈ L ∞ (X) with |||f 1 ||| k = 0, one has lim N →∞ 1 N 2 0≤m,n<N T p1(m,n) f 1 • • • T pr(m,n) f r L 2 (X) = 0.
(ii) We use the notation of Section 2.2. It follows from (3) that |||f 1 ||| k+1 = 0 implies |||f 1 ⊗ f1 ||| k,s = 0 for m-almost every s. Using the previous result to the ergodic system (X × X, µ s , T × T ), one gets: if

|||f 1 ||| k+1 = 0, then lim N →∞ 1 N 2 0≤m,n<N T p1(m,n) f 1 • • • T pr(m,n) f r dµ 2 = lim N →∞ 1 N 2 0≤m,n<N Z X×X (T × T ) p1(m,n) f 1 ⊗ f1 • • • (T × T ) pr(m,n) f r ⊗ fr dµ s dm(s) ≤ Z lim N →∞ X×X 1 N 2 0≤m,n<N (T × T ) p1(m,n) f 1 ⊗ f1 • • • (T × T ) pr(m,n) f r ⊗ fr dµ s dm(s) ≤ Z lim N →∞ 1 N 2 0≤m,n<N (T × T ) p1(m,n) f 1 ⊗ f1 • • • (T × T ) pr(m,n) f r ⊗ fr L 2 (µs)
dm(s) = 0.

(iii) Expanding the square, one sees that

1 N N -1 n=0 c n T p1(n) f 1 • • • T pr(n) f r 2 L 2 is equal to 1 N 2 0≤m,n<N c n cm T p1(n) f 1 • • • T pr(n) f r T p1(m) f1 • • • T pr(m)
fr dµ, which is less than a constant times

1 N 2 0≤m,n<N T p1(n) f 1 • • • T pr(n) f r T p1(m) f1 • • • T pr(m) fr dµ , (5)
since by our assumption the sequence c n is bounded. Notice that p 1 (n), . . . , p r (n), p 1 (m), . . . , p r (m) is a family of 2r essentially distinct polynomials of degree at most b.

(iv) Combining the previous parts we see that if |||f 1 ||| k(2r,b)+1 = 0, where k(r, b) is defined in part (i), then for every family of nonconstant essentially distinct polynomials p 1 , . . . , p r : Z → Z of degree ≤ b, we have that the average (5) converges to 0. This proves Theorem 2.2.

Acknowledgement. We thank the referee for very useful remarks and the simplification of the proof.


J Friedrich

J Peinke

A Pumir

R Grauer

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 [START_REF] Mandelbrot | The fractal geometry of nature[END_REF]. 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 [START_REF] Barenblatt | Scaling, self-similarity, and intermediate asymptotics: dimensional analysis and intermediate asymptotics[END_REF][START_REF] Barenblatt | [END_REF][START_REF] Goldenfeld | Lectures on phase transitions and the renormalization group[END_REF]. Typical examples include subsurface hydrology [START_REF] Molz | [END_REF]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 [START_REF] Lévy | Processus stochastiques et mouvement brownien[END_REF][START_REF] Mandelbrot | [END_REF][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 [START_REF] Frisch | Turbulence[END_REF]. 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 [START_REF] Ghashghaie | [END_REF], hydraulic conductivity measurements in geophysics [15], or urban rent price fields [16,[START_REF] Lengyel | Multiscale urban analysis and modeling: Trends in the Ruhr Area[END_REF]. 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 [START_REF] Juneja | [END_REF][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 [START_REF] Frisch | Turbulence[END_REF]. 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 [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF]. 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)

f n (u 1 , x 1 ; . . . ; u n , x n ) = n i=1 δ(u i -u(x i )) , (1)
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 [START_REF] Nawroth | [END_REF][31][32][START_REF] Friedrich | Closure of the Lundgren-Monin-Novikov hierarchy in turbulence via a Markov property of velocity increments in scale[END_REF][START_REF] Sinhuber | [END_REF], 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 [START_REF] Friedrich | Closure of the Lundgren-Monin-Novikov hierarchy in turbulence via a Markov property of velocity increments in scale[END_REF][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 [START_REF] Mandelbrot | The fractal geometry of nature[END_REF] 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

ϕ[α] = e i dxα(x)u(x) , (2)
in terms of a superposition of characteristic functionals ϕ ξ [α] belonging to Gaussian sub-ensembles

ϕ[α] = dξg(ξ)ϕ ξ [α] . (3)
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 [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF]59,[START_REF] Lukassen | Progress in Turbulence VII[END_REF]) as

ϕ ξ [α] = e -1 2 dx dx α(x)C ξ (x-x )α(x ) , (4)
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

dw(x) = - 1 L w(x)dx + dB H (x) , (5)
where dB H (x) denotes the increment of fractional Gaussian noise with covariance

B H (x)B H (x ) = σ 2 2 |x| 2H + |x | 2H -|x -x | 2H .
(6) 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. [START_REF] Mardoukhi | [END_REF] 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

C ξ (r) = - σ 2 2 ρ ξ (r) 2H + σ 2 L 2H 2 Γ(2H + 1) cosh ρ ξ (r) L + σ 2 ρ ξ (r) 2H+1 4L(2H + 1) e -ρ ξ (r) L K H ρ ξ (r) L -e ρ ξ (r) L K H - ρ ξ (r) L , (7)
where we imply that r ≥ 0. Furthermore, we introduce Kummer's confluent hypergeometric function

1 F 1 (a, b, z) = Γ(b) Γ(b -a)Γ(a) 1 0 dte zt t a-1 (1 -t) b-a-1 , (8) 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

ρ ξ (r) = ξ A+µ ln+ x r+x ≺ r + x ≺ x µ 2 r , (9)
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

ln + x = ln x for x ≥ 1 , 0 for 0 < x < 1 . (10)
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 [START_REF] Mardoukhi | [END_REF]. Furthermore, in the limit r → ∞, the asymptotic behavior of Kummer's confluent hypergeometric function suggests

1 F 1 (a, b, z) ∼ Γ(b) e z z a-b Γ(a) + (-z) -a Γ(b-a)
and we obtain lim r→∞ C ξ (r) = 0. In the final step, we demand that the parameter ξ follows a lognormal distribution

g(ξ) = 1 √ 2πξ e -[ln ξ] 2 2 , (11)
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

PDF h(v, r) = δ(v -u(x + r) + u(x)) = 1 2π dwe -iwv e iw[u(x+r)-u(x)] = 1 2π dwe -iwv ϕ [α(x ) = w{δ(x -x -r) -δ(x -x)}] = 1 2π dwe -iwv ∞ 0 dξg(ξ)e -w 2 2 [2C ξ (0)-2C ξ (r)] . (12)
Here, we can introduce the second-order structure function

S 2,ξ (r) = 2C ξ (0) -2C ξ (r) (13)
which yields

h(v, r) = ∞ 0 dξg(ξ) 1 2πS 2,ξ (r) exp - v 2 2S 2,ξ (r) . (14)
Taking the moments of h(v, r) yields

[u(x + r) -u(x)] n = dvv n h(v, r) = ∞ 0 dξg(ξ)(n -1)!![S 2,ξ (r)] n/2 , (15)
and for large L, we can approximate S 2,ξ (r) ≈ σ 2 ρ ξ (r) 2H and obtain

[u(x + r) -u(x)] n =(n -1)!!σ n r + x ≺ x µ 2 nH r nH dξg(ξ)ξ nH A+µ ln + x r+x ≺ =e -1 2 n 2 H 2 A+µ ln + x r+x ≺ =C n r + x ≺ x µ 2 (nH-n 2 H 2 ) r nH , (16)
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

h(v, r) = ∞ 0 dξ e -ξ 2 2 exp       -   ln v σξ r+x ≺ x µH 2 r H    2 2H 2 A+µ ln+ x r+x ≺       2πvH A + µ ln + x r+x ≺ . ( 17
)
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 [START_REF] Mandelbrot | [END_REF]. Hence, by virtue of these three step process, we introduce strong correlations which culminate in the non-Gaussian statistics of field increments [START_REF] Lengyel | Multiscale urban analysis and modeling: Trends in the Ruhr Area[END_REF]. 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

f n (u 1 , x 1 ; . . . ; u n , x n ) = n i=1 δ(u i -u(x i )) . ( 18
)
From Eq. (3), we obtain the n-point PDF (19) 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 [START_REF] Mandelbrot | [END_REF]. Hence, the n-point PDF (19) belongs to the class of compound probability distributions [START_REF] Feller | An introduction to probability theory and its applications[END_REF]. 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 [START_REF] Dillon | [END_REF]. 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.

f n (u 1 , x 1 ; . . . ; u n , x n ) = dξg(ξ)f n,ξ (u 1 , x 1 ; . . . ; u n , x n ),
f n,ξ (u 1 , x 1 ; . . . ; u n , x n ) = 1 (2π) n det σ ξ e -1 2 u T σ -1 ξ u , (20)

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

f 1 (u 1 , x 1 ) = dξg(ξ) 1 2πC ξ (0) e - u 2 1 2C ξ (0) . (21)
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

f 1 (u 1 ) = 1 πΓ(2H + 1)σL H e - u 2 1 σ 2 L 2H Γ(2H+1) . (22)
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

C(r) = dξg(ξ)C ξ (r) , (23)
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-

ments δ r v(x) = u(x + r) -u(x) S n (r) = (δ r v) n = dξg(ξ)(n -1)!! [S 2,ξ (r)] n/2 , ( 24
)
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 [START_REF] Frisch | Turbulence[END_REF]. Here, ε is defined as the average local energy dissipation rate, which for the present one-dimensional velocity field is defined according to

ε = 2ν ∂u(x) ∂x 2 . ( 25
)
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 and4 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 [START_REF] Eyink | [END_REF]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 p(w, r, r

) = δ w - δ r v(x) δ r v(x) , (26)
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

p f Bm (w, r, r ) = 1 π γ(r, r ) [w -w 0 (r, r )] 2 + γ(r, r ) 2 , (27)
with scale and location parameter given by

γ(r, r ) = 1 - [r 2H + r 2H -(r -r) 2H ] 2 4r 2H r 2H r H r H , ( 28
)
w 0 (r, r ) = r 2H + r 2H -(r -r) 2H 2r H . ( 29
)
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 [START_REF] Rényi | Foundations of probability[END_REF]) and thus might not be observable in turbulence experiments or numerical simulations. Similar arguments have already been put forth by Siefert and Peinke [START_REF] Siefert | [END_REF] who reproduced the Cauchy distribution by an ordinary Ornstein-Uhlenbeck process.

The multiplier PDF ( 26) can be obtained from the twoincrement PDF

h 2 (v, r; v , r ) = δ(v -δ r v(x))δ(v -δ r v(x)) , (30) according to p(w, r, r ) = dv dv δ w - v v h 2 (v, r; v , r ) . ( 31
)
We thus derive the two-increment PDF from the threepoint PDF (Eq. ( 19) for n = 3) according to

h 2 (v, r; v , r ) = du 1 du 2 du 3 δ(v -u 3 + u 1 )δ(v -u 2 + u 1 ) dx 1 dx 2 dx 3 δ(r -x 3 + x 1 )δ(r -x 2 + x 1 )δ(x -x 1 ) f 3 (u 1 , x 1 ; u 2 , x 2 ; u 3 , x 3 ) , (32)
which results in a compound bivariate normal distribution

h 2 (v, r; v , r ) = dξg(ξ) 1 2π S 2,ξ (r)S 2,ξ (r ) 1 -ρ(r, r ) 2 exp   - v 2 S2(r) -2ρ(r,r )vv √ S 2,ξ (r)S 2,ξ (r ) + v 2 S 2,ξ (r) 2(1 -ρ(r, r ) 2 )    , (33)
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 ) . ( 34
)
Inserting the two-increment PDF into Eq. ( 31) yields

p(w, r, r ) = dξg(ξ) 1 π γ ξ (r, r ) [w -w 0,ξ (r, r )] 2 + γ ξ (r, r ) 2 . ( 35
)
Hence, the multiplier PDF is a compound Cauchy distribution with location and scale parameter defined as

γ ξ (r, r ) 2 = 1 -ρ(r, r ) 2 S 2,ξ (r) S 2,ξ (r ) , (36)
and w 0,ξ (r, r ) = S 2,ξ (r) S 2,ξ (r ) ρ(r, 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 rr, 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,[START_REF] Siefert | [END_REF]. 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 [START_REF] Eyink | [END_REF]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.

µ = 0.227 3µ/4 µ/2 µ/4 r, r ) (c)
FIG. 5. (a) Multiplier PDF p(w, r, r ) determined from numerical integration of Eq. (35). Dashed lines correspond to the prediction of the fractional Brownian motion model (27). (b) Cumulative PDF of the multipliers P (w, r, r ), dashed lines correspond to the fBm prediction P f Bm (w, r, r

) = 1 2 + 1 π arctan w-w 0 (r,r ) γ(r,r )
with scale and location parameter defined in Eqs. ( 28) and ( 29). (c) Kullback-Leibler divergence (38) as a function of the scale separation rr for different values of the intermittency coefficient µ. Deviations between the distributions become visible at small and large scale separations and are increasing with increasing intermittency coefficient µ.

The n-point PDF (19) can readily be generalized to the multipoint statistics of a three-dimensional velocity field according to

f n (u 1 , x 1 ; . . . ; u n , x n ) = dξg(ξ) 1 (2π) 3 det σ ξ e -1 2 ui,ασ -1 ξ,i,α;j,β u j,β , (39)
where greek indices indicate spatial components α = 1, 2, 3, and where we introduce the covariance matrix σ ξ,i,α;j,β = C ξ,αβ (x ix 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:

dξg(ξ) C ξ,tt (r) - 1 2r ∂ ∂r r 2 C ξ,rr (r) = 0 , (41)
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 [START_REF] Mandelbrot | [END_REF], 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 [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF].

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].

Due to its straightforward generalizations to multiple dimensions outlined in Sec. IV, the model is directly applicable to the specific tailoring of turbulent wind fields. A central problem in the field of wind energy is the interpolation of sparse wind field measurements [55]. It would thus be desirable to incorporate multiscaling features into a recently developed stochastic interpolation scheme of sparse measurements by fractional Brownian motion [86]. Moreover, a multiscale refinement method [START_REF] Sinhuber | [END_REF], which is based on the empirical knowledge of the three-point PDF, could be modeled by the compound distribution (19) for n = 3. Instead of drawing velocity field realizations similar to the one in Fig. 1 from the full n-point statistics, latter refinement method could thus be used to generate surrogate data at a relatively low computational cost (although deviations from the K62 scaling at small scales might be expected). On a more basic level, the multipoint statistics (39) could also be used as a closure method for the multipoint hierarchy derived from the Navier-Stokes equation [40,[START_REF] Friedrich | Non-perturbative Methods in Statistical Descriptions of Turbulence[END_REF]. Other potential applications include the modeling of financial data [START_REF] Ghashghaie | [END_REF], geophysical settings (non-Gaussian features in the measurements of hydraulic conductivity [15] or subsurface hydrology [START_REF] Molz | [END_REF]), as well as the synthesis of astrophysical magnetic fields and their impact on cosmic particle transport [START_REF] Reichherzer | [END_REF][START_REF] Giacalone | [END_REF][90].

from financial support through the Project IDEXLYON of the University of Lyon in the framework of the French program "Programme Investissements d'Avenir" (ANR-16-IDEX-0005).

Appendix A: Derivation of the covariance of the fractional Ornstein-Uhlenbeck process

In this appendix we want to derive the stationary limit of the covariance of the fractional Ornstein-Uhlenbeck process (5), i.e., Eq. ( 7) with ρ ξ (r) set to r. The stationary solution of Eq. ( 5) reads

w(x) = x -∞ dB H (x )e -(x-x )/L , (A1)
which can be integrated by parts according to

w(x) = B H (x) - 1 L x -∞ dx e -(x-x )/L B H (x ) . (A2)
The covariance of the fractional Ornstein-Uhlenbeck process can thus be calculated from w(x 1 )w(x 2 ) = B H (x 1 )B H (x 2 )

- 1 L x1 -∞ dx 1 e -(x1-x 1 )/L B H (x 1 )B H (x 2 ) - 1 L x2 -∞
dx 2 e -(x2-x 2 )/L B H (x 1 )B H (x 2 )

+ 1 L 2 x1 -∞ dx 1 x2 -∞
dx 2 e -(x1-x 1 )/L e -(x2-x 2 )/L B H (x 1 )B H (x 2 ) . (A3)

We start with the second integral in Eq. (A3) and obtain where we assumed that x 1 ≤ x 2 .

FIG. 1 .

 1 FIG.1. Typical realization of the velocity field u(x) drawn from the compound n-point PDF(19) with spatial resolution n = 4096.

FIG. 2 . 6 FIG. 3 .

 263 FIG.2. Single-point PDF from 5000 realizations of velocity fields similar to the one depicted in Fig.1Dashed lines correspond to the prediction(22).

FIG. 4 .

 4 FIG.[START_REF] Goldenfeld | Lectures on phase transitions and the renormalization group[END_REF]. Evolution of the velocity increment PDF determined from 5000 realizations similar to Fig.1. Dashed lines correspond to the prediction[START_REF] Ghashghaie | [END_REF].

dx 2 0 dx 2 e 1 0 1 F 1 2H + 1 , 2H + 2, x 2 L 1 0

 202111121 e x 2 /L |x 1 | 2H + |x 2 | 2H -|x 1x 2 | 2H =Le x2/L |x 1 | 2H + 0 -∞ dx 2 e x 2 /L |x 2 | 2H + x2 0 dx 2 e x 2 /L |x 2 | 2H -x1 -∞ dx 2 e x 2 /L (x 1x 2 ) 2H -0 x1 dx 2 e x 2 /L (x 2x 1 ) 2H =Le x2/L |x 1 | 2H + ∞ -x 2 /L (x 2 ) 2H + |x 2 | 2H+1 dte x2t/L t 2H -∞ 0 dre (x1-r)/L r 2H -x2-x1 0 dre (x1+r)/L r 2H =Le x2/L |x 1 | 2H + Γ(2H + 1)L 2H+1 + |x 2 | 2H+1 2H + 1 e x1/L Γ(2H + 1)L 2H+1 e x1/L |x 2x 1 | 2H+1 dte (x2-x1)t t 2H =Le x2/L |x 1 | 2H + Γ(2H + 1)L 2H+1

dx 1 0 dx 2 e x 2 /L x 2H 2 -dx 2 e 1 dx 2 e 1 F 1 2H + 1 , 2H + 2, x 1 L+ 1 F 1 2H + 1 , 1 x2-x 1 0∞ dx 1 x2-x 1 0

 102221211111111111 e x 1 /L (-x 1 ) 2H + x1 0 dx 1 e x 1 /L (x 1 ) 2H + Le x1/L 0 -∞ dx 2 e x 2 /L (-x 2 ) 2H + x2 (x 1 +x 2 )/L (x 1x 2 ) 2H -(x 1 +x 2 )/L (x 2 x 1 ) 2H =Le x2/L Γ(2H + 1)L 2H+1 + |x 1 | 2H+1 2H + 1 Le x1/L Γ(2H + 1)L 2H+1 + |x 2 | 2H+1 2H + 1 dre (2x 1 +r)/L r 2H . (A6)Here, the remaining double integral can be evaluated by a partial integration which yields x1 -dre (2x 1 +r)/L r 2H = L 2

  dx 1 e x 1 /L |x 1 | 2H + |x 2 | 2H -|x 1x 2 | 2H =Le x1/L |x 2 | 2H + Γ(2H + 1)L 2H+1 1e x2/L + |x 1 | 2H+1 2H + 1 1 F 1 2H + 1, 2H + 2, dx 2 e (x 1 +x 2 )/L |x 1 | 2H + |x 2 | 2H -|x 1x 2 | 2H

				x1
				-∞
							x 1
							L
			+	|x 2 -x 1 | 2H+1 2H + 1	1 F 1 2H + 1, 2H + 2, -	x 2 -x 1 L	.	(A5)
							The double integral in Eq. (A3) is treated according to
		x1	x2
		dx 1
	-∞	-∞
						1 -e x1/L +	|x 2 | 2H+1 2H + 1	1 F 1 2H + 1, 2H + 2,	x 2 L
	-	|x 2 -x 1 | 2H+1 2H + 1	1 F 1 2H + 1, 2H + 2,	x 2 -x 1 L	.	(A4)
							Similarly, the first integral reads

  1 )/L (x 2x 1 ) 2H

	w(x 1 )w(x 2 ) = -	σ 2 2	(x 2 -x 1 ) 2H +	σ 2 L 2H 2	Γ(2H + 1) cosh	x 2 -x 1 L
	+	σ 2 (x 2 -x 1 ) 2H+1 4L(2H + 1)	e -x 2 -x 1 L	1 F 1 2H + 1, 2H + 2,	x 2 -x 1 L	) -e	x 2 -x 1 L	1 F 1 2H + 1, 2H + 2, -	x 2 -x 1 L	. (A7)
	=	L 2	e 2x1/L |x 2 -x 1 | 2H+1 2H + 1	1 F 1 2H + 1, 2H + 2,	x 2 -x 1 L	-	L 2	e 2x2/L	0 ∞	dre -r/L r 2H +	0	x2-x1	dre -r/L r 2H
	=	L 2	e 2x1/L |x 2 -x 1 | 2H+1 2H + 1	1 F 1 2H + 1, 2H + 2,	x 2 -x 1 L
			+	L 2H+2 2	e 2x2/L Γ(2H + 1) -	L 2	e 2x2/L |x 2 -x 1 | 2H+1 2H + 1	1 F 1 2H + 1, 2H + 2,	x 1 -x 2 L	.
	Inserting the terms (A4)-(A7) into the correlation func-	tion (A3) yields

ACKNOWLEDGMENTS

J.F. acknowledges funding from the Humboldt Foundation within a Feodor-Lynen fellowship and benefited