# Space-time statistics of a linear dynamical energy cascade model

Gabriel B. Apolinário<EMAIL_ADDRESS>Univ Lyon, ENS de
Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique, 46 allée d’Italie
F-69342 Lyon, France  Laurent Chevillard  Univ Lyon, ENS de Lyon, Univ Claude
Bernard, CNRS, Laboratoire de Physique, 46 allée d’Italie F-69342 Lyon, France

###### Abstract

A linear dynamical model for the development of the turbulent energy cascade
was introduced in Apolinário _et al_ [arxiv 2107.03309]. This partial
differential equation, randomly stirred by a forcing term which is smooth in
space and delta-correlated in time, was shown to converge at infinite time
towards a state of finite variance, without the aid of viscosity. Furthermore,
the spatial profile of its solution gets rough, with the same regularity as a
fractional Gaussian field. We here focus on the temporal behavior and derive
explicit asymptotic predictions for the correlation function in time of this
solution and observe that their regularity is not influenced by the spatial
regularity of the problem, only by the correlation in time of the stirring
contribution. We also show that the correlation in time of the solution
depends on the position, contrary to its correlation in space at fixed times.
We then investigate the influence of a forcing which is correlated in time on
the spatial and time statistics of this equation. In this situation, while for
small correlation times the homogeneous spatial statistics of the white-in-
time case are recovered, for large correlation times homogeneity is broken,
and a concentration around the origin of the system is observed in the
velocity profiles. In other words, this fractional velocity field is a
representation in one-dimension, through a linear dynamical model, of the
self-similar velocity fields proposed by Kolmogorov in 1941, but only at fixed
times, for a delta-correlated forcing, in which case the spatial statistics is
homogeneous and rough, as expected of a turbulent velocity field. The
regularity in time of turbulence, however, is not captured by this model.

Keywords— fractional Gaussian fields, statistical theory of turbulence,
stochastic partial differential equations, homogeneous operators of degree 0,
energy cascade, pseudo-spectral numerical simulation, Ornstein-Uhlenbeck
process

##  1 Introduction

Since the important work of Taylor and Richardson in the first half of the XX
century, a direct energy cascade has been invoked to interpret the rough
nature of velocity fluctuations that are observed in Navier-Stokes turbulence
in three dimensions [1]. This picture has been further supplemented by the
statistical description provided by Kolmogorov in 1941 [2], in which
fractional stochastic processes play a paradigmatic role in modeling the rough
behavior of incompressible velocity fields. The energy cascade can then be
seen as the dynamical process through which energy is transported from the
large scale smooth external force to the small scales of the flow, where it
generates rough random fields. Furthermore, it is conjectured that this rough
stationary state is reached even in the absence of viscous dissipation [3, 4,
5, 6].

Consider that $u_{i}(x,t)$, $i\in\\{1,2,3\\}$, $x\in\mathbb{R}^{3}$ is a
divergence-free velocity field, solution of the incompressible Navier-Stokes
equations. The statistical behavior of turbulence in [2], based on two
propositions regarding symmetries and self-similarity of these solutions,
describes a fluid of kinematic viscosity $\nu$, stirred at a large scale $L$
by an external force, and determines a power-law scaling for the statistical
moments of increments of the velocity field (also called structure functions
of order $q$),

$\lim_{\nu\to 0}\mathbb{E}(\delta_{\ell}u_{i})^{q}\underset{\ell\to
0}{\sim}c_{q}\left(\frac{\ell}{L}\right)^{qH}\ ,$ (1.1)

where the velocity increments are defined as
$\delta_{\ell}u_{i}(t,x)=u_{i}(t,x+\ell)-u_{i}(t,x)$, and the constants
$c_{q}$ do not depend on the exact setup of the flow, but only on the large
scale $L$ and on the mean energy dissipation rate. The exponent $H$ is
approximately $1/3$, as predicted by dimensional analysis and verified by
numerical and experimental measurements [1, 7, 8, 9]. It should be remarked
that the characterization of Eq. (1.1) has been later supplemented by the
intermittency hypothesis, in which the structure function exponents do not
grow linearly with the order of the statistical moments, meaning that the
underlying velocity field is not exactly self-similar, but rather a
multifractal [1].

In this setting, the fractional Gaussian field of exponent $H=1/3$ can be seen
as the simplest model of a turbulent velocity field, with statistics given by
Eq. (1.1). These mathematical objects, introduced by Kolmogorov (in the
framework we henceforth refer to as K41) and named by Mandelbrot and Van Ness
[10, 11, 12], are uniquely characterized by the Hölder exponent $H$, which
describes their local regularity, coinciding with the behavior of Brownian
motion for $H=1/2$, while being rougher than that for $H<1/2$ and smoother for
$H>1/2$. Nevertheless, being Gaussian and of mean zero, all of the odd order
increment moments of fractional Gaussian fields vanish, in contrast with the
phenomenology of turbulence, where strong skewness has been observed in the
statistical distribution of velocity increments, and where an important exact
result on the third order structure function has been obtained, given by
$\mathbb{E}(\delta_{\ell}u_{i})^{3}=-\frac{4}{5}\varepsilon\ell$, with
$\varepsilon$ being the mean energy dissipation rate [1, 13]. This discrepancy
has lead to several proposals for synthetic random fields with a local
regularity prescribed either by Eq. (1.1) or by its intermittent extensions
[14, 15, 16, 17].

Nevertheless, these models are static, and do not address the fundamental
dynamics of the formation of small scales from the random advection of large
scale flows, established by the external forcing. Some success in the
heuristics of turbulent statistics from a dynamical description has been
obtained in the framework of shell models [18, 19, 20], which are systems of
nonlinear differential equations for the variables $u_{k}$, called shells,
inspired by the Navier-Stokes equations in Fourier space. Each shell has a
characteristic wavelength $k$ and interacts only with its closest neighbors in
a discretized, often logarithmic, $k$-space, as a way to model the local
transport of energy in Fourier space observed in turbulence. With constraints
such as energy and helicity conservation, and a with large scale forcing (i.e.
a forcing concentrated around shells of wavelength $k_{0}$, where $k_{0}$ is
small), a scaling behavior very similar to that of Navier-Stokes is observed,
albeit only in a numerical context.

The present work addresses the modeling of self-similar random velocity fields
in a dynamical framework, building upon the results of [21], in which a
nonlinear partial differential equation (PDE) for a complex one-dimensional
velocity field was introduced. The stationary solution of this PDE, when
subject to a random force on the large scales, is multifractal and fractional,
as controlled respectively by the intermittency parameter $\gamma$ and by the
Hölder parameter $H$, and as a step to its construction, a linear PDE whose
stationary state is a fractional Gaussian field of parameter $H$ is presented
as well. Small scales are built dynamically by an homogeneous operator of
degree 0 [22, 23, 24], which has been identified in the phenomenon of
concentration of internal waves around geometrical attractors in stably
stratified flows [25, 26, 27, 28], and whose discretized version has been used
to describe a linear cascade of energy in the context of shell models [29].

For the linear PDE introduced in [21], consider a complex velocity field
$u_{H,\nu}(t,x)$, defined on $x\in\mathbb{R}$ and characterized by a positive
kinematic viscosity $\nu$ and a Hölder exponent $H\in]0,1[$. This velocity
field is governed by

$\partial_{t}u_{H,\nu}(t,x)=P_{H}\mathcal{L}P_{H}^{-1}u_{H,\nu}(t,x)+\nu\partial_{x}^{2}u_{H,\nu}(t,x)+f(t,x)\
,$ (1.2)

where temporal and spatial derivatives are denoted respectively by
$\partial_{t}$ and $\partial_{x}$. A viscous dissipation term, proportional to
the viscosity $\nu$ and an external large scale force $f$ are present. The
operator $\mathcal{L}$ on the right-hand side of Eq. (1.2) is identified with
the homogeneous operator of degree 0 described in [22, 23, 24] and its action
on a general velocity field $u(t,x)$ is defined as

$\mathcal{L}u(t,x)\equiv 2i\pi cxu(t,x)\ ,$ (1.3)

where $c$ is a constant. This term is responsible for the transport of energy
from the large scale $L$ to the small scales, where viscous dissipation acts
strongly in regularizing the velocity field. It can be seen that the dynamics
induced by $\mathcal{L}$ is a linear transport in Fourier space, towards
increasing wavelengths in the case of $c>0$. In this paper, we always consider
$c$ to be positive, but the case of $c<0$ is equivalent, with the energy
cascade flowing to large negative wavelengths. Ultimately, the operator
$P_{H}$, responsible for the local regularity of order $H$ of the stationary
state solution, is defined by

$P_{H}u(t,x)=\int_{\mathbb{R}}e^{2\pi
ikx}\frac{1}{|k|_{1/L}^{H+1/2}}\widehat{u}(t,k)\,dk\ ,$ (1.4)

where $|k|_{1/L}$ is a regularization of the norm $|k|$ on small wavelengths,
whose exact description is not crucial, and $\widehat{u}$ represents a Fourier
transform of $u$, with the chosen definition

$\mathcal{F}[u](t,k)\equiv\widehat{u}(t,k)=\int_{\mathbb{R}}\,e^{-2i\pi
kx}u(t,x)\,dx\ .$ (1.5)

The combined action of $\mathcal{L}$ in generating a delta-correlated velocity
field and of $P_{H}$ in building power-law correlations from this delta-
correlated field lead to the formation of a stationary state which resembles a
fractional Gaussian field, as explained in [21].

The forcing contribution of Eq. (1.2) is taken to be a complex Gaussian random
field of mean zero and covariance

$\mathbb{E}[f(t,x)f^{*}(t^{\prime},y)]=\chi(t-t^{\prime})C_{f}(x-y)\ ,$ (1.6)

where $*$ represents complex conjugation. Moreover, the real and imaginary
components of $f$ are independent, meaning that
$\mathbb{E}\left[f(t,x)f(t^{\prime},y)\right]=0$. $\chi$ and $C_{f}$ are even
functions of their arguments, representing the correlations in time and in
space of the external force, respectively. These functions have finite
characteristic scales, a length $L$ for $C_{f}$ and a time $T$ for $\chi$,
beyond which they rapidly vanish. The scale $L$ is furthermore identified in
the turbulence phenomenology as the integral length scale, at which energy is
injected in the flow.

The case of white in time forcing (that is, $\chi(t)$ is a Dirac delta
function) has been considered in [21] and is revisited in this work, where we
investigate the properties in time and in space of the solution $u_{H,\nu}$.
Additionally, we examine the case of a finite correlation time $T$ and
establish a connection with the previous results in the limit $T\to 0$. We
choose the initial condition $u_{H,\nu}(0,x)=0$, and, for the correlation
functions, we specify $\chi(t)=e^{-|t|/T}/2T$ and $C_{f}(x)=e^{-x^{2}/2L^{2}}$
as particular cases which have been applied to the numerical simulations.

This paper is organized as follows. In Sec. 2, we discuss the statistical
properties in space and in time of the solution of Eq. (1.2) driven by a
delta-correlated external force, with analytical results and numerical
simulations. In Sec. 3, these statistical properties are revisited for a force
of finite correlation time, with numerical simulations and analytical results.
Concluding remarks and future perspectives are drawn in Sec. 4.

##  2 Dynamical fractional Gaussian fields

In this section we discuss the statistical properties of the stationary
solution of Eq. (1.2) under delta-correlated forcing in time. We begin with
recalling the spatial statistics that were first analyzed in [21], and we then
introduce statistical properties in time as well. This is also an opportunity
to present the details of our numerical simulations and directly compare them
with the analytical results.

###  2.1 Statistical properties for a delta-correlated in time forcing term

####  2.1.1 Spatial statistics

As discussed in [21], the solution of Eq. (1.2), $u_{H,\nu}$, reaches a
statistically stationary state of finite variance even for vanishing
viscosity, and spatial slices of this stationary solution resemble fractional
Gaussian fields. At zero viscosity, the system only reaches this exactly self-
similar description at infinite time, since the external energy injected at
scale $L$ by the forcing term is linearly transported in $k$-space to higher
frequencies, and an infinite amount of time is necessary for arbitrarily small
wavelengths to be manifest. In the presence of a finite viscosity, though,
wavelengths higher than some typical dissipative wavelength are strongly
suppressed, and the system reaches its stationary state in a finite time.
Nevertheless, since both the solutions at finite or zero viscosity are
statistically stationary, their statistical features are very similar,
differing only in the dissipative range of wavelengths, where viscosity
regularizes the rough local behavior of the solution. Therefore, we discuss
the solution of the inviscid fractional equation

$\partial_{t}u_{H,0}(t,x)=P_{H}\mathcal{L}P_{H}^{-1}u_{H,0}(t,x)+f(t,x)\ ,$
(2.1)

with initial condition $u_{H,0}(t=0,x)=0$ and $f$ a large scale Gaussian
forcing, smooth in space, exactly as in Eq. (1.6). We first consider the
delta-correlated case, $\chi(t)=\delta(t)$. This equation admits an analytical
solution in its Fourier space formulation,

$\widehat{u}_{H,0}(t,k)=|k|_{1/L}^{-H-\frac{1}{2}}\int_{0}^{t}|k-cs|_{1/L}^{H+\frac{1}{2}}\widehat{f}(t-s,k-cs)ds\
,$ (2.2)

for any forcing scheme. The small wavelength regularization might be defined
as $|k|_{1/L}=(k^{2}+L^{-2})^{1/2}$, a choice which is important in the
numerical simulations. The solution in Eq. (2.2) is Gaussian since it is
defined as a linear operation on the forcing $f$, which is itself chosen
Gaussian. Hence, to characterize the solution from a statistical point of
view, it suffices to investigate second order observables, namely the power
spectrum and the structure function, which elucidate the nature of the energy
cascade and of the local regularity of the velocity field.

In its stationary state, obtained as the limit of $t\to\infty$, Eq. (2.2) is a
rough velocity field, correlated over the large scale $L$. Its local
regularity can be described by the power-law decay of the power spectrum
density, defined as

$\widehat{C}_{H,0}(t,k)=\mathcal{F}\,\mathbb{E}[u_{H,0}(t,x+\ell)u_{H,0}^{*}(t,x)](t,k)\
,$ (2.3)

where the Fourier transform is performed over the separation variable $\ell$,
since the correlation function in this equation is homogeneous, as can be seen
from the analytical solution. The exact expression for the power spectrum of
$u_{H,0}$ is

$\widehat{C}_{H,0}(t,k)=|k|_{1/L}^{-(2H+1)}\int_{0}^{t}|k-cs|_{1/L}^{2H+1}\widehat{C}_{f}(k-cs)ds\
,$ (2.4)

which in the limit of large wavelengths, in the stationary state, becomes

$\lim_{t\to\infty}\widehat{C}_{H,0}(t,k)\underset{k\to\infty}{\sim}b_{H}\,|k|^{-(2H+1)}$
(2.5)

with a constant given by

$b_{H}=\frac{1}{c}\int_{\mathbb{R}}|s|_{1/L}^{2H+1}\widehat{C}_{f}(s)ds\ ,$
(2.6)

an expression which can be explicitly evaluated in terms of special functions
in the case of a Gaussian correlation function $\widehat{C}_{f}$, and can also
be easily numerically integrated in a general setting. Eq. (2.5) illustrates
the direct energy cascade, in which higher wavelengths are progressively
excited, linearly in time. As previously stated, the variance of the solution
$u_{H,0}$ is finite, even with a vanishing viscosity. It is given by

$\mathbb{E}|u_{H,0}|^{2}=\lim_{t\to\infty}\int_{\mathbb{R}}\widehat{C}_{H,0}(t,k)\,dk\
,$ (2.7)

and its finiteness is a consequence of the integrability of the power
spectrum.

The local regularity of the solution is also measured by the small scale
asymptotic statistics of the increments $\delta_{\ell}u_{H,0}(t,x)$, defined
as

$\delta_{\ell}u_{H,0}(t,x)=u_{H,0}(t,x+\ell)-u_{H,0}(t,x)\ .$ (2.8)

The increment variance, or second order structure function, is asymptotically
given by the power-law relation

$\lim_{t\to\infty}\mathbb{E}|\delta_{\ell}u_{H,0}|^{2}\underset{\ell\to
0^{+}}{\sim}c_{H}\,\ell^{2H}\ ,$ (2.9)

where the constant is

$c_{H}=\frac{(2\pi)^{2H+1}}{2c\sin(\pi
H)\Gamma(1+2H)}\int_{\mathbb{R}}|s|_{1/L}^{2H+1}\widehat{C}_{f}(s)\,ds$ (2.10)

and $\Gamma$ represents the Gamma function. Higher order expectations of the
increments, $\mathbb{E}(\delta_{\ell}u_{H,0})^{q}$, can be obtained directly
from simple combinatorial arguments, since the increments are Gaussian as
well. Furthermore, $q$-order expectation values display an exponent linear in
$H$, and the expectation value themselves vanish for all odd $q$, as a
consequence of the vanishing mean of the velocity field.

This representation of a turbulent velocity field is thus compatible with the
portrait of Eq. (1.1). With the Hölder parameter chosen as $H=1/3$, the power-
law decay in Eq. (2.5) has an exponent $-5/3$ and the increment variance has
an exponent $2/3$, two values which are widely regarded as universal in the
statistical description of turbulence, and which are given by the K41
framework. Nonetheless, we remark that the multifractal description of
turbulence induces a small deviation to the value $1/3$, which requires the
introduction of an intermittency parameter. The nonlinear equation addressed
in [21], where the intermittency parameter governs the intensity of the
nonlinearity, induces non-Gaussian statistics on its solution, which capture
both the multifractal behavior and the nontrivial nature of the odd order
structure functions of real turbulent velocity fields.

####  2.1.2 Temporal statistics

We now turn to the statistical properties in time of the inviscid solution. We
observe that the $P_{H}$ operator does not significantly influence the
temporal behavior of $u_{H,0}$, and that its regularity in time is solely
governed by the shape of the forcing correlation in time. The regularity in
time is again investigated through increments $\delta_{\tau}u_{H,0}(t,x)$,
analogously defined as

$\delta_{\tau}u_{H,0}(t,x)=u_{H,0}(t+\tau,x)-u_{H,0}(t,x)\ .$ (2.11)

It is possible to write the following expression for the temporal second-order
structure function,

$\lim_{t\to\infty}\mathbb{E}|\delta_{\tau}u_{H,0}(t,x)|^{2}=\lim_{t\to\infty}2\mathbb{E}|u_{H,0}(t,x)|^{2}-2\mathcal{R}\,\mathbb{E}[u_{H,0}(t+\tau,x)u_{H,0}^{*}(t,x)]\
,$ (2.12)

in terms of the variance and the real part of the two-point correlation-
function in time, equal to

$\begin{split}\lim_{t\to\infty}\mathbb{E}[u_{H,0}(t+\tau,x)u_{H,0}^{*}(t,x)]&=\frac{1}{c}\int_{k\in\mathbb{R}}e^{2\pi
ic\tau x}|k|_{1/L}^{-H-1/2}|k+c\tau|_{1/L}^{-H-1/2}\\\
&\times\int_{s=-\infty}^{k}|s|_{1/L}^{2H+1}\widehat{C}_{f}(s)\,dsdk\
.\end{split}$ (2.13)

This equation shows an explicit oscillatory behavior with the position $x$,
which is inherited by the structure function. To see the behavior of the
temporal second-order structure function at small scales $\tau\to 0$, we
perform a Taylor expansion of the integrand and get

$\lim_{t\to\infty}\mathbb{E}|\delta_{\tau}u_{H,0}(t,x)|^{2}=d_{1,H}\,|\tau|+d_{2,H}(x)\,\tau^{2}+\mathcal{O}(\tau^{3})\
,$ (2.14)

where the constants are given by

$\begin{split}&d_{1,H}=(2H+1)\int_{k\in\mathbb{R}}k|k|_{1/L}^{-2H-3}\int_{s=-\infty}^{k}|s|_{1/L}^{2H+1}\widehat{C}_{f}(s)\,dsdk\
\mbox{, and}\\\
&d_{2,H}(x)=4c\pi^{2}x^{2}\int_{k\in\mathbb{R}}|k|_{1/L}^{-2H-1}\int_{s=-\infty}^{k}|s|_{1/L}^{2H+1}\widehat{C}_{f}(s)\,dsdk\\\
&+\frac{c}{4}(2H+1)\int_{k\in\mathbb{R}}\left[2/L^{2}-(2H+3)k^{2}\right]|k|_{1/L}^{-2H-5}\int_{s=-\infty}^{k}|s|_{1/L}^{2H+1}\widehat{C}_{f}(s)\,dsdk\
.\end{split}$ (2.15)

Eq. (2.14) says that the spatio-temporal field $u_{H,0}(t,x)$ has the same
local regularity in time as a Brownian motion, as given by the linear term in
the expansion. Remark that, contrary to its spatial behavior, where it has
been pinpointed a regularity similar to a fractional Brownian motion of
parameter $H$ (Eq. 2.9), the regularity in time is independent of $H$. It can
be observed that the integrals in Eqs. (2.14) and (2.15) are finite given that
the function $\widehat{C}_{f}$ decays fast enough. A Gaussian decay for
instance, which is employed in the numerical simulations, is sufficient.
Furthermore, $d_{1,H}$ is positive, since it can be written as an integral
over positive terms only, and the first contribution of $d_{2,H}$ is already
expressed as an integral of positive terms. The second term of $d_{2,H}$,
which does not depend on $x$, can be shown to be positive for the specific
regularization we employ through an integration by parts, even though its
exact value is hard to obtain. We reinforce that the $d_{1,H}$ term is
independent of $x$ and of $c$, which play a role only in the second order
contribution.

An important quantity in the temporal dynamics of $u_{H,0}$ is its integral
correlation time, which is usually defined as a time integral of the two-point
correlation function in time (Eq. 2.13), normalized by the variance.
Nevertheless, it can be observed that Eq. (2.13) is not integrable at
$\tau\to\infty$ for $H<1/2$. For this reason, in order to obtain a finite
value of the integral correlation time, we propose an estimation based on the
Taylor expansion of the second-order structure function in time. In the limit
of large time intervals, Eq. (2.14) is no longer valid and the temporal
second-order structure function evolves to the homogeneous value of twice the
variance of $u_{H,0}$, as can be seen in the exact expression, Eq. (2.12). We
then define the position-dependent correlation time $T_{H}(x)$ as the value in
which the quadratic expansion equals the asymptotic value at large times. The
positive solution of this quadratic equation is

$T_{H}(x)=\frac{-d_{1,H}+\sqrt{d_{1,H}^{2}+8d_{2,H}(x)\mathbb{E}|u_{H,0}|^{2}}}{2d_{2,H}(x)}\
,$ (2.16)

which elicits that the asymptotic behavior at large $x$ of $T_{H}(x)$ is a
decay inversely proportional to $x$. Since all terms in the above equation are
finite and positive, $T_{H}(x)$ is positive for all $x$. Furthermore, a
numerical evaluation of $T_{H}(0)$ confirms the expectation from dimensional
analysis that this quantity is approximately inversely proportional to the
integral length scale $L$.

These results show that the statistical properties in space at fixed time of
the fractional linear cascade model, Eq. (1.2), reproduce well the fractional
Gaussian field. On the other hand, statistical properties in time are
inhomogeneous, in contrast with the phenomenology of stochastic Navier-Stokes
turbulence, in which homogeneous statistics is observed both for time and
spatial increments. In turbulent velocity fields, the Eulerian structure
function in time (that is, increments in time measured at some fixed point in
space) show the same local regularity as their spatial equivalent, with a
power-law structure function approximately given by (for the second order)
$\mathbb{E}(\delta_{\tau}u)^{2}\propto\tau^{2/3}$ in the limit of vanishing
viscosity, for a small time interval $\tau$. This effect, which has been
measured in direct numerical simulations [30, 31], can be phenomenologically
explained with the sweeping effect [32], the random advection of the small
scales by the large scale motion of the flow, and reproduced with spatio-
temporal random fields [33, 34].

We briefly mention that the aforementioned properties hold as well for the
viscous fractional equation, Eq. (1.2), and can be derived from its analytical
solution [21]. With a finite viscosity, the solution develops a finite
inertial range where the behavior described by the spatial statistics of Eqs.
(2.5) and (2.9) is valid, but which is interrupted by the dissipative length
introduced by the viscosity, of the order of $(\nu/c)^{1/3}$. With the action
of dissipation, the system furthermore reaches its stationary state in a
finite time of the order of $1/(c^{4/3}\nu^{1/3})$. One other consequence of
the dissipative term is regularizing the solution at small spatial scales,
creating a dissipative range in which the second order structure function in
space is proportional to $\ell^{2}$, as expected for a smooth function. This
regularization in space however is not transmitted to the time dependence,
which remains with the same roughness as a Brownian motion.

###  2.2 Numerical results

Numerical simulations of Eq. (1.2) have been performed with a pseudo-spectral
method in space and a first order predictor-corrector scheme for the time
evolution, as described in Ref. [21], where the same numerical method has been
applied to the study of the equivalent nonlinear equation. An important
difference in the method is that a spatial truncation of the force is applied
in [21], through a smooth bump function which vanishes at the boundaries of
the domain, $\pm L_{\mathrm{tot}}/2$. This truncation is fundamental in the
nonlinear setting, since the simulation domain is finite and thus the operator
$\mathcal{L}$ is discontinuous at these boundaries. The interaction of the
nonlinearities with this discontinuous operator leads to a numerical
divergence which is avoided by the truncation. The simulations performed in
this work do not employ such truncation, since the linear equation is much
more well behaved numerically. An advantage of this approach, without the
truncation, is that the velocity field is homogeneous, as expected from the
theory, and entire time slices of it can be used for the evaluation of
statistical observables, easily generating large statistical ensembles. The
effects of the discontinuous operator are small, but still noticeable, and
will be discussed in the next paragraphs.

The chosen parameters for our simulations are the following: The spatial
dimension, with a total length of $L_{\mathrm{tot}}=1.0$ is composed of
$N=2^{12}$ points, with a spacing of $\Delta x=N/L_{\mathrm{tot}}$, which
leads to a maximum wavelength of $k_{\mathrm{max}}=N/2L_{\mathrm{tot}}$. For
the external force, the functional form $C_{f}(x)=\exp(-x^{2}/2L^{2})$ was
chosen, with three different values of the correlation length,
$L=(L_{\mathrm{tot}}/8,L_{\mathrm{tot}}/16,L_{\mathrm{tot}}/32)$. Three values
were chosen for the Hölder exponent as well, $H=(1/3,1/2,2/3)$, where $H=1/3$
is adopted for most figures because of its prominence in the statistics of
turbulence. $H=1/2$ represents a standard Brownian motion and $H=2/3$ is a
random field which is more regular than Brownian motion. In order to show the
action of the dissipation in regularizing small scales, the following values
of viscosity are chosen: $\nu=(2\times 10^{-7},3\times 10^{-8},2\times
10^{-9},2\times 10^{-10})$. We set the initial condition to $u_{H,\nu}(0,x)=0$
and the simulation first runs for a transient time large enough for the system
to reach its statistically stationary state. This transient time is estimated
as $T_{\mathrm{stat}}=2^{12}$, since the transport rate can be set to $c=1.0$
without loss of generality. Since the $\mathcal{L}$ operator induces a linear
transport in Fourier space, this transient time is linearly proportional to
the dissipative wavelength of the system and the chosen value is enough for it
to reach a statistically stationary state. After this transient regime, the
simulation runs for more $T_{\mathrm{tot}}=2^{12}$ units of time, and during
this stage we store the data for later analysis as several velocity profiles
in space at different instants in time and as profiles in time at different
positions. The time steps were taken as $\Delta t=\Delta x$, as in [21], a
value which is seen to produce well-resolved velocity fields, since the PDE we
investigate is linear with a very small viscosity. For larger values of the
viscosity, a choice of $\Delta t\propto(\Delta x)^{2}$ would have been
necessary, as expected from the stability analysis of the heat equation (the
dissipative term) [35].

Figure 1: Spatial profile of the real part of the numerical solution of Eq.
(1.2) at a fixed time $t$ in the statistically stationary regime for different
values of the parameter $H$ (see legend), for $N=2^{12}$,
$L=L_{\mathrm{tot}}/8$ and $c=1.0$. In (a), the viscosity is $\nu=2\times
10^{-7}$ and in (b) it is $\nu=2\times 10^{-10}$. The same random space-time
realization of the forcing was used for the three values of $H$ for each value
of viscosity.

In Fig. 1, the real part of the spatial profile of the velocity field
$u_{H,\nu}$ at a fixed time is shown for two different values of the
viscosity, $\nu=2\times 10^{-7}$ (a) and $\nu=2\times 10^{-10}$ (b), for all
three values of $H$, shown with different colors. The same random instance of
the external forcing was applied to all three velocity fields for each fixed
viscosity, which explains the fact that all three velocity fields at either
(a) or (b) look very similar, since Eq. (1.2) is linear. Nevertheless, the
curves display different roughness based on their values of the parameter $H$,
with the blue curve ($H=1/3$) being the one with highest fluctuations, while
the red one ($H=2/3$) is the smoothest. The velocity fields look homogeneous
over the whole system, since the forcing is not truncated, and a discontinuity
at the boundary can be noticed, due to the discontinuity in the $\mathcal{L}$
operator. Furthermore, the effect of the large difference in viscosity between
(a) and (b) is noticeable in the fact that Fig. 1(b) develops more roughness
and smaller scales than at (a). Yet, both velocity fields share the same large
scale correlation length $L=L_{\mathrm{tot}}/8$, noticeable in the large scale
behavior of all curves. The chosen value for $L$ as a few fractions of
$L_{\mathrm{tot}}$ allows for the easier visualization of the spatial
correlations in $u_{H,\nu}$, but this value must be small enough such that the
correlation function $C_{f}(x)$ is well resolved and sufficiently small at the
boundaries of the system, otherwise an extra numerical discontinuity would be
introduced. We have only shown the real part of the velocity fields since the
imaginary part looks similar.

Figure 2: Numerical estimations of the power spectrum (a) and second order
structure function in space (b) of the solution $u_{H,\nu}$ with the same
parameters as in Fig. 1 but only one value of the Hölder exponent, $H=1/3$.
Different colors represent different viscosities, following the legend. The
dashed black lines show the asymptotic analytical expressions, for the
inviscid equation, of the power spectrum (Eq. 2.4) in (a) and of the increment
variance (Eq. 2.9) in (b), together with their corresponding power-law
dependence indicated. The inset shows the negative part of the power spectrum
and its analytical result.

The second order statistical properties in space of the velocity field
$u_{H,\nu}(t,x)$ are depicted in Fig. 2, namely the power spectrum (a) and the
second order structure function in space (b). The estimation of the power
spectrum is obtained from the periodogram, as the square norm of the discrete
Fourier transform of the numerical solution divided by $N$, averaged over
several independent instants in time, obtained from a single long realization
of Eq. (1.2). In (a), the periodogram is shown for different values of the
viscosity, $\nu=(2\times 10^{-7},3\times 10^{-8},2\times 10^{-9},2\times
10^{-10})$, with different colors, along with the analytical expression for
the power spectrum of the inviscid velocity field (Eq. 2.4), in a black dashed
curve. It can be observed that the numerical curves exactly follow the
analytical result up to some maximum wavelength where the numerical power
spectrum rapidly decays. Furthermore, this maximum wavelength grows as the
viscosity is reduced, and for the smallest value of viscosity a large inertial
range with a power-law decay of $k^{-(2H+1)}$ can be observed. We mention
that, since the viscous fractional equation also possesses an explicit
analytical solution, we can also plot its analytical power spectrum, and
observe that it matches the numerical curve not only on the energy containing
and inertial range, but also on the dissipative range (data not shown), up to
a new maximum wavelength where boundary effects in the simulations can be
noticed. These boundary effects can be seen in Fig. 2(a) as a change of slope
in the dissipative range. This is a consequence of the discontinuity of the
operator $\mathcal{L}$ due to the finite size of the system in the numerical
simulations, but we remark that it only differs from the analytical results in
the far dissipative range. The inset of this figure shows the negative part of
the spectrum, which is not symmetric in $k$, a consequence of the complex
nature of the velocity field $u_{H,\nu}$. The negative part of the spectrum is
dominated solely by the forcing correlation function and decays fast, as can
be seen in the numerical or in the analytical result. Nevertheless, the effect
of the discontinuity can be more clearly noticed as the accumulation of energy
at negative frequencies, as seen in the inset. This deviation from the
analytical results leads to a numerical variance which is slightly bigger than
the analytical result, since the variance is an integral of the power spectrum
(Eq. 2.7), this difference is however numerically quite small. Additionally,
the inverse cascade due to the discontinuity of $\mathcal{L}$ has been
analogously observed in Ref. [23], where this operator is replaced by a
periodic version, as $2i\pi cx\to iL_{\mathrm{tot}}\sin(2\pi
cx/L_{\mathrm{tot}})$, but the regions of negative slope in the sine produce a
strong accumulation of energy at the boundaries. For these reasons, the
original operator is directly employed in our simulations without
periodization.

In Fig. 2(b) the second order structure function in space is shown, as a
function of the separation distance $\ell$, rescaled by the integral length
scale $L$, in log-log scale. To estimate this quantity, we average the square
norm of the increment,
$\delta_{\ell}u_{H,\nu}(t,x)=u_{H,\nu}(t,x+\ell)-u_{H,\nu}(t,x)$, over
independent instants in time and in space as well. This velocity field is
homogeneous, hence a large number of samples can be obtained. An inertial
range with an $\ell^{2H}$ power law behavior, which agrees with the analytical
predictions at vanishing viscosity (black dashed line) can be noticed,
although it is smaller than the inertial range for the power spectrum. At
small separation, the dissipative range manifests as an $\ell^{2}$ behavior in
the structure function, signaling a smooth velocity field. It can be observed
as well that the transition from the inertial range to the dissipative range
happens at smaller scales for the smallest values of viscosity. At large
separations, all curves converge to twice the value of the variance of the
velocity field, where the small difference in the variance of each curve is
due to dissipation.

Figure 3: Profile in time of the real part of the numerical solution at its
stationary state at two different positions, with the same parameters as in
the previous figures and $H=1/3$. In (a), the evolution in time is shown at
the origin of the system, and in (b), at $x=0.1L_{\mathrm{tot}}$. Each plot
shows two values of viscosity, a large one and a small one. The plots have
been rescaled in terms of the position-dependent integral correlation time of
the velocity field, as given by the estimation in Eq. (2.16), and show
equivalent intervals after this rescaling. The inset of (b) shows the profile
in time at $x=0.1L_{\mathrm{tot}}$ but for the same total time as in (a), in
the units of the simulation.

We now turn to a discussion of the statistical properties in time of
$u_{H,\nu}$ when driven by a delta-correlated force. In Fig. 3, the profile in
time of the real part of the velocity field is shown, at two different
positions, $x=0$ (a) and $x=L_{\mathrm{tot}}/10$ (b), with two different
values of viscosity in each plot, as indicated in the legend. The regularity
in time of this random process is that of a Brownian motion, independently of
the spatial Hölder parameter, as shown by Eq. (2.14). For this reason, only
$H=1/3$ is shown, with a similar behavior observed in other cases. We also
observe that viscous dissipation does not affect the small time scales, since
both curves at each plot are roughly indistinguishable in a statistical sense,
showing the same correlation time, amplitude of fluctuations and level of
roughness. The horizontal axes in (a) and (b) have been rescaled in terms of
the local integral correlation time, according to Eq. (2.16), and it can be
seen that the large scale patterns are similar in both figures. Since
$T_{H}(x)$ at (b) is a few times smaller than $T_{H}(0)$, we also show, in the
inset of (b), the profile in time of the numerical solution for the same total
time as in (a), in units of the simulation. In this inset, many multiples of
the integral time behavior can be noticed. This reinforces the relevance of a
position-dependent integral time scale, since the large-scale behavior in (a)
and the inset of (b) show strong qualitative differences.

Figure 4: Estimations of the second order structure function in time of
$u_{H,\nu}$ at its stationary state for different values of viscosity $\nu$
(see legend). In (a), the structure function is measured at the origin of the
system and in (b) at $x=0.1L_{\mathrm{tot}}$. Black dashed lines indicate
asymptotic analytical results, given by Eq. (2.14), and each plot is rescaled
by its respective position-dependent integral correlation time $T_{H}(x)$. In
the inset, the second order structure function in time, at the origin, for
$\nu=2\times 10^{-10}$ and different forcing correlation lengths $L$ is shown,
with the following values: $L_{\mathrm{tot}}/8$ (blue), $L_{\mathrm{tot}}/16$
(green) and $L_{\mathrm{tot}}/32)$ (red). The inset has not been rescaled by
$T_{H}(0)$, since it depends on $L$, but by $T_{\mathrm{stat}}$, which depends
only on the viscosity.

The nonhomogeneities of the velocity profile in time are inherited by the
temporal second order structure function, which is shown in Fig. 4 for two
different positions, the same as in the previous figure. In Fig. 4(a), the
analytical prediction at small time separations given by Eq. (2.14) is plotted
in a black dashed curve, with the coefficients $d_{1,H}$ and $d_{2,H}$
evaluated numerically. It can be seen that, at small times, the analytical
curve matches all numerical observations with different values of viscosity,
confirming the claim that the local regularity in time is determined only by
the correlations in time of the force and not by viscosity. The analytical
curve, furthermore, has a slope of one for small times, corresponding to the
local regularity of a Brownian motion. The concavity of the black curve,
nonetheless, is positive, since it was obtained from a quadratic
approximation, and to capture the negative concavity of the numerical curves a
higher order expansion in $\tau$ would be required. The inset of (a) shows the
same observable for the smallest viscosity and three different values of $L$,
along with the analytical prediction which does not depend on $\nu$ or $L$.
The correlation length influences the value of the variance, leading thus to
different correlation times, as given by Eq. (2.16). Moreover, the reasoning
behind the approximation taken to obtain Eq. (2.16) is visible in the inset,
where we observe a rapid transition from the linear behavior to the large
separation constant behavior. In Fig. 4(b), the same analytical prediction is
seen (Eq. 2.14), calculated at $x=0.1L_{\mathrm{tot}}$. The numerical curves
are seen to agree on the positive concavity of the analytical approximation,
for small values of $\tau$. Just as in Fig. 3, a different rescaling is used
at (a) and (b), corresponding to each local integral time scale, calculated
numerically. In the inset, the chosen rescaling is $T_{\mathrm{stat}}$, which
does not depend on $L$, but only on the viscosity, unlike the correlation time
$T_{H}(0)$ which depends on $L$ but not on the viscosity.

##  3 Correlated in time forcing

The natural limit to contrast the results of a delta-correlated forcing would
be that of a frozen force, a random force correlated in space on a length $L$
but constant in time. Nevertheless, it can be shown that the variance of
$u_{H,0}$ at the origin, that is,
$\lim_{t\to\infty}\mathbb{E}|u_{H,0}(t,x=0)|^{2}$, diverges under a frozen
force if the Hölder parameter is smaller than $1/2$. This means the dynamical
evolution of Eq. (1.2) is not able to regularize a constant forcing scheme,
producing a state of infinite variance at the origin which is unphysical. As
we have said, the range $H<1/2$ is of great interest due to its connection
with the phenomenology of turbulence, and for this reason we concentrate in
this section on a forcing scheme with finite correlation time $T$. We choose
the correlation function in time as

$\chi(t)=\frac{1}{2T}e^{-|t|/T}\ ,$ (3.1)

such that the results of the previous section can be reobtained in the limit
of $T\to 0$, in which the above correlation function is reduced to a Dirac
delta. In the next section, we revisit analytical and numerical results of Eq.
(1.2), now driven by an exponentially correlated in time forcing. As will be
discussed, the statistical homogeneity of the spatial observables is broken, a
strong qualitative difference in comparison to the delta-correlated forcing
case. We discuss the time statistics as well, which were already not
homogeneous under a white-in-time forcing.

###  3.1 Numerical results

Under a correlated forcing in time, a numerical scheme for deterministic
differential equations must be employed to obtain $u_{H,\nu}$, since the
external force is now continuous in time and space. For this reason, we
generate the forcing with a numerical scheme for stochastic differential
equations, and for the time evolution of the system, given the forcing, we use
a fourth-order Runge-Kutta scheme. The stochastic scheme for the forcing is
the same as in the previous section, a first-order predictor-corrector
algorithm [36]. A forcing with the specific correlation profile of Eq. (1.6)
is provided by an Ornstein-Uhlenbeck process, defined by the solution of the
following dynamics:

$df(t,x)=-\frac{1}{T}f(t,x)dt+\frac{1}{T}d\eta(t,x)\ .$ (3.2)

To completely specify the random force $d\eta$, we furthermore say it is of
zero mean, smooth in space and delta-correlated in time, with independent real
and imaginary parts. These conditions can be stated precisely in terms of the
integration against a smooth test function $h(t)$ as

$\mathbb{E}\int h(t)d\eta(t,x)=0,$ (3.3) $\mathbb{E}\int
h(t_{1})h(t_{2})d\eta(t_{1},x)d\eta(t_{2},y)=0\ \mathrm{and}$ (3.4)
$\mathbb{E}\int
h(t_{1})h(t_{2})d\eta(t_{1},x)d\eta^{*}(t_{2},y)=C_{f}(x-y)\int h^{2}(t)dt,$
(3.5)

where $C_{f}(x)$ is chosen the same as in Eq. (1.6). The drift and diffusion
terms entering in Eq. 3.2 have been chosen such that the correlation function
of the solution, $f$, in its statistically stationary state, is

$\mathbb{E}[f(t+\tau,x)f^{*}(t,y)]=\frac{1}{2T}e^{-|\tau|/T}C_{f}(x-y)\ ,$
(3.6)

as stated before. The correlations in space of the force are obtained from a
circulant embedding algorithm, that is, the convolution of independent complex
variables (of mean zero and variance $\sqrt{dx/2}$ for their real and
imaginary parts) with a Gaussian kernel, which produces a random series with
Gaussian autocorrelation function. To begin the simulations, the force must be
at its stationary state, a process that takes a time of the order of $T$.
After this first transient phase, the force is then used in the numerical
solution of Eq. (1.2) in the same manner a deterministic external contribution
would be. Despite the different algorithm, at this stage the simulation occurs
as described in Sec. 2.2: There is a second transient phase determined by the
linear transport in Fourier space, which takes a time of the order of
$T_{\mathrm{stat}}\propto\nu^{-1/3}$, and after that a stationary phase, from
which we collect random series in space and in time. The parameters chosen for
these numerical simulations are similar to the ones already used: $N=2^{12}$,
$L_{\mathrm{tot}}=1.0$, $L=L{\mathrm{tot}}/8$, $c=1.0$, $\nu=2\times 10^{-10}$
and $\Delta t=\Delta x$. The same three values of the Hölder parameter were
used, $H=1/3$, $H=1/2$ and $H=2/3$, and for the forcing correlation time the
following values were chosen: $T=(0.1,0.2,0.5,2.0,5.0,10.0)$.

Figure 5: Spatial profile of the real part of the numerical solution of the
fractional equation (Eq. 1.2) in the statistically stationary state, under an
exponentially correlated in time forcing (Eq. 3.2), for different correlation
times. The parameters in the figures are $N=2^{12}$, $L=L_{\mathrm{tot}}/8$,
$c=1.0$ and $\nu=2\times 10^{-10}$. Each figure shows three different values
of the $H$ parameter for the same space-time realization of the force, these
values are $H=1/3$ (blue), $H=1/2$ (green) and $H=2/3$ (red). Each figure
corresponds to a different forcing correlation time, respectively, $T=0.5$
(a), $T=2.0$ (b) and $T=5.0$ (c).

In Fig. 5, the spatial profile of some of these numerical simulations is shown
at a fixed time, in the stationary state. Analogously to Fig. 1, the same
space-time realization of the forcing was applied to three different values of
the Hölder parameter, and a different roughness can be observed for each
curve. On the other hand, all of the profiles shown display the same value of
viscosity, $\nu=2\times 10^{-10}$, with each figure corresponding to a
different forcing correlation time, where in (a) it is the smallest, equal to
$T=0.5$, in (b) it is $T=2.0$ and in (c), $T=5.0$. The velocity fields in (a)
look homogeneous and similar to those shown in Fig. 1(b), while in Fig. 5(b),
as the correlation time grows, a tendency of accumulation at the center of the
domain can be noticed, an effect that is confirmed in (c), where the
nonhomogeneity of the system is clear and the correlation time is the largest.

To understand these observations, let us consider the spatio-temporal
correlation of the solution $u_{H,0}(t,x)$. After some calculations, we obtain

$\mathbb{E}\left[u_{H,0}(t+\tau,x+\ell)u^{*}_{H,0}(t,x)\right]=\frac{1}{2T}\int_{s_{1}=0}^{t+\tau}\int_{s_{2}=0}^{t}e^{2i\pi
c[s_{1}\ell+x(s_{1}-s_{2})]}\mathcal{D}_{\ell}(s_{1},s_{2})e^{-\frac{|\tau+s_{2}-s_{1}|}{T}}ds_{1}ds_{2},$
(3.7)

where enters the symmetric quantity

$\begin{split}\mathcal{D}_{\ell}(s_{1},s_{2})&=\int_{k\in\mathbb{R}}e^{2i\pi
k\ell}|k+cs_{1}|^{-H-1/2}_{1/L}|k+cs_{2}|^{-H-1/2}_{1/L}|k|^{2H+1}_{1/L}\widehat{C}_{f}(k)dk\\\
&=\mathcal{D}_{\ell}(s_{2},s_{1}).\end{split}$ (3.8)

Concerning the variance at large times, consider the expression provided in
Eq. (3.7) for $\ell=0$ and $\tau=0$. Using the symmetry of the function
$\mathcal{D}_{\ell}$, we get

$\lim_{t\to\infty}\mathbb{E}\left|u_{H,0}(t,x)\right|^{2}=\frac{1}{2T}\int_{(s_{1},s_{2})\in(\mathbb{R}^{+})^{2}}\cos(2\pi
cx(s_{1}-s_{2}))\mathcal{D}_{0}(s_{1},s_{2})e^{-\frac{|s_{2}-s_{1}|}{T}}ds_{1}ds_{2},$
(3.9)

which depends explicitly on the position $x$ and is bounded by its value at
the origin $x=0$. First, remark again that, for a frozen force, obtained by
replacing $e^{-|s_{1}-s_{2}|/T}/2T\rightarrow 1$, this expression diverges at
$x=0$ and $H<1/2$. Going back to the correlated in time forcing (Eq. 3.9), we
can see its finiteness for $H>0$ while noticing that, along the diagonal
$s_{1}=s_{2}$, using the expression of $\mathcal{D}_{\ell}$ (Eq. 3.8) and the
rapid decay of $\widehat{C}_{f}(k)$, all integrals converge at large
arguments. This explains the observations made in Fig. 5.

Figure 6: Estimation of the spatial statistics of the fractional equation for
different correlation times of the forcing, in comparison to a solution with
delta-correlated forcing (in purple, indicated by $T=0$ in the legend). The
parameters $H=1/3$ and $L=L_{\mathrm{tot}}/8$ are used and all correlation
times are indicated in the legend. In (a), the power spectrum shown is
estimated from the periodogram, as $\mathbb{E}|\widehat{u}_{H,\nu}|^{2}/N$,
and the black dashed line corresponds to the analytical expression of the
power spectrum in the inviscid and white-in-time case (Eq. 2.4). The inset
shows the corresponding negative frequencies of the power spectrum, with the
same analytical expression in a black dashed line. In (b), the second order
structure function, calculated with respect to the origin (see text) is shown,
along with the analytical prediction for the inviscid delta-correlated
equation (Eq. 2.9).

The spatial statistics of this velocity field is shown in Fig. 6, for
different correlation times of the forcing. Similarly to Fig. 2(a), we begin
displaying in Fig. 6(a) the periodogram. For a direct comparison, we reproduce
the white-in-time forcing case, indicated by $T=0.0$ and plotted with a dark
purple curve. For small values of the correlation time, $T=0.1$ and $T=0.2$,
the curves can be hardly distinguished from the white-in-time case, both of
them also matching the analytical prediction of Eq. (2.4), in a black dashed
line. As the correlation time grows, the general shape of the curve is
preserved, with an inertial range displaying power-law decay with an exponent
of $-(2H+1)$, but at significantly smaller values, due to the $1/T$ prefactor
in the time correlation of the force. In the inset, the power spectrum at
negative frequencies is shown for the same parameters, and the agreement with
the original simulation is observed for small correlation times, while smaller
values with the same shape are observed for larger correlation times.

As we have seen, in particular while considering the variance in presence of a
forcing term which is correlated in time (Eq. 3.9), the solution is not
statistically homogeneous. This is why we calculate the second order structure
function in space of the nonhomogeneous velocity fields with respect to the
origin. That is, all velocity increments at scale $\ell$ are defined as

$\delta_{\ell}u_{H,\nu}=u_{H,\nu}(t,\ell)-u_{H,\nu}(t,0)\ ,$ (3.10)

and averaged only over independent instants in time. Our results are displayed
in Fig. 6(b). Being calculated only with respect to the origin, the statistics
available is much smaller and the numerical curves are not as well resolved.
Nevertheless, a tendency to follow the power-law behavior with a $2H$ exponent
is still observed, and even the variance of the system at the origin, as
measured by the large separation results, has a similar value to that of the
$T=0$ curve. For larger correlation times, this agreement is lost, since the
horizontal asymptote is
$\mathbb{E}|u_{H,\nu}(t,\ell)|^{2}+\mathbb{E}|u_{H,\nu}(t,0)|^{2}$, which
differs significantly from $2\mathbb{E}|u_{H,\nu}(t,0)|^{2}$, as already shown
in Eq. (3.9). An analytical expression for the increment variance can be
obtained from Eq. (3.7):

$\begin{split}\mathbb{E}|\delta_{\ell}u_{H,0}|^{2}&=\frac{1}{T}\int_{(s_{1},s_{2})\in(\mathbb{R}^{+})^{2}}\left(\mathcal{D}_{0}(s_{1},s_{2})-\mathcal{R}e^{2i\pi
cs_{1}\ell}\mathcal{D}_{\ell}(s_{1},s_{2})\right)e^{-|s_{1}-s_{2}|/T}ds_{1}ds_{2}\
.\end{split}$ (3.11)

This expression is positive and real, as expected, and in the limit $T\to 0$,
the delta-correlated result, Eq. (2.9), is recovered.

Figure 7: Profile in time of the velocity field $u_{H,\nu}$ at the origin (a)
and at $x=0.1L_{\mathrm{tot}}$ (b), each rescaled by the local integral
correlation time calculated for the white-in-time forcing case in Eq. (2.16).
Each figure shows two values of the correlation time: $T=0.1$ (blue) and
$T=10.0$ (red). The inset of (b) shows the same velocity profile as in (b),
for a total time equal to the one of (a), in units of the simulation.

The statistical properties in time of Eq. (1.2) under correlated forcing in
time are investigated in the next two figures. For these observables, we
remind that their statistics under delta-correlated force was already
nonhomogeneous, as shown in Eq. (2.14) and discussed in Figs. 3 and 4. An
important difference, though, is that the solution $u_{H,\nu}$ under
correlated forcing is continuous and differentiable in time as well as in
space, as can be seen in Fig. 7, where a profile in time is shown at different
positions, the origin (a) and at $x=0.1L_{\mathrm{tot}}$ (b). In each figure,
two values of the forcing correlation time are shown, $T=0.1$ (blue) and
$T=10.0$ (red), and it can be noticed that, for the small correlation time,
the curve looks rough exactly as the one in Fig. 3(a), its regularization in
time happening only at very small scales, while the red curve, with a larger
correlation time, is noticeably smoother. Outside the origin, as shown in Fig.
7(b), the horizontal axis has been rescaled by the corresponding integral time
(Eq. 2.16), and an equivalent interval as in (a) is shown. In the inset, a
total time equal to the one in (a), in units of the simulation, is shown, and
many multiples of the integral time can be seen. An important difference with
the case of white-in-time forcing, though, is that the variance is almost
homogeneous when the correlation time is small, while it decays fast for
larger $T$, reason for which the red curve has a smaller amplitude than the
blue curve, as can be clearly noticed in the inset.

Figure 8: Numerical estimation of the second order structure function in time
of Eq. (1.2) under correlated forcing in time. The colors follow the same
pattern as in Fig. 6, where the purple curve represents a delta-correlated
simulation and yellow represents the largest correlation time ($T=10.0$). In
(a), this structure function is evaluated at the origin and in (b), at
$x=0.1L_{\mathrm{tot}}$. Each horizontal axis has been rescaled according to
the local integral correlation time, as given by the estimate of Eq. (2.16)
and the black dashed curves are the same as in Fig. 4.

The corresponding statistical properties in time are reported in Fig. 8, where
we plot the temporal second order structure function for different correlation
times of the forcing and at two different positions, the origin (a) and at
$x=0.1L_{\mathrm{tot}}$ (b). In (a), the black dashed curve, corresponding to
the small-time analytical result for the delta-correlated forcing (Eq. 2.14
with $x=0$) can be seen to match only the corresponding delta-correlated
numerical simulation (dark purple) at very small time intervals. All other
curves, with a finite correlation time, show a dissipative range in time,
where this structure function is proportional to $\tau^{2}$ instead of $\tau$.
This can be understood by looking at the derivatives of the correlation
function in time in the stationary state. Consider the Taylor expansion of Eq.
(3.7) around $\tau$, with $\ell=0$. The first order contribution is
proportional to

$\begin{split}&\frac{d\mathbb{E}[u_{H,0}(t+\tau,x)u_{H,0}^{*}(t,x)]}{d\tau}=\frac{1}{2T}\int_{s_{2}=0}^{t}e^{2\pi
icx(t+\tau-s_{2})}\mathcal{D}_{0}(t+\tau,s_{2})e^{-|t-s_{2}|/T}ds_{2}\\\
&+\frac{1}{2T^{2}}\int_{s_{1}=0}^{t+\tau}\int_{s_{2}=0}^{t}e^{2\pi
icx(s_{1}-s_{2})}\mathcal{D}_{0}(s_{1},s_{2})e^{-|\tau+s_{2}-s_{1}|/T}\frac{\tau+s_{2}-s_{1}}{|\tau+s_{2}-s_{1}|}ds_{1}ds_{2}\
.\end{split}$ (3.12)

The first term in this equation vanishes in the limit $t\to\infty$ and the
second term is identically zero at $\tau=0$ because it is odd under the
exchange of $s_{1}$ and $s_{2}$. Therefore, the contribution proportional to
$\tau$ vanishes. Deriving Eq. (LABEL:eq:corrf-time-taylor-first) with respect
to $\tau$, we obtain a single nonvanishing contribution in the stationary
limit, equal to

$\begin{split}\lim_{t\to\infty}\frac{d^{2}\mathbb{E}[u_{H,0}(t+\tau,x)u_{H,0}^{*}(t,x)]}{d\tau^{2}}\Bigg{\lvert}_{\tau=0}&=-\frac{1}{2T^{3}}\int_{(s_{1},s_{2})\in(\mathbb{R}^{+})^{2}}\cos(2\pi
cx(s_{1}-s_{2}))\\\
&\times\mathcal{D}_{0}(s_{1},s_{2})e^{-|s_{1}-s_{2}|/T}ds_{1}ds_{2}\
,\end{split}$ (3.13)

which is proportional to the variance, and is thus finite for any $H>0$.
Moreover, this contribution is negative, since the autocorrelation function at
$\tau$ must necessarily be smaller than the variance. This confirms the
numerical observation in Fig. 8 that the introduction of a finite time
correlation in the forcing induces a smooth behavior of the velocity field at
small time intervals.

This new dissipative range depends only on the correlation time and not on
viscosity, but it can still be seen that, for small values of $T$, there is an
inertial range where the linear behavior is recovered, in agreement with Eq.
(2.14) at $x=0$. Furthermore, we observe that the variance of all of these
solutions is approximately the same, as can be seen in the large time
separation end of the graph. Hence, at the origin the velocity field
$u_{H,\nu}$ retains a variance close to that of the white-in-time scenario, as
already observed in Fig. 6 for small correlation times.

In Fig. 8(b), the second order structure function in time is shown at the
position $x=0.1L_{\mathrm{tot}}$. Again, it can be noticed that, for small
values of the correlation time, the darker curves closely follow the behavior
of the delta-correlated solution, but differ only at very small time
increments, where the finitely correlated solutions show smooth behavior,
proportional to $\tau^{2}$, while the white-in-time forcing solution resembles
a Brownian motion in time, as given by Eq. (2.14) at very small values of
$\tau$. For larger correlation times, the same general behavior is observed,
with a dissipative range at small time separations, oscillations with the same
pattern as those of the delta-correlated forcing, but significantly smaller
values of the variance, as given by the asymptotic values of these curves at
large time. This is again a consequence of the loss of homogeneity of the
solution $u_{H,\nu}$ as the correlation time of the force is increased.

##  4 Concluding remarks

We have investigated the statistical properties in space and in time of a
linear dynamical model of the turbulent energy cascade, introduced in Ref.
[21]. This is a linear stochastic partial differential equation with a smooth
large-scale forcing in space and general correlation in time. The spatial
statistics under a white-in-time forcing were first investigated in [21],
where it is observed that the solution, at fixed times, is homogeneous and
resembles a fractional Gaussian field with Hölder parameter $H$, a stationary
state which is reached even in the absence of viscosity. Hence, with the
specific choice of $H=1/3$, the spatial profiles of the solution offer a
representation of the K41 framework of turbulence as the result of a
stochastic PDE. Nevertheless, the statistics in time under delta-correlated
forcing, first investigated in the present work (Sec. 2), are not homogeneous,
and are not influenced by the regularity in space, displaying the same local
regularity in time as a Brownian process. At the origin of the system, the
correlation time of the solution is larger, but it rapidly decays as one moves
away from this point. One can contrast this picture with turbulent velocity
fields in the Eulerian setting, in which the same local roughness in time and
in space is observed [30, 31].

The effects of a finitely correlated in time forcing were investigated in Sec.
3. For very small correlation times $T$, many of the previous results for a
delta-correlated in time force are reproduced, but as $T$ is increased, the
solution becomes clearly nonhomogeneous, as can be seen in its spatial
profiles (Fig. 5). The local roughness of $H$ is preserved for small spatial
scales, but these nonhomogeneous velocity fields are no longer a good portrait
of turbulent flows. Since the forcing is now smooth in space and time, the
solution is smooth in time as well, as shown by its second order structure
functions, which display a dissipative range (Figs. 6 and 8). This dissipative
regime in time is again not influenced by viscous dissipation or by the Hölder
exponent of the solution.

This simple cascade model allows for the construction of fractional Gaussian
fields in a dynamical setting, starting from a smooth force in space. In [21],
a nonlinearity which induces non-Gaussian statistics on this velocity field
was investigated, but only for a delta-correlated in time forcing. Therefore,
the combined effects of nonlinearity and correlations in time of the forcing
in this model remain to be understood. In particular, we notice the contrast
between the model and hydrodynamic turbulence, in which correlations in time
of the stirring do not crucially change statistical properties of the
solution.

## Acknowledgments

We thank Jules Guioth and Malo Tarpin for enriching discussions. The authors
are partially supported by the Simons Foundation Award ID 651475.

## Conflict of interest

The authors declare no conflict of interest.

## References

  * [1] U. Frisch.  Turbulence: the legacy of A. N. Kolmogorov.  Cambridge university press, 1995.
  * [2] A. N. Kolmogorov.  The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers.  Cr. Acad. Sci. URSS, 30:301–305, 1941.
  * [3] L. Onsager.  Statistical hydrodynamics.  Il Nuovo Cimento, 6(2):279–287, 1949.
  * [4] P. Constantin, W. E, and E. S. Titi.  Onsager’s conjecture on the energy conservation for solutions of Euler’s equation.  Commun. Math. Phys., 165(1):207–209, 1994.
  * [5] J. Duchon and R. Robert.  Inertial energy dissipation for weak solutions of incompressible Euler and Navier-Stokes equations.  Nonlinearity, 13(1):249, 2000.
  * [6] G. Eyink and K. R. Sreenivasan.  Onsager and the theory of hydrodynamic turbulence.  Rev. Mod. Phys., 78(1):87, 2006.
  * [7] K. P. Iyer, K. R. Sreenivasan, and P. K. Yeung.  Reynolds number scaling of velocity increments in isotropic turbulence.  Phys. Rev. E, 95:021101, 2017.
  * [8] P. Debue, D. Kuzzay, E.-W. Saw, F. Daviaud, B. Dubrulle, L. Canet, V. Rossetto, and N. Wschebor.  Experimental test of the crossover between the inertial and the dissipative range in a turbulent swirling flow.  Phys. Rev. Fluids, 3:024602, 2018.
  * [9] B. Dubrulle.  Beyond Kolmogorov cascades.  J. Fluid Mech., 867, 2019.
  * [10] B. B. Mandelbrot and J. W. Van Ness.  Fractional brownian motions, fractional noises and applications.  SIAM review, 10:422–437, 1968.
  * [11] R. Robert and V. Vargas.  Hydrodynamic turbulence and intermittent random fields.  Commun. Math. Phys., 284:649–673, 2008.
  * [12] A. Lodhia, S. Sheffield, X. Sun, and S. S. Watson.  Fractional Gaussian fields: a survey.  Probab. Surv., 13:1–56, 2016.
  * [13] A. N. Kolmogorov.  Dissipation of energy in the locally isotropic turbulence.  Dokl. Akad. Nauk SSSR A, 32:16–18, 1941.
  * [14] F. Schmitt and D. Marsan.  Stochastic equations generating continuous multiplicative cascades.  Eur. Phys. J. B, 20:3–6, 2001.
  * [15] R. M. Pereira, C. Garban, and L. Chevillard.  A dissipative random velocity field for fully developed fluid turbulence.  J. Fluid Mech., 794:369–408, 2016.
  * [16] 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.  Ann. Henri Poincaré, 20:3693–3741, 2019.
  * [17] J. Friedrich, J. Peinke, A. Pumir, and R. Grauer.  Explicit construction of joint multipoint statistics in complex systems.  arXiv preprint arXiv:2105.03223, 2021.
  * [18] E. B. Gledzer.  System of hydrodynamic type admitting two quadratic integrals of motion.  Sov. Phys. Dokl., 18:216–217, 1973.
  * [19] K. Ohkitani and M. Yamada.  Temporal intermittency in the energy cascade process and local lyapunov analysis in fully-developed model turbulence.  Prog. Theor. Phys., 81(2):329–341, 1989.
  * [20] L. Biferale.  Shell models of energy cascade in turbulence.  Annu. Rev. Fluid Mech., 35:441–468, 2003.
  * [21] G. B. Apolinário, L. Chevillard, and J.-C. Mourrat.  Dynamical fractional and multifractal fields.  arXiv preprint arXiv:2107.03309, 2021.
  * [22] Y. Colin de Verdière and L. Saint-Raymond.  Attractors for two-dimensional waves with homogeneous Hamiltonians of degree 0.  Commun. Pur. Appl. Math., 73(2):421–462, 2020.
  * [23] S. Dyatlov and M. Zworski.  Microlocal analysis of forced waves.  Pure and Applied Analysis, 1:359–384, 2019.
  * [24] Y. Colin de Verdière.  Spectral theory of pseudodifferential operators of degree 0 and an application to forced linear waves.  Analysis & PDE, 13(5):1521–1537, 2020.
  * [25] L. R. M. Maas, D. Benielli, J. Sommeria, and F.-P. A. Lam.  Observation of an internal wave attractor in a confined, stably stratified fluid.  Nature, 388:557–561, 1997.
  * [26] M. Rieutord and L. Valdettaro.  Inertial waves in a rotating spherical shell.  J. Fluid Mech., 341:77–99, 1997.
  * [27] H. Scolan, E. Ermanyuk, and T. Dauxois.  Nonlinear fate of internal wave attractors.  Phys. Rev. Lett., 110:234501, 2013.
  * [28] C. Brouzet, E. V. Ermanyuk, S. Joubaud, I. Sibgatullin, and T. Dauxois.  Energy cascade in internal-wave attractors.  Europhys. Lett., 113:44001, 2016.
  * [29] J. C. Mattingly, T. Suidan, and E. Vanden-Eijnden.  Simple systems with anomalous dissipation and energy cascade.  Commun. Math. Phys., 276:189–220, 2007.
  * [30] L. Chevillard, S. G. Roux, E. Lévêque, N. Mordant, J.-F. Pinton, and A. Arnéodo.  Intermittency of velocity time increments in turbulence.  Phys. Rev. Lett., 95:064501, 2005.
  * [31] A. Gorbunova, G. Balarac, L. Canet, G. Eyink, and V. Rossetto.  Spatio-temporal correlations in three-dimensional homogeneous and isotropic turbulence.  Phys. Fluids, 33:045114, 2021.
  * [32] H. Tennekes and J. L. Lumley.  A First Course in Turbulence.  MIT Press, 1972.
  * [33] M. Chaves, K. Gawedzki, P. Horvai, A. Kupiainen, and M. Vergassola.  Lagrangian dispersion in Gaussian self-similar velocity ensembles.  J. Stat. Phys., 113:643–692, 2003.
  * [34] J. Reneuve and L. Chevillard.  Flow of spatiotemporal turbulentlike random fields.  Phys. Rev. Lett., 125:014502, 2020.
  * [35] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang.  Spectral methods: fundamentals in single domains.  Springer Science & Business Media, 2007.
  * [36] P. E. Kloeden and E. Platen.  Numerical solution of stochastic differential equations, volume 23.  Springer Science & Business Media, 2013.