# A Linear Stochastic Model of Turbulent Cascades and Fractional Fields

Gabriel B. Apolinário, Geoffrey Beck, Laurent Chevillard, Isabelle Gallagher, Ricardo Grande

► To cite this version:

Gabriel B. Apolinário, Geoffrey Beck, Laurent Chevillard, Isabelle Gallagher, Ricardo Grande. A Linear Stochastic Model of Turbulent Cascades and Fractional Fields. 2023. [https://hal.science/hal-03919233v2](https://hal.science/hal-03919233v2)

HAL Id: hal-03919233

https://hal.science/hal-03919233v2

Preprint submitted on 4 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A LINEAR STOCHASTIC MODEL OF TURBULENT CASCADES AND FRACTIONAL FIELDS

GABRIEL B. APOLINÁRIO, GEOFFREY BECK, LAURENT CHEVILLARD, ISABELLE GALLAGHER, AND RICARDO GRANDE

**ABSTRACT.** Turbulent cascades characterize the transfer of energy injected by a random force at large scales towards the small scales. In hydrodynamic turbulence, when the Reynolds number is large, the velocity field of the fluid becomes irregular and the rate of energy dissipation remains bounded from below even if the fluid viscosity tends to zero. A mathematical description of the turbulent cascade is a very active research topic since the pioneering work of Kolmogorov in hydrodynamic turbulence and that of Zakharov in wave turbulence. In both cases, these turbulent cascade mechanisms imply power-law behaviors of several statistical quantities such as power spectral densities. For a long time, these cascades were believed to be associated with nonlinear interactions, but recent works have shown that they can also take place in a dynamics governed by a linear equation with a pseudo-differential operator of degree 0. In this spirit, we construct a linear equation that mimics the phenomenology of energy cascades when the external force is a statistically homogeneous and stationary stochastic process. In the Fourier variable, this equation can be seen as a linear transport equation, which corresponds to an operator of degree 0 in physical space. Our results give a complete characterization of the solution: it is smooth at any finite time, and, up to smaller order corrections, it converges to a fractional Gaussian field at infinite time.

## 1. INTRODUCTION

1.1. **Background and motivation.** This work is mainly motivated by some important aspects of the phenomenology of three-dimensional homogenous and isotropic fluid turbulence [35, 42, 20], of which several aspects have been also observed and formalized for waves in various situations when they are weakly interacting [43]. As has been repeatedly observed in geophysical and laboratory flows, and in numerical simulations of the incompressible Navier-Stokes equations, a fluid that is stirred by a statistically stationary random force $f(t, x)$, assumed to be smooth in space, will eventually reach a statistically stationary state in which the velocity variance is finite. To dissipate all the energy that is constantly injected into the system in such an efficient way, the velocity field of that fluid will develop a complex multiscale structure ending up with high values of velocity gradients such that viscosity can easily transform mechanical energy into heat. In other words, the fluid has transferred the energy pumped at large scales by the forcing towards small scales, at which viscous diffusion efficiently acts. This picture is known as the cascading process of energy.

The purpose of this article is to model and reproduce this phenomenon of transfer of energy as a cascading process through the scales. We propose a partial differential equation, which is of course much simpler than the nonlinear Navier-Stokes equations, stochastically forced by an additive random force $f(t, x)$ that we take to be smooth in space and correlated over a typical large lengthscale (known in the turbulence literature as the integral lengthscale), whose solution develops roughness as time goes on. More precisely, our goal is to generate rough fractional Gaussian Hölder continuous random fields of parameter $H$ (see for instance the textbook [16]) from smooth forcing through a dynamical evolution, which can be seen as a simple stochastic representation of the phenomenology mainly developed by Kolmogorov [25].
As mentioned earlier, a striking feature of three-dimensional turbulent motion is its ability to efficiently dissipate the energy that is injected at large scales in a statistically stationary and homogeneous manner. To be more precise, let us consider a solution of the incompressible Navier-Stokes equation, i.e. a divergence-free velocity field $u(t, x) \in \mathbb{R}^3$ with periodic boundary conditions. This dynamics is stirred by a divergence-free vector forcing term $f(t, x)$, that we take delta-correlated in time and smooth in space, say Gaussian, of zero-average and of covariance

$$ \mathbb{E}[f(t, x) \otimes f(s, y)] = \delta_{t-s} C_f(x - y), $$

where $\otimes$ stands for the matrix product, and the matrix $C_f(x)$ is made of a linear combination of the matrix $x \otimes x$ and the identity, with multiplicative coefficients depending only on $|x|$, i.e. a typical covariance matrix of a statistically homogeneous and isotropic vector field [6, 36]. We furthermore require that these scalar functions of $|x|$ are smooth and compactly supported over a range of the size order of the aforementioned large length scale, so as to mimic the energy injection at the so-called integral length scale. As time goes on, it has been repeatedly observed that the velocity field $u$ reaches a statistically stationary state, which is furthermore statistically homogeneous, of finite variance, and with the additional striking property that it becomes independent of the viscosity $\nu$ as $\nu$ goes to zero, i.e.

$$ (1.1) \qquad \lim_{\nu \to 0} \lim_{t \to \infty} \mathbb{E} [|u(t, x)|^2] < +\infty \quad \text{for all } x. $$

The former asymptotic behavior of the velocity variance illustrates clearly how a turbulent fluid can dissipate energy with high efficiency. For instance, in the same setup but considering the heat equation instead of the Navier-Stokes equations, a statistically stationary regime would also be reached at $t \to \infty$. However, the variance of the solution is then inversely proportional to the viscosity $\nu$, see [14]. Instead, turbulent motion dissipates energy in a way that the velocity variance is eventually independent of viscosity, which is a far more efficient way of dissipating energy. In order to ensure the independence of said variance on viscosity, (1.1), the fluid develops a rough behavior of Hölder-type at small scales, in such a way that the variance of the velocity increments asymptotically behaves as follows:

$$ (1.2) \qquad \lim_{\nu \to 0} \lim_{t \to \infty} \mathbb{E} [|u(t, x + \ell) - u(t, x)|^2] \underset{|\ell| \to 0}{\propto} |\ell|^{2H} \quad \text{for all } x, $$

where the power-law exponent is determined by Kolmogorov's prediction $H \approx 1/3$ [20]. Much more could be said on a more precise characterization of the distribution of the increments than only its variance, (1.2), such as its higher order moments that quantify its non Gaussian, skewed and intermittent nature [20]. In this article, we will focus on a second-order modeling of these fluctuations, (1.1) and (1.2); we leave finer descriptions for future research.

A first precise formalization of the cascade phenomenon could be built by imposing a particular dynamical relation between the coefficients of a decomposition of the velocity field, such as a continuous wavelet transform or a discrete (dyadic) decomposition on a tree [18]. This has been explored in the literature [5, 4, 13] leading to precise statements on Hölder regularity and its relationship with scaling behaviors of the coefficients. Although great progress has been made in the understanding of such models and their formalization, which usually exploits a typical quadratic interaction between neighboring coefficients, these approaches often avoid the important question of the relation of these coefficients in space. This is necessary in order to design a model that leads to statistically homogeneous velocity fields, as observed in nature and in numerical simulations. Nonetheless, these models can be seen as a sophistication of the so-called shell models<fnref n="1" />. In this spirit, we believe an important step was made in [33], where the authors investigate a simple linear

<footnote n="1">¹See for instance [10, 9] which consist in exploring quadratic interactions between shells, that share some behaviors with velocity Fourier modes and wavelet coefficients, along a single branch of a tree decomposition, lacking thus a discussion of the spatial relationships between coefficients.</footnote>
relation between shells, which is shown to be able to transfer energy from large to small scales. Let us also mention [31] where some ideas to build a PDE from these shell models are proposed.

From a somewhat different side of fluid mechanics, more focused on the implications of global rotation [40, 39] or stratification of the density field [30, 41] on a flow, it has been evidenced a phenomenon of focusing of waves onto attractors, whose precise shape are determined by the boundaries. Based on a linearization of the fluid equations, this phenomenon has been interpreted as a cascading process through scales. These ideas have been then formalized and rigorously studied from a mathematical viewpoint in a series of recent articles [17, 19], which underline the importance of operators of degree 0 as a deterministic mechanism able to transfer energy through scales.

## 1.2. Main results.
The rough and disordered nature of a turbulent velocity field $u(t, x)$ has been repeatedly observed in laboratory and numerical flows, and in geophysical situations [42, 20]. From this signal, considering for instance a component of the velocity vector field as a function of space, depending on the experimental possibilities and the large-scale geometry of the flows, one can construct the energy spectrum $|k| \mapsto \mathbb{E}|\hat{u}(t, k)|^2$ where $\hat{u}$ stands for the spatial Fourier transform. According to the standard phenomenology of fluid turbulence, which has been multiply confirmed by observations in very different situations, the energy-spectrum resembles a curve [42, 20] that can be schematically decomposed as follows:

* (injection range) for small $|k|$ of the order of the characteristic wavelength of energy injection, the energy-spectrum is mainly determined by the forcing and the associated large-scale geometry of the flow,

* (inertial range) for intermediate $|k|$, the energy-spectrum develops a power-law behavior whose exponent is found universal, i.e. independent of viscosity and of the nature of the flow, and can be interpreted as the generation of small scales by the internal motion of the fluid following a transfer of energy from small wave-numbers to large wave-numbers,

* (dissipative range) for large $|k|$, the energy-spectrum is governed by dissipation processes which damp efficiently all the energy coming from the large scales, making the spatial velocity profile a smooth function.

The intermediate range of scales, called the inertial range in the turbulence literature [42, 20], is where this mechanism of transport of energy takes place. The universally observed power-law exponent of the energy-spectrum can be written as $-(2H + d)$, i.e. $\mathbb{E}|\hat{u}(t, k)|^2 \sim |k|^{-(2H+d)}$, where we have introduced for the sake of generality the space dimension $d$, and the parameter $H$ that will be eventually interpreted as a Hurst, or Hölder, exponent, in a statistically averaged sense. In real situations, for $d = 3$, it is indeed observed that $H \approx 1/3$, as predicted by dimensional arguments mainly attributed to Kolmogorov [25, 28].

The main goal of this paper is to propose a family of partial differential equations, such that, when stirred by a statistically stationary and spatially homogenous, smooth in space forcing term, its solution $u(t, x)$ reaches at long times a statistically stationary state which displays the typical spectral behavior detailed above. We will achieve this with the following transport equation in Fourier space:

$$ (1.3) \qquad \begin{cases} \partial_t \hat{u}(t, k) + \text{div}_k \left( \frac{ck}{|k|} \hat{u}(t, k) \right) + c \frac{H + \frac{1}{2}}{|k|} \hat{u}(t, k) = \hat{f}(t, k) & t > 0, k \in \mathbb{R}^d, |k| > \kappa > 0, \\ \hat{u}(t, k) = 0 & t > 0, k \in \mathbb{R}^d, |k| \le \kappa, \\ \hat{u}(0, k) = 0. & \end{cases} $$

Here $c, \kappa > 0$ and $H \in \mathbb{R}$ are fixed, and the source $f$ satisfies

$$ \mathbb{E}[f(t, x)f(s, y)] = \delta_{t-s} C_f(x - y), $$
where $C_f$ is smooth and satisfies some additional assumptions detailed below. Our main result is the following:

**Theorem 1.1.** Let $H \in (0, 1)$ and let the forcing $f$ be

$$ (1.4) \qquad f(t, x) = \int_{\mathbb{R}_y^d} \varphi(x - y) dW(t, y), $$

where $dW$ is a space-time Gaussian real white noise and $\varphi \in S(\mathbb{R}_x^d)$ is a radial function such that $\hat{\varphi}(k) = 0$ for all $|k| < \kappa$.

(i) The transport equation in wavenumber space (1.3) with source (1.4) can be rigorously formulated in physical space as an a.s. well-posed PDE. Moreover, at any $t > 0$, the solution $u(t, x)$ has finite variance and a.s. smooth paths with respect to $x$.

(ii) As $t \to \infty$, $u(t, x)$ converges in law to a zero-mean Gaussian field $u_\infty(x)$ which has a.s. $\alpha$-Hölder continuous paths for any $0 < \alpha < H$.

(iii) The correlations are given by

$$ \mathbb{E}[u_{\infty}(x_1)u_{\infty}(x_2)] = C(d, H)K_H(x_1 - x_2) - (K_H * \mathcal{J}_H)(x_1 - x_2), $$

where

$$ K_H := \mathcal{F}^{-1} \left[ \chi_{|k|>\kappa} |k|^{-(2H+d)} \right], $$

while $C(d, H)$ is an explicit constant and the function $\mathcal{J}_H \in S(\mathbb{R}_x^d)$ depends explicitly on $\varphi$ in (1.4).

A more detailed version of this result is presented in Theorem 4.4, page 23.

**Remark 1.2.** The parameter $\kappa$ can be chosen as the smallest non-vanishing wavenumber in the support of the Fourier transform of the forcing. Following the analogy with the Navier-Stokes equations presented in the introduction, $\kappa$ may be interpreted as a quantity linked to the inverse of the integral lengthscale<fnref n="2" />, i.e the typical lengthscale of the correlations of the forcing. The fact that our force acts at large but finite scales means that $\kappa$ is small but non-zero. The kernel $K_H$ is a function when $H \in (0, 1)$ and $\kappa > 0$. However, when $\kappa = 0$ we have the following operator:

$$ K_H \xrightarrow{\kappa \to 0} (-\Delta)^{-(H+\frac{d}{2})}. $$

**Remark 1.3.** Note that the limiting Gaussian field $u_\infty$ shares some statistical properties, such as roughness, with statistical homogeneous fractional gaussian fields [16, 29] defined by

$$ (-\Delta)^{-\frac{H+d}{2}} dW, $$

that are classically encountered in the turbulence literature [26, 12, 15]. Indeed, for $H \in (0, 1)$ both have a.s. $\alpha$-Hölder continuous paths for any $0 < \alpha < H$ and one can show that

$$ (1.5) \qquad u_{\infty} = C(d, H) \mathcal{F}^{-1} \left[ \chi_{|k|>\kappa} \right] * (-\Delta)^{-\frac{H+d}{2}} dW - u_{\text{reg}}. $$

Here, $u_{\text{reg}}$ is a smooth zero-mean Gaussian field with correlations

$$ \mathbb{E}[u_{\text{reg}}(x_1)u_{\text{reg}}(x_2)] = (K_H * \mathcal{J}_H)(x_1 - x_2), $$

which is a smooth function with respect to $(x_1 - x_2)$ even if $H \le 0$. The case $H \in [-d/2, 0]$ will be discussed in Section 4.

<footnote n="2">²If $L_f$ is the integral lengthscale, then there exists two real positive numbers $a < b$ such that the support of $\hat{f}$ is contained in the annulus of inner radius $a/L_f$ and outer radius $b/L_f$. Thus one may set $\kappa = a/L_f$.</footnote>
**Remark 1.4.** When $H \in [-d/2, 0]$, as $t \to \infty$, $u(t)$ still converges to a zero-mean Gaussian field $u_\infty$, but this field is not necessarily Hölder continuous with respect to $x$ anymore. In this case, one may view $u_\infty$ as a distribution living in the dual of an appropriate test function space $T$, see Theorem 4.5 for more details. The correlation structure of the limiting Gaussian measure is given by:

$$ (1.6) \qquad \mathbb{E}[\langle u_\infty, g_1 \rangle \langle u_\infty, g_2 \rangle] = \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} (C(d, H) - \hat{\mathcal{J}}_H(k)) \hat{g}_1(k) \hat{g}_2(k) dk $$

for any test functions $g_1, g_2 \in T$, where $\langle \cdot, \cdot \rangle$ stands for the duality product in $T$. In Theorem 4.5, we show that, for any test functions in $T \cap S(\mathbb{R}_d^x)$, (1.6) yields a rate of convergence proportional to $(ct)^{-(2H+d+2n)}$ for $n$ as large as desired. The expression (1.6) corresponds to the energy-spectrum picture described above. Indeed, $|k|^{-(2H+d)}$ corresponds precisely to the inertial range previously described, while $\hat{\mathcal{J}}_H(k)$ captures the contribution from the forcing, which is a correction in the injection range. In fact, if the source is spectrally supported in small wavenumbers, then $\hat{\mathcal{J}}_H(k)$ vanishes in the inertial range.

Finally, let us highlight the difference between the properties of the solution at finite and infinite time. At finite time, the solution is smooth with respect to $x$, whereas at infinite time the solution is only Hölder continuous (or even rougher if $H \le 0$, as explained in Remark 1.4). This loss of regularity at infinite time is what is expected in linear turbulence. Turbulence is usually associated to a nonlinear equation. For example, in the case of wave turbulence, nonlinearities create wave interactions which allow the transfer of energy to higher and higher wavenumbers. Such transfers of energy typically result in a loss of regularity. However, nonlinearities might not be the only way in which such loss of regularity can occur. Indeed, Y. Colin de Verdière and L. Saint-Raymond [17] have shown that, in the context of internal waves, a loss of regularity can also take place in the case of a linear equation with an operator of degree 0. Linear equations with operators of degree 0 are also common whenever one introduces a dispersive perturbation in a hyperbolic system. In such cases, these operators of degree 0 are used to model wave propagation under strong dispersive effects and they are responsible for memory effects. For example, in the context of wave-energies, the second author and D. Lannes show that the waves generated by a moving floating object are governed in the linear regime by a non-local transport equation of degree 0, see [8]. In the context of electrical circuits, there are cases in which 1D models of electromagnetic waves propagating along a coaxial cable are governed by operators of degree 0, see for instance [7, Chapter 5].

One issue of our model is that it only features a single Hölder exponent. The velocity field of a concrete turbulent fluid consists of many Hölder exponents, i.e. the Hölder-regularity of the velocity field $u(t, x)$ around a point $x \in \mathbb{R}^d$ depends on the point itself. This is known as the multifractal formalism [14, 16]. The term multifractal refers to the fact that the sets of points with same regularity are often fractal. Moreover, our model does not capture finer descriptions (beyond the variance) of the distribution of the increments of the velocity field. Such descriptions should quantify its non-Gaussian and intermittent nature [15], and therefore our linear model does not suffice. It is known that one can construct a multifractal and intermittent field with the theory of Gaussian multiplicative chaos [38], however our actual goal is to obtain a multifractal and intermittent field dynamically, i.e. as the solution to a non-linear equation forced by a white-noise in time, but that admits a rigorous mathematical treatment. We can also consider a forcing which is not a white-noise in time whose temporal correlation function is given by an oscillating function in order to make a comparison with [17, 21, 1] and [11]. Finally, we could investigate other linear models of cascade such as in the case of a compact operator plus a potential of degree 0 as in [32]. These issues will be tackled in future papers.
1.3. **Outline.** The article is organized as follows. In Section 2, we present a simple transport equation that converges to a complex white noise (up to lower order terms). A small tweak to this model allows us to construct a model that gives rise to a real white noise. In Section 3, we explain how to generalize the latter model to higher dimensions and give a heuristic proof of the main results in this paper. In Section 4, we provide a mathematically rigorous study of our model: we introduce the right functional setting, we develop a global well-posedness theory and give a complete description of the asymptotic behavior of the solution, as well as its properties. This constitutes the proof of Theorem 1.1. Finally, in Section 5 we propose a numerical method and conduct numerical simulations in dimensions 1, 2 and 3 to validate our model.

1.4. Acknowledgements. We thank Oliver Bühler for his interesting suggestion about replacing the pure transport term ∂_k in (2.1) by sign(k)∂_k as a way of fixing the physically undesirable behavior of the solution. We thank Charles-Edouard Bréhier and the referee for their careful reading and their corrections of a previous version of this manuscript. All five authors were funded by the Simons Collaboration Grant on Wave Turbulence, Simons Award ID: 651475 and 651675.

1.5. Notation. For an integrable function $f: \mathbb{R}^d \to \mathbb{C}$, we denote by $\hat{f}$ its **Fourier transform**, namely

$$ \forall k \in \mathbb{R}_k^d, \quad \hat{f}(k) := \mathcal{F}f(k) := \int_{\mathbb{R}_k^d} e^{-2\pi i x \cdot k} f(x) dx. $$

Whenever defined, the inverse Fourier transform is

$$ \forall x \in \mathbb{R}_x^d, \quad f(x) = \mathcal{F}^{-1}\hat{f}(x) = \int_{\mathbb{R}_k^d} e^{2\pi i x \cdot k} \hat{f}(k) dk. $$

It is well known that the Fourier transform is an isometry from $L^2(\mathbb{R}_x^d)$ to $L^2(\mathbb{R}_k^d)$, from $S(\mathbb{R}_x^d)$ to $S(\mathbb{R}_k^d)$, where $S(\mathbb{R}_x^d)$ denote the space of **Schwartz functions** (i.e. smooth functions whose derivatives are rapidly decreasing), and from $S'(\mathbb{R}_x^d)$ to $S'(\mathbb{R}_k^d)$ where $S'(\mathbb{R}_x^d)$ denote the space of tempered distribution (i.e. the dual space of $S(\mathbb{R}_x^d)$). We will denote by $\langle \cdot, \cdot \rangle$ the duality product between $S'$ and $S$.

We will also need some spaces that quantify the regularity of functions more precisely. For a fixed integer $n$, the **Sobolev space** of order $n$ is defined by

$$ H^n(\mathbb{R}^d) := \{u \in L^2(\mathbb{R}^d) \mid \partial_{x_i}^j u \in L^2(\mathbb{R}^d) \text{ with } 1 \le i \le d \text{ and } 0 \le j \le n\}, $$

and their dual spaces are denoted by $H^{-n}(\mathbb{R}^d)$. We will denote by $\langle \cdot, \cdot \rangle_{H^{-n}, H^n}$ the duality product between $H^n$ and $H^{-n}$. For $\alpha \in (0, 1)$ the **Hölder space** $C^{0,\alpha}(\mathbb{R}^d)$ is defined by

$$ C^{0,\alpha}(\mathbb{R}^d) := \{u \text{ continuous and bounded } |\exists C > 0, \forall x, \ell \in \mathbb{R}^d, |\ell| \le 1, |\delta_\ell u(x)| \le C|\ell|^\alpha\}, $$

where $\delta_\ell$ denotes the **increment** defined by

$$ \delta_{\ell}u(x) := u(x + \ell) - u(x) \quad \text{for } x, \ell \in \mathbb{R}^d. $$

We will denote by $\chi_A$ the **characteristic function**$^3$ of the set $A$.

Let $(\Omega, \sigma(\Omega), \mathbb{P})$ be a probability space. A **Gaussian field** $u: \mathbb{R}^d \to L^2(\Omega)$ is a field such that for all $n \ge 1$ and for all $(x_1, x_2, \dots, x_n) \in (\mathbb{R}^d)^n$, the random vector $(u(x_1), u(x_2), \dots, u(x_n))$ is a Gaussian random vector. When $d = 1$, a 1D Gaussian field is usually called a Gaussian process. A **Gaussian random measure** $\mu$ acting on $S(\mathbb{R}^d)$ is a random tempered distribution such that for every $g \in S(\mathbb{R}^d)$, the random variable $\langle \mu, g \rangle$ is a centered Gaussian. A **white noise**

<footnote>$^3$This means that $\chi_A(k) = 1$ if $k \in A$ and $\chi_A(k) = 0$ if $k \notin A$.</footnote>
$dW(x)$ is a Gaussian random measure acting on $L^2(\mathbb{R}^d)$ that satisfies the following for any functions $f, g \in L^2(\mathbb{R}^d)$

$$ \mathbb{E} \left[ \left( \int_{\mathbb{R}^d} f(x) dW(x) \right) \left( \int_{\mathbb{R}^d} g(x) dW(x) \right) \right] = \int_{\mathbb{R}^d} f(x) g(x) dx. $$

Since we always integrate a deterministic function against $dW(x)$, the choice between Itô and Stratonovich integrals is unimportant, since Wiener integration suffices (see, for instance, [22]).

## 2. ONE-DIMENSIONAL TRANSPORT IN WAVENUMBER SPACE

### 2.1. Building one-dimensional white noise: real vs complex.
In order to mimic the transport of energy from large scales to small scales, the authors in [3, 2] proposed a simple transport equation in Fourier space. To present these ideas, we first consider a one-dimensional model for a velocity field $u(t, x)$, whose spatial Fourier transform aims to solve the linear evolution

$$ (2.1) \qquad \begin{cases} \partial_t \hat{u}(t, k) + c \partial_k \hat{u}(t, k) = \hat{f}(t, k), & (t, k) \in (0, \infty) \times \mathbb{R}, \\ \hat{u}(t, k)|_{t=0} = 0. \end{cases} $$

Here $c > 0$ is fixed and can be viewed as a transport rate in wavenumber space. On the right-hand side of (2.1), we have included an additive term $\hat{f}$ which is the Fourier transform of a spatial forcing term. The support of $\hat{f}$ is localized at small wavenumbers, which is consistent with the assumption that the forcing term in physical space acts at large scales.

As we can see, the dynamical evolution proposed in (2.1) is a genuine transport equation, and only the presence of a forcing makes it inhomogeneous. Adopting such a setup immediately imposes the complex nature of the velocity field in physical space, as it can be seen when formally taking the inverse Fourier transform of (2.1), and obtaining the following evolution in physical space:

$$ (2.2) \qquad \begin{cases} \partial_t u(t, x) - 2\pi i c x u(t, x) = f(t, x), & (t, x) \in (0, \infty) \times \mathbb{R}, \\ u(t, x)|_{t=0} = 0. \end{cases} $$

Note that the operator in (2.2) corresponds to multiplication by the space variable $2\pi i c x$. In [3], it is shown that when the forcing $f$ is a white noise in time and statistically homogeneous in space, then the solution to (2.2), $u: (0, \infty)_t \times \mathbb{R}_x \to \mathbb{C}$, converges<fnref n="4" /> to a complex white noise in space as $t \to \infty$. In other words, the evolution that has been proposed, expressed in Fourier space as a transport equation (2.1) and in physical space as an equation (2.2) involving an operator of degree 0 (multiplication by $-2\pi i c x$), is able to transfer energy through scales. Moreover, the solution is statistically homogeneous at any time and it develops the regularity of a white noise as time goes on (technically it's a sudden drop in regularity at $t = \infty$). To complete the program suggested by the phenomenology of turbulence, the additional linear action of a fractional operator allows, in a similar setup, to generate a solution with asymptotic Hölder-type regularity of parameter $H \in (0, 1)$ instead of the one of the white noise, as explained in [3].

If some energy is introduced by the forcing at a negative wavelength $k < 0$, the transport equation (2.1) will move it to smaller negative wavenumbers, going through $k = 0$, and then to infinitely large positive wavenumbers $k \to \infty$. In order to avoid this pathological behavior, one could replace $\partial_k$ by $\partial_{|k|} = \text{sign}(k)\partial_k$ which leads to a transport in the direction of $|k|$ instead of $k$. Note however that $\text{sign}(k)\partial_k$ is not properly defined at $k = 0$, and so one needs to be careful in order to propose a well-posed mathematical problem. With this in mind, the heart of this article

<footnote n="4">⁴Up to lower order terms, see Theorem 2.4 for a full asymptotic expansion.</footnote>
will be the theoretical and numerical study of the following formal evolution:

$$ (2.3) \qquad \begin{cases} \partial_t \hat{u}(t, k) + c \partial_{|k|} \hat{u}(t, k) = \hat{f}(t, k), \quad (t, k) \in (0, \infty) \times \mathbb{R}, \\ \hat{u}(t, k)|_{t=0} = 0, \\ \hat{u}(t, k)|_{|k|=0} = 0. \end{cases} $$

Note that it is necessary to introduce a transmission condition between negative and positive $k$. In (2.3), we have decided to add the boundary condition $\hat{u}(t, k)|_{|k|=0} = 0$ to decouple negative from positive wavenumbers, so that no energy crosses $k = 0$. In particular, this means that the integral over space of $u$ is zero for all times. This new dynamics proposed in (2.3) can be written in physical space after formally applying the inverse Fourier transform:

$$ (2.4) \qquad \left\{ \begin{aligned} & \partial_t u + 2\pi c x \mathcal{H}(u) = f, \quad (t, x) \in (0, \infty) \times \mathbb{R}, \\ & u|_{t=0} = 0, \\ & \int_{\mathbb{R}_x} u \, dx = 0, \end{aligned} \right. $$

where $\mathcal{H}$ denotes the Hilbert transform defined in the usual way:

$$ (2.5) \qquad \mathcal{H}f(x) := \frac{1}{\pi} \text{p.v.} \int_{\mathbb{R}_y} \frac{f(y)}{x-y} dy = -\frac{1}{\pi} \lim_{\epsilon \to 0^+} \int_{\epsilon}^{\infty} \frac{f(x+y) - f(x-y)}{y} dy. $$

Notice that we have used the fact the integral of $u(t, x)$ over space vanishes to get the expression of (2.4). Notice also that, despite the fact that the spectral evolutions (2.1) and (2.3) (whose equivalent expressions in physical space are provided respectively in (2.2) and (2.4)), look very similar, the solution to the new dynamics (2.4) is now real-valued. Equivalently, the dynamics in Fourier space (2.3) conserves the Hermitian symmetry of an appropriate initial condition, here assumed to be zero. Moreover, the solution $u : (0, \infty)_t \times \mathbb{R}_x \to \mathbb{R}$ of (2.4) can be shown to asymptotically converge to a real white noise in space. As we will explain in the sequel, the additional linear action of a fractional operator will allow the generation, from smooth forcing, of a real fractional Gaussian field.

Finally, it is tempting to generalize (2.3) to higher dimensions by replacing $\partial_{|k|}$ by $\frac{k}{|k|} \nabla_k$. As we will develop in Section 4, this eventually generates a statistically homogeneous and isotropic solution that will converge to a real $d$-dimensional Gaussian random measure which is rougher than a white noise (in space) whenever $d > 1$. As we will explain in the sequel, the additional linear action of a fractional operator will allow us to generate a real $d$-dimensional Gaussian random measure with the desired Hölder regularity. Interestingly, it is not obvious to generalize (2.1) to space dimension $d \ge 1$ which would generate a similar statistically homogeneous and isotropic solution in physical space. We provide at the end of the section some additional discussions on this matter.

However, for the time being we focus on developing a good understanding in the one-dimensional setting. In the case of (2.4), we have the following result:

**Theorem 2.1.** Let the forcing $f$ in (2.4) be

$$ (2.6) \qquad f(t, x) = \int_{\mathbb{R}_y} \varphi(x - y) dW(t, y), $$

where $dW$ is a space-time Gaussian real white noise, and $\varphi \in S(\mathbb{R})$ is a non-negative, non-identically null, even function with null average. Then:

(i) Equation (2.4) admits a global (in time) solution $u(t, x)$, which is a Gaussian process with a.s. smooth paths in $x$, and with $\alpha$-Hölder continuous paths in $t$ for any $0 < \alpha < 1/2$.
(ii) As $t \to \infty$, the solution $u(t)$ converges in $S'(\mathbb{R})$ to a random Gaussian measure $u_\infty$ acting on $S(\mathbb{R})$ with zero-mean, i.e. for every $g \in S(\mathbb{R})$, $\langle u(t), g \rangle$ converges in law to a complex Gaussian random variable $\langle u_\infty, g \rangle$ with $\mathbb{E}[\langle u_\infty, g \rangle] = 0$.

(iii) We have the following asymptotic behavior:

$$ (2.7) \qquad \begin{aligned} \mathbb{E}[\langle u_\infty, g_1 \rangle \langle u_\infty, g_2 \rangle] &= \lim_{t \to \infty} \mathbb{E}[\langle u(t), g_1 \rangle \langle u(t), g_2 \rangle] \\ &= C \int_{\mathbb{R}_x} g_1(x) g_2(x) dx - \int_{\mathbb{R}_x \times \mathbb{R}_y} \mathcal{I}(x-y) g_1(x) g_2(y) dx dy. \end{aligned} $$

for any $g_1, g_2 \in S(\mathbb{R})$. Here $C > 0$ is a constant and $\mathcal{I}$ is an explicit continuous, even function that depends on $\varphi$ in (2.6).

A more detailed version of this result is presented in Theorem 4.5, see also Section 4.4. Note that the first term on the right-hand side of (2.7) corresponds to a delta function, while the second term given by $\mathcal{I}$ is a smooth lower order term.

**Remark 2.2.** The forcing introduced in (2.6) is indeed a Gaussian white noise in time and statistically homogeneous in space, i.e.

$$ (2.8) \qquad \mathbb{E}[f(s, x)f(t, y)] = C_f(x - y) \delta_{s-t}, $$

where the spatial correlation function $C_f = \varphi * \varphi$ is a convolution. Given that $\varphi \in S(\mathbb{R})$, $C_f \in S(\mathbb{R})$. The constant $C$ in (2.7) is precisely

$$ C = \frac{C_f(0)}{2c} = \frac{1}{2c} \int_{\mathbb{R}_x} |\varphi(x)|^2 dx > 0. $$

**Remark 2.3.** The solution to the stochastic PDE (2.4) is an explicit Gaussian Itô process. A precise formula will be given in Section 4. Even if this solution is continuous in time and space, one can lose regularity at $t = +\infty$, which is why one needs to consider $u_\infty$ on the left-hand side of (2.7) as a distribution.

The proof of this theorem is posponed to the next section where a more general case, i.e. multidimesional white noise, will be tackled. Before we develop the techniques needed to prove Theorem 2.1, it is important to understand the asymptotic behavior of solutions to (2.2) in the complex setting, which is less technical, but informative.

In this setting, we have the following result:

**Theorem 2.4.** Let the forcing $f$ in (2.2) be

$$ (2.9) \qquad f(t, x) = \int_{\mathbb{R}_y} \varphi(x - y) dW(t, y), $$

where $dW$ is a space-time Gaussian complex white noise, and $\varphi \in S(\mathbb{R})$ is a complex, non-identically null, even function. Then:

(i) Equation (2.2) admits a global (in time) solution $u(t, x)$, which is a Gaussian process with a.s. continuous paths in time and space.

(ii) As $t \to \infty$, the solution $u(t)$ converges (in $S'(\mathbb{R}^d)$) to a random Gaussian measure $u_\infty$ acting on $S(\mathbb{R}^d)$, i.e for every $g \in S(\mathbb{R})$, $\langle u(t), g \rangle$ converges in law to $\langle u_\infty, g \rangle$ where $u_\infty$ is a random Gaussian measure with $\mathbb{E}[\langle u_\infty, g \rangle] = 0$.
(iii) We have the following asymptotic behavior (in the sense of distributions). For any $g_1, g_2 \in S(\mathbb{R})$,

$$ (2.10) \qquad \begin{aligned} \mathbb{E}[\langle u_{\infty}, g_1 \rangle \overline{\langle u_{\infty}, g_2 \rangle}] &= \lim_{t \to \infty} \mathbb{E}[\langle u(t), g_1 \rangle \overline{\langle u(t), g_2 \rangle}] \\ &= \frac{1}{2c} C_f(0) \int_{\mathbb{R}_z} g_1(z)g_2(z)dz \\ &\quad + \frac{1}{2\pi ic} \text{p.v.} \int_{\mathbb{R}_z} \frac{C_f(z)}{z} \left( \int_{\mathbb{R}_y} g_1(z+y)g_2(y)dy \right) dz. \end{aligned} $$

The function $C_f = \varphi * \bar{\varphi}$ is the spatial correlation function given by

$$ (2.11) \qquad \mathbb{E}[f(s, x) \overline{f(t, y)}] = \delta_{s-t} C_f(x - y), $$

and p.v. $\frac{C_f(z)}{z}$ is the principal value of the distribution $C_f(z)/z$. That is, for any test function $g \in S(\mathbb{R})$,

$$ (2.12) \qquad \left\langle \text{p.v.} \frac{C_f(z)}{z}, g \right\rangle := \text{p.v.} \int_{\mathbb{R}} \frac{C_f(z)g(z)}{z} dz = \int_{0}^{\infty} C_f(z) \frac{g(z) - g(-z)}{z} dz. $$

As we mentioned in Remark 2.2, the forcing in (2.9) is a complex Gaussian white noise in time and statistically homogeneous in space. Moreover, the solution admits an explicit formula:

$$ (2.13) \qquad u(t, x) = \int_{0}^{t} \int_{\mathbb{R}_{y}} e^{2\pi i c x (t-s)} \varphi(x - y) dW(s, y). $$

The same comments as in Remark 2.3 apply in this case.

**Remark 2.5.** The asymptotic expansion (2.10) remains valid when testing against functions with a finite number of derivatives. Indeed, $u_{\infty}$ in (2.10) can also be interpreted as a Gaussian random measure in $H^{-n}(\mathbb{R})$ for any integer $n \ge 2$. More precisely, we will show that for any test functions $g_1, g_2 \in H^n(\mathbb{R})$ one gets

$$ (2.14) \qquad \begin{aligned} \mathbb{E}[\langle u(t), g_1 \rangle_{H^{-n}, H^n} \overline{\langle u(t), g_2 \rangle_{H^{-n}, H^n}}] &\sim \frac{C_f(0)}{2c} \int_{\mathbb{R}} g_1(z)g_2(z) dz \\ &\quad + \frac{1}{2\pi ic} \text{p.v.} \int_{\mathbb{R}_z} \frac{C_f(z)}{z} \left( \int_{\mathbb{R}_y} g_1(z+y)g_2(y)dy \right) dz \\ &\quad + \text{d}(n) \|g_1\|_{H^n} \|g_2\|_{H^n} \left( \frac{1}{ct} \right)^{n-1}. \end{aligned} $$

where $\text{d}(n)$ depends only on $n$. This characterizes the rate of convergence.

**Remark 2.6.** One interpretation of (2.7) (resp. (2.10)) is that the correlation function (resp. the real part of it) asymptotically behaves like a white noise. However, this theorem gives a lower-order correction, in the sense that the regularity of the correction is higher than that of the white noise. In (2.7), such a regular correction is given by a Schwartz function, whereas in (2.10), the regular correction is a purely imaginary principal value which has no singularity at zero: by (2.12), the principal value is “controlled” by $C'_f(0)$ near zero. It is important to note that in both cases the regular correction is fast-decaying.

**Remark 2.7.** Note that we recover the result in proposition 2.1 in [3], i.e.

$$ (2.15) \qquad \lim_{t \to \infty} \mathbb{E}[u(t, x) \overline{u(t, y)}] = \frac{1}{2c} C_f(0) \delta_{x-y} $$
as long as one only tests against even functions with respect to the variable $x - y$, as is easily seen from the right-hand side of (2.12).

One way to recover a result similar to that in [3] that holds for all test functions is to define the function $v(t, x) = e^{-\pi i c t} u(t, x)$. This function now satisfies

$$ \mathbb{E}[v(t, x) \overline{v(t, y)}] = t \sinh(ct(x - y)) C_f(x - y) $$

where $\sinh(x) := \frac{\sin(\pi x)}{\pi x}$ denotes the normalized sinc function. It immediately follows that

$$ (2.16) \qquad \lim_{t \to \infty} \mathbb{E}[v(t, x) \overline{v(t, y)}] = \frac{C_f(0)}{c} \delta_{x-y}. $$

However, it is unclear whether this transformation of $u$ is an interesting object from the physical viewpoint. Assuming one can take the Fourier transform, (2.16) can be rewritten in wavenumber space as

$$ (2.17) \qquad \lim_{t \to \infty} \mathbb{E}[\hat{v}(t, k) \overline{\hat{v}(t, k')}] = \lim_{t \to \infty} \mathbb{E} \left[ \hat{u}(t, k + \pi ct) \overline{\hat{u}(t, k' + \pi ct)} \right] = \delta_{k-k'} C_f(0). $$

The transformation given by $v$ is therefore equivalent to computing the correlation between the $k+\pi ct$ and $k' + \pi ct$ Fourier modes as $t \to \infty$.

2.2. **Proof of Theorem 2.4.** First of all, note that equation (2.2) admits the explicit solution (2.13) thanks to the Duhamel formula.

**Step 1:** The solution $u$ is a well defined Gaussian field whose limit at $t \to \infty$ is a Gaussian random measure.

Clearly, $u(t, x)$ is a well defined Itô process with zero average and variance

$$ (2.18) \qquad \mathbb{E}[|u(t, x)|^2] = \int_0^t \int_{\mathbb{R}_y} |\varphi(x - y)|^2 dy ds = t \|\varphi\|_{L^2}^2. $$

For any test function $g$ (we will soon see that actually $g \in H^1(\mathbb{R})$ suffices), one gets

$$ \langle u(t), g \rangle = \int_0^t \int_{\mathbb{R}_y} \left( \int_{\mathbb{R}_x} e^{2\pi i c x (t-s)} \varphi(x-y) g(x) dx \right) dW(s, y). $$

For fixed $t$, this random variable has the same distribution as

$$ (2.19) \qquad \langle v(t), g \rangle := \int_0^t \int_{\mathbb{R}_y} \left( \int_{\mathbb{R}_x} e^{-2\pi i c x s} \varphi(x-y) g(x) dx \right) dW(s, y). $$

As $t \to \infty$, one has the mean-square convergence property

$$ (2.20) \qquad \lim_{t \to \infty} \langle v(t), g \rangle = \langle v_\infty, g \rangle = \int_0^\infty G(s, y) dW(s, y) $$

where

$$ G(s, y) = \int_{\mathbb{R}_x} e^{-2\pi i c x s} \varphi(x - y) g(x) dx. $$

In order to justify the convergence result above, it is sufficient to assume that $G \in L^2(\mathbb{R}_s^+ \times \mathbb{R}_y)$, which we set out to prove next. Firstly, the Young convolution inequality immediately yields

$$ \|G\|_{L^2([0,1]_s \times \mathbb{R}_y)} \le \|\varphi\|_{L^1} \|g\|_{L^2}. $$

Next, to handle the case $|s| \ge 1$, we assume that $g \in H^1(\mathbb{R})$ and we integrate by parts:

$$ 2\pi i c s G(s, y) = - \int_{\mathbb{R}_x} e^{-2\pi i c x s} \partial_x(\varphi(x - y) g(x)) dx. $$
By the Young convolution inequality,

$$ \|2\pi i c s G(s, y)\|_{L_y^2} \le \|g'\|_{L^2} \|\varphi\|_{L^1} + \|g\|_{L^2} \|\varphi'\|_{L^1}. $$

Thus one easily finds that:

$$ \|G(s, y)\|_{L^2(\mathbb{R}_s^+ \times \mathbb{R}_y)} \lesssim \frac{1}{c} \|g\|_{H^1} \|\varphi\|_{W^{1,1}}, $$

which justifies (2.20).

Since the mean and the variance of $\{\langle v(t), g \rangle\}_{t>0}$ converge and since they coincide with those of $\{\langle u(t), g \rangle\}_{t>0}$ for each fixed $t$, we deduce that $\langle u(t), g \rangle$ converges in law to a Gaussian random variable

$$ (2.21) \qquad \langle u_{\infty}, g \rangle_{\text{law}} = \int_{0}^{\infty} G(s, y) dW(s, y). $$

Moreover, all the moments of $\{\langle u(t), g \rangle\}_{t>0}$ converge to those of $\langle u_{\infty}, g \rangle$.

**Step 2:** We would like to show that

$$ (2.22) \qquad \mathbb{E}[\langle u_{\infty}, g_1 \rangle \overline{\langle u_{\infty}, g_2 \rangle}] = \lim_{t \to \infty} \mathbb{E}[\langle u(t), g_1 \rangle \overline{\langle u(t), g_2 \rangle}]. $$

By (2.20) and (2.21),

$$ \mathbb{E}[\langle u_{\infty}, g_1 \rangle \overline{\langle u_{\infty}, g_2 \rangle}] = \mathbb{E}[\langle v_{\infty}, g_1 \rangle \overline{\langle v_{\infty}, g_2 \rangle}] $$

for $v$ defined in (2.19).

Since $\mathbb{E}[\langle v(t), g_1 \rangle \overline{\langle v(t), g_2 \rangle}] = \mathbb{E}[\langle u(t), g_1 \rangle \overline{\langle u(t), g_2 \rangle}]$ for each fixed $t$, our objective (2.22) follows from the identity

$$ (2.23) \qquad \mathbb{E}[\langle v_{\infty}, g_1 \rangle \overline{\langle v_{\infty}, g_2 \rangle}] = \lim_{t \to \infty} \mathbb{E}[\langle v(t), g_1 \rangle \overline{\langle v(t), g_2 \rangle}]. $$

In order to justify (2.23), we may simply use the mean-square convergence of $\langle v(t), g \rangle$ to $\langle v_{\infty}, g \rangle$ for any $g \in S(\mathbb{R}^d)$, together with the following “polarization” identity:

$$ \langle v(t), g_1 \rangle \overline{\langle v(t), g_2 \rangle} = \frac{1}{4} \left( |\langle v(t), g_1 + g_2 \rangle|^2 - |\langle v(t), g_1 - g_2 \rangle|^2 \right) \\ + \frac{i}{4} \left( |\langle v(t), g_1 + ig_2 \rangle|^2 - |\langle v(t), g_1 - ig_2 \rangle|^2 \right). $$

**Step 3:** Calculation of the correlations as $t \to \infty$.

We start by computing the correlations for a finite time $t > 0$.

$$ \mathbb{E}[u(t, x_1) \overline{u(t, x_2)}] = \int_0^t e^{2\pi i c(x_1 - x_2)s} C_f(x_1 - x_2) ds = \frac{e^{2\pi i c(x_1 - x_2)} - 1}{2\pi i c(x_1 - x_2)} C_f(x_1 - x_2). $$

This is a well-defined function for each finite $t$, but we must treat it as a distribution if we want to take the limit $t \to \infty$. To do so, we test it against some $g_1, g_2 \in S(\mathbb{R})$:

$$ (2.24) \qquad \begin{aligned} \mathbb{E}[\langle u(t), g_1 \rangle \overline{\langle u(t), g_2 \rangle}] &= \int_{\mathbb{R}_{x_1} \times \mathbb{R}_{x_2}} \mathbb{E}[u(t, x_1) \overline{u(t, x_2)}] g_1(x_1) g_2(x_2) dx_1 dx_2 \\ &= \int_{\mathbb{R}_z} \frac{e^{2\pi i c t z} - 1}{2\pi i c z} \psi(z) dz \end{aligned} $$

with $\psi(z) = C_f(z) \int_{\mathbb{R}_y} g_1(z + y) g_2(y) dy$, $z = x_1 - x_2$ and $y = x_2$. Our goal is to study the last integral as $t \to \infty$.

We start with a simple identity that gives us a way to integrate the function $(e^{2\pi i c t z} - 1)/(2\pi i c z)$. Note that this function is not absolutely integrable in $\mathbb{R}$. However, its integral does converge
conditionally, i.e. the final result might depend on how we integrate it. More precisely, we recall that for all $t > 0$

$$ (2.25) \qquad \int_{\mathbb{R}} \frac{e^{itz} - 1}{iz} dz := \lim_{R \to \infty} \int_{-R}^{R} \frac{e^{itz} - 1}{iz} dz = \pi. $$

As a result of (2.25), we have that

$$ \lim_{R \to \infty} \int_{-R}^{R} \frac{e^{-2\pi i \mathit{ctz}} - 1}{2\pi i \mathit{cz}} \psi(z) dz - \frac{\psi(0)}{2c} = \lim_{R \to \infty} \int_{-R}^{R} \frac{e^{-2\pi i \mathit{ctz}} - 1}{2\pi i \mathit{cz}} [\psi(z) - \psi(0)] dz. $$

In view of

$$ \psi(0) = C_f(0) \int_{\mathbb{R}_y} g_1(y)g_2(y)dy, $$

it suffices to prove the following in order to obtain (2.10):

$$ (2.26) \qquad \lim_{t \to \infty} \lim_{R \to \infty} \int_{-R}^{R} \frac{e^{-2\pi i \mathit{ctz}} - 1}{2\pi i \mathit{cz}} [\psi(z) - \psi(0)] dz dy = \frac{1}{2\pi c} \int_{0}^{\infty} \frac{\psi(z) - \psi(-z)}{iz} dz \\ = \frac{1}{2\pi i c} \text{p.v.} \int_{\mathbb{R}_z} \frac{\psi(z)}{z} dz. $$

To prove (2.26), we rewrite the left-hand side as follows:

$$ (2.27) \qquad \lim_{R \to \infty} \int_{-R}^{R} \frac{e^{-2\pi i \mathit{ctz}} - 1}{2\pi c iz} \frac{\psi(z) - \psi(0)}{iz} dz = \frac{1}{2\pi i c} \left( \lim_{R \to \infty} \int_{-R}^{R} e^{-2\pi i \mathit{ctz}} \frac{\psi(z) - \psi(0)}{z} dz - \lim_{R \to \infty} \int_{-R}^{R} \frac{\psi(z) - \psi(0)}{z} dz \right) $$

Note that the last term gives the desired limit after using the fact that

$$ \int_{-R}^{R} \frac{\psi(z) - \psi(0)}{z} dz = \int_{0}^{R} \frac{\psi(z) - \psi(-z)}{z} dz $$

and taking $R \to \infty$.

The final step is to show that the first term on the right-hand side of (2.27) tends to zero as $t \to \infty$. Note that

$$ F(z) := \frac{\psi(z) - \psi(0)}{z} $$

is not integrable in $\mathbb{R}_z$. However the following lemma shows that its derivatives have better properties. Its proof is postponed to the end of the proof of Theorem 2.4. Recall that $\psi$ was defined in terms of $g_1, g_2$ right after (2.24).

**Lemma 2.8.** Let $n \ge 1$. If $g_1 \in H^{n+1}(\mathbb{R})$ and $g_2 \in L^2(\mathbb{R})$ then $F \in W^{n,1}(\mathbb{R}_z)$ and

$$ (2.28) \qquad \lim_{|z| \to \infty} (\partial_z^m F)(z) = 0, \text{ for any } 0 \le m \le n, $$

$$ (2.29) \qquad \|\partial_z^m F\|_{L^1(\mathbb{R}_z)} \lesssim \|g_1\|_{H^{m+1}} \|g_2\|_{L^2} \text{ for any } 1 \le m \le n. $$

We integrate by parts the first term on the right-hand side of (2.27)

$$ \int_{-R}^{R} e^{-ictz} \frac{\psi(z) - \psi(0)}{z} dz = - \int_{-R}^{R} \frac{1}{ict} \partial_z e^{-ictz} F(z) dz \\ = - \frac{1}{ict} e^{-ictz} F(z) \Big|_{z=-R}^{R} + \frac{1}{ict} \int_{-R}^{R} e^{-ictz} \partial_z F(z) dz. $$
Using Lemma 2.8, we are able to take the limit $R \to \infty$:

$$ \lim_{R \to \infty} \int_{-R}^{R} e^{-ictz} \frac{\psi(z) - \psi(0)}{z} dz = \frac{1}{ict} \int_{\mathbb{R}_z} e^{-ictz} \partial_z F(z) dz. $$

Next, we continue to integrate by parts using (2.28) (which gets rid of the boundary terms)

$$ \lim_{R \to \infty} \int_{-R}^{R} e^{-ictz} \frac{\psi(z) - \psi(0)}{z} dz = \left(\frac{1}{ict}\right)^n \int_{\mathbb{R}_z} e^{-ictz} \partial_z^n F(z) dz, $$

thus we have

$$ \left| \lim_{R \to \infty} \int_{-R}^{R} e^{-ictz} \frac{\psi(z) - \psi(0)}{z} dz \right| \le \left(\frac{1}{ct}\right)^n \|\partial_z^n F\|_{L^1(\mathbb{R}^2)}. $$

We use (2.29) to finish the proof of Theorem 2.4. Now, we need to prove the technical Lemma 2.8.

*Proof of Lemma 2.8.* First we can easily show that

$$ (2.30) \qquad \partial_z^n F(z) = \int_0^1 s^n \partial_z^{n+1} \psi(zs) ds $$

$$ (2.31) \qquad = \int_0^z \frac{s^n}{z^{n+1}} \partial_z^{n+1} \psi(s) ds. $$

Then we show that for all $0 \le m \le n + 1$ and all $p \in [1, \infty]$

$$ (2.32) \qquad \|\partial_z^m \psi\|_{L^p(\mathbb{R})} \lesssim \|C_f\|_{W^{m,p}(\mathbb{R}_z)} \|g_1\|_{H^m} \|g_2\|_{L^2}. $$

Indeed, the Leibniz rule yields

$$ \partial_z^m \psi(z) = \sum_{j=0}^{m} \binom{m}{j} \partial_z^j C_f(z) \int_{\mathbb{R}_y} \partial_z^{m-j} g_1(y+z) g_2(y) dy. $$

By Cauchy-Schwarz inequality,

$$ |\partial_z^m \psi(z)| \le \left( \sum_{j=0}^m \binom{m}{j} |\partial_z^j C_f(z)| \right) \|g_1\|_{H^m} \|g_2\|_{L^2} $$

thus one gets (2.32).

Then we show that for all $0 \le m \le n$,

$$ \lim_{|z| \to \infty} \partial_z^m F(z) = 0. $$

Indeed, by the Cauchy-Schwarz inequality, (2.31) and (2.32),

$$ |\partial_z^m F(z)| \le \|\partial_z^{m+1} \psi\|_{L^2(\mathbb{R}_z)} \left( \int_0^{|z|} \frac{s^{2m}}{|z|^{2(m+1)}} ds \right)^{\frac{1}{2}} = \frac{\|\partial_z^{m+1} \psi\|_{L^2(\mathbb{R}_z)}}{|z|^{\frac{1}{2}}} \underset{|z| \to \infty}{\longrightarrow} 0. $$

Finally, for all $1 \le m \le n$, (2.30) implies

$$ \|\partial_z^m F(z)\|_{L^1(\mathbb{R}^2)} \le \int_0^1 s^m \|\partial_z^{m+1} \psi(zs)\|_{L^1(\mathbb{R}^2)} ds \le \|\partial_z^{m+1} \psi\|_{L^1(\mathbb{R}^2)} \left( \int_0^1 s^{m-1} ds \right) $$

and thus we conclude with (2.32).

□
Natural generalizations of the dynamics in (2.1)-(2.2) to higher dimensions could be obtained in two ways. Firstly, the product $c \partial_k \hat{u}$ in the transport equation (2.1) could be generalized to a scalar product of a given unit vector $e \in \mathbb{R}^d$ with the gradient $\nabla_k \hat{u}(t, k)$. Unfortunately, this puts too much weight on the constant vector $e$ and results in an obvious statistical anisotropy in physical space. For applications to turbulence, we must require that statistical laws are not only invariant by translation, but also under rotation (i.e. statistical isotropy), as commonly observed in laboratory and numerical experiments, and as expected from a physical point of view. Another option would be to replace the multiplication by $ix$ in the physical space formulation (2.2) by the multiplication by $i|x|$, where $|x|$ is the modulus of $x \in \mathbb{R}^d$. Once again, this would introduce anisotropy in the system, and more importantly, it would break statistical homogeneity even in dimension $d=1$. One may check these claims directly using the exact solution (2.2) to compute the covariance function at a given time $t$ and any two positions $x, y$ (see [3] for details). Indeed, this covariance function eventually depends on the difference $|x| - |y|$, and not on $x - y$ as would be desirable. Beyond these issues, none of these propositions would ensure that a given real-valued initial condition $u(0, x) \in \mathbb{R}$ gives rise to a real-valued solution $u(t, x) \in \mathbb{R}$ at all future times. In other words, in order to construct a real-valued solution in physical space, one needs to propose a dynamical picture able to preserve the Hermitian symmetry of the Fourier transform $\hat{u}(t, k)$, as does the dynamics in (2.3).

### 3. HIGHER DIMENSIONAL REAL FRACTIONAL GAUSSIAN FIELDS: HEURISTIC

In the previous section we gave a rigorous proof of the construction of a dynamical complex white noise. This proof was carried out in physical space. For the dynamical real white noise of (2.3), on the other hand, it is more convenient to think of its wavenumber formulation. However, even in the case of a dynamical complex white noise, the solution $u(t, x)$ is not in $L^p(\mathbb{R}_x)$ for any $1 \le p < \infty$ (see (2.18)), hence its Fourier transform is not defined pointwise. In conclusion, it is difficult to make sense of equation (2.1) in wavenumber space. However, we will not concern ourselves with such difficulties in this section, and we will work as if the solution to (2.1) were well-defined pointwise: we refer to Section 4 for a rigorous analysis. As we will later see, working in wavenumber space is very convenient to formally show that (2.3) builds a dynamical real white noise, as well as to extend this construction to $d$-dimensional fractional Gaussian fields.

**3.1. A transport equation in wavenumber space.** We propose the following initial value problem as a generalization of (2.3)

$$ (3.1) \quad \begin{cases} \partial_t \hat{u}(t, k) + \operatorname{div}_k \left( \frac{ck}{|k|} \hat{u}(t, k) \right) + c \frac{H + \frac{1}{2}}{|k|} \hat{u}(t, k) = \hat{f}(t, k) & t > 0, k \in \mathbb{R}^d, |k| > \kappa > 0, \\ \hat{u}(t, k) = 0 & t > 0, k \in \mathbb{R}^d, |k| \le \kappa, \\ \hat{u}(0, k) = 0. & \end{cases} $$

where $c, \kappa > 0$ are fixed, $H$ is a real constant (which will be eventually connected with the Hölder exponent of the solution), and $\operatorname{div}_k$ stands for the usual divergence operator in wavenumber space. As one can easily verify with a few vector calculus identities,

$$ (3.2) \quad \operatorname{div}_k \left( \frac{ck}{|k|} \hat{u}(k) \right) + c \frac{H + \frac{1}{2}}{|k|} \hat{u}(k) = c \partial_{|k|} \hat{u}(k) + c \frac{H + d - \frac{1}{2}}{|k|} \hat{u}(k) \quad \text{with} \quad \partial_{|k|} := \frac{k}{|k|} \cdot \nabla_k. $$

When $d=1$ and $H = -1/2$, one recovers (2.3) from the above. The initial value problem (3.1) can be regarded as a conservation law in wavenumber space with a source term $\hat{f}(t, k)$, a damping term $(H + 1/2) \hat{u}(t, k)/|k|$ and Dirichlet boundary conditions at the sphere $|k| = \kappa$.

Equation (3.1) has been thoroughly studied when $\hat{f}(t, k)$ is regular enough (see [27], [34]), but in our case $\hat{f}(t, k)$ is too rough for such classical results to be applicable. In particular, we will
assume that the forcing term $\hat{f}$ satisfies

$$ (3.3) \qquad \mathbb{E}[\hat{f}(t, k_1) \overline{\hat{f}(s, k_2)}] = \hat{C}_f(k_1) \delta_{t-s} \delta_{k_1-k_2}, $$

where $\hat{C}_f(k)$ is radial and null in the ball $|k| < \kappa$. Condition (3.3) formally follows from considering the Fourier transform of a white noise in time satisfying

$$ (3.4) \qquad \mathbb{E}[f(t, x)f(s, y)] = \delta_{t-s} C_f(x - y). $$

**Remark 3.1.** As part of equation (3.1) one has the following technical condition

$$ (3.5) \qquad \hat{u}(t, k) = 0 \quad t > 0, \quad k \in \mathbb{R}^d, \quad |k| \le \kappa, $$

for $\kappa > 0$. Indeed, if $|k| = 0$, then $\text{div}_k \left( \frac{k}{|k|} \cdot \right)$ and $\frac{H+1}{|k|}$ are not well-defined. This condition (3.5) implies in particular that $\mathcal{F}^{-1}\hat{u}$ has null spatial average, whenever defined. It might be possible to make sense of the problem (3.1) for $\kappa = 0$ by adequately changing condition (3.5), imposing $\hat{f}(t, k = 0) = 0$ and an appropriate behaviour near $k = 0$, but this is outside the scope of this paper.

3.2. **Asymptotic behavior: power-law.** In this section, we show that the two-point correlation of the solution to (3.1) displays a power-law behavior. In our first result (Theorem 3.2), we obtain this asymptotic behavior as $t \to \infty$ under fairly mild assumptions on the forcing $\hat{f}$. Under stronger assumptions on $\hat{f}$, we derive a second result (Proposition 3.5) showing this power law behavior in finite time. Among other things, such power laws are important because their exponent determines the Hölder regularity of the solution in physical space (should it be possible to take the inverse Fourier transform). The main idea of the heuristic proof of our desired results is to perform a change of variables to rewrite (3.1) as a 1D transport equation with respect to $|k|$ and parametrized by the “angular variable” $\frac{k}{|k|}$. Such an equation admits an explicit solution that we will exploit. We will further discuss such consequences in Section 4.

**Theorem 3.2 (Heuristic version).** Let the forcing $\hat{f}$ satisfy (3.3) in such a way that $\hat{C}_f(k) = \psi(|k|)$ is radial, non-negative, non identically null, with $s^{2H+d} \psi(s) \in L^1(\mathbb{R}_s^+)$.
Furthermore, we assume that $\psi$ is null when $|k| < \kappa$. Then, (3.1) admits a solution that satisfies the following asymptotic behavior

$$ \lim_{t \to \infty} \mathbb{E}[\hat{u}(t, k) \overline{\hat{u}(t, k')}] = |k|^{-(2H+d)} (C(d, H) - \Psi_{d,H}(|k|)) \delta_{k-k'} $$

where

$$ (3.6) \qquad \Psi_{d,H}(|k|) := \frac{1}{c} \int_{|k|}^{\infty} s^{2H+d} \psi(s) ds. $$

is a positive non-increasing absolutely continuous function and

$$ (3.7) \qquad C(d, H) := \Psi_{d,H}(0) > 0. $$

As we have already pointed out in Remark 2.6, the function $\Psi_{d,H}$ can be seen as lower-order correction in comparison with the Dirac distribution. Moreover, if we assume that $\psi$ is fast-decaying, then $\Psi_{d,H}$ will also be fast-decaying.

**Heuristic proof.** In order to find a formal solution to (3.1), we let

$$ \hat{v}(t, k) := |k|^{H+d-\frac{1}{2}} \hat{u}(t, k), \quad \hat{g}(t, k) := |k|^{H+d-\frac{1}{2}} \hat{f}(t, k). $$
Next we rewrite (3.1) in terms of $\hat{v}$, namely

$$ (3.8) \qquad \begin{cases} \partial_t \hat{v}(t, k) + c \partial_{|k|} \hat{v}(t, k) = \hat{g}(t, k) & t > 0, |k| > \kappa, \\ \hat{v}(t, k) = 0 & t > 0, |k| \le \kappa, \\ \hat{v}(0, k) = 0, & t < 0, |k| > \kappa. \end{cases} $$

where we use the fact that $\partial_{|k|} = \frac{k}{|k|} \cdot \nabla_k$.

Equation (3.8) can be regarded as a 1D transport equation with respect to $|k|$ and parametrized by the “angular variable” $k/|k|$. Thus it is easy to give an explicit solution:

$$ (3.9) \qquad \hat{v}(t, k) = \int_{(t - \frac{|k| - \kappa}{c})_+}^{t} \hat{g} \left( s, (|k| - ct + cs) \frac{k}{|k|} \right) ds. $$

Next we compute the correlations of $\hat{v}(t)$ using those of $\hat{f}$. We have

$$ \mathbb{E}[\hat{v}(t, k_1) \overline{\hat{v}(t, k_2)}] = \int_{(t - \frac{|k_1| - \kappa}{c})_+}^{t} \int_{(t - \frac{|k_2| - \kappa}{c})_+}^{t} \delta_{cs_1 - cs_2} \delta_{\frac{|k_1| - ct + cs_1}{|k_1|}} k_1 - \delta_{\frac{|k_2| - ct + cs_2}{|k_2|}} k_2 \hat{C}_g \left( \frac{|k_1| - ct + cs_1}{|k_1|} k_1 \right) ds_1 ds_2. $$

We assume that we can write

$$ \delta_{cs_1 - cs_2} \delta_{\frac{|k_1| - ct + cs_1}{|k_1|}} k_1 - \delta_{\frac{|k_2| - ct + cs_2}{|k_2|}} k_2 = \delta_{cs_1 - cs_2} \delta_{\frac{|k_1| - ct + cs_1}{|k_1|}} k_1 - \delta_{\frac{|k_2| - ct + cs_1}{|k_2|}} k_2 $$

even if this not mathematically rigorous. Moreover, by the change of variables from cartesian to polar coordinates

$$ \delta_{\frac{|k_1| - ct + cs}{|k_1|}} k_1 - \delta_{\frac{|k_2| - ct + cs}{|k_2|}} k_2 = \left( \frac{|k_1| - ct + cs}{|k_1|} \right)^{-(d-1)} \delta_{k_1 - k_2}. $$

This is possible since the Jacobian $(\frac{|k_1| - ct + cs}{|k_1|})^{-(d-1)}$ has no singularities in the region of integration thanks to $\kappa > 0$. As a result, one obtains

$$ \mathbb{E}[\hat{v}(t, k) \overline{\hat{v}(t, k')}] = \left( \int_{(t - \frac{|k| - \kappa}{c})_+}^{t} \left( \frac{|k| - ct + cs}{|k|} \right)^{-(d-1)} \hat{C}_g \left( \frac{|k| - ct + cs}{|k|} k \right) ds \right) \delta_{k-k'}. $$

This immediately implies

$$ \mathbb{E}[\hat{u}(t, k) \overline{\hat{u}(t, k')}] = |k|^{-(2H+d)} \left( \int_{(t - \frac{|k| - \kappa}{c})_+}^{t} (|k| - ct + cs)^{2H+d} \hat{C}_f \left( \frac{|k| - ct + cs}{|k|} k \right) ds \right) \delta_{k-k'}. $$

The change of variables $s \mapsto |k_1| - ct + cs$ yields

$$ \mathbb{E}[\hat{u}(t, k) \overline{\hat{u}(t, k')}] = |k|^{-(2H+d)} \left( \chi_{|k|>ct+\kappa} \int_{|k|-ct}^{|k|} s^{2H+d} \hat{C}_f \left( \frac{s}{|k|} k \right) \frac{ds}{c} \right) \delta_{k-k'} \\ + |k|^{-(2H+d)} \left( \chi_{|k| \le ct+\kappa} \int_{\kappa}^{|k|} s^{2H+d} \hat{C}_f \left( \frac{s}{|k|} k \right) \frac{ds}{c} \right) \delta_{k-k'} $$

where $\chi_A$ denotes the characteristic function of the set A. Remember that $\psi(|k|) := \hat{C}_f(k)$ since $C_f$ is radial. Using (3.6), we may rewrite

$$ (3.10) \qquad \mathbb{E}[\hat{u}(t, k) \overline{\hat{u}(t, k')}] = |k|^{-(2H+d)} \left( \chi_{|k|>ct+\kappa} [\Psi_{d,H}(|k| - ct) - \Psi_{d,H}(|k|)] \right. \\ \qquad \qquad \qquad \qquad \qquad \left. + \chi_{|k| \le ct+\kappa} [\Psi_{d,H}(\kappa) - \Psi_{d,H}(|k|)] \right) \delta_{k-k_2}. $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. The functions $|k| \mapsto F(t_1, |k|)$ and $|k| \mapsto F(t_2, |k|)$ for $t_1 < t_2$ and $\psi$ supported in $[\kappa, r_*]$.</figcaption></figure>

We conclude by taking the limit $t \to \infty$ and noting that $\Psi_{d,H}(\kappa) = \Psi_{d,H}(0)$. $\square$

**Remark 3.3.** In the previous proof, we have shown that (see (3.10)) for any fixed $t > 0$

$$ \mathbb{E}[\hat{u}(t, k)\overline{\hat{u}(t, k')}] = |k|^{-(2H+d)} F(t, |k|) \delta_{k-k'} $$

where $F$ is a non-negative absolutely continuous function in $t$ and $|k|$ that admits the following expression

$$ (3.11) \qquad F(t, |k|) = \begin{cases} \Psi_{d,H}(\kappa) - \Psi_{d,H}(|k|), & |k| < ct + \kappa, \\ \Psi_{d,H}(|k| - ct) - \Psi_{d,H}(|k|), & |k| > ct + \kappa. \end{cases} $$

In particular, for any time $t_1 \le t_2$, $F(t_1, |k|) = F(t_2, |k|)$ for all $|k| \le ct_1 + \kappa$. As $t$ grows the window of $|k|$ where $F(t, |k|)$ is stationary (and displays a power-law behavior) grows too (see figure 1). As a result, it is not necessary to wait until $t = \infty$ in order to observe the power-law. This is helpful when performing numerical simulations.

**Remark 3.4.** The limits $t \to \infty$ and $|k| \to \infty$ in (3.11) don't commute. Indeed, one gets

$$ \lim_{t \to \infty} \mathbb{E}[\hat{u}(t, k)\overline{\hat{u}(t, k')}] \underset{|k| \to \infty}{\sim} |k|^{-(2H+d)} C(d, H) \delta_{k-k'}, $$

whereas for any fixed time $t > 0$

$$ \lim_{|k| \to \infty} \mathbb{E}[\hat{u}(t, k)\overline{\hat{u}(t, k')}] = 0. $$

In simulations, it is usual to look for the power-law for large $|k|$. But, as pointed out in the previous remark, if we don't wait long enough then we will see nothing. The next proposition shows that if the two-point correlation of the source is spectrally located in the ball $|k| \le r_*$, then it is enough to wait until $t \ge r_*$ to observe the desired power-law.

**Proposition 3.5.** Let $\psi$ supported in $[\kappa, r_*]$ for $r_* > \kappa > 0$. For any $t > r_*$ fixed, the function $F$ given in Remark 3.3 satisfies

$$ F(t, |k|) = \begin{cases} 0 & |k| \in (0, \kappa) \\ \text{stationary in } t \text{ and non-decreasing in } |k| & |k| \in (\kappa, r_*), \\ \text{constant} & |k| \in (r_*, ct + \kappa), \\ \text{non-decreasing in } t \text{ and non-increasing in } |k| & |k| \in (ct + \kappa, ct + r_*), \\ 0 & |k| \in (ct + r_*, \infty). \end{cases} $$
**Remark 3.6.** Since $\psi$ vanishes in the ball $|k| < \kappa$,

$$ \Psi(0) = \Psi(\kappa) = \Psi(|k| - ct) \quad \text{for } ct \le |k| \le ct + \kappa, $$

and thus the expression

$$ F(t, |k|) = \begin{cases} \Psi_{d,H}(0) - \Psi_{d,H}(|k|), & |k| < ct, \\ \Psi_{d,H}(|k| - ct) - \Psi_{d,H}(|k|), & |k| > ct. \end{cases} $$

is independent of $\kappa$. On this note, see Remark 3.1.

In order that the solution of (3.12) be closer to the experimental energy-spectrum picture (see beginning of section 1.2), one can add some viscosity in the equation. More precisely, one can consider $f - \nu\Delta u$ instead of a sole forcing term $f$.

**Theorem 3.7 (Heuristic version).** Let $\nu \ge 0$ and let the forcing $\hat{f}$ satisfy (3.3) such that $\hat{C}_f(k) = \psi(|k|)$ is radial, non-negative, non identically null, with $s^{2H+d}e^{\frac{8\pi^2\nu}{3c}s^3} \psi(s) \in L^1(\mathbb{R}_s^+)$ and $\psi$ null in the ball $|k| < \kappa$. Then (3.12)

$$ \left\{ \begin{array}{ll} \partial_t \hat{u}_\nu(t,k) + \text{div}_k \left( \frac{ck}{|k|} \hat{u}_\nu(t,k) \right) + \left( c \frac{H + \frac{1}{2}}{|k|} + \nu |2\pi k|^2 \right) \hat{u}_\nu(t,k) = \hat{f}(t,k) & t > 0, k \in \mathbb{R}^d, |k| > \kappa > 0, \\ \\ \hat{u}_\nu(t,k) = 0 & t > 0, k \in \mathbb{R}^d, |k| \le \kappa, \\ \\ \hat{u}_\nu(0,k) = 0 & \end{array} \right. $$

admits a solution that satisfies for any fixed $t > 0$

$$ \mathbb{E}[\hat{u}_\nu(t,k) \overline{\hat{u}_\nu(t,k')}] = |k|^{-(2H+d)} e^{-\frac{8\pi^2\nu}{3c} |k|^3} F_\nu(t, |k|) \delta_{k-k'} $$

where $F_\nu$ is a non-negative absolutely continuous function in $t$ and $|k|$ that admits the following expression

$$ (3.13) \qquad F_\nu(t, |k|) = \begin{cases} \Psi_{d,H,\nu}(\kappa) - \Psi_{d,H,\nu}(|k|), & |k| < ct + \kappa, \\ \Psi_{d,H,\nu}(|k| - ct) - \Psi_{d,H,\nu}(|k|), & |k| > ct + \kappa. \end{cases} $$

and where

$$ (3.14) \qquad \Psi_{d,H,\nu}(|k|) := \frac{1}{c} \int_{|k|}^{\infty} s^{2H+d} e^{\frac{8\pi^2\nu}{3c}s^3} \psi(s) ds. $$

is a positive non-increasing absolutely continuous function.

*Proof.* In order to find a formal solution to (3.12), we let

$$ \hat{v}(t,k) := |k|^{H+d-\frac{1}{2}} e^{\frac{4\pi^2\nu}{3c} |k|^3} \hat{u}(t,k), \quad \hat{g}(t,k) := |k|^{H+d-\frac{1}{2}} e^{\frac{4\pi^2\nu}{3c} |k|^3} \hat{f}(t,k) $$

and rewriting (3.12) for $\hat{v}$, we find that $\hat{v}$ solves (3.8). Thus the proof is similar to the one of Theorem 3.2. $\square$

Note that in the presence of viscosity, the limits $\nu \to 0$ and $t \to \infty$ commute.
3.3. Main issues with the heuristic proof. The calculations in this section provide invaluable intuition regarding the construction of an equation whose solution displays a power-law behavior as well as the right parameters involved. Nevertheless, there are several steps in our calculations that are difficult to justify from a rigorous mathematical viewpoint (if possible at all). More precisely:

(1) Fourier-transform. As we already mentioned, the Fourier transform of a rough non-decreasing source (and therefore the solution of equation (3.1)) is unlikely to be defined pointwise.

(2) Transport equation in wavenumber space. We have never showed that (3.1) is well-posed. Even if it were in some weak sense, it is unclear whether one can define the trace of the solution at the boundary $|k| = \kappa$ or to define the domain of the operator $\text{div}_k (\frac{k}{|k|})$.

(3) Solution of a stochastic transport equation. In the simplest case, (3.1) becomes (3.8) (for $H = \frac{1}{2} - d$) whose solution we gave in (3.9), namely

$$ (3.15) \qquad \hat{u}(t,k) = \int_{(ct-(|k|-\kappa))_+}^{t} \hat{f} \left( s, (|k|-ct+cs) \frac{k}{|k|} \right) ds. $$

In the case where $f$ is a white noise in time of the form (2.6) (in $\mathbb{R}^d$ instead of $\mathbb{R}$), the integral on the right-hand side of (3.15) is not a Riemann integral. Worse than that, one would need to define the object $\hat{f}(t,k)$ and evaluate it along the characteristics of the transport equation. To sum up, one would need to make sense of the “stochastic integral”:

$$ \int_{(t-(|k|-\kappa))_+}^{t} \hat{\varphi} \left( (|k|-ct+cs) \frac{k}{|k|} \right) \widehat{dW} \left( s, (|k|-ct+cs) \frac{k}{|k|} \right). $$

Moreover, without a rigorous functional setup it is hard to study the important problem of uniqueness of solutions. In fact, it is unclear whether the solution (3.15) is unique.

(4) Two-point correlations. Even if we can make sense of (3.1) weakly, e.g. as a tempered distribution, the two-point correlation is not necessary well defined. Indeed, it involves a product of distributions whose singular supports might overlap.

(5) Mixing coordinates in Dirac distributions. In the proof, we use the following identity

$$ \delta_{s_1-s_2} \delta \frac{|k_1|-ct+cs_1}{|k_1|} k_1 - \delta \frac{|k_2|-ct+cs_2}{|k_2|} k_2 = \delta_{s_1-s_2} \delta \frac{|k_1|-ct+cs_1}{|k_1|} k_1 - \delta \frac{|k_2|-ct+cs_1}{|k_2|} k_2. $$

In other words, we first apply $\delta_{s_1-s_2}$ to $\delta \frac{|k_1|-ct+cs_1}{|k_1|} k_1 - \delta \frac{|k_2|-ct+cs_2}{|k_2|} k_2$. Can this be justified?

This is in particular possible if the inner delta were a function, which is of course not the case. Later, in the proof, we perform a change of coordinates to obtain

$$ \delta_{s_1-s_2} \delta \frac{|k_1|-ct+cs_1}{|k_1|} k_1 - \delta \frac{|k_2|-ct+cs_2}{|k_2|} k_2 = \left( \frac{|k_1|-ct+cs_1}{|k_1|} \right)^{-(d-1)} \delta_{s_1-s_2} \delta_{k_1-k_2}. $$

Does the same final result hold if we exchanged the order in which the deltas on the left-hand side are applied?

Ultimately, we would also like to write equation (3.1) in physical space and to interpret the power-law in physical space. One first idea is to perform an inverse Fourier transform of (3.1) but the difficulties mentioned above make it complicated.
4. HIGHER DIMENSIONAL REAL FRACTIONAL GAUSSIAN FIELDS: RESULTS

As mentioned in the previous subsection, it is difficult to justify the Fourier transform of a rough
source such as the one that interests us. Similarly, the solution to such a problem (3.1) may not
admit a Fourier transform either. For that reason, our goal is to try to formulate a weaker notion
of (3.1) which is entirely given in physical space, thereby avoiding such problems. To do so, we
first introduce the Hilbert space

$$
(4.1) \qquad X := \{u \in L^2(\mathbb{R}_x^d) \mid \hat{u}(k) = 0 \text{ for all } |k| < \kappa \}
$$

together with the inner product

$$
(u, v) := \int_{\mathbb{R}_k^d \setminus B(0, \kappa)} \hat{u}(k) \overline{\hat{v}(k)} dk.
$$

Next we introduce the unbounded operator

$$
(4.2) \qquad A : u \mapsto \mathcal{F}^{-1} \left[ c \operatorname{div}_k \left( \frac{k}{|k|} \hat{u} \right) + c \frac{H + \frac{1}{2}}{|k|} \hat{u} \right]
$$

with domain

$$
(4.3) \qquad
\begin{aligned}
D(A) &:= \{u \in X \mid Au \in X \text{ and } \hat{u}|_{|k|=\kappa} = 0\} \\
&= \{u \in X \mid \partial_{|k|} \hat{u} \in L^2(\mathbb{R}_k^d \setminus B(0, \kappa)) \text{ and } \hat{u}|_{|k|=\kappa} = 0\}.
\end{aligned}
$$

For any function $u \in X$ such that $\partial_{|k|} \hat{u} \in L^2(\mathbb{R}_k^d \setminus B(0, \kappa))$, the trace in Fourier space $\hat{u}|_{|k|=\kappa}$ is well-defined thanks to the Sobolev embedding theorem. This allows us to impose $\hat{u}|_{|k|=\kappa} = 0$. Consequently, the transport equation in wavenumber space (3.1) can be written in physical space after formally applying the inverse Fourier transform:

$$
(4.4) \qquad \begin{cases} \partial_t u + Au = f & \text{for } t > 0, \\ u|_{t=0} = 0. \end{cases}
$$

At this stage, note that for any regular enough forcing $f \in L_{\text{loc}}^1((0, +\infty)_t, X)$, the problem (4.4)
is rigorously defined. Its mild solution is given by:

$$
(4.5) \qquad u(t) = \int_0^t e^{-(t-s)A} f(s) ds
$$

where $e^{-tA} : X \to X$ is defined by

$$
(4.6) \qquad e^{-tA}u_0 := \mathcal{F}^{-1} \left[ \chi_{|k|>ct+\kappa} \left( \frac{|k|-ct}{|k|} \right)^{H+d-\frac{1}{2}} \hat{u}_0 \left( \frac{|k|-ct}{|k|}k \right) \right].
$$

Indeed, it is easy to check that (4.5) solves (4.4). Moreover, it is also easy to check that if $u_0 \in X \cap S(\mathbb{R}_x^d)$ then $e^{-tA}u_0 \in X \cap S(\mathbb{R}_x^d)$.

Unfortunately, our source $f$ is a white noise in time (and colored in space) and therefore the above considerations do not apply. Rigorously speaking, $f$ is a Gaussian random measure acting on $L^2((0, +\infty)_t)$, i.e.

$$
(4.7) \qquad \langle f, g \rangle(x) := \int_{(0, +\infty)_t \times \mathbb{R}_y^d} \tau_y \varphi(x) g(t) dW(t, y),
$$

for all test functions $g \in L^2((0, +\infty)_t)$ where $dW$ is a space-time Gaussian real white noise, $\tau_y$ is a translation by $y$ (i.e. $\tau_y g(x) = g(x - y)$) and $\varphi \in X \cap S(\mathbb{R}_x^d)$ is a real, non-identically null, radial function.
In order to give a meaning to solutions of (4.4) with rough source, we need a weak formulation.
To propose one, one first needs to introduce the adjoint of A.

**Proposition 4.1.** *The adjoint of A in (4.2) is defined by*

$$
(4.8) \qquad A^* : u \mapsto \mathcal{F}^{-1} \left[ -c \operatorname{div}_k \left( \frac{k}{|k|} \hat{u} \right) + c \frac{H+d-\frac{1}{2}}{|k|} \hat{u} \right]
$$

and its domain is

$$
D(A^*) := \{u \in X \mid \partial_{|k|} \hat{u} \in L^2(\mathbb{R}_k^d \setminus B(0, \kappa)) \} \supseteq D(A).
$$

*Proof.* Let $u \in D(A)$ and $v \in X$. Using (3.2), one gets

$$
c^{-1}(Au, v)_{L^2} = \int_{\mathbb{R}_k^d \setminus B(0, \kappa)} \partial_{|k|} \hat{u}(k) \bar{\hat{v}}(k) dk + \int_{\mathbb{R}_k^d \setminus B(0, \kappa)} \frac{H+d-\frac{1}{2}}{|k|} \hat{u}(k) \bar{\hat{v}}(k) dk.
$$

We perform a change of variables from cartesian to polar coordinates in wavenumber space:

$$
\int_{\mathbb{R}_k^d} \partial_{|k|} \hat{u}(k) \overline{\hat{v}(k)} dk = \int_{\mathbb{R}^+} \int_{S_\theta} \partial_{|k|} \hat{u}(k) \overline{\hat{v}(k)} |k|^{d-1} d\sigma_\theta d|k|.
$$

We now assume that $\partial_{|k|} \hat{v} \in L^2(\mathbb{R}_k^d \setminus B(0, \kappa))$ so that we can integrate by parts in variable $|k|$ and immediately return to cartesian coordinates:

$$
\int_{\mathbb{R}_k^d} \partial_{|k|} \hat{u} \bar{\hat{v}} dk = - \int_{\mathbb{R}_k^d} \hat{u} \partial_{|k|} \bar{\hat{v}} dk - \int_{\mathbb{R}_k^d} \frac{(d-1)}{|k|} \hat{u} \bar{\hat{v}} dk - \int_{S_\theta} \kappa^{d-1} \hat{u} |_{|k|=\kappa} \bar{\hat{v}} |_{|k|=\kappa} d\sigma_\theta.
$$

Since $\hat{u}|_{|k|=\kappa} = 0$, the previous equality yields

$$
c^{-1}(Au, v)_{L^2} = - \int_{\mathbb{R}_k^d} \hat{u} \partial_{|k|} \bar{\hat{v}} + \int_{\mathbb{R}_k^d} \frac{H + \frac{1}{2}}{|k|} \hat{u} \bar{\hat{v}} dk.
$$

Using (3.2) once again, one readily gets

$$
c^{-1}(Au, v)_{L^2} = \int_{\mathbb{R}_k^d \setminus B(0, \kappa)} \hat{u}(k) \left[ -\operatorname{div}_k \left( \frac{k}{|k|} \bar{\hat{v}}(k) \right) + \frac{H+d-\frac{1}{2}}{|k|} \bar{\hat{v}}(k) \right] dk
$$

for all $u \in D(A)$ and $v \in X$ such that $\partial_{|k|} \hat{v} \in L^2(\mathbb{R}_k^d \setminus B(0, \kappa))$.

Note that it is not necessary to impose $\hat{v}|_{|k|=\kappa} = 0$ and thus $D(A^*)$ is strictly larger than $D(A)$. $\square$

**Remark 4.2.** In the definition of $D(A)$ the condition $\hat{u}|_{|k|=\kappa} = 0$ can be interpreted as a boundary condition in wavenumber space. This boundary condition is lost in $D(A^*)$, which is strictly larger than $D(A)$. It might be possible to take $\kappa = 0$ by adequately exchanging $\hat{u}|_{|k|=\kappa} = 0$ by $\hat{u}(t, k = 0) = 0$, and by imposing an appropriate behaviour near $k = 0$, but this is outside the scope of this paper.

We are ready to give the weak formulation of (4.4).

**Definition 4.3.** We say that a stochastic process $u$ is a weak solution to (4.4) with $f$ given by (4.7), if $u \in L_{\text{loc}}^1([0, \infty)_t \times \mathbb{R}^d)$ almost surely and

$$
(4.9) \quad \int_{\mathbb{R}_x^d} u(t, x)g(t, x)dx + \int_0^t \int_{\mathbb{R}_x^d} u(s, x)(-\partial_s + A^*)g(s, x)dsdx
$$
$$= \int_0^t \int_{\mathbb{R}_y^d} \left( \int_{\mathbb{R}_x^d} \varphi(x-y)g(s,x)dx \right) dW(s,y)$$

for all $g \in C^\infty((0, \infty), D(A^*) \cap S(\mathbb{R}_x^d))$, and such that

$$(4.10) \qquad \int_0^t \int_{\mathbb{R}_x^d} u(s,x)g(s,x)dsdx = 0$$

for all $g \in C^\infty((0, \infty), S(\mathbb{R}_x^d))$ with $\hat{g}(t,k) = 0$ for all $|k| > \kappa$.

Note that for a test function $g \in X \cap S(\mathbb{R}_x^d)$, $g \in D(A^*)$ and $A^*g \in D(A^*) \cap S(\mathbb{R}_x^d)$ and thus the terms on the left-hand side of (4.9) are well-defined.

We are ready to give the main result of this paper for positive $H \in (0, 1)$.

**Theorem 4.4.** Let $H \in (0, 1)$, and let the forcing $f$ in (4.4) be

$$(4.11) \qquad f(t,x) = \int_{\mathbb{R}_y^d} \tau_y \varphi(x) dW(t,y),$$

where $dW$ is a space-time Gaussian real white noise, $\tau_y$ is a translation by $y$ (i.e. $\tau_y g(x) = g(x-y)$) and $\varphi \in S(\mathbb{R}_x^d)$ is a real, non-identically null, radial function such that $\hat{\varphi}(k) = 0$ for all $|k| < \kappa$. Then:

(i) Equation (4.4) admits a unique global (in time) weak solution $u(t,x)$ in $C((0, \infty), (X \cap S(\mathbb{R}_x^d))')$ given by

$$(4.12) \qquad u(t,x) = \int_0^t \int_{\mathbb{R}_y^d} e^{-(t-s)A}[\tau_y \varphi](x) dW(s,y)$$

where $e^{-tA}: X \cap S(\mathbb{R}_x^d) \to X \cap S(\mathbb{R}_x^d)$ is defined in (4.6).

(ii) For each $t > 0$ fixed, $u(t,x)$ is a Gaussian process with a.s. smooth paths in $x$, and with $\alpha$-Hölder continuous paths in $t$ for any $0 < \alpha < 1/2$.

(iii) As $t \to \infty$, the solution $u(t,x)$ converges in law to a zero-mean Gaussian field $u_\infty(x)$ with $\alpha$-Hölder continuous paths in $x$ for any $0 < \alpha < H$.

(iv) The limiting process $u_\infty(x)$ is characterized by the correlations:

$$(4.13) \qquad \begin{aligned} \mathbb{E}[u_\infty(x_1)u_\infty(x_2)] &= \lim_{t \to \infty} \mathbb{E}[u(t, x_1)u(t, x_2)] \\ &= C(d, H) \mathcal{K}_H(x_1 - x_2) - (\mathcal{K}_H * \mathcal{J}_H)(x_1 - x_2), \end{aligned}$$

where

$$\mathcal{K}_H := \mathcal{F}^{-1}[\chi_{|k|>\kappa}|k|^{-(2H+d)}] \quad \text{and} \quad \mathcal{J}_H := \mathcal{F}^{-1}[\chi_{|k|>\kappa} \Psi_{d,H}(|k|)]$$

and where $C(d, H)$ and $\Psi_{d,H}$ are given in Theorem 3.2.

In the case of non-positive $H$, one needs to be more careful since the solution converges to a very rough object. In this direction, we have the following result:

**Theorem 4.5.** Let $H \in [-d/2, 0]$, and let the forcing $f$ as in (4.11). Then points (i) and (ii) in Theorem 4.4 still hold. Moreover, we have that:

(iii) As $t \to \infty$, the solution $u(t)$ converges in law to a random Gaussian measure $u_\infty$ acting on $X \cap S(\mathbb{R}_x^d)$ with zero-mean, i.e. for any $g \in X \cap S(\mathbb{R}^d)$, $\langle u(t), g \rangle$ converges in law to a Gaussian random variable $\langle u_\infty, g \rangle$ with $\mathbb{E}[\langle u_\infty, g \rangle] = 0$.
(iv) We have the following asymptotic behavior:

$$ (4.14) \qquad \begin{aligned} \lim_{t \to \infty} \mathbb{E}[\langle u(t), g_1 \rangle \langle u(t), g_2 \rangle] &= \mathbb{E}[\lim_{t \to \infty} \langle u(t), g_1 \rangle \langle u(t), g_2 \rangle] \\ &= C(d, H) \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \\ &\quad - \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} \Psi_{d,H}(|k|) \hat{g}_1(k) \overline{\hat{g}_2(k)} dk. \end{aligned} $$

for any $g_1, g_2 \in X \cap S(\mathbb{R}^d)$, where $\Psi_{d,H}$ are given in Theorem 3.2.

Before we prove these theorems, let us make a final comment about the operator $A$. This operator can be regarded as an operator of degree 0 plus a bounded operator of degree $-1$ as shown in the following proposition.

**Proposition 4.6.** For all $u \in D(A)$ and $v \in D(A^*)$, the following equalities hold:

$$ (4.15) \qquad \text{Au}(x) = \int_{\mathbb{R}_k^d / B_\kappa} a(x, k) e^{ik \cdot x} \hat{u}(k) dk \quad \text{and} \quad \text{A}^*v(x) = \int_{\mathbb{R}_k^d / B_\kappa} a^*(x, k) e^{ik \cdot x} \hat{v}(k) dk. $$

where

$$ a(x, k) := -i \frac{k}{|k|} \cdot x + \frac{H + \frac{1}{2}}{|k|} \quad \text{and} \quad a^*(x, k) := i \frac{k}{|k|} \cdot x + \frac{H + d - \frac{1}{2}}{|k|}. $$

For any multi-index $(\alpha, \beta)$ with $|\alpha| \ge 1$, one has

$$ |x||k|^{-|\beta|} \lesssim |\partial_k^\beta a(x, k)| + |\partial_k^\beta a^*(x, k)| \lesssim |x||k|^{-|\beta|} + |k|^{-1-|\beta|} $$

and

$$ |\partial_x^\alpha \partial_k^\beta a(x, k)| + |\partial_x^\alpha \partial_k^\beta a^*(x, k)| \lesssim |k|^{-|\beta|}, $$

for any $(x, k) \in \mathbb{R}^{2d}$ such that $|k| \ge \kappa$. The implicit constant doesn't depend on $\kappa$, $x$ or $k$.

**Proof.** Let $u \in D(A)$. Using the vector calculus identities (3.2), one gets

$$ Au = \int_{\mathbb{R}_k^d / B_\kappa} \left[ \frac{k}{|k|} \cdot \nabla_k \hat{u} + \frac{H + d - \frac{1}{2}}{|k|} \hat{u} \right] e^{ik \cdot x} dk. $$

Integrating by parts

$$ Au = - \int_{\mathbb{R}_k^d / B_\kappa} \mathrm{div}_k \left( \frac{k}{|k|} e^{ik \cdot x} \right) \hat{u} + \int_{\mathbb{R}_k^d / B_\kappa} \frac{H + d - \frac{1}{2}}{|k|} \hat{u} e^{ik \cdot x} dk. $$

As one can easily verify with a few vector calculus identities

$$ (4.16) \qquad \mathrm{div}_k \left( \frac{k}{|k|} e^{ik \cdot x} \right) = \frac{ik \cdot x + (d-1)}{|k|} e^{ik \cdot x} $$

such that

$$ Au = - \int_{\mathbb{R}_k^d / B_\kappa} \frac{ik \cdot x + (d-1)}{|k|} \hat{u} e^{ik \cdot x} dk + \int_{\mathbb{R}_k^d / B_\kappa} \frac{H + d - \frac{1}{2}}{|k|} \hat{u} e^{ik \cdot x} dk $$

which gives the expression in (4.15). Similar computations hold for $A^*$.

<figure><img src="image_3.png" /></figure>

Several recent works [17, 11, 21] show that linear equations involving operators of degree 0 may be useful in showcasing behavior related to turbulence, such as loss of regularity at long time. In [17], Yves Colin de Verdière and Laure Saint-Raymond in [17] consider an equation $\partial_t u + iAu = f$ where $f$ is an $L^2$ deterministic monochromatic forcing and $A$ is an homogeneous self-adjoint operator of degree 0 that satisfies some technical (but general) assumptions. The authors
show that, as $t \to \infty$, the solution resembles more and more the generalized eigenfunctions of the operator $A$, which do not have $L^2$ finite energy. The fact that $A$ is of order 0 is essential to ensure that its generalized eigenfunctions live in a space akin to $H^{-\frac{1}{2}-}$, which is less regular than $L^2$. R. Carles and C. Cheverry in [11] have also introduced an operator of degree 0 in the context of nuclear magnetic resonance. More precisely, they have shown that for some highly oscillatory source, at long times (i-e diffractive time), the solution produces constructive and destructive interferences which are interpreted as turbulent effects.

It is therefore interesting to compare our model with this existing literature. A first difference is that these other works are developed in deterministic setting. But even if one were to consider our model with a monochromatic deterministic source, it is hard to use similar techniques as those in [17]. Firstly, because our operator satisfies $D(A) \subsetneq D(A^*)$ and it cannot thus be self-adjoint or skew-adjoint. As a result, one cannot use Mourre's commutator theory as in [17] to prove a limiting amplitude principle which allows a nice spectral representation of the solution at infinite time. Secondly, the principal symbol of our operator $A$ (see (4.15)) grows with respect to $x$. This is not a standard pseudo-differential operator<fnref n="5" /> of degree 0. A more profound study of this special operator (4.2) is therefore interesting, and it is postponed for a future work. Nevertheless one can mention some conserved quantities in equation (4.4) when the source is null and initial datum is considered:

* conservation of volume for all $H$, i-e

$$ \partial_t \int_{\mathbb{R}_x^d} u(t, x) \, dx = 0, $$

* conservation of $L^2$-norm for $H = -d/2$, i-e

$$ \partial_t \int_{\mathbb{R}_x^d} |\hat{u}(t, x)|^2 \, dx = 0, $$

* conservation of trace at origin $H = -1/2$, i-e

$$ \partial_t u|_{x=0} = 0. $$

## 4.1. Well-posedness.

4.1.1. Uniqueness. Suppose that we have two solutions which live a.s. in $L_{\text{loc}}^1([0, \infty)_t \times \mathbb{R}^d)$ satisfying (4.9). When tested against test functions whose Fourier transform is supported in $B(0, \kappa)$, both solutions are null and thus they agree. When tested against $g \in C^\infty((0, \infty), D(A^*) \cap S(\mathbb{R}_x^d))$, their difference $v$ a.s. satisfies

$$ (4.17) \qquad \int_{\mathbb{R}_x^d} v(t, x)g(t, x)dx + \int_0^t \int_{\mathbb{R}_x^d} v(s, x)(-\partial_s + A^*)g(s, x)dsdx = 0 \quad \text{for all } t > 0. $$

If $g_0 \in D(A^*) \cap S(\mathbb{R}_x^d)$ then

$$ g_t(s, \cdot) = \mathcal{F}^{-1} \left[ |k|^{-H-\frac{1}{2}} \hat{g}_0 \left( \frac{|k|-c(s-t)}{|k|}k \right) \right] \in C^\infty((0, \infty), D(A^*) \cap S(\mathbb{R}_x^d)) $$

satisfies

$$ (-\partial_s + A^*)g_t = 0 \quad \text{and} \quad g_t|_{s=t} = g_0. $$

$$ |\partial_x^\alpha \partial_k^\beta a(x, k)| \lesssim (1 + |k|)^{m-|\beta|} $$

for all multi-index $\alpha$ and $\beta$.

<footnote n="5">⁵A standard operator of degree $m$ is an operator whose symbol $a \in C^\infty(\mathbb{R}^{2d})$ satisfies</footnote>
Using this choice of test function, (4.17) yields that a.s. for all $t > 0$

$$\int_{\mathbb{R}_x^d} v(t, x)g_0(x)dx = 0.$$

As a result, almost surely, for all $t$, $v(t, \cdot) = 0$ almost everywhere. More precisely, for any two solutions $u_1, u_2$ of our problem, we have:

$$\mathbb{P}[\forall t > 0, u_1(t, x) = u_2(t, x) \text{ almost everywhere in } x] = 1.$$

**Remark 4.7.** Note that one could obtain a stronger version of uniqueness by slightly changing the definition of weak solution. For instance, one could add the condition that $\mathbb{E}|u(t, x)|^2 \in L_{\text{loc}}^1((0, \infty)_t \times \mathbb{R}_x^d)$. Under this new definition one could show that any two solutions $u_1$ and $u_2$ must satisfy

$$\mathbb{E}|u_1(t, x) - u_2(t, x)| = 0$$

$(t, x)$-almost everywhere. This follows directly from the above argument to prove uniqueness together with the fact that the variance is a.e. finite.

4.1.2. Existence of a weak solution. Let us start by giving an expression for $e^{-tA}[\tau_y\varphi](x)$, which will enable us to prove that (4.12) does provide a weak solution to (4.4).

**Proposition 4.8.** For all $\varphi \in X \cap S(\mathbb{R}_x^d)$, the function $G(t, x, y) := e^{-tA}[\tau_y\varphi](x)$ given by

$$ (4.18) \qquad \begin{aligned} G(t, x, y) &= \int_{\mathbb{R}_k^d} \chi_{|k|>ct+\kappa} \left( \frac{|k|-ct}{|k|} \right)^{H+d-\frac{1}{2}} \hat{\varphi} \left( \frac{|k|-ct}{|k|}k \right) e^{-i2\pi \left( \frac{|k|-ct}{|k|} \right)k \cdot y + i2\pi k \cdot x} dk \\ &= \int_{\mathbb{R}_k^d} \left( \frac{|k|}{|k|+ct} \right)^{H+\frac{1}{2}} \hat{\varphi}(k) e^{i2\pi \left( \frac{|k|+ct}{|k|} \right)k \cdot x - i2\pi k \cdot y} dk. \end{aligned} $$

is a smooth function in all variables. Moreover:

* For each multi-index $\beta$ and each fixed $y \in \mathbb{R}^d$, $\partial_y^\beta G(\cdot, \cdot, y) \in C^1([0, \infty)_t, X \cap S(\mathbb{R}_x^d))$. Additionally, any Schwartz seminorm of $\partial_y^\beta G(\cdot, \cdot, y)$ (with respect to $x$) is locally uniform in $y$.

* For each $x \in \mathbb{R}^d$, $\partial_x^\beta G(\cdot, x, \cdot) \in C^1([0, \infty)_t, S(\mathbb{R}_y^d))$. Additionally, any Schwartz seminorm of $\partial_x^\beta G(\cdot, x, \cdot)$ (with respect to $y$) is locally uniform in $x$.

Finally, its Fourier transform in $y$ is given by

$$(\mathcal{F}_2 G)(t, x, k) = \left( \frac{|k|}{|k| + ct} \right)^{H + \frac{1}{2}} \overline{\hat{\varphi}(k)} e^{-i2\pi \left( \frac{|k| + ct}{|k|} \right)k \cdot x}.$$

*Proof.* The expression (4.18) comes from the definition of $e^{-tA}$ in (4.6), and the change of variables $k \mapsto (|k| + ct)/|k|$ $k$. The second expression directly yields the formula for $(\mathcal{F}_2 G)(t, x, k)$.

The fact that $G$ is smooth in all three variables follows easily from the integral expression and the fact that we can interchange differentiation and the integral. The latter step will be justified next.

Let us show that for any multi-index $\beta$ and any fixed $x$, $\partial_x^\beta G(\cdot, x, \cdot) \in C^1([0, \infty)_t, S(\mathbb{R}_y^d))$ (the case of fixed $y$ follows a similar proof). Note that any derivative $\partial_y^n \partial_x^\beta G(t, x, y)$ results in a factor
of $(-i2\pi k)^n$, which can be absorbed by $\hat{\varphi}(k)$, and thus smoothness is no problem. Let us focus on decay. For $|y| \le 1$, note that

$$|\partial_x^\beta G(t, x, y)| \le (2\pi)^{|\beta|} \int_{\mathbb{R}_k^d} \left( \frac{|k| + ct}{|k|} \right)^{|\beta|} |\hat{\varphi}(k)| dk.$$

Recall that the support of $\varphi$ guarantees that there are not problems of integration around $k = 0$.

For $|y| \ge 1$ and for any $j \in \{1, \dots, d\}$,

$$ \begin{aligned} |(-iy_j)^n \partial_x^\beta G(t, x, y)| &= \left| \int_{\mathbb{R}_k^d} \left( \frac{|k|}{|k| + ct} \right)^{H + \frac{1}{2} - |\beta|} \hat{\varphi}(k) e^{i2\pi \left( \frac{|k| + ct}{|k|} \right) k \cdot x} \partial_{k_j}^n e^{-i2\pi k \cdot y} dk \right| \\ &\le (2\pi)^{|\beta|} \int_{\mathbb{R}_k^d} \left| \partial_{k_j}^n \left[ \left( \frac{|k|}{|k| + ct} \right)^{H + \frac{1}{2} - |\beta|} \hat{\varphi}(k) e^{i2\pi \left( \frac{|k| + ct}{|k|} \right) k \cdot x} \right] \right| dk \end{aligned} $$

The integrand of the last expression can be shown to be absolutely integrable on the support of $\hat{\varphi}$. Note, however, that one needs to pay some powers of $x$ (which is fixed) to control this term. The implicit constant will therefore depend on such powers, but it is locally uniform in $x$ (i.e. the same constant for all $x$ in a compact set).

This concludes the proof that for each fixed $y \in \mathbb{R}^d$, $G(\cdot, \cdot, y) \in C([0, \infty)_t, X \cap S(\mathbb{R}_x^d))$. In order to improve this space to $C^1([0, \infty)_t, X \cap S(\mathbb{R}_x^d))$, it suffices to differentiate $G$ in time and repeat the above procedure.

An analogous argument based on integration by parts shows that for each fixed $y \in \mathbb{R}^d$ and each multi-index $\beta$, $\partial_y^\beta G(\cdot, \cdot, y) \in C^1([0, \infty)_t, X \cap S(\mathbb{R}_x^d))$. $\square$

Using these nice properties of $G$, one obtains the following corollary.

**Corollary 4.9.** The function $u(t, x)$ defined in (4.12) is a well-defined Gaussian process. For any $0 < \alpha < 1/2$, it has a.s. $\alpha$-Holder continuous paths with respect to $t$ and a.s. smooth paths with respect to $x$.

In particular, for any multi-index $\beta$, any $t_1, t_2 > 0$ and any $x_1, x_2 \in \mathbb{R}^d$ there exists a (locally uniform) constant $C > 0$ such that

$$ (4.19) \qquad \mathbb{E}|\partial_x^\beta u(t_1, x_1) - \partial_x^\beta u(t_2, x_2)|^2 \le C (|t_1 - t_2| + |x_1 - x_2|^2). $$

*Proof.* Let $G$ be as in (4.18), then $u$ can be written as

$$ u(t, x) = \int_0^t \int_{\mathbb{R}_y^d} G(t - s, x, y) \, dW(s, y) $$

By Proposition 4.8, for any fixed $(t, x)$, $G(\cdot, x, \cdot) \in C^1([0, \infty)_s, X \cap S(\mathbb{R}_x^d)) \subset L^2([0, t]_s \times \mathbb{R}_y^d)$ and therefore $u$ is a well-defined Itô process with zero mean and variance

$$ \mathbb{E}|u(t, x)|^2 = \int_0^t \int_{\mathbb{R}_y^d} |G(s, x, y)|^2 \, dy \, ds. $$

Next let us study the regularity. Firstly, note that for any multi-index $\beta$

$$ \partial_x^\beta u(t, x) = \int_0^t \int_{\mathbb{R}_y^d} \partial_x^\beta G(t - s, x, y) \, dW(s, y). $$
A simple calculation shows that for $t_1 < t_2$,

$$ (4.20) \qquad \begin{aligned} \mathbb{E}|\partial_x^\beta u(t_1, x) - \partial_x^\beta u(t_2, x)|^2 &= \int_0^{t_1} \int_{\mathbb{R}_y^d} |\partial_x^\beta G(t_2 - s, x, y) - \partial_x^\beta G(t_1 - s, x, y)|^2 dy ds \\ &\quad + \int_{t_1}^{t_2} \int_{\mathbb{R}_y^d} |\partial_x^\beta G(t_2 - s, x, y)|^2 dy ds \\ &\le \int_0^{t_1} \int_{\mathbb{R}_y^d} |\partial_x^\beta G(t_2 - t_1 + s, x, y) - \partial_x^\beta G(s, x, y)|^2 dy ds \\ &\quad + |t_1 - t_2| \left\| \partial_x^\beta G(\cdot, x, \cdot) \right\|_{L^\infty([0,t_2], L^2(\mathbb{R}_y^d))}^2 \\ &\lesssim t_1 |t_1 - t_2|^2 \left\| \partial_x^\beta \partial_s G(\cdot, x, \cdot) \right\|_{L^\infty([0,t_2], L^2(\mathbb{R}_y^d))}^2 \\ &\quad + |t_1 - t_2| \left\| \partial_x^\beta G(\cdot, x, \cdot) \right\|_{L^\infty([0,t_2], L^2(\mathbb{R}_y^d))}^2 \end{aligned} $$

using the fact that $\partial_x^\beta G(\cdot, x, \cdot) \in C^1([0, \infty)_s, X \cap S(\mathbb{R}_x^d))$.

Similarly, for any $x_1, x_2$, we have that

$$ (4.21) \qquad \begin{aligned} \mathbb{E}|\partial_x^\beta u(t, x_1) - \partial_x^\beta u(t, x_2)|^2 &= \int_0^t \int_{\mathbb{R}_y^d} |G(s, x_1, y) - G(s, x_2, y)|^2 dy ds \\ &= \int_0^t \int_{\mathbb{R}_y^d} \left| \int_0^1 (x_2 - x_1) \cdot \nabla_x G(s, (1-\lambda)x_1 + \lambda x_2, y) d\lambda \right|^2 dy ds \\ &\le |x_2 - x_1|^2 \int_0^t \int_{\mathbb{R}_y^d} \int_0^1 |\nabla_x G(s, (1-\lambda)x_1 + \lambda x_2, y)|^2 d\lambda dy ds \\ &\le |x_2 - x_1|^2 \int_0^1 \|\nabla_x G(s, (1-\lambda)x_1 + \lambda x_2, y)\|_{L^2([0,t] \times \mathbb{R}_y^d)}^2 d\lambda \end{aligned} $$

The latter factor is finite since $\partial_x^\beta G(\cdot, x, \cdot) \in C^1([0, \infty)_t, S(\mathbb{R}_y^d))$ locally uniformly in $x$.

One may use (4.20) and (4.21) to obtain (4.19). Given that $\partial_x^\beta u(t, x)$ is a Gaussian process, (4.19) immediately implies that for any $m \in \mathbb{N}$,

$$ \mathbb{E}|\partial_x^\beta u(t_1, x_1) - \partial_x^\beta u(t_2, x_2)|^{2m} \le C_m (|t_1 - t_2|^m + |x_1 - x_2|^{2m}). $$

The Kolmogorov continuity theorem guarantees that $\partial_x^\beta u(t, x)$ is therefore a.s. $\alpha$-Hölder continuous in time for any $0 < \alpha < (m-d)/(2m)$. $\square$

We are ready to show that the function $u$ defined in (4.12) is indeed a weak solution to (4.4) in the sense given by Definition 4.3.

*Proof that $u$ is a weak solution. Step 1.* We start by checking (4.10). Choose any $g \in C^\infty((0, \infty), S(\mathbb{R}_x^d))$ with $\hat{g}(t, k) = 0$ for all $|k| > \kappa$. Then

$$ (4.22) \qquad \int_0^t \int_{\mathbb{R}_x^d} \left( \int_0^s \int_{\mathbb{R}_y^d} G(s-s', x, y) dW(s', y) \right) g(s,x) ds dx = \\ \qquad \int_0^t \int_0^s \int_{\mathbb{R}_y^d} \left( \int_{\mathbb{R}_x^d} G(s-s', x, y) g(s,x) dx \right) dW(s', y) ds, $$
In order to justify this, it suffices to show that $G(s - s', x, y) g(s, x)$ is absolutely integrable with respect to the Lebesgue measure in $x \in \mathbb{R}^d$ and $s \in [0, t]$, and with respect to $d\tilde{W}(s', y)$ for $s' \in [0, t]$ and $y \in \mathbb{R}^d$.

As proved in Proposition 4.8, one may show that $|G(s - s', x, y)| \le C(s - s', x) \langle y \rangle^{-d-1}$ where $C(s - s', x)$ depends on $s - s'$ and $x$ in a polynomial way. To compensate for this growth in $x$, recall that for any $N$, $|g(s, x)| \lesssim_N \langle x \rangle^{-N}$. As a result, $|G(s - s', x, y) g(s, x)| \le C(t, N) \langle y \rangle^{-d-1} \langle x \rangle^{-N}$. This guarantees that we may integrate in whichever order we prefer.

In order to check (4.10), it suffices to show that for all $s' \in [0, s]$, $s \in [0, t]$ and $y \in \mathbb{R}^d$,

$$ (4.23) \qquad \int_{\mathbb{R}_x^d} G(s - s', x, y) g(s, x) dx = 0. $$

This follows from the Plancherel theorem, which allows us to rewrite the integral in (4.23) as

$$ \int_{\mathbb{R}_k^d} \hat{G}(s - s', k, y) \hat{g}(s, k) dk = \int_{\mathbb{R}_k^d} \chi_{|k|>c(s-s')+\kappa} \left( \frac{|k|-c(s-s')}{|k|} \right)^{H+d-\frac{1}{2}} \hat{\varphi} \left( \frac{|k|-c(s-s')}{|k|} k \right) e^{-i2\pi \left( \frac{|k|-c(s-s')}{|k|} \right) k \cdot y} \hat{g}(s, k) dk. $$

Note that the support of $\hat{g}(s, \cdot)$ and that of $\hat{\varphi}$ do not overlap. Therefore this integral is zero.

**Step 2.** Next we check (4.9). Fix some $g \in C^\infty((0, \infty), D(A^*) \cap S(\mathbb{R}_x^d))$. After exchanging the order of integration (which may be justified using Proposition 4.8 as in Step 1), (4.9) is equivalent to:

$$ (4.24) \qquad \int_{\mathbb{R}_x^d} G(t-s', x, y) g(t, x) dx + \int_{s'}^t \int_{\mathbb{R}_x^d} G(s-s', x, y) (-\partial_s + A^*) g(s, x) dx ds \\ = \int_{\mathbb{R}_x^d} \varphi(x-y) g(s', x) dx $$

for all $s' \in [0, t]$ and all $y \in \mathbb{R}^d$. Next we may use the pre-dual of $-\partial_s + A^*$ to rewrite the second term on the left-hand side. After using the fact that $G(0, x, y) = \varphi(x-y)$, we find that (4.24) holds if and only if:

$$ (4.25) \qquad \int_{s'}^t \int_{\mathbb{R}_x^d} (\partial_s + A) G(s - s', x, y) g(s, x) dx ds = 0 \quad \forall s' \in [0, t], y \in \mathbb{R}^d. $$

Therefore, it suffices to prove that $(\partial_t + A)G(t, x, y) = 0$. Taking the Fourier transform in $x$, which is allowed by Proposition 4.8, it suffices to check that:

$$ (4.26) \qquad \partial_t \hat{G}(t, k, y) + c \operatorname{div}_k \left( \frac{k}{|k|} \hat{G}(t, k, y) \right) + c \frac{H+d-\frac{1}{2}}{|k|} \hat{G}(t, k, y) = 0 $$

One may easily check that (4.26) does indeed hold by direct calculation using the fact that:

$$ \hat{G}(t, k, y) = \chi_{|k|>ct+\kappa} \left( \frac{|k|-ct}{|k|} \right)^{H+d-\frac{1}{2}} \hat{\varphi} \left( \frac{|k|-ct}{|k|} k \right) e^{-i2\pi \left( \frac{|k|-ct}{|k|} \right) k \cdot y}. $$

□
4.2. **Correlations for $H \in (0, 1)$**. Our first result is an explicit formula for the correlations of our solution at any finite time. Note that this result holds regardless of the value of $H$, because $t > 0$ is finite. However, we will only be able to take the limit (as a function) when $H \in (0, 1)$.

**Proposition 4.10.** For any finite time, the correlation function satisfies

$$ \mathbb{E}[u(t, x_1)u(t, x_2)] = \int_{\mathbb{R}_k^d} |k|^{-(2H+d)} F(t, |k|) e^{i2\pi k(x_1-x_2)} dk $$

where $F$ is the function given in Remark 3.3.

*Proof.* For any finite time $t > 0$, one has

$$ \mathbb{E}[u(t, x_1)u(t, x_2)] = \int_0^t \int_{\mathbb{R}^d} e^{-(t-s)A}[\tau_y \varphi](x_1)e^{-(t-s)A}[\tau_y \varphi](x_2) dy ds. $$

The change of variables $s \mapsto t-s$ yields

$$ \mathbb{E}[u(t, x_1)u(t, x_2)] = \int_0^t \int_{\mathbb{R}^d} e^{-sA}[\tau_y \varphi](x_1)e^{-sA}[\tau_y \varphi](x_2) dy ds. $$

Using Proposition 4.8 together with the Plancherel theorem (in the $y$ variable), we obtain

$$ \mathbb{E}[u(t, x_1)u(t, x_2)] = \int_0^t \int_{\mathbb{R}^d} \left( \frac{|k|}{|k| + cs} \right)^{2H+1} |\hat{\varphi}(k)|^2 e^{-i2\pi \left( \frac{|k|+cs}{|k|} \right) k \cdot (x_1 - x_2)} dk \, ds. $$

Next we perform the change of variables $k \mapsto \frac{|k|}{|k|+cs} k$, which yields

$$ \mathbb{E}[u(t, x_1)u(t, x_2)] = \int_0^t \int_{\mathbb{R}^d} \chi_{|k|>cs+\kappa} \left( \frac{|k|-cs}{|k|} \right)^{2H+d} \left| \hat{\varphi} \left( \frac{|k|-cs}{|k|} k \right) \right|^2 e^{-i2\pi k \cdot (x_1-x_2)} dk \, ds $$

$$ = \int_{\mathbb{R}^d} e^{-i2\pi k \cdot (x_1 - x_2)} |k|^{-2H-d} \left( \int_0^t \chi_{|k|>cs+\kappa} (|k|-cs)^{2H+d} \left| \hat{\varphi} \left( \frac{|k|-cs}{|k|} k \right) \right|^2 ds \right) dk. $$

This concludes the proof.

<figure><img src="image_4.png" /></figure>

As a corollary, we take the limit as $t \to \infty$ and show points (ii) and (iv) of Theorem 4.4.

**Corollary 4.11.** Suppose that $H \in (0, 1)$. Then $u(t, x)$ converges in law to a zero-mean Gaussian field $u_\infty(x)$ with the following correlation structure:

$$ \mathbb{E}[u_\infty(x_1)u_\infty(x_2)] = \lim_{t \to \infty} \mathbb{E}[u(t, x_1)u(t, x_2)] $$

Moreover

$$ \mathbb{E}[u(t, x_1)u(t, x_2)] = C(d, H) K_H(x_1 - x_2) - (K_H * J_H)(x_1 - x_2) + O((ct)^{-2H}) $$

where

$$ (4.27) \qquad K_H := \mathcal{F}^{-1}\left[\chi_{|k|>\kappa} |k|^{-(2H+d)}\right] \quad \text{and} \quad \mathcal{J}_H := \mathcal{F}^{-1}\left[\chi_{|k|>\kappa} \Psi_{d,H}(|k|)\right] $$

and where $C(d, H)$ and $\Psi_{d,H}$ are given in Theorem 3.2.

*Proof.* Defining

$$ v(t, x) = \int_0^t \int_{\mathbb{R}_y^d} G(s, x, y) \, dW(s, y), $$

we note that for fixed $t$, $v(t, x)$ and $u(t, x)$ are Gaussian random variables with the same law, since the mean and the correlations coincide.
First we show that

$$v_{\infty}(x) := \int_{0}^{\infty} \int_{\mathbb{R}_{y}^{d}} G(s, x, y) \, dW(s, y)$$

is a well-defined Gaussian field. This is due to the fact that for fixed $x \in \mathbb{R}^d$, the function $G(\cdot, x, \cdot) \in L^2([0, \infty]_s \times \mathbb{R}_y^d)$. Indeed, following the arguments in Proposition 4.10 one easily finds that

$$\int_{0}^{\infty} \int_{\mathbb{R}_{y}^{d}} |G(s, x, y)|^{2} \, dy \, ds = \int_{0}^{\infty} \int_{\mathbb{R}_{d}} \left( \frac{|k|}{|k| + cs} \right)^{2H+1} |\hat{\varphi}(k)|^{2} \, dk \, ds.$$

The convergence of $v(t)$ to $v_{\infty}$ is now standard. Note that:

$$v_{\infty}(x) - v(t, x) = \int_{t}^{\infty} \int_{\mathbb{R}_{y}^{d}} G(s, x, y) \, dW(s, y),$$

and therefore

$$\begin{align*}
\mathbb{E}|v_{\infty}(x) - v(t, x)|^2 &= \int_t^\infty \int_{\mathbb{R}_y^d} |G(s, x, y)|^2 dy \, ds \\
&= \int_t^\infty \int_{\mathbb{R}^d} \left( \frac{|k|}{|k| + cs} \right)^{2H+1} |\hat{\varphi}(k)|^2 dk \, ds \\
&= \frac{1}{2cH} \int_{\mathbb{R}^d} (|k| + ct)^{-2H} |k|^{2H+1} |\hat{\varphi}(k)|^2 dk
\end{align*}$$

The same argument allows us to take the limit as $t \to \infty$ for the two-point correlations of $v(t)$, which coincide with the two-point correlations of $u(t)$ in Proposition 4.10, and which yield the correlations of $u_{\infty}$. The fact that we can exchange the limit and the expectation admits the same argument as in Step 2 in the proof of Theorem 2.4. Indeed, note that the variance $\mathbb{E}|u_{\infty}(x)|^2$ is finite and independent of $x \in \mathbb{R}^d$.

The rate of convergence of the correlations as $t \to \infty$ depend directly on the function $F$ in Proposition 4.10. A careful analysis shows that the rate of convergence is given by a multiple of $(ct)^{-2H}$. $\square$

Now that we know that, whenever $H \in (0, 1)$, $u_{\infty}$ is a well-defined Gaussian field indexed by $x \in \mathbb{R}^d$, we may wonder about its continuity and its regularity. In this direction, we have the following result:

**Corollary 4.12.** Suppose that $H \in (0, 1)$. The limit $u_{\infty} = \lim_{t \to \infty} u(t)$ is a Gaussian field in $x \in \mathbb{R}^d$ which has a modification with a.s. $\alpha$-Hölder continuous paths for any $0 < \alpha < H$. More precisely, for all $x, \ell \in \mathbb{R}^d$, the variance of the increment satisfies

$$(4.28) \qquad \mathbb{E}|\delta_{\ell}u_{\infty}(x)|^2 \le \frac{C_K}{H(1-H)}|\ell|^{2H} + C_J|\ell|^2,$$

where $C_K$ and $C_J$ only depend of $H$ and $d$ without blowing up when $H$ tends to 0 or 1.

*Proof.* By (4.13), we have that

$$(4.29) \qquad \mathbb{E}|\delta_{\ell}u_{\infty}(x)|^2 = 2C(d, H)[\mathcal{K}_H(0) - \mathcal{K}_H(\ell)] - 2[(\mathcal{K}_H * \mathcal{J}_H)(0) - (\mathcal{K}_H * \mathcal{J}_H)(\ell)].$$

**Step 1.** The function $(\mathcal{K}_H * \mathcal{J}_H)$ in (4.27) is smooth, since $\Psi_{d,H}(|k|)$ decays as rapidly as desired in $k$. Moreover, we may write:

$$(4.30) \qquad 1 - e^{2\pi ik \cdot \ell} = 1 - \cos(2\pi k \cdot \ell) - i \sin(2\pi k \cdot \ell).$$
Then note that

$$ \int_{\mathbb{R}_k^d} \sin(2\pi k \cdot \ell) \chi_{|k|>\kappa} |k|^{2-(2H+d)} |\Psi_{d,H}(|k|)| dk = 0 $$

thanks to the fact that $\chi_{|k|>\kappa} |k|^{-(2H+d)} \Psi_{d,H}(|k|)$ is rotationally invariant together with the change of variables $k \mapsto -k$.

As a result,

$$ \begin{aligned} |(\mathcal{K}_H * \mathcal{J}_H)(0) - (\mathcal{K}_H * \mathcal{J}_H)(\ell)| &= \int_{\mathbb{R}_k^d} [1 - \cos(2\pi k \cdot \ell)] \chi_{|k|>\kappa} |k|^{-(2H+d)} |\Psi_{d,H}(|k|)| dk \\ &\le \int_{\mathbb{R}_k^d} \frac{(2\pi |k|\ell)^2}{2} \chi_{|k|>\kappa} |k|^{-(2H+d)} |\Psi_{d,H}(|k|)| dk \\ &\le 2\pi^2 |\ell|^2 \int_{\mathbb{R}_k^d} |k|^{2-(2H+d)} |\Psi_{d,H}(|k|)| dk. \end{aligned} $$

The latter integral is finite thanks to the rapid decay of $\Psi_{d,H}(|k|)$. This takes care of the second term on the right-hand side of (4.29).

**Step 2.** The leading term, for the purpose of studying the Hölder regularity, is the first term on the right-hand side of (4.29). We separate the integral into two parts:

$$ \begin{aligned} \mathcal{K}_H(0) - \mathcal{K}_H(\ell) &= \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} (1 - e^{2\pi i k \cdot \ell}) |k|^{-(2H+d)} dk \\ (4.31) \qquad \qquad \qquad \qquad &= \int_{\mathbb{R}_k^d} \chi_{|k|<r} \chi_{|k|>\kappa} (1 - e^{2\pi i k \cdot \ell}) |k|^{-(2H+d)} dk \\ &\quad + \int_{\mathbb{R}_k^d} \chi_{|k|>r} \chi_{|k|>\kappa} (1 - e^{2\pi i k \cdot \ell}) |k|^{-(2H+d)} dk. \end{aligned} $$

We choose $r = (2\pi|\ell|)^{-1}$; one can carefully check that this choice is optimal. Let us further assume that $|\ell| < \kappa$, the alternative scenario is easier and can be studied separately.

The second term on the right-hand side of (4.31) yields:

$$ \left| \int_{\mathbb{R}_k^d} \chi_{|k|>r} \chi_{|k|>\kappa} (1 - e^{2\pi i k \cdot \ell}) |k|^{-(2H+d)} dk \right| \le \left| \int_{\mathbb{R}_k^d} \chi_{|k|>r} |k|^{-(2H+d)} dk \right|. $$

Then, changing to spherical coordinates in the variable $k$,

$$ \left| \int_{\mathbb{R}_k^d} \chi_{|k|>r} |k|^{-(2H+d)} dk \right| \lesssim \left| \int_{|k|>r} |k|^{-(2H+1)} d|k| \right| \lesssim \frac{1}{H} r^{-2H} \lesssim \frac{1}{H} |\ell|^{2H}. $$

The first term on the right-hand side of (4.31) requires more care. Using (4.30) as before, note that

$$ \int_{|k|<r} \sin(2\pi k \cdot \ell) \chi_{|k|>\kappa} |k|^{-(2H+d)} dk = 0, $$

as is easily seen from the change of variables $k \mapsto -k$. Therefore it suffices to bound

$$ \begin{aligned} \left| \int_{|k|<r} [1 - \cos(2\pi k \cdot \ell)] \chi_{|k|>\kappa} |k|^{-(2H+d)} dk \right| &\le \int_{|k|<r} \frac{(2\pi |k||\ell|)^2}{2} \chi_{|k|>\kappa} |k|^{-(2H+d)} dk \\ &\lesssim 2\pi^2 |\ell|^2 \int_{\kappa}^{r} |k|^{1-2H} d|k| \\ &\lesssim \frac{1}{2-2H} |\ell|^{2H}. \end{aligned} $$
This finishes the proof of (4.28). Using this bound on the variance, a standard argument based on the Gaussianity of $u_\infty$ and the Kolmogorov continuity theorem show that $u_\infty$ has a.s. $\alpha$-Hölder continuous paths for any $0 < \alpha < H$. $\square$

4.3. **Correlations for $H \in [-d/2, 0]$**. For values of $H$ in this range, we will still show that $u_\infty = \lim_{t \to \infty} u(t)$ exists. However, $u_\infty$ cannot be interpreted as a function, in fact its variance is not finite. Instead, $u_\infty$ is a random distribution acting on a space of test functions $X \cap S(\mathbb{R}^d)$. Our first result is the analogue of Proposition 4.10 in the context of distributions.

The rate of convergence depends on some integrability and differentiability properties of the test functions chosen. If one tests against Schwartz functions in $X \cap S(\mathbb{R}^d)$, convergence will happen faster than $t^{-n}$ for any $n \in \mathbb{N}$ (as we will show below). However, one may want to test agains less regular or less decaying functions, in which case this rate of convergence can be quantified.

Before doing so, one needs to identify a good space of test functions for which “testing” is still well-defined. Given that $u(t, x)$ does not generally decay in the $x$-variable, one needs to impose a certain decay on the test functions. For $n \in \mathbb{N}$, let us define the space$^6$ $H^{n,n}(\mathbb{R}^d)$, which is the closure of $S(\mathbb{R}^d)$ with respect to the norm

$$ \|g\|_{H^{n,n}(\mathbb{R}^d)} = \sum_{j=0}^{n} \|\langle x \rangle^j g\|_{H^{n-j}(\mathbb{R}^d)} $$

It is not hard to show that $H^{n,n}(\mathbb{R}^d)$ is a Banach space algebra living inside $L^2(\mathbb{R}_x^d)$, and thus the Fourier transform is well-defined for such functions. Moreover, $g \in H^{n,n}(\mathbb{R}_x^d)$ if and only if $\hat{g} \in H^{n,n}(\mathbb{R}_k^d)$ (and both norms are comparable).

**Proposition 4.13.** Let $H \in [-d/2, 0]$. For any $t > 0$ and any test functions $g_1, g_2 \in H^{n,n}(\mathbb{R}^d)$ with $3(d+1)/2 < n$, the correlations satisfy

$$ \begin{align} \mathbb{E}[\langle u(t), g_1 \rangle \langle u(t), g_2 \rangle] &= C(d, H) \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \nonumber \\ &\quad - \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} \Psi_{d,H}(|k|) \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \nonumber \\ &\quad + \int_{\mathbb{R}_k^d} \chi_{|k|>ct+\kappa} |k|^{-(2H+d)} [\Psi_{d,H}(|k| - ct) - \Psi_{d,H}(0)] \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \tag{4.32} \\ &= C(d, H) \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \nonumber \\ &\quad - \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} \Psi_{d,H}(|k|) \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \nonumber \\ &\quad + C(d, H) \left(\frac{1}{ct}\right)^{-(2H+d)-2n} \|g_1\|_{H^{n,n}} \|g_2\|_{H^{n,n}} \nonumber \end{align} $$

where $C(d, H)$ and $\Psi_{d,H}$ are given in Theorem 3.2, and $F$ is given by (3.3).

**Remark 4.14.** If one chooses test functions $g_1, g_2 \in X \cap S(\mathbb{R}^d)$, then the above result yields a rate of convergence proportional to $(ct)^{-(2H+d)-2n}$ for $n$ as large as desired.

<footnote>$^6$Let us mention that a better space might be possible, in the sense that one might require less decay in physical space. However, this requires a careful study of the right functional space for $u(t, x)$ and $G(s, x, y)$. We leave this study for future research.</footnote>
**Proof. Step 1.** Firstly, note that $\langle u(t), g_1 \rangle$ is well-defined for each $t > 0$. This follows after a similar argument to the one that justifies (4.22), i.e. the fact that $G(s, x, y)g_1(x)$ is absolutely integrable with respect to the Lebesgue measure in $x$ and with respect $dW(s, y)$. Indeed, we showed that $|G(s, x, y)| \lesssim_t \langle x \rangle^{d+1} \langle y \rangle^{-d-1}$, so one needs only show that $\langle x \rangle^{d+1} |g_1(x)| \in L^1(\mathbb{R}_x^d)$. But this follows from the choice of $n$ together with the Cauchy-Schwarz inequality:

$$ \int_{\mathbb{R}^d} \langle x \rangle^{d+1} |g_1(x)| dx \lesssim \|g_1\|_{H^{n,n}} \left( \int_{\mathbb{R}^d} \langle x \rangle^{-d-1} dx \right)^{1/2}. $$

By Proposition 4.10, we have:

$$ \begin{aligned} \mathbb{E}[\langle u(t), g_1 \rangle \langle u(t), g_2 \rangle] &= \int_{\mathbb{R}_{x_1}^d \times \mathbb{R}_{x_2}^d} \left( \int_{\mathbb{R}_k^d} |k|^{-2H-d} F(t, |k|) e^{i2\pi k(x_1-x_2)} dk \right) g_1(x_1) g_2(x_2) dx_1 dx_2 \\ &= \int_{\mathbb{R}_k^d} |k|^{-2H-d} F(t, |k|) \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \end{aligned} $$

Finally, we use the expression for $F(t, |k|)$ obtained in Remark 3.3. This yields (4.32).

**Step 2.** In order to obtain an asymptotic expansion in terms of $t$, note that it suffices to study the last term on the right-hand side of (4.32) which, using the definition of $\Psi_{d,H}$, can be rewritten as follows:

$$ - \int_{\mathbb{R}_k^d} \chi_{|k|>ct+\kappa} |k|^{-(2H+d)} \left( \int_{\kappa}^{|k|-ct} s^{2H+d} \psi(s) ds \right) \hat{g}_1(k) \overline{\hat{g}_2(k)} dk. $$

Note that generally, we cannot expect a better bound than

$$ \left| \chi_{|k|>ct+\kappa} \left( \int_{\kappa}^{|k|-ct} s^{2H+d} \psi(s) ds \right) \right| \le C(d, H) \lesssim 1. $$

Certainly, better bounds could be obtained should $\psi$ vanish on a large ball around zero, but this is not the physical setting we are interested in.

Therefore, by the Cauchy-Schwarz inequality,

$$ \begin{aligned} & \left| \int_{\mathbb{R}_k^d} \chi_{|k|>ct+\kappa} |k|^{-(2H+d)} \left( \int_{\kappa}^{|k|-ct} s^{2H+d} \psi(s) ds \right) \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \right| \\ & \le C(d, H) \int_{\mathbb{R}_k^d} \chi_{|k|>ct+\kappa} |k|^{-(2H+d)} |\hat{g}_1(k)| |\hat{g}_2(k)| dk \\ & \lesssim C(d, H) (ct)^{-(2H+d)-2n} \|g_1\|_{H^{n,n}} \|g_2\|_{H^{n,n}}. \end{aligned} $$

□

Thanks to the formula for the correlations obtained in Proposition 4.13, we may now study the asymptotic behavior of $u(t)$ as $t \to \infty$. That is the content of the next result.

**Corollary 4.15.** Let $H \in [-d/2, 0]$ and $n > 3(d+1)/2$. The solution $u(t)$ in (4.12) converges in law, as $t \to \infty$, to a random Gaussian measure acting on $X \cap H^{n,n}(\mathbb{R}^d)$. Moreover, the correlation
structure of the limiting Gaussian measure is given by:

$$ (4.33) \qquad \begin{aligned} \mathbb{E}[\langle u_{\infty}, g_1 \rangle \langle u_{\infty}, g_2 \rangle] &= \lim_{t \to \infty} \mathbb{E}[\langle u(t), g_1 \rangle \langle u(t), g_2 \rangle] \\ &= C(d, H) \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} \hat{g}_1(k) \overline{\hat{g}_2(k)} dk \\ &\quad - \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} |k|^{-(2H+d)} \Psi_{d,H}(|k|) \hat{g}_1(k) \overline{\hat{g}_2(k)} dk. \end{aligned} $$

for any $g_1, g_2 \in X \cap H^{n,n}(\mathbb{R}^d)$.

*Proof.* The convergence of $\langle u(t), g \rangle$ as $t \to \infty$ follows from a similar argument to the one in the proof of Corollary 4.11, so let us focus on proving (4.33). The second equality in (4.33) follows by taking $t \to \infty$ in (4.32). Finally, the first equality admits the same argument as Step 2 in the proof of Theorem 2.4, so we omit it. □

4.4. The case $H = -d/2$. In this case, note that the limiting Gaussian field $u_{\infty}$ obtained in Corollary 4.15 seems to differ from the usual white noise measure in two ways:

(i) the space of test functions $X \cap S(\mathbb{R}^d)$, which only admits functions whose Fourier transform vanishes in $|k| < \kappa$; and

(ii) the top order term in (4.33), which involves a cut-off to wavenumbers $|k| > \kappa$.

Point (ii) can be fixed by defining a continuous zero-mean Gaussian field $u_{\mathcal{I}}$ with correlations

$$ (4.34) \qquad \begin{aligned} \mathbb{E}[u_{\mathcal{I}}(x_1)u_{\mathcal{I}}(x_2)] &= \mathcal{J}_{-d/2}(x_1 - x_2) + \mathcal{F}^{-1}[\chi_{|k|\le\kappa}](x_1 - x_2) \\ &= \int_{\mathbb{R}_k^d} \chi_{|k|>\kappa} \Psi_{d,H}(|k|) e^{2\pi i k \cdot (x_1 - x_2)} dk + \mathcal{F}^{-1}[\chi_{|k|\le\kappa}](x_1 - x_2) \\ &=: \mathcal{I}(x_1 - x_2), \end{aligned} $$

where $\mathcal{J}$ was defined in (4.27). Note that $\mathcal{I}(x)$ is a continuous, radial, square-integrable function, and therefore they are much more regular than a delta function. For instance, in the one-dimensional case $d=1$, the second term on the right-hand side of (4.34) corresponds to the sinc function.

Using (4.34), we may rewrite the correlations of our solution $u_{\infty}$ in (4.33) as follows:

$$ (4.35) \qquad \begin{aligned} \mathbb{E}[\langle u_{\infty}, g_1 \rangle \langle u_{\infty}, g_2 \rangle] &= \mathbb{E}[\langle dW, g_1 \rangle \langle dW, g_2 \rangle] - \int_{\mathbb{R}^{2d}} \mathcal{I}(x_1 - x_2) g_1(x_1) g_2(x_2) dx_1 dx_2 \\ &= \int_{\mathbb{R}_x^d} g_1(x) g_2(x) dx - \int_{\mathbb{R}^{2d}} \mathcal{I}(x_1 - x_2) g_1(x_1) g_2(x_2) dx_1 dx_2 \end{aligned} $$

for any $g_1, g_2 \in X \cap S(\mathbb{R}_x^d)$. This solves point (ii).

Let us next focus on point (i), regarding the space of test functions. There is a natural way to regard $u(t)$ as a distribution in $S(\mathbb{R}^d)'$ (as opposed to $(X \cap S(\mathbb{R}^d))'$). Fix a test function $h \in S(\mathbb{R}^d)$ such that

$$ (4.36) \qquad \begin{cases} \hat{h}(k) = 0 & \text{if } |k| \ge \kappa, \\ \hat{h}(k) & \text{is radial and takes values in } [0, 1] \text{ for } |k| < \kappa. \end{cases} $$

For any test function $g \in S(\mathbb{R}^d)$, one may decompose

$$ g = \mathcal{L}_1 g + \mathcal{L}_2 g := \mathcal{F}^{-1}[\hat{g}\hat{h}] + \mathcal{F}^{-1}[\hat{g}(1 - \hat{h})]. $$
Then one may make a small modification to the notion of solution in Definition 4.3. Indeed, one could define a weak solution to (4.4) to be a stochastic process $u$ such that $u \in L_{\text{loc}}^1([0, \infty)_t \times \mathbb{R}^d)$ almost surely and such that for all $g \in C^\infty((0, \infty), \mathcal{S}(\mathbb{R}_x^d))$

$$ \int_{\mathbb{R}_x^d} u(t, x) \mathcal{L}_2 g(t, x) dx + \int_0^t \int_{\mathbb{R}_x^d} u(s, x) (-\partial_s + A^*) \mathcal{L}_2 g(s, x) ds dx \\ = \int_0^t \int_{\mathbb{R}_y^d} \left( \int_{\mathbb{R}_x^d} \varphi(x-y) \mathcal{L}_2 g(s, x) dx \right) dW(s, y) $$

and

$$ \int_0^t \int_{\mathbb{R}_x^d} u(s, x) \mathcal{L}_1 g(s, x) ds dx = 0. $$

Note that the solution found in (4.12) is still a solution with this new definition. Moreover, the solution does not depend on the choice of function $h$ as long as it satisfies (4.36).

This allows us to view $u(t)$ and $u_\infty$ as distributions in $\mathcal{S}(\mathbb{R}^d)'$, and thus to extend (4.35) to test functions in $\mathcal{S}(\mathbb{R}^d)$. Then the top order of (4.35) does indeed correspond to a white noise in the classical sense.

## 5. NUMERICAL SIMULATIONS

Our goal in this section is to perform numerical simulations to illustrate the theory presented in the previous sections. We recall the continuous problem (3.12):

$$ (5.1) \qquad \begin{cases} \partial_t \hat{u}(t, k) + \mathcal{L}(\hat{u})(t, k) + \mathcal{D}(\hat{u})(t, k) = \hat{f}(t, k) & \text{for } t > 0, k \in \mathbb{R}^d, |k| > \kappa \\ \hat{u}(t, k) = 0 & \text{for } t > 0, k \in \mathbb{R}^d, |k| \le \kappa, \\ \hat{u}(0, k) = 0 & \text{for } k \in \mathbb{R}^d, \end{cases} $$

where

$$ (5.2) \qquad \mathcal{L}(\hat{u})(t, k) := c \operatorname{div}_k \left( \frac{k}{|k|} \hat{u}(t, k) \right) \quad \text{and} \quad \mathcal{D}(\hat{u})(t, k) := \left( c \frac{H + \frac{1}{2}}{|k|} + (2\pi)^2 \nu |k|^2 \right) \hat{u}(t, k) $$

Here, $H$ is eventually the Hurst exponent of the solution (at infinite time), and $\nu > 0$ denotes viscosity. The introduction of viscosity is necessary in order to reach a statistically stationary state at a finite time, state in which a statistical analysis is possible. In the inviscid problem proposed in (3.1), on the other hand, scales as small as $1/ct$ are populated at time $t$, leading to numerical instabilities and a blow up when simulated in a finite periodic box with a finite resolution.

Our numerical method combines ideas from time predictor-corrector schemes [24] and pseudo-spectral methods [36]. In particular, this implies using the Discrete Fourier Transform (DFT), which forces certain choices regarding the discretization.

We will denote by $\hat{u}[t, k]$ the (discrete) vector whose continuous counterpart is $\hat{u}(t, k)$. We denote by $\Delta t > 0$ the time stepping and by $\Delta x > 0$ the mesh size. The mesh size is also the smallest accessible length scale.

### 5.1. Numerical method.
5.1.1. Discretization. We discretize our physical space by considering a discrete periodic box $x \in (\mathbb{Z}/N\mathbb{Z})^d$ of unit length $L_{tot} = 1$, using $N = 2^n$ collocation points in each direction, with $n \in \mathbb{N}$ where we adopt the convention that 0 is not a natural number. Therefore, the mesh size is

$$ \Delta x = L_{tot}/N. $$

The wave vector $k = (k_i)_{1\le i\le d}$ is discretized as

$$ k_i = [0, 1, \dots, N/2, -N/2 + 1, -N/2 + 2, \dots, -1] \Delta k, $$

where the spectral resolution is given by

$$ \Delta k = 1/L_{tot}. $$

This discretization is standard when using the Discrete Fourier Transform (DFT). The choice of starting from $k_i = 0$ as the first element of the array is dictated by the convention that is used to define the DFT and its inverse.

In order to discretize the divergence term $\mathcal{L}(\hat{u})$ in (5.2), we introduce the following discretization of derivatives in each component $k_j$ of the wavenumber:

$$ (5.3) \qquad \partial_{k_j} \hat{g}[t, k] = \text{DFT} [-2\pi i \tilde{x}_j \text{DFT}^{-1}[\hat{g}[t, k]]], $$

where the component $\tilde{x}_j = [0, 1, \dots, N/2 - 1, 0, -N/2 + 1, -N/2 + 2, \dots, -1] \Delta x$, considered as a vector and $\text{DFT}^{-1}[\hat{g}[k]]$ as a scalar, corresponds to the component $j$ of the position $x_j$ where the “Nyquist” mode has been set to 0 in order to respect the parity of the differentiation.

5.1.2. Discretization of the forcing. In order to produce the forcing term $\hat{f}[t, k]$, we generate an instance of $N^d$ independent, zero-average and unit-variance Gaussian random variables at each time step, which we store in a vector $g[t, x]$. Then we weigh them by the appropriate factor $(\Delta x)^{d/2}$, we take the DFT and multiply by the indicator function $\chi_{3\le k|L_{tot}\le 5}$, which ensures in particular that no energy is injected at the mode $k = 0$, and mostly at large scales. Mathematically, one gets

$$ (5.4) \qquad \hat{f}[t, k] = \chi_{3\le k|L_{tot}\le 5} \text{DFT} \left[ (\Delta x)^{d/2} g[t, x] \right]. $$

Finally, we choose $\kappa = 1/L_{tot}$ which is smaller than the smallest non vanishing wavenumber of the source $\hat{f}[t, k]$ (see Remark 1.2).

5.1.3. Algorithm.

* We pre-compute the forcing $\hat{f}[t, k]$ as in (5.4).

* Initialization step: $\hat{u}[0, k] = 0$.

* Induction step: given $\hat{u}[t, k]$, we compute $\hat{u}[t + \Delta t, k]$ via the following procedure:

(1) (Prediction step - spatial part) Compute $\mathcal{L}(\hat{u})[t, k]$ according to (5.3), and compute the numerical damping $\mathfrak{D}(\hat{u})[t, k]$.

(2) (Prediction step - temporal part) For each $k$ such that $|k| \ge \kappa$, compute the predictor $\hat{u}^*[t, k]$ according to

$$ \frac{\hat{u}^*[t, k] - \hat{u}[t, k]}{\Delta t} + \mathcal{L}(\hat{u})[t, k] + \mathfrak{D}(\hat{u})[t, k] = \hat{f}[t, k](\Delta t)^{-\frac{1}{2}}. $$

(3) (Correction step - spatial part) Compute $\mathcal{L}(\hat{u}^*)[t, k]$ with (5.3), and the numerical damping $\mathfrak{D}(\hat{u}^*)[t, k]$.
(4) (Correction step - temporal part) For each $k$ such that $|k| \ge \kappa$, compute the corrector $\hat{u}[t + \Delta t, k]$ according to

$$ \frac{\hat{u}[t + \Delta t, k] - \hat{u}[t, k]}{\Delta t} + \frac{\mathcal{L}(\hat{u}^*)[t, k] + \mathcal{D}(\hat{u}^*)[t, k] + \mathcal{L}(\hat{u})[t, k] + \mathcal{D}(\hat{u})[t, k]}{2} = \hat{f}[t, k](\Delta t)^{-\frac{1}{2}}. $$

## 5.2. Discussion of the numerical method.

5.2.1. Discretization of the operators $\mathcal{L}$. For any choice made for the estimation of the derivatives entering in $\mathcal{L}(\hat{u})$ (see (5.2)), the way the wavevector is discretized in a Cartesian fashion is not well adapted to this numerical problem which has natural spherical symmetry. Besides the forcing $\hat{f}$, which is only statistically isotropic, the deterministic part of the evolution is spherically symmetric$^7$. Nonetheless, the way it is presently discretized in a Cartesian form allows an easy and standard implementation of the DFT (using the Fast Fourier Transform algorithm) such that the solution $u[t, x]$ in physical space, over the discrete set of positions $x = (x_i)_{1 \le i \le d}$ with $x_i = [0, 1, ..., N/2, -N/2 + 1, -N/2 + 2, ..., -1]\Delta x$, is obtained with an inverse DFT of $\hat{u}[t, k]$.

As we have already mentioned, we cannot make sense of $\hat{u}(t, k)$ as a pointwise function of $k$. To this regard, the formulation of the problem in Fourier space needs to be justified and a proper meaning has to be given to the divergence of a rough field, see Sections 3 and 4 for full details. For the purpose of numerical simulations we have thus proposed (5.3) as a numerical realization of the derivatives, with the use of back and forth DFTs and an appropriate multiplication by a space-dependent factor.

5.2.2. Discretization of the operators $\partial_t$. Concerning time integration for a given time stepping $\Delta t$, we use an explicit predictor-corrector method with a single independent realization of the forcing, consisting in predicting the solution $\hat{u}[t + \Delta t, k]$ using an explicit Euler discretization scheme taking previous time step $\hat{u}[t, k]$ as an initial condition, and correcting by another Euler scheme which weighs the initial condition and the prediction equally.

5.2.3. Choice of discretization parameters. Let us give the last relevant parameters that we use for our simulations, keeping in mind that the spatial resolution $\Delta x = L_{tot}/N$ is determined by the physical length of the side of the periodic box (henceforth, we choose $L_{tot} = 1$ without loss of generality), and the number of collocation points $N$ in each direction. As we will see, the viscosity $\nu$ is chosen accordingly.

Firstly, the time step $\Delta t$ has to be chosen. At this stage, a CFL criterion for the dynamics proposed in our algorithm is not clear. Thus for the time being, we assume that a reasonable stability criterion may come from the time derivative and the viscous term. In that spirit, it seems reasonable to consider the numerical stability of the heat equation, which is imposed by the viscous term and which requires$^8$ that $\nu\Delta t < (\Delta x)^2/2$ (see for instance Ref. [37]).

Forthcoming simulations will be done for very high values of the number of collocation points $N$ in order to consider small values for the viscosity $\nu$. Even in the most comfortable situation where $\nu\Delta t$ is of order $\Delta x$, the numerical cost will eventually become prohibitive. For this reason, we choose a time step independent of resolution and dimension $d$, and given by $\Delta t = 5 \times 10^{-3}$. Additional simulations (data not shown) for $d = 1$ and $d = 2$ using $\Delta t = (\Delta x)^2/2$ and for moderate numbers of collocation points $N$ have given similar numerical results as the ones obtained with

<footnote>$^7$Indeed, all quantities involve the vector amplitude $|k|$. Moreover, $\mathcal{L}(\hat{u})$ can be written as a radial derivative with respect to $|k|$ (see (3.2)).</footnote>

<footnote>$^8$A more sophisticated scheme could be written including explicitly the exact solution of the underlying heat equation in the time stepping; this is known as an exponential scheme and would allow for $\nu\Delta t$ of order $\Delta x$.</footnote>
$\Delta t = 5 \times 10^{-3}$. It is unclear at this stage why this chosen value of the time step does not lead to numerical instabilities, although we could invoke the fact that viscosities will be chosen very small.

### 5.3. Simulations.

5.3.1. Determination of viscosity and averaging procedure. All predictions that have been made in the theoretical sections concern the statistical behavior of the solution $u(t, x)$ of the continuous problem recalled in (3.12), with or without viscosity. From a numerical point of view, we need to define a time $T_*$ at which the system has reached a statistically stationary regime, in which mathematical expectations will be approximated by an empirical average in time. Because of the cascading process of energy will eventually populate modes at higher and higher wavenumber $k$ as time goes on, we introduce viscosity in order to damp all the energy once it has reached the highest accessible wave-number

$$k_{\max} := (N/2)\Delta k.$$

In an inviscid ($\nu = 0$) regime, energy injected at low wavenumbers eventually reaches $k_{\max}$ at a time scale of the order of

$$T_* := k_{\max}/c.$$

We choose $\nu$ smaller than $ck_{\max}^{-3}$ to ensure a strong decrease of the spectral density at high wavenumbers (see Theorem 3.7). In practice, exploratory simulations will be carried out with various values of $\nu$, all of them satisfying this above criterion. This criterion is essential to ensure that no waves are reflected at the boundary of the artificial wavenumber periodic box. We will consider that $T_*$ is the time when the statistically stationary regime is reached.

We let then the simulation run after $T_*$, and we will average statistical estimators that we will define later over a set made of 100 instances every $(10^3\Delta t)$ to ensure statistical independence$^9$. More explicitly, forthcoming time averages $\langle g[t] \rangle_t$ of a given function of time $g[t]$ will be taken as

$$(5.5) \qquad \langle g[t] \rangle_t := \frac{1}{100} \sum_{n=1}^{100} g[T_* + n \times (10^3 \Delta t)].$$

Accordingly, we stop the integration in time at

$$T^* = T_* + 100 \times (10^3 \Delta t).$$

With these choices, simulations are longer and longer as the number of collocation points $N$ increases, but it allows for smaller and smaller viscosities, which is necessary to develop an extended inertial range.

5.3.2. Definition and estimation of key statistical quantities. First of all, the $L^2$-norm $\sigma_u^2[t]$ of the solution $u[t, x]$ at a given instant $t$ is defined by

$$(5.6) \qquad \sigma_u^2[t] = \sum_k |\hat{u}[t, k]|^2 \Delta k,$$

where the sum is taken over all possible values of the discrete wave vector $k$. Initially, we have $\sigma_u^2[0] = 0$, it will then quickly grow and end up fluctuating around a certain average value way before $T_*$ (data not shown). This aforementioned mean value $\langle \sigma_u^2[t] \rangle_t$, where the time average procedure is defined in (5.5), could be considered as the variance of the solution, and we have checked that it is independent of viscosity as expected if $\nu$ is chosen small enough. From the physical point of view, this is consistent with the observation that the variance of the solution of the Navier-Stokes equations get independent of viscosity at large Reynolds numbers (1.1).

<footnote>$^9$The number of instances and the time between each samples are empirically chosen.</footnote>
To characterize more precisely the statistical behavior of the solution when the $L^2$-norm starts fluctuating around a mean value (5.6), we define the energy spectral density estimated as a periodogram, i.e. the norm square of the Fourier mode, that is

$$ (5.7) \qquad \hat{C}_u(t, k) = |\hat{u}[t, k]|^2, $$

and its averaged version

$$ (5.8) \qquad \hat{C}_u(k) = \langle |\hat{u}[t, k]|^2 \rangle_t $$

where the time-average is defined in (5.5).

Another quantity of great importance is the so-called second-order structure function, i.e. the variance of the increment over a scale $\ell \in \mathbb{R}^d$, and given by

$$ (5.9) \qquad S_2(\ell) = \langle (u[t, x + \ell] - u[t, x])^2 \rangle_{t,x}, $$

where the time average is defined in (5.5) and the additional spatial-average is taken over $x \in \mathbb{R}^d$ for a given function of space $g[x]$:

$$ \langle g[x] \rangle_x := \frac{1}{L_{\text{tot}}^d} \sum_x g[x] (\Delta x)^d. $$

where the sum is taken over all possible values of the discrete wave vector $x$. We recall that concerning solutions of the Navier-Stokes equations, it is observed that $S_2(\ell)$ behaves as $|\ell|^{2/3}$ at infinite Reynolds numbers, corresponding to a regularity of Hölder type with $H = 1/3$ (c.f. (1.2)).

### 5.3.3. Results and comments in dimension $d=1$
We display in Figure 2 the results of our simulations for space dimension $d=1$. We have used $c=1$ and $H=1/3$, and the time step $\Delta t = 5 \times 10^{-3}$, whose value is motivated in section 5.2.3. As explained in Section 5.3.1, we run the simulation until $T_* = k_{\max}/c$ at which all accessible length scales and wave lengths have been populated. After this transient, we then average various estimators such as the energy spectral density (5.8) and the second order structure function (5.9) every 2 units of time. We indeed observe (data not shown) much before $T_*$ that the $L^2$-norm of the solution (5.6) fluctuates around a mean value. Typical snapshots of the solution $u[\cdot, x]$ are displayed in Figs. 2(a) and (c) at a time pertaining to the statistically stationary regime. A moderate viscosity $\nu = 10^{-8}$ has been used in (a), whereas a much smaller one $\nu = 10^{-12}$ is used in (c). We can see that as $\nu$ gets smaller, the velocity profile develops smaller scales and for the smallest viscosity that has been considered, it looks rougher. The averaged periodograms (5.8) obtained for all considered viscosities are represented in Figure 2(b). All estimated power spectral densities coincide at low wave lengths since the same forcing term has been used in all simulations. Then, at higher wave lengths, the spectra develop a universal power-law, independent of the characteristics of the forcing, with an exponent $-5/3$ which coincides with the expected $-(2H+d)$ exponent (see Eq. (3.10)) when choosing $H=1/3$. Then, at a characteristic wave length determined by the value of viscosity, spectra undergo a strong decrease, which is reminiscent of viscous damping. Similar behaviors are observed on the second order structure function $S_2(\ell)$ (5.9). At large length scales $\ell$ of the order of the scales where energy is injected into the system, determined by the spectral support of the forcing, $S_2(\ell)$ is independent of viscosity. At these scales, larger than the correlation length of the spatial profile of $u[t, x]$, $S_2(\ell)$ is approximatively equal to two times the variance of the solution. The fact that all curves superimposed at these scales show that this variance is indeed independent of viscosity. In the so-called inertial range of scales, as it is recalled in Section 1.2, $S_2(\ell)$ develops a power-law behavior, whose exponent $2H$ governs the Hölder regularity of the asymptotic solution. At smaller length scales, viscosity smooths out any irregular variations such that $S_2(\ell)$ can be Taylor expanded and becomes proportional to $\ell^2$.
<figure><img src="image_5.png" /><figcaption>FIGURE 2. Local and statistical behaviors of the solution $u[t, x]$ to the evolution provided in paragraph 5.1.3 for space dimension $d = 1$, for a given viscosity $\nu$ in the statistically stationary range. All simulations have been done with $c = 1$ and $H = 1/3$ (a): spatial profiles of $u[\cdot, x]$ at a given time $t$ in the statistically stationary regime for $\nu = 10^{-8}$ obtained using $N = 2^{10}$ collocation points. (b): estimations of the energy spectral density based on the averaged periodograms (see Section 5.3.2 and Equation (5.8)) of the solution for various values of viscosity $\nu = 10^{-8}, 10^{-9}, 10^{-10}, 10^{-11}$ and $10^{-12}$ (from left to right), using respectively $N = 2^{10}, 2^{11}, 2^{12}, 2^{13}$ and $2^{14}$ collocation points. We superimpose with a dashed line the asymptotic prediction $|k|^{-5/3}$. (c): same plot as in (a), but for a lower value of viscosity $\nu = 10^{-12}$. (d): Similar plot as for (b) but for the second order structure function $S_2(\ell)$ (5.9), i.e. the variance of the increments, following an averaging procedure detailed in the text. We superimpose the expected asymptotic power-law behavior $\ell^{2/3}$.</figcaption></figure>

5.3.4. *Results and comments in dimension d = 2.* We now display in Figure 3 the results of our simulations for space dimension $d = 2$. We have used again $c = 1$ and $H = 1/3$, and the time step $\Delta t = 5 \times 10^{-3}$. Similarly to the $d = 1$ case, we propagate in time the discrete evolution provided in paragraph 5.1.3 using the predictor-corrector method until time $T_*$ at which all accessible scales have been populated. Before then, viscosity has damped all the energy at small scales such that the Fourier mode at $k_{max}$ is exponentially small. Hence, we will consider that starting from time $T_*$, the statistically stationary regime has also been reached. We display in Figure 3(a) the amplitude of the Fourier modes $\hat{u}[t, k]$ at a time $t$ lying in the statistically stationary regime, as a function of the two components $k_x$ and $k_y$ of the wave vector $k$, in a logarithmic representation. We can see that a lot of energy is concentrated along the two lines ($k_x$, 0) and (0, $k_y$), whereas elsewhere in the plane, up to fluctuations, energy is distributed in a rotation invariant (i.e. isotropic) way. Recall
<figure><img src="image_6.png" /><figcaption>FIGURE 3. Local and statistical behaviors of the solution $u[t, x]$ to the evolution provided in paragraph 5.1.3 for space dimension $d = 2$, for a given viscosity $\nu$ in the statistically stationary range. All simulations have been done with $c = 1$ and $H = 1/3$ (a): Representation of the logarithm of the absolute value of $\hat{u}[\cdot, k]$ at a given time in the statistically stationary regime for $\nu = 10^{-9}$ obtained using $N = 2^{11}$ collocation points. (b): estimations of the angle-averaged energy spectral density based on the averaged periodograms (see Section 5.3.2 and (5.8)) of the solution for various values of viscosity $\nu = 10^{-5}, 10^{-6}, 10^{-7}, 10^{-8}$ and $10^{-9}$ (from left to right), using respectively $N = 2^7, 2^8, 2^9, 2^{10}$ and $2^{11}$ collocation points. We superimpose with a dashed line the asymptotic prediction $|k|^{-2H-d}$, with $H = 1/3$ and $d = 2$. (c): similar plot as in (a) but for the corresponding spatial profiles of $u[\cdot, x]$. (d) Similar plot as for (b) but for the second order structure function $S_2(\ell)$ (5.9), i.e. the variance of the increments, following an averaging procedure detailed in the text. We superimpose the expected asymptotic power-law behavior $\ell^{2H}$, with $H = 1/3$.</figcaption></figure>

that in a continuous framework, the statistical properties of Fourier modes $\hat{u}(t, k)$ at any time $t$ are expected to depend only on the amplitude $|k|$. This shows that our numerical representation $\hat{u}[t, k]$ is intrinsically anisotropic. We interpret this spurious anisotropy as the consequence of the finiteness of our simulation domain, with the implied finiteness of the resolution $\Delta k = 1/L_{tot}$, but also the fact that the Fourier modes $\hat{u}[t, k]$ are distributed on a Cartesian grid, whereas the continuous solution $\hat{u}(t, k)$ is expected in average to be spherically symmetric. Another limitation of our numerical approach is related to the rough nature of the expected solution. It is clear from the inspection of Figure 3(a) that Fourier modes are correlated and smoother than what is expected from the Fourier transform of a statistically isotropic random field. This is very possibly related to the estimation of the divergence operator entering in the evolution (3.12) with finite differences, as
it is implicitly done using back and forth DFTs (see Equation (5.3)). A specifically devoted article aimed at exploring the numerical representation of the continuous formulation provided in (3.12) would be needed, designing for instance some finite volume algorithms able to deal with the rough nature of the underlying fields. It will be the subject of future publications.

Accordingly, the representation of the solution in physical space, that we display in Figure 3(c), exhibits two types of anisotropies. The first type of anisotropy can be observed along the two directions $x$ and $y$ of the Cartesian frame as straight lines. This anisotropy is consistent with what is observed on the Fourier transform displayed in Figure 3(a) along $k_x$ and $k_y$. From the inspection of Figure 3(c), another type of anisotropy can be evidenced around the origin, i.e. around $x = y = 0$. Recall that for a statistically homogenous field, probability laws are expected to be invariant by translation, and thus the origin of the domain should not play a particular role. In our numerical solution, this is clearly not the case, and we believe that this anisotropy is related to the aforementioned correlated nature of the Fourier modes displayed in Figure 3(a). Once again, more work is needed to design proper numerical schemes able to get rid of these spurious anisotropies.

Nonetheless, despite these anisotropies, that are not present in the solution of the continuous framework, our numerical solution behaves in a statistical sense in the expected way. For instance, we display in Figure 3(b) the estimation of the power spectral densities, obtained while averaging in time the square of the amplitude of Fourier modes. We underline that these spectral densities are furthermore averaged over the angles that the wave vector $k$ makes with the axes of the Cartesian frame, such that displayed spectral densities depend only on the wave length amplitude $|k|$. We indeed observe that, as viscosity gets smaller and smaller, spectral densities develop a power-law behavior in the inertial range of scales, with the expected exponent $-(2H + d)$ which is derived in (3.10). At higher wavenumbers, Fourier modes are exponentially damped by viscosity. Similar conclusions could be drawn from inspection of the scale dependence of the second order structure function $S_2(\ell)$ that is shown in Figure 3(d). Let us mention here that spatial averages have been obtained over the full spatial domain, averaging furthermore over two structure functions obtained while considering spatial lags $\ell = |\ell|e_x$ and $\ell = |\ell|e_y$, where $e_x$ and $e_y$ are the two orthonormal unit vectors of the Cartesian frame. At the largest scales, above the characteristic ones of the forcing, we can see that $S_2(\ell)$ reaches a plateau, which gets independent of viscosity as $\nu \to 0$. Similarly to the $d=1$ case, above the correlation length of $u[t, x]$, $S_2(\ell)$ coincides with two times the variance of $u[t, x]$, which says in other words that the variance gets itself independent of viscosity if $\nu$ is chosen small enough. At lower scales, i.e. in the inertial range, $S_2(\ell)$ develops a power-law behavior of exponent $2H$, as it is expected from the predictions made in the continuous framework (Corollary 4.12), over a range of scales which grows as viscosity gets smaller and smaller. Finally, at even smaller scales, viscous effects dominate and smooth out the spatial profiles, such that $S_2(\ell)$ gets proportional to $\ell^2$.

5.3.5. *Results and comments in dimension d = 3.* Let us now finish this Section devoted to numerical simulations by presenting in Figure 4 the results for space dimension $d=3$. Let us mention that in a three-dimensional setup, simulations are much more demanding from a computational perspective, and the cost of performing back and forth FFTs gets much higher, because derivatives along the 3 directions have to be considered, and also because the number of discretization points, $N^3$, increase tremendously as $N$ increases. For these reasons, and because we are propagating the integration in time from a vanishing initial condition towards the statistically stationary regime before taking averages, we could not go above $N = 2^9 = 512$ collocation points along each direction. Consequently, we have not been able to perform simulations for viscosities smaller than $\nu = 10^{-7}$. Nonetheless, we observe (data not shown) that fluctuations as quantified by the $L^2$-norm (5.6) get independent of viscosity in a good approximation starting from $\nu = 10^{-6}$, as it is expected from
<figure><img src="image_7.png" /><figcaption>FIGURE 4. Local and statistical behaviors of the solution $u[t, x]$ to the evolution provided in paragraph 5.1.3 for space dimension $d = 3$, for a given viscosity $\nu$ in the statistically stationary range. All simulations have been done with $c = 1$ and $H = 1/3$ (a): Representation of some slices of $u[\cdot, x]$ at a given time in the statistically stationary regime for $\nu = 10^{-7}$ obtained using $N = 2^9$ collocation points in each directions. (b): estimations of the angle-averaged energy spectral density based on the averaged periodograms (see Section 5.3.2 and (5.8)) of the solution for various values of viscosity $\nu = 10^{-3}, 10^{-4}, 10^{-5}, 10^{-6}$ and $10^{-7}$ (from left to right), using respectively $N = 2^5, 2^6, 2^7, 2^8$ and $2^9$ collocation points. We superimpose with a dashed line the asymptotic prediction $|k|^{-2H-d}$, with $H = 1/3$ and $d = 3$.</figcaption></figure>

the behavior of the solution in a continuous framework. As mentioned, and similarly to the $d = 1$ and $d = 2$ cases, we go through the transient while integrating the solution until time $T_*$, and only then we start taking averages.

We display in Figure 4(a) a rendering of our three-dimensional simulations in physical space in the Cartesian frame. For the sake of clarity, and because visualizations gets more complicated, we only show three slices along the planes $(x, y, 0)$, $(x, 0, z)$ and $(0, y, z)$. Bright and dark colors correspond respectively to large positive and large negative fluctuations of the solution, similarly to what has been observed for the $d = 2$ case, which is displayed in Figure 3(c). Once again, we observe the two types of anisotropies that we evidenced in the $d = 2$ case, one along the three directions of the Cartesian frame, and one around the origin.

Nonetheless, as it is displayed in Figure 4(b), the simulations behave as expected in a statistical manner, as it can be observed on the power spectral densities (5.8). Once again, these densities are not only averaged in time once the statistically stationary regime has been reached, but they are also averaged over the two angles that the wave vector $k$ does with the Cartesian axes, such that $\hat{C}_u(k)$ is a function of the norm $|k|$ only. We can see that at large scales, i.e. at low wavenumbers $|k|$, fluctuations are independent of viscosity, even for the highest value $\nu = 10^{-3}$ which is used for the lowest number of collocation points $N = 2^5$. It is nonetheless crucial to consider smaller
values of viscosity, necessitating thus higher values of $N$, up to $N = 2^9$, in order to develop an extended inertial range. It is clear for the smallest value of viscosity ($\nu = 10^{-7}$ for $N = 2^9$) that the spectral density has developed a power-law behavior in the inertial range, with the expected exponent $-(2H + d)$, in a consistent manner with our theoretical predictions (3.10).

We thus see that we are able to give an appropriate numerical representation of the continuous framework using the DFT to define the divergence operator entering in the spectral evolution provided in (3.12). In particular, we are able to reproduce in a discrete setup the statistical behaviors of the spectral densities (Eq. 3.10) and second-order structure functions (Corollary 4.12), and their related power-law behaviors. Nonetheless, more work is needed to get rid of the anisotropies that are clearly observed for the $d=2$ and $d=3$ cases. To do so, a promising direction could be given while designing a finite volume scheme able to respect the inherent spherical symmetry of the deterministic part of the evolution (3.12). This is required in order to propose a realistic model of fully developed fluid turbulence that could be used in an efficient way in various applications, in which spatial fluctuations of the velocity field have crucial consequences on the evolutions of dynamical quantities of interest. We keep these important developments for future investigations.

## REFERENCES

[1] J. ALMONACID AND N. NIGAM, Characterization of singular flows of zeroth-order pseudo-differential operators via elliptic eigenfunctions: a numerical study, 2022. arXiv preprint.

[2] G. B. APOLINÁRIO AND L. CHEVILLARD, Space-time statistics of a linear dynamical energy cascade model, Mathematics in Engineering, 5 (2023), pp. 1–23.

[3] G. B. APOLINÁRIO, L. CHEVILLARD, AND J.-C. MOURRAT, Dynamical fractional and multifractal fields, J Stat Phys, 186 (2022).

[4] D. BARBATO, L. A. BIANCHI, F. FLANDOLI, AND F. MORANDIN, A dyadic model on a tree, Journal of Mathematical Physics, 54 (2013), p. 021507.

[5] D. BARBATO, F. MORANDIN, AND M. ROMITO, Smooth solutions for the dyadic model, Nonlinearity, 24 (2011), p. 3083.

[6] G. K. BATCHELOR, The theory of homogeneous turbulence, Cambridge University Press, 1953.

[7] G. BECK, Modélisation et étude mathématique de réseaux de câbles électriques, thesis, Université Paris-Saclay, Mar. 2016.

[8] G. BECK AND D. LANNES, Freely floating objects on a fluid governed by the boussinesq equations, Annales de l'Institut Henri Poincaré C, Analyse non linéaire, 39 (2022).

[9] L. BIFERALE, Shell models of energy cascade in turbulence, Annual review of fluid mechanics, 35 (2003), pp. 441–468.

[10] T. BOHR, M. H. JENSEN, G. PALADIN, AND A. VULPIANI, Dynamical systems approach to turbulence, Cambridge University Press, 1998.

[11] R. CARLES AND C. CHEVERRY, Constructive and destructive interferences in nonlinear hyperbolic equations, Mémoires de la Société Mathématique de France, 174 (2022).

[12] M. CHAVES, K. GAWEDZKI, P. HORVAI, A. KUPIAINEN, AND M. VERGASSOLA, Lagrangian dispersion in gaussian self-similar velocity ensembles, Journal of statistical physics, 113 (2003), pp. 643–692.

[13] A. CHESKIDOV, M. DAI, AND S. FRIEDLANDER, Dyadic models for fluid equations: a survey, arXiv preprint arXiv:2209.10203, (2022).

[14] L. CHEVILLARD, Une peinture aléatoire de la turbulence des fluides, habilitation à diriger des recherches, ENS Lyon, Oct. 2015.

[15] L. CHEVILLARD, C. GARBAN, R. RHODES, AND V. VARGAS, On a skewed and multifractal unidimensional random field, as a probabilistic representation of kolmogorov's views on turbulence, in Annales Henri Poincaré, vol. 20, Springer, 2019, pp. 3693–3741.

[16] S. COHEN AND J. ISTAS, Fractional fields and applications, vol. 2, Springer, 2013.

[17] Y. COLIN DE VERDIÈRE AND L. SAINT-RAYMOND, Attractors for two-dimensional waves with homogeneous hamiltonians of degree 0, Communications on Pure and Applied Mathematics, 73 (2020), pp. 421–462.

[18] I. DAUBECHIES, Ten lectures on wavelets, SIAM, 1992.

[19] S. DYATLOV AND M. ZWORSKI, Microlocal analysis of forced waves, Pure and Applied Analysis, 1 (2019), pp. 359–384.

[20] U. FRISCH, Turbulence: the legacy of AN Kolmogorov, Cambridge University Press, 1995.
[21] J. GALKOWSKI AND M. ZWORSKI, *Viscosity limits for zeroth-order pseudodifferential operators*, Communications on Pure and Applied Mathematics, 75 (2022), pp. 1798–1869.

[22] S. JANSON, *Gaussian Hilbert Spaces*, Cambridge Tracts in Mathematics, Cambridge University Press, 1997.

[23] I. KARATZAS AND S. SHREVE, *Brownian Motion and Stochastic Calculus*, vol. 113 of Graduate Texts in Mathematics, Springer-Verlag New York, 1998.

[24] P. KLOEDEN AND E. PLATEN, *Numerical Solution of Stochastic Differential Equations*, Stochastic Modelling and Applied Probability, Springer Berlin Heidelberg, 2011.

[25] A. N. KOLMOGOROV, *Dissipation of Energy in Locally Isotropic Turbulence*, Akademiia Nauk SSSR Doklady, 32 (1941), p. 16.

[26] R. H. KRAICHNAN, *Small-scale structure of a scalar field convected by turbulence*, The Physics of Fluids, 11 (1968), pp. 945–953.

[27] D. LANNES, *Lecture notes volume of the winter school on fluid dynamics, dispersive equations and quantum fluids (Bressanone, December 2018)*. To appear, private communication.

[28] R. LEWANDOWSKI AND B. PINIER, *The Kolmogorov Law of Turbulence What Can Rigorously Be Proved? Part II*, The Foundations of Chaos Revisited: From Poincaré to Recent Advancements. Understanding Complex Systems. Springer, Cham, (2016).

[29] A. LODHIA, S. SHEFFIELD, X. SUN, AND S. S. WATSON, *Fractional Gaussian fields: A survey*, Probability Surveys, 13 (2016), pp. 1–56.

[30] L. R. MAAS, D. BENIELLI, J. SOMMERIA, AND F.-P. A. LAM, *Observation of an internal wave attractor in a confined, stably stratified fluid*, Nature, 388 (1997), pp. 557–561.

[31] A. A. MAILYBAEV, *Continuous representation for shell models of turbulence*, Nonlinearity, 28 (2015), p. 2497.

[32] A. MASPERO, *Generic transporters for the linear time dependent quantum harmonic oscillator on $\mathbb{R}$*, 2022. arXiv preprint.

[33] J. C. MATTINGLY, T. SUIDAN, AND E. VANDEN-EIJNDEN, *Simple systems with anomalous dissipation and energy cascade*, Communications in Mathematical Physics, 276 (2007), pp. 189–220.

[34] G. MÉTIVIER, *On the $L^2$ well posedness of Hyperbolic Initial Boundary Value Problems*, Annales de l'Institut Fourier, 67 (2017), pp. 1809–1863.

[35] A. MONIN AND A. YAGLOM, *Statistical fluid dynamics*, Vol. I and II MIT Press, Cambridge, (1971).

[36] S. B. POPE, *Turbulent flows*, Cambridge University Press, 2000.

[37] W. H. PRESS, S. A. TEUKOLSKY, W. T. VETTERLING, AND B. P. FLANNERY, *Numerical recipes 3rd edition: The art of scientific computing*, Cambridge University Press, 2007.

[38] R. RHODES AND V. VARGAS, *Gaussian multiplicative chaos and applications: a review*, Probab. Surv., 11 (2014), pp. 315–392, 1305.6221.

[39] M. RIEUTORD, B. GEORGEOT, AND L. VALDETTARO, *Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum*, Journal of Fluid Mechanics, 435 (2001), pp. 103–144.

[40] M. RIEUTORD AND L. VALDETTARO, *Inertial waves in a rotating spherical shell*, Journal of Fluid Mechanics, 341 (1997), pp. 77–99.

[41] H. SCOLAN, E. ERMANYUK, AND T. DAUXOIS, *Nonlinear fate of internal wave attractors*, Physical review letters, 110 (2013), p. 234501.

[42] H. TENNEKES AND J. L. LUMLEY, *A first course in turbulence*, MIT press, 1972.

[43] V. E. ZAKHAROV, V. S. L'VOV, AND G. FALKOVICH, *Kolmogorov spectra of turbulence I: Wave turbulence*, Springer Science & Business Media, 2012.
THEORETICAL PHYSICS I, UNIVERSITY OF BAYREUTH, UNIVERSITÄTSSTR. 30, 95447 BAYREUTH, GERMANY

Email address: gapolinario@uni-bayreuth.de

UNIV RENNES, IRMAR UMR 6625 & CENTRE INRIA DE L'UNIVERSITÉ DE RENNES (MINGUS) & ENS RENNES, FRANCE

Email address: geoffrey.a.beck@inria.fr

UNIV LYON, ENS DE LYON, UNIV. CLAUDE BERNARD, CNRS, LABORATOIRE DE PHYSIQUE, 46 ALLÉE D'ITALIE, 69342 LYON, FRANCE & UNIV LYON, UNIVERSITÉ CLAUDE BERNARD LYON 1, CNRS UMR 5208, INSTITUT CAMILLE JORDAN, 43 BOULEVARD DU 11 NOVEMBRE 1918, F-69622 VILLEURBANNE, FRANCE

Email address: laurent.chevillard@cnrs.fr

DÉPARTEMENT DE MATHÉMATIQUES ET APPLICATIONS, ÉCOLE NORMALE SUPÉRIEURE, CNRS, PSL UNIVERSITY, AND UNIVERSITÉ PARIS CITÉ, 75005 PARIS, FRANCE

Email address: isabelle.gallagher@ens.fr

DÉPARTEMENT DE MATHÉMATIQUES ET APPLICATIONS, ÉCOLE NORMALE SUPÉRIEURE, CNRS, PSL UNIVERSITY, 75005 PARIS, FRANCE

Email address: ricardo.grande@ens.fr