Regularized fractional Ornstein-Uhlenbeck processes, and their relevance to the modeling of fluid turbulence

TL;DR

This paper introduces a regularized fractional Ornstein-Uhlenbeck process to model turbulent fluid velocity, capturing roughness and increment behavior similar to fractional Brownian motion, with implications for simulation and analysis.

Contribution

It proposes a novel linear stochastic differential equation with fractional Gaussian noise, analyzing its properties and relevance to turbulence modeling, especially for rough cases where H<1/2.

Findings

The process reaches a stationary regime with bounded marginals as regularization vanishes.

Increment variance at small scales matches that of fractional Brownian motion for any H in (0,1).

The approach is suitable for simulating very rough stochastic processes.

Abstract

Motivated by the modeling of the temporal structure of the velocity field in a highly turbulent flow, we propose and study a linear stochastic differential equation that involves the ingredients of a Ornstein-Uhlenbeck process, supplemented by a fractional Gaussian noise, of parameter $H$, regularized over a (small) time scale $\epsilon>0$. A peculiar correlation between these two plays a key role in the establishment of the statistical properties of its solution. We show that this solution reaches a stationary regime, which marginals, including variance and increment variance, remain bounded when $\epsilon \to 0$. In particular, in this limit, for any $H\in ]0,1[$, we show that the increment variance behaves at small scales as the one of a fractional Brownian motion of same parameter $H$. From the theoretical side, this approach appears especially well suited to deal with the (very) rough case $H<1/2$, including the boundary value $H=0$, and to design simple and efficient numerical simulations.

Equations

$$dv(t)=-\frac{1}{T}v(t)dt+\sigma \sqrt{\frac{2}{T}}dW(t),$$

$$du(t)=-\frac{1}{T}u(t)dt+‘‘d{W}_H(t)",$$

$$d{X}_{ϵ,H}(t)=-\alpha X_{ϵ,H}(t)dt+d{W}_{ϵ,H}(t),$$

$$d{W}_{ϵ,H}(t)={\omega }_{ϵ,H}(t)dt+{ϵ}^{H-\frac{1}{2}}dW(t),$$

$${\omega }_{ϵ,H}(t)=\left(H-\frac{1}{2}\right){\int }_{-\infty }^t\frac{1}{(t-s+ϵ{)}^{\frac{3}{2}-H}}dW(s),$$

$$X_{ϵ,H}(t)=X_{ϵ,H}^{(1)}(t)+X_{ϵ,H}^{(2)}(t),$$

$$X_{ϵ,H}^{(1)}(t)={\int }_{t_0}^t{e}^{-\alpha (t-s)}{\omega }_{ϵ,H}(s)ds$$

$$X_{ϵ,H}^{(2)}(t)={ϵ}^{H-\frac{1}{2}}{\int }_{t_0}^t{e}^{-\alpha (t-s)}dW(s)$$

$$\mathbb{E}{X}_H^2$$

$$=\frac{{\alpha }^{-2H}{\left[\Gamma \left(H+\frac{1}{2}\right)\right]}^2}{2\sin (\pi H)}<\infty ,$$

$$\mathbb{E}{\left({\delta }_{\tau }{X}_H\right)}^2=\underset{ϵ\to 0}{\lim }\underset{t\to \infty }{\lim }\mathbb{E}\left[{\left(X_{ϵ,H}(t+\tau )-X_{ϵ,H}(t)\right)}^2\right].$$

$$\mathbb{E}{\left({\delta }_{\tau }{X}_H\right)}^2\underset{\tau \to 0}{\sim }\frac{1}{\sin (\pi H)}\frac{{\left[\Gamma \left(H+\frac{1}{2}\right)\right]}^2}{\Gamma (2H+1)}|\tau {|}^{2H}.$$

$$W_{ϵ,H}(t)$$

$${\int }_{-\infty }^0\left[(t-s+ϵ{)}^{H-1/2}-(-s+ϵ{)}^{H-1/2}\right]dW(s)$$

$$+{\int }_0^t(t-s+ϵ{)}^{H-1/2}dW(s).$$

$$d{W}_{ϵ,H}(t)\equiv W_{ϵ,H}(t+dt)-W_{ϵ,H}(t)=$$

$${\int }_{-\infty }^t\left[(t+dt-s+ϵ{)}^{H-1/2}-(t-s+ϵ{)}^{H-1/2}\right]dW(s)$$

$$+{\int }_t^{t+dt}(t+dt-s+ϵ{)}^{H-1/2}dW(s).$$

$$\underset{t\to \infty }{\lim }\mathbb{E}\left[X_{ϵ,0}^2(t)\right]\underset{ϵ\to 0}{\sim }\ln \frac{1}{ϵ},$$

$$\underset{ϵ\to 0}{\lim }\underset{t\to \infty }{\lim }\mathbb{E}\left[X_{ϵ,0}(t)X_{ϵ,0}(t+\tau )\right]\underset{\tau \to 0}{\sim }\ln \frac{1}{\tau }.$$

$$d{W}_{ϵ,H}(t)={\omega }_{ϵ,H}(t)dt+{ϵ}^{H-1/2}dW(t),$$

$${\mathcal{C}}_{ϵ,H}^{f,g}(t_2-t_1)$$

$$=\mathbb{E}\left[{\int }_{\mathbb{R}}f(t_1-s_1)d{W}_{ϵ,H}(s_1){\int }_{\mathbb{R}}g(t_2-s_2)d{W}_{ϵ,H}(s_2)\right]$$

$$\displaystyle=\int_{\mathbb{R}^{2}}f(t_{1}-s_{1})g(t_{2}-s_{2})\mathbb{E}\Big{[}dW_{\epsilon,H}(s_{1})dW_{\epsilon,H}(s_{2})\Big{]}.$$

$$\mathcal{C}_{H}^{f,g}(\tau)=\lim_{\epsilon\to 0}\mathcal{C}_{\epsilon,H}^{f,g}(\tau)=-\frac{1}{2}\frac{1}{\sin(\pi H)}\frac{\left[\Gamma\left(H+\frac{1}{2}\right)\right]^{2}}{\Gamma(2H)}\int_{0}^{+\infty}\Big{[}(f\star g)^{\prime}(\tau+h)-(f\star g)^{\prime}(\tau-h)\Big{]}h^{2H-1}dh,$$

$$(f\star g)(h)={\int }_{\mathbb{R}}f(s)g(h+s)ds.$$

$${\mathcal{C}}_{ϵ,1/2}^{f,g}(\tau )=(f\star g)(\tau ).$$

$$\displaystyle\lim_{\epsilon\to 0}\mathbb{E}\Big{[}dW_{\epsilon,H}(s_{1})dW_{\epsilon,H}(s_{2})\Big{]}$$

$$\displaystyle=\lim_{\epsilon\to 0}\mathbb{E}\Big{[}\omega_{\epsilon,H}(s_{1})\omega_{\epsilon,H}(s_{2})\Big{]}ds_{1}ds_{2}$$

$$={\left(H-\frac{1}{2}\right)}^2{\int }_{u=0}^{\infty }\frac{1}{u^{3/2-H}}\frac{1}{(u+|s_1-s_2|{)}^{3/2-H}}dud{s}_{1}d{s}_2$$

$$=\left(H-\frac{1}{2}\right)\frac{1}{\sin (\pi H)}\frac{{\left[\Gamma \left(H+\frac{1}{2}\right)\right]}^2}{\Gamma (2H)}|s_1-s_2{|}^{2H-2}d{s}_{1}d{s}_2.$$

$${\mathcal{C}}_{ϵ,H>\frac{1}{2}}^{f,g}(\tau )=$$

$$\times {\int }_{\mathbb{R}}(f\star g)(\tau +h)|h{|}^{2H-2}dh,$$

$$\displaystyle\mathbb{E}\Big{[}dW_{\epsilon,H}(s_{1})dW_{\epsilon,H}(s_{2})\Big{]}=$$

$$\displaystyle\left(\mathbb{E}\Big{[}\omega_{\epsilon,H}(s_{1})\omega_{\epsilon,H}(s_{2})\Big{]}+\epsilon^{2H-1}\delta(s_{1}-s_{2})\right)ds_{1}ds_{2}+$$

$$\displaystyle\epsilon^{H-\frac{1}{2}}\left(\mathbb{E}\Big{[}\omega_{\epsilon,H}(s_{1})dW(s_{2})\Big{]}ds_{1}+\mathbb{E}\Big{[}\omega_{\epsilon,H}(s_{2})dW(s_{1})\Big{]}ds_{2}\right),$$

$$\displaystyle\mathbb{E}\Big{[}\omega_{\epsilon,H}(s_{1})\omega_{\epsilon,H}(s_{2})\Big{]}=\left(H-\frac{1}{2}\right)^{2}$$

Full text

I Introduction

The motivation of this present work comes from the stochastic modeling of certain aspects of fluid turbulence Fri95 ; Pop00 . In three dimensional space, a fluid, when stirred at large scale $L$ by a force, reacts while developing fluctuations: the fluid velocity $\textbf{u}(\textbf{x},t)$ reaches a statistically stationary regime characterized by a standard deviation, say $\sigma$. The more intense is the forcing, the bigger is $\sigma$. If it happens that the product $\sigma L$ is much bigger than the kinematic viscosity $\nu$ of the fluid, such that the so-called Reynolds number $\mathcal{R}_{e}=\sigma L/\nu$ is $\gg 1$, the fluid motions become turbulent. In this regime, velocity fluctuations, characterized by $\sigma$, are observed to be independent on the very nature of the dissipative mechanisms that are taking place, in particular of the precise value of $\nu$. In order to achieve such a efficient way to dissipate energy, the fluid will develop small spatial scales such that asymptotically, as $\mathcal{R}_{e}\to\infty$, the velocity field will develop infinite gradients, and becomes rough Kol41 ; Fri95 .

Rephrased in terms of regularity of functions, this phenomenology, mainly due to Kolmogorov Kol41 ; Fri95 , and that we depicted schematically in the former paragraph, has a precise formulation if we assume underlying Gaussian statistics: the velocity field $\textbf{u}(\textbf{x},t)$ is a finite-variance (zero-average under the assumptions of homogeneity and isotropy) Gaussian random field, which spatial increments over $\ell$ along any directions behave as those of a fractional Brownian motion ManVan68 of parameter $H=1/3$ (see Refs RobVar08 ; CheRob10 ; PerGar16 for recent developments on this matter).

Until now, we focused on the statistical characterization of the scale-invariance properties of the velocity field in space. We would like now to focus on the stochastic modeling of the temporal structure of turbulence. A similar phenomenology can be developed for the Lagrangian velocity $\textbf{v}(t)=\textbf{u}(\textbf{r}(t),t)$ of a fluid particle, initially at the position say $\textbf{r}_{0}$, along its trajectory $\textbf{r}(t)$ defined through the Eulerian velocity field u as $\partial\textbf{r}(t)/\partial t=\textbf{u}(\textbf{r}(t),t)$. In this case, experimental and numerical observations (see the reviews Yeu02 ; TosBod09 ; PinSaw12 and references therein) suggest strongly, as expected from a dimensional analysis, that Lagrangian velocity $v$ is also a finite-variance process with the same standard deviation $\sigma$, typically correlated over the large time scale $T=L/\sigma$, and which increments behave, up to the variance, at small scale as those of the Brownian motion. In a Gaussian approximation, which is a good starting point in the context of stochastic modeling, but which is too simplistic to reproduce all the observed statistics (in particular the intermittency phenomenon Yeu02 ; TosBod09 ; PinSaw12 ; CheCas12 ), a natural stochastic model for the dynamics of $v$ would be a Ornstein-Uhlenbeck process that reads

[TABLE]

where $W$ is a Wiener process. The statistical properties of the solution $v$ to this stochastic differential equation (Eq. 1) are well known, and reproduce adequately at this level of phenomenology the fluctuations of the turbulent Lagrangian velocity. In this view, Lagrangian velocity is Hölder continuous of parameter $H=1/2$.

If now we are interested in the temporal description of the fluctuations of the full Eulerian velocity field $\textbf{u}(x,y,z,t)$, or let us say at the description at a fixed position of one component of the velocity field, for example $u(t)\equiv u_{x}(x_{0},y_{0},z_{0},t)$, then dimensional arguments TenLum72 ask for a finite variance process which is Hölder continuous as those of a fractional Brownian motion of parameter $H=1/3$, thus less regular than the Lagrangian counterpart $v$. This extrapolated non-differential behavior at infinite Reynolds number can be observed in numerical simulations, as it has been done in Ref. CheRou05 , and recently revisited in the context of renormalization group theory CanRos17 . In this case, a natural stochastic model for the dynamics of $u$ would be given by a generalization of the Ornstein-Uhlenbeck process (Eq. 1) of the form

[TABLE]

that would lead to a stationary regime, in which $u$ is of finite variance, and behaves at small scales as a fractional Brownian motion of parameter $H\in]0,1[$. All the purpose of this article is devoted to give a precise mathematical (and numerical) meaning to the random measure $``dW_{H}"$ that enters the dynamics of $u$ (Eq. 2). Let us mention that this level of roughness ($H<1/2$) is also a hallmark of subdiffusive dynamics, as it is studied in MagWer09 .

A natural way to generalize the dynamics of the Ornstein-Uhlenbeck process (Eq. 1) to the fractional case, i.e. a fractional Ornstein-Uhlenbeck process (Eq. 2), that would lead to a stationary finite variance process with a appropriate rough behavior at small scales, is to consider a fractional Brownian motion (fBm) $W_{H}(t)$ of parameter $H$. The fBm is a well defined probabilistic object ManVan68 , and allows to define accordingly the integrated dynamics as $u(t)-u(0)=-\frac{1}{T}\int_{0}^{t}u(s)ds+W_{H}(t)-W_{H}(0)$, as it is studied in CheKaw03 ; KaaSal11 ; HuNua10 . This procedure is clearly well defined, can be extended to a more general framework allowing to elaborate a stochastic calculus with respect to fractional Brownian motion PipTaq00 ; BiaHu08 ; Mis08 . Such a dynamics for $u$ indeed leads to the statistical properties that have been listed. But it eludes the question regarding the meaning of the infinitesimal increment $``dW_{H}(t)"$.

Nonetheless, as already proposed in ManVan68 , it is tempting to use the so-called fractional Gaussian noise, that we will specify later, to give a meaning to this infinitesimal increment. Indeed, for $H>1/2$, this fractional Gaussian noise has a well-behaved covariance, which is bounded for non vanishing argument. The purpose of this article is to include in this picture the (very) rough case $H<1/2$, as demanded by the physics of turbulence ($H=1/3$), that requires a different method of construction due to the pathological nature of the fractional Gaussian noise (its covariance is no more bounded) at these levels of roughness. This will be achieved using a regularized form of this noise over a small time scale, say $\epsilon$, supplemented by a Gaussian white noise weighted by a factor that diverges with $\epsilon$ (see the following Eq. 4). As we will see, an additional correlation between these two plays a key role. Once injected in the dynamics (Eq. 2), we will then study the limiting behavior of the marginals of the Gaussian process $u$ when $\epsilon\to 0$. In this sense, the regularization procedure brings to light the underlying mechanisms at play, and allows to give a precise mathematical meaning (up to $\epsilon>0$) to the random measure $dW_{H}(t)$ that enters in the dynamics of $u$ (Eq. 2). This overall picture clarifies that the expected dynamics leading to stationary and rough processes has to be non Markovian. In particular, we recover standard interpretations pertaining to fractional Brownian motions, namely that the infinitesimal increment is positively correlated (or persistent) for $H>1/2$, and negatively correlated (or anti-persistent) for $H<1/2$. Incidentally, it also gives a way to build up a realistic numerical approximation of the trajectories of $u$ at a given regularization scale $\epsilon$. Finally, this approach allows to deal with the boundary value $H=0$ in a non ambiguous manner, a case that is tricky to consider using the standard approach consisting in working directly with a fractional Brownian motion, as it is proposed in Refs. CheKaw03 ; KaaSal11 ; HuNua10 ; PipTaq00 ; BiaHu08 ; Mis08 .

Going back to the physics of turbulence, let us mention that the fractional Ornstein-Uhlenbeck process as a model of the temporal structure of the velocity field has been already proposed and studied in the literature Sha95 ; Sch12 in the usual sense developed in the mathematical references CheKaw03 ; KaaSal11 ; HuNua10 ; PipTaq00 ; BiaHu08 ; Mis08 . The novelty of the present work, besides the theoretical and numerical aspects mentioned earlier, is the introduction of a new characteristic time scale $\epsilon$ that eventually depends on viscosity, or equivalently on the Reynolds number $\mathcal{R}_{e}$. Following dimensional arguments developed in TenLum72 , we expect $\epsilon$ to coincide with the Kolmogorov time scale $\propto T\mathcal{R}_{e}^{-1/2}$ in a Lagrangian context (Eq. 1) and with the sweeping time scale $\propto T\mathcal{R}_{e}^{-3/4}$ in the Eulerian reference frame (Eq. 2). From a physical point of view, in both cases, it is expected that temporal velocity profiles are smooth below $\epsilon$. In this article, we will mainly focus on the limit $\epsilon\to 0$, but a further modeling step consisting in filtering the Gaussian white noise entering in the construction over the time scale $\epsilon$ could be performed to impose this smooth behavior imposed by viscosity. We keep this aspect for future investigations.

We give in Section II a proper meaning of this stochastic differential equation (Eq. 2) and define the random measure $dW_{H}$ that enters in the dynamics, more precisely a regularized version of it $dW_{\epsilon,H}$ over a small time scale $\epsilon$, and set our notations. The statistical properties in the stationary regime and their limit when $\epsilon\to 0$ of the unique solution of such a differential equation are exposed in Section III. Section IV is devoted to the proofs of the propositions made in Section III in a general framework. We gather all the numerical experiments in Section V and conclude our work in Section VI.

II Setup and notations

We are interested here in studying the statistical properties of the solution $X_{\epsilon,H}(t)$ of the following linear stochastic differential equation

[TABLE]

where the random measure $dW_{\epsilon,H}$ is defined by

[TABLE]

with $H\in]0,1[$, $\alpha>0$ the inverse of a characteristic time scale, $W$ an instance of the Wiener process, and

[TABLE]

a regularized version of the fractional Gaussian noise ManVan68 over the time scale $\epsilon>0$. Remark that the noise $\omega_{\epsilon,H}$ (Eq. 5) is a zero-average and finite-variance Gaussian stationary process for any $\epsilon>0$. Remark also that the very same instance of the Wiener process $W$ enters in both the dynamics (Eq. 3) and in the definition of $\omega_{\epsilon,H}$ (Eq. 5), making these two correlated.

The unique solution $X_{\epsilon,H}(t)$ of Eq. 3 with initial condition, for instance, $X_{\epsilon,H}(t_{0})=0$ can be conveniently written as

[TABLE]

with

[TABLE]

which coincides with the fractional Ornstein-Uhlenbeck (FOU) process CheKaw03 ; KaaSal11 ; HuNua10 in the asymptotic limit $\epsilon\to 0$ for $H>1/2$, and

[TABLE]

a standard Ornstein-Uhlenbeck (OU) process. The very origin of taking the sum of these two correlated Gaussian processes will become clear later when we will present a heuristics, proposed in Ref. Sch03 , that we adapt to our context. Let us now look at the statistical properties of the Gaussian process $X_{\epsilon,H}(t)$ in the stationary regime, if any, and asymptotically when the regularizing scale $\epsilon$ tends to zero.

III Convergence of the covariance in the stationary regime as the regularization scale $\epsilon$ goes to zero

III.1 Statement of the results

For $\alpha>0$ and a given $H\in]0,1[$, the Gaussian process $X_{\epsilon,H}(t)$ (Eq. 6) reaches a stationary regime. Furthermore, it is zero-average and the variance remains bounded when $\epsilon\to 0$. We note

[TABLE]

where enters the Gamma function $\Gamma(z)=\int_{0}^{\infty}x^{z-1}e^{-x}dx$ defined $\forall z>0$. Let us call then $\delta_{\tau}X_{H}$ the corresponding increment over $\tau$, note its variance as

[TABLE]

For $H\in]0,1[$, we have the following behavior of the increment variance at small scales

[TABLE]

III.2 Remarks

Given the statements of Section III.1, we can see that the proposed dynamics (Eq. 3), for a given $\epsilon>0$ and $0<H<1$, converges when $t\to\infty$ towards a stationary process $X_{H}$ which variance remains bounded. Its precise value can be calculated (see the section devoted to proofs) and is strictly positive. Furthermore, the increments of this process $X_{H}$ behave as those of the fractional Brownian motion (Eq. 10).

III.2.1 The case $H=1/2$

The case $H=1/2$ corresponds to a standard Ornstein-Uhlenbeck process $X_{\epsilon,1/2}=X_{1/2}$ since $\omega_{\epsilon,1/2}=0$ (Eq. 5) at any time, of variance $1/(2\alpha)$ and increments behaving as in Eq. 10 (for $H=1/2$).

III.2.2 The case $H\in]1/2,1[$

The case $H\in]1/2,1[$ has already been studied in Ref. CheKaw03 and can be understood easily. For $H>1/2$, the increment of the Wiener process entering in the dynamics (i.e. the third term in the RHS of Eq. 3) will have no contribution when $\epsilon\to 0$ since it is multiplied by a factor $\epsilon^{H-1/2}$. In this case, both variance and increment variance of $X_{\epsilon,H}$ are given by those of $X_{\epsilon,H}^{(1)}$ (Eq. 7) that can be shown to remain bounded when $\epsilon\to 0$ with a proper scaling (Eq. 10) at small scales. The proof relies on the fact that the noise $\omega_{\epsilon,H}$ (Eq. 5) entering in the dynamics, even if its variance diverges when $\epsilon\to 0$ (as expected from a fractional Gaussian noise), has a bounded covariance structure in this limit, such that the variance of $X_{\epsilon,H}^{(1)}$ remains also bounded. We develop these ideas in Section IV.2.2.

III.2.3 The case $H\in]0,1/2[$

The case $H\in]0,1/2[$ is more surprising since both variances of $X_{\epsilon,H}^{(1)}$ (Eq. 7) and $X_{\epsilon,H}^{(2)}$ (Eq. 8) diverge when $\epsilon\to 0$ whereas the variance of $X_{\epsilon,H}=X_{\epsilon,H}^{(1)}+X_{\epsilon,H}^{(2)}$ will remain bounded in this limit. Cancelations in the variance will take place because of the negative correlation existing in between the processes $X_{\epsilon,H}^{(1)}$ and $X_{\epsilon,H}^{(2)}$. This negative correlation originates from the fact that the Gaussian noise $\omega_{\epsilon,H}$ (Eq. 5) is made up of the very same instance of the Wiener process $W$ that enters in the dynamics of $X_{\epsilon,H}$ (Eq. 3).

To justify the form of the infinitesimal increment $dW_{\epsilon,H}$ (Eq. 4) and make a connection with fBm ManVan68 , let us here rephrase some arguments developed in Sch03 . To do so, consider a regularized version $W_{\epsilon,H}(t)$, over $\epsilon>0$, of a fBm of parameter $H$, that is

[TABLE]

Remark that the regularization procedure entering in Eq. III.2.3 is not unique. For instance, we could also define a regularized version of the fBm as its convolution with a mollifier (properly rescaled over $\epsilon$) as it is done in ManVan68 . In all cases, these regularized versions (including Eq. III.2.3) can easily be shown to converge toward the canonical fBm when $\epsilon\to 0$. The regularization procedure that we propose (Eq. III.2.3) allows to compute in a simple way the infinitesimal increment of the process $W_{\epsilon,H}$ that will eventually coincide with the noise $dW_{\epsilon,H}$ that we defined in Eq. 4. Indeed, regrouping terms in a convenient way, we get from Eq. III.2.3

[TABLE]

Performing then a Taylor development (as $dt\to 0$) inside the first integral entering in the RHS of Eq. III.2.3, we recover the contribution of order $dt$, proportional to $\omega_{\epsilon,H}$ (Eq. 5) entering in our initial proposition for the noise $dW_{\epsilon,H}$ (Eq. 4). The second integral entering in the RHS of Eq. III.2.3 justifies the second term $\epsilon^{H-1/2}dW(t)$ entering in Eq. 4. This is indeed true in average and for the variance. More mathematical developments would be needed to fully justify this locally (pathwise). Nonetheless, we can see that the proposition that we made for the noise $dW_{\epsilon,H}$ (Eq. 4) can be justified in a convincing, but not fully rigorous, manner while considering the infinitesimal increment of a regularized version of the fBm (Eq. III.2.3).

Let us also mention recent works FyoKho16 that use a similar, although different, type of regularization procedure to define and use a fractional Gaussian noise for the very rough case $H<1/2$ (see also Ref. Unt09 for similar developments). Their process (see Eq. 1.9 of Ref. FyoKho16 ) shares similar features as the fractional Brownian motion of Ref. ManVan68 when their regularization parameter tends to zero. More work is needed to characterize precisely the differences between the construction of Ref. FyoKho16 and ours, and this is beyond the scope of the present article.

III.2.4 The case $H=0$

The case $H=0$ is of special interest and is fully treated in Ref. PerMor17 in view of applications in theory of fluid turbulence. It can be shown, and this is not purpose of the article (see Ref. PerMor17 for developments on this matter), that indeed the process $X_{\epsilon,0}(t)$ reaches a stationary regime, but the variance diverges when $\epsilon\to 0$ as

[TABLE]

and, furthermore, for $\tau>0$,

[TABLE]

As we can see, $X_{\epsilon,0}$ converges towards a Gaussian process which is logarithmically correlated in time, and thus of infinite variance. This type of random distributions have been studied for some time (see the review Ref. DupRho14 and citations therein), find many applications in the theory of multiplicative chaos RhoVar14 , which is used as a model of the intermittency phenomenon in turbulence (see Fri95 ; PerGar16 ), explaining why it was considered in Refs. Sch03 ; PerMor17 . Once again, the construction of Ref. FyoKho16 is similar to the one of the present work, although different, and more work is needed to compare them precisely.

IV A more general framework and proofs

IV.1 Stochastic integration against fractional Gaussian noise

In order to show the statistical properties of the Gaussian process $X_{\epsilon,H}$ in the limit $\epsilon\to 0$ as announced in Section III.1 (Eqs. III.1 and 10), let us rephrase former considerations in a more general way. We are here interested in calculating the covariance of the Gaussian processes obtained as a linear operation on the Gaussian random measure $dW_{\epsilon,H}(t)$ (Eq. 4) that we recall here for convenience

[TABLE]

where $W(t)$ is a instance of the Wiener process over $t\in\mathbb{R}$, and $\omega_{\epsilon,H}(t)$ the respective causal fractional Gaussian that we consider in Eq. 5. Henceforth, we will only consider convolutions as linear operations, and so only consider stationary processes (of possibly infinite variance), for the sake of simplicity. Following developments could be adapted to non stationary processes, we keep them for future investigations. As far as we are concerned, in particular while being interested by the statistics of the fractional Ornstein-Uhlenbeck process $X_{\epsilon,H}(t)$ (Eq. 6), we can choose as a initial condition $X_{\epsilon,H}(t_{0}=-\infty)=0$, such as $X_{\epsilon,H}(t)$ is directly stationary. With these given precisions, for any suitable test functions $f$ and $g$ that we will specify latterly, let us define the covariance function

[TABLE]

Let us then rephrase in this more general framework the results of Section III.1.

We have, for any $H\in]0,1[$ and a appropriate set of test function $f$ and $g$, the following explicit expression of the correlation function $\mathcal{C}_{\epsilon,H}^{f,g}(\tau)$ when the regularization scale $\epsilon$ tends to 0:

[TABLE]

where the symbol $\star$ stands for the correlation product, i.e.

[TABLE]

and where ′ stands for the derivative. As we will see more precisely in the following, the correlation function $\mathcal{C}_{H}^{f,g}$ (Eq. 15) makes sense for any $H\in]0,1[$ since the function $h^{2H-1}$ is locally integrable everywhere. As far as the test functions $f$ and $g$ are concerned, we are asking them to be such that the derivatives of their correlation product decreases fast enough at large arguments such that the integral entering in Eq. 15 exists. More precise constraints on test functions in terms of functional spaces are developed in particular in Refs. PipTaq00 ; BiaHu08 ; Mis08 ; Jol07 . We will see in Section IV.3 how to apply Eq. 15 to the case of the fractional Ornstein-Uhlenbeck $X_{\epsilon,H}(t)$ in order to prove the results of Section III.1. Let us discuss now and demonstrate the behavior of the correlation function $\mathcal{C}_{H}^{f,g}$ (Eq. 15) given a value of the parameter $H$, in the spirit of the remarks made in paragraph III.2.

IV.2 Dependence of the covariance on $H$

IV.2.1 Integration over the Wiener process

The case $H=1/2$ corresponds to the integration over the Wiener process, known as the Itô integral Nua00 . In this case, the random measure $dW_{\epsilon,1/2}=dW$ (Eq. 13) is the increment of the Wiener process, and is independent on $\epsilon$, and we get the following simple expression for the covariance function (Eq. IV.1)

[TABLE]

Remark that in this case, the expression of the correlation function Eq. 17 corresponds to the more general expression given in Eq. 15 when $H=1/2$.

IV.2.2 Integration over the fractional Gaussian noise when $H>1/2$

The case $H>1/2$ is also well understood in this framework PipTaq00 ; Jol07 . For this range of $H$, we can give a meaning of the limiting value of the correlation function $\mathcal{C}_{\epsilon,H}^{f,g}(\tau)$ (Eq. IV.1) when $\epsilon\to 0$. In other words, besides giving a way to perform numerical simulations of the random measure $dW_{\epsilon,H}(t)$ (Eq. 13), as we will see in Section V, there is no need theoretically to introduce a regularization at the scale $\epsilon$ in the construction, and we can safely take the limit $\epsilon\to 0$ pointwise (see for instance Refs. PipTaq00 ; Jol07 for mathematical developments in the framework of random distributions). Taking a pointwise limit, we have formally

[TABLE]

Inserting then Eq. IV.2.2 into the expression of the correlation function (Eq. IV.1), making the change of variable $h=s_{1}-s_{2}$ and integrating over $s_{2}$, we get

[TABLE]

which eventually coincides with the given general expression for this correlation function (Eq. 15) after splitting the integral in two and performing a integration by parts over the dummy variable $h$. Remark that the intermediate expression of the correlation function (Eq. IV.2.2) makes perfectly sense since the singularity $|h|^{2H-2}$ is locally integrable in the neighborhood of the origin for $1/2<H<1$.

IV.2.3 Integration over the fractional Gaussian noise when $H<1/2$

This higher level of roughness requires some more work since the correlation function of the fractional Gaussian noise $\omega_{\epsilon,H}$ entering in Eq. IV.2.2 has no meaning in the limit $\epsilon\to 0$. The mechanism of regularization over $\epsilon$ will here play a key role, and will allow several key cancellations of diverging quantities, such that the correlation function $\mathcal{C}_{\epsilon,H}^{f,g}(\tau)$ (Eq. IV.1) remains a bounded function of $\epsilon$. Remark first that, formally, we can write

[TABLE]

where the contribution of the Wiener process corresponding to $\mathbb{E}[dW(s_{1})dW(s_{2})]=\delta(s_{1}-s_{2})ds_{1}ds_{2}$ can be conveniently noted with a Dirac function $\delta$. Let us work out first the contribution coming from the fractional Gaussian noise $\omega_{\epsilon,H}$. It reads

[TABLE]

which is indeed a function of $s_{1}-s_{2}$, and is not bounded in $\epsilon$. To extract diverging quantities, perform a integration by parts and obtain

[TABLE]

Noticing that we have formally,

[TABLE]

we can see that, once inserted in Eq. 20, the first diverging term entering in the RHS of Eq. IV.2.3 will cancel out with the contribution of Eq. IV.2.3. It will only remain in Eq. 20 the second term of he RHS of Eq. IV.2.3 and the $\epsilon^{2H-1}\delta(s_{1}-s_{2})ds_{1}ds_{2}$ term. Once inserted in the expression of the correlation function (Eq. IV.1), we end up with

[TABLE]

Perform then a integration by parts over $h$ in the remaining integral, we can see that the first term of the RHS of the former equation (of order $\epsilon^{2H-1}$) will be compensated, and we obtain

[TABLE]

which coincides with the expression given in Eq. 15 once is performed the integral over $u$.

IV.3 Fractional Ornstein-Uhlenbeck processes

As a application of the formula given in Eq. 15 to our center of interest, namely the fractional Ornstein-Uhlenbeck processes $X_{\epsilon,H}(t)$ (Eq. 6), with initial condition $X_{\epsilon,H}(t_{0}=-\infty)=0$, consider the test functions

[TABLE]

such that indeed

[TABLE]

Making use of

[TABLE]

once inserted in Eq. 15, we obtain

[TABLE]

which shows once simplified the proposition made in Eq. III.1.

Similarly, we obtain for the increment variance, assume for instance $\tau>0$,

[TABLE]

Regrouping terms in a convenient way, we obtain

[TABLE]

As $\tau\to 0$, the first term decreases toward 0 as $\tau^{2}$, and thus will be negligible in front of the second term that will behave as $\tau^{2H}$. To see this, rescale in the second term the dummy variable $h$ by $\tau$, then take safely the limit $\tau\to 0$ inside the integral such that to get

[TABLE]

which coincides with the power-law announced in Eq. 10 once simplified.

To finish this Section, let us focus on the boundary case $H=0$. This case, fully developed in Ref. PerMor17 , deserves more care since we are dealing with infinite variance processes as mentioned in Section III.2.4 in the limit $\epsilon\to 0$. This logarithmic divergence with $\epsilon$ of the variance $\mathcal{C}_{\epsilon,0}^{f,f}(0)$ can be seen in Eq. IV.2.3. Nonetheless, as we already explained, even if the variance of such processes diverge with $\epsilon$ (logarithmically), the correlation function $\mathcal{C}_{0}^{f,f}(\tau)$ remains bounded for $\tau>0$. Using the general expression given in Eq. 15, applied to the kernel $f$ of the fractional Ornstein-Uhlenbeck process (Eq. 25), we obtain formally, noticing that the prefactor $\frac{1}{2}\frac{1}{\sin(\pi H)}\frac{\left[\Gamma\left(H+\frac{1}{2}\right)\right]^{2}}{\Gamma(2H)}$ tends to 1 when $H\to 0$,

[TABLE]

Since the first term entering in the RHS of Eq. IV.3 remains bounded when $\tau$ gets smaller and smaler, this demonstrates the logarithmic diverging behavior of the correlation function as $\tau\to 0$, as claimed in Section III.2.4. This can be readily seen while performing a further integration by parts over the dummy variable $h$ entering in the second term of the RHS side of the equation.

V Numerical simulations

V.1 A periodized approximation of the fractional Gaussian noise

We here propose a numerical method that allows to estimate the trajectories of the solution $X_{\epsilon,H}(t)$ for $H\in]0,1[$ and $\epsilon>0$. To do so, we need to come up with a approximation of the fractional Gaussian noise $\omega_{\epsilon,H}(t)$ (Eq. 5) entering in the dynamics (Eq. 3). The first idea would be to truncate, say over a large time scale $T^{\prime}$, the Wiener integral entering in (Eq. 5), and define accordingly the corresponding estimator $\hat{\omega}_{\epsilon,H,n^{\prime}}(t)$ as

[TABLE]

where $\Delta t$ is the resolution time scale of this numerical problem, $n^{\prime}=\left\lfloor T^{\prime}/\Delta t\right\rfloor$ the integer corresponding to the large time scale $T^{\prime}$ of the truncation, and $\Delta W(t_{i})=\sqrt{\Delta t}\mathcal{N}(0,1)$ is a discrete collection of the increment at time $t_{i}=i\Delta t$ of the underlying wiener process, and is made up of independent zero-average Gaussian random variable of variance $\Delta t$. As we will see in the following, the regularizing time scale $\epsilon$ is chosen as a multiple of the resolution time scale, typically $\epsilon=10\Delta t$. The large-scale truncation $T^{\prime}$ should be chosen, depending of the values of $\Delta t$, $\epsilon$ and $H$, large “enough”. Recall that the fractional noise $\omega_{\epsilon,H}(t)$ (Eq. 5) is a well defined random process for a finite $\epsilon>0$, so its estimator $\hat{\omega}_{\epsilon,H,n^{\prime}}$ (Eq. V.1) should become independent on $n^{\prime}$ as $n^{\prime}\to\infty$. The reason that explains this independence on $T^{\prime}$ is connected to the fact that the kernel $(t+\epsilon)^{H-3/2}$ decreases fast enough, such that its square (entering in the variance of $\omega_{\epsilon,H}$) is integrable at $t\to\infty$. Remark also that the estimator $\hat{\omega}_{\epsilon,H,n^{\prime}}$ (Eq. V.1) requires of the order $n^{\prime}$ operations at each time step, which is numerically demanding. This numerical strategy has been nonetheless followed in Ref. PerMor17 in order to simulate and explore a more complex process involving a tensorial and non Gaussian generalization of the noise $\omega_{\epsilon,H}(t)$.

Thus, instead of using Eq. V.1, in order to minimize the error made while truncating the integral entering in Eq. 5 and perform numerical simulations in a efficient way, we will rely on the discrete Fourier transform and approximate $\omega_{\epsilon,H}(t)$ (Eq. 5) by a periodical estimator $\widetilde{\omega}_{\epsilon,H}$. For full benefit of the fast Fourier transform algorithm, we consider $N=2^{k}$ with $k\in\mathbb{N}^{*}$. More precisely, call $N$ the number of collocation points of your numerical approximation $\widetilde{\omega}_{\epsilon,H}$ and set $T_{0}$ the physical time duration of the trajectory, such that $\Delta t=T_{0}/N$. Define $t_{n}=n\Delta t$ for $n\in[0,N-1]$ and $t_{n}^{(N)}$ its periodized version, i.e. $t_{n}^{(N)}=t_{n}$ for $n\leq N/2$ and $t_{n}^{(N)}=t_{n}-T_{0}$ for $N/2+1\leq n\leq N-1$ . Define then the regularized over $\epsilon$, periodized and causal kernel $\varphi(t_{n})=(H-1/2)(t^{(N)}_{n}+\epsilon)^{H-3/2}1_{t_{n}^{(N)}\geq 0}$. Note by $\mathcal{F}$ the discrete Fourier transform. We thus get a periodized approximation $\widetilde{\omega}_{\epsilon,H}$ of the fractional Gaussian noise $\omega_{\epsilon,H}(t)$ (Eq. 5) taking

[TABLE]

where again, $\Delta W(t_{n})$ are $N$ independent realizations of a zero-average normal random variable of variance $\Delta t$. Trajectories of $X_{\epsilon,H}(t)$ are finally obtained while integrating, using a Euler discretization scheme, their dynamics (Eq. 3) as, starting for example with $X_{\epsilon,H}(0)=0$,

[TABLE]

We recall here that the very same instance of the white noise $\Delta W$ enters both at the level of the stochastic differential equation (Eq. V.1) and in the definition of the noise $\widetilde{\omega}_{\epsilon,H}$ (Eq. 28). The additional implied correlation between $X_{\epsilon,H}$ and $\Delta W$ is crucial and plays a key role in the statistical properties of $X_{\epsilon,H}$ in the stationary regime.

V.2 Numerical results

We consider the simulation of trajectories $X_{\epsilon,H}(t)$ of the stochastic differential equation Eq. 3 under the approximations developed in the former Section. We choose $T_{0}=1$ and time is nondimensionalized accordingly. In order to minimize any statistical effects of the transitory regime, we choose $\alpha=50/T_{0}$, since we expect $X_{\epsilon,H}(t)$ to be correlated over typically the time scale $1/\alpha$. Present simulations are performed using $N=2^{28}$ collocation points, such that $\Delta t=T_{0}/N\approx 6.10^{-8}T_{0}$. We use for the regularization scale the value $\epsilon=10\Delta t$, and generate and analyze the statistical properties of $150$ independent trajectories of $X_{\epsilon,H}(t)$, for various values of $H$.

We represent in Fig. 1 the estimation of the variance of the simulated trajectories of the process $X_{\epsilon,H}(t)$ (Eq. V.1) as a function of the parameter $H$. We indeed observe a stationary regime, and compute the variance from $150$ independent trajectories, from which we estimate the error bars. We furthermore compare with our analytical asymptotic prediction (Eq. III.1). Let us first mention that this comparison is more and more demanding as $H$ increases, since we are rescaling the estimated variance by a factor $\alpha^{-2H}$ that may become very small as $H$ approaches unity. This being said, we indeed observe that our prediction (Eq. III.1) is compatible with the variance estimated on our trajectories for $0.3\lesssim H\lesssim 0.8$. For $H\lesssim 0.3$, our prediction does not reproduce the observed variance which is characterized by a strong variability. This could be explain by several facts: (i) as a general remark, the statistical convergence of such a large-scale quantity as the variance requires many realizations and we may miss some of them, (ii) $H\lesssim 0.3$ corresponds to the very rough case, this may require to take the small scale regularization $\epsilon$ to be taken larger than what we chose (recall that here $\epsilon=10\Delta t$), at the cost of missing the scaling properties in the asymptotics, (iii) the smoother cases $0.8\lesssim H$ have a strong statistical variability as shown by extended error bars, this might be due to a lack of statistical convergence, or a slow convergence towards the asymptotic $\epsilon\to 0$ regime. Overall, given the aforementioned limitations, our predictions seem to reproduce in a acceptable manner the variance of the simulated trajectories over a extended range of values of $H$. Let us add that numerical tests have been performed over half of the samples (i.e. 75 trajectories) without a quantitative change of the amplitudes of the error bars (data not shown), showing that discrepancies between estimated variances and theoretical prediction can be barely minimized while working on a larger set of realizations.

We represent in Fig. 2 the scaling behavior of the so-called structure function of second order, namely the variance $\mathbb{E}\left(\delta_{\tau}X_{\epsilon,H}\right)^{2}$ of the increments of the process $X_{\epsilon,H}(t)$ as a function of the scale $\tau$, for various values of the parameter $H$, from the roughest case $H=0.1$ to the smoothest case $H=0.9$. We indeed observe a power-law behavior $\tau^{2H}$ for any of the values of the parameter $H$. We superimpose on this representation the predicted behavior in the asymptotic limit $\epsilon\to 0$ (Eq. 10) and observe that indeed this prediction reproduces both the scale-dependence, and the prefactors. We can see also that the comparison between predictions and estimated variances deteriorates for the smallest and largest values of the parameter $H$, as it is also observed in Fig. 1. Since we are studying here the small scales of the process, that benefit from a large statistical sampling, we infer that we might not have reached the asymptotic regime of vanishing regularization scale $\epsilon\to 0$.

In order to quantify precisely the differences between the observed power-laws of the increment variance and our asymptotical prediction (Eq. 10) as depicted in Fig. 2, we represent in Fig. 3 the compensated variance $\mathbb{E}\left(\delta_{\tau}X_{\epsilon,H}\right)^{2}$ by our analytical prediction (Eq. 10), for 7 different values of the parameter $H$. For the sake of clarity, we superimpose also with vertical dashed lines the two characteristic time scales $\epsilon$ and $1/\alpha$ in between which we expect a power law behavior of exponent $2H$. We indeed observe that over almost three decades in scale, increment variance obeys a extended power-law behavior correctly captured by our analytical prediction (Eq. 10) when $H>1/2$. We can also observe that our prediction deteriorates as $H$ gets smaller and smaller compared to $1/2$. Similarly, we can see two reasons for this: (i) we did not reach yet the asymptotic regime $\epsilon\to 0$ and we are observing a slow convergence towards it, and (ii) some consequences are expected in the very rough case if we chose the regularizing scale $\epsilon$ not big enough compared to the resolution scale $\Delta t$. We keep for future investigations a more developed numerical study of this process. Nonetheless, we can see that the fractional Ornstein-Uhlenbeck process can be easily simulated and our theoretical predictions (Eqs. III.1 and 10) compare reasonably well with our numerical results.

VI Conclusions

We have proposed a generalization (Eq. 3) of the Ornstein-Uhlenbeck process (Eq. 1) which is of finite-variance in the stationary regime, and which regularity is governed by the parameter $H\in]0,1[$ when the regularizing scale tends to zero, in the same fashion as the fractional Brownian motion ManVan68 . This article includes a theoretical study aimed at getting exact expressions for the variance and increment variance of the process, and their asymptotical behavior when $\epsilon\to 0$ and at vanishing scale. We propose then a numerical study showing that the fractional Gaussian noise entering in the dynamics can be approximated in a accurate way, and observe to some extend for any $H\in]0,1[$ the convergence toward the asymptotical regime developed in the theoretical section.

Further efforts in this spirit will be devoted to the modeling of fluid turbulence, that asks for further developments in order to give a realistic picture of the intermittency phenomenon and energy transfers in scale. We keep these perspectives for future investigations.

I would like to thank P. Abry, P. Borgnat, G. Didier, P. Flandrin, C. Garban, K. Gawedzki, L. Moriconi, R.M. Pereira, R. Rhodes, S.G. Roux and V. Vargas for many discussions concerning this matter, and J.F. Palierne for a critical proofreading of a first version of the manuscript. The author is supported by ANR grants Liouville ANR-15-CE40-0013 and MultiFracs ANR-16-CE33-0020.