# Multiscale turbulence synthesis: Validation in 2D hydrodynamics

P. Lesaffre, J.-B. Durrive, J. Goossaert, S. Poirier, S. Colombi, P. Richard, E.

Allys, W. Béthune

► To cite this version:

P. Lesaffre, J.-B. Durrive, J. Goossaert, S. Poirier, S. Colombi, et al. Multiscale turbulence synthesis: Validation in 2D hydrodynamics. Astronomy & Astrophysics - A&A, 2025, 703, pp.A291. 10.1051/0004-6361/202453264. hal-05410977

HAL Id: hal-05410977

https://hal.science/hal-05410977v1

Submitted on 11 Dec 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Multiscale Turbulence Synthesis: validation in 2D hydrodynamics

P. Lesaffre<fnref n="1" />,<fnref n="2" />, J.-B. Durrive<fnref n="3" />, J. Goossaert<fnref n="1" />, S. Poirier<fnref n="1" />, S. Colombi<fnref n="4" />, P. Richard<fnref n="1" />,<fnref n="2" />, E. Allys<fnref n="1" />,<fnref n="2" />, W. Béthune<fnref n="5" />

?? / ??

## ABSTRACT

**Context.** Numerical simulations can follow the evolution of fluid motions through the intricacies of developed turbulence. However, they are rather costly to run, especially in 3D. In the past two decades, generative models have emerged which produce synthetic random flows at a computational cost equivalent to no more than a few time-steps of a simulation. These simplified models qualitatively bear some characteristics of turbulent flows in specific contexts (incompressible 3D hydrodynamics or magnetohydrodynamics), but generally struggle with the synthesis of coherent structures.

**Aims.** We aim at generating random fields (e.g. velocity, density, magnetic fields, etc.) with realistic physical properties for a large variety of governing partial differential equations and at a small cost relative to time-resolved simulations.

**Methods.** We propose a set of simple steps applied to given sets of partial differential equations: filtering from large to small scales, first order time evolution approximation from Gaussian random initial conditions during a prescribed coherence time, with a final summation over scales to generate the whole fields.

**Results.** We test the validity of our method in the simplest framework: 2D decaying incompressible hydrodynamical turbulence. We compare results of 2D decaying simulations with snapshots of our synthetic turbulence. We assess quantitatively the difference first with standard statistical tools: power spectra, increments and structure functions. These indicators can be reproduced by our method during up to about a third of the turnover time scale. We also consider recently developed scattering transforms statistics, able to efficiently characterise non-Gaussian structures. This reveals more significant discrepancy, which can however be bridged by bootstrapping. Finally, the number of Fourier transforms necessary for one synthesis scales logarithmically in the resolution, compared to linearly for time-resolved simulations.

**Conclusions.** We have designed a multiscale turbulence synthesis (MuScaTS) method to efficiently short-circuit costly numerical simulations to produce realistic instantaneous fields.

**Key words.** hydrodynamics, magnetohydrodynamics (MHD), turbulence, ISM:structure

## 1. Introduction and context

Turbulence has distinctive features, which are hard to characterise quantitatively. Nevertheless, our subjective experience tells us that a lot of information is encoded even in a single 2D projection of a tracer in a 3D turbulent flow. Indeed, when we give a vigorous push to a puff of smoke, it soon develops convoluted motions, diverging trajectories blur the initial impulse, but nevertheless some characteristic signatures remain. When we look at a still picture of smoke trailing in the wake of a candle or an incense stick, our brain not only immediately recognises the context, it also grasps a sense for the fluid shearing and spiralling motions in vortices and gets a feeling of the large scale wind direction. Astronomers who study the interstellar medium are very much in the same situation: the sky before their eyes is a still projection of tracers evolving through the fluid

dynamics at play in these dilute media, which they attempt to infer.

Quantitative probes of the nature of the underlying fluid dynamics are however quite elusive. Turbulence research so far has managed to isolate only a few generic statistical properties shared by turbulent flows. In the case of incompressible hydrodynamics (HD), power-law scalings for the velocity power spectrum were first obtained by Kolmogorov (1941) using conservation arguments across scales. But scaling laws for the velocity and magnetic fields power spectra are still a matter of debate even in the relatively simple case of incompressible magnetohydrodynamic (MHD) turbulence (Iroshnikov 1963; Kraichnan 1965; Goldreich & Sridhar 1995; Schekochihin 2022). Then, anomalous scalings for the structure functions were discovered (Kolmogorov 1962), which led to a large legacy of intermittency models and observations of statistical properties on the increments of velocity, energy dissipation rates,

<footnote n="1">¹ Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005, Paris, France</footnote>

<footnote n="2">² Observatoire de Paris, PSL University, Sorbonne Université, LERMA, 75014, Paris, France</footnote>

<footnote n="3">³ Institut de recherche en astrophysique et planétologie Université Toulouse III - Paul Sabatier, Observatoire Mid-Pyrénées, Centre National de la Recherche Scientifique, UMR5277 France</footnote>

<footnote n="4">⁴ Sorbonne Université, CNRS, UMR7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, 75014 Paris, France</footnote>

<footnote n="5">⁵ DPHY, ONERA, Université Paris-Saclay, 91120, Palaiseau, France</footnote>
vorticity or magnetic fields (Frisch 1996; She & Leveque 1994; Grauer et al. 1994; Politano & Pouquet 1995). Finally, few exact analytical results were predicted such as the Kármán-Howarth-Monin 4/5th relation (de Karman & Howarth 1938; Monin 1959; Frisch 1996) in 3D incompressible HD and its analogues in compressible or MHD turbulence (Galtier & Banerjee 2011; Banerjee & Galtier 2013, 2014, 2016).

In contrast, time-resolved simulations can accurately predict the evolution of fluid motions from random initial conditions and we could try to assess what statistical laws control their outcomes. But such simulations are rather costly to run, especially in 3D and in a developed turbulent regime, which limits the number of possible independent realisations for these experiments.

In the past two decades, several authors have developed efficient techniques to produce random fields that bear some of the known characteristics of turbulence. A natural starting step is simply to produce a Gaussian velocity field with a prescribed power spectrum (Mandelbrot & van Ness 1968). Application of a non-linear transformation to this fractional Gaussian field can then yield some degree of non-Gaussianity, such as intermittency. However, these non-Gaussianities may have nothing to do with actual physics. For instance, it is challenging to generate asymmetries in the longitudinal velocity increments which generate non-zero energy transfer, let alone transfer in the right scale direction.

Some techniques focus on advection. Rosales & Meneveau (2006) (multiscale minimal Lagrangian map, MMLM) and Rosales & Meneveau (2008) (multiscale turnover minimal Lagrangian map, MTML) transport ballistically an initial fractional Gaussian velocity field under its own influence in a hierarchy of embedded smoothing scales. They show this procedure generates anomalous scalings and energy transfer although they do not compare it quantitatively to experiments. Subedi et al. (2014) later extended the same approach in MHD. More recently Lübke et al. (2024) showed they could initiate the formation of some current sheets if they applied the same technique to an already non-Gaussian initial field.

Other techniques focus on deformation, neglecting advection. Chevillard & Meneveau (2006) and Chevillard et al. (2011) deform a Gaussian field with a stretching matrix exponential to produce a tunable random field with given power spectrum and intermittency properties. Durrive et al. (2020) later generalised the same ideas to incompressible MHD. However, these models fail at reproducing coherent structures seen in direct numerical simulations (DNS, see Momferratos et al. 2014; Richard et al. 2022), so Durrive et al. (2022) extended their model to include prescribed random dissipative sheets.

In this paper, we present a set of generic ideas that allow to produce synthetic models for many types of turbulent flows (HD incompressible or isothermal, self-gravitating or not, incompressible MHD, etc...). We present our multiscale turbulence synthesis (MuScaTS) technique in section 2 and compare it to previous models in section 3. In section 4, we carefully assess its validity in the simplest framework of 2D incompressible hydrodynamical turbulence. To this aim, on one hand we compute the results of DNS evolving an initial Gaussian field after some finite time, and, on the other hand, we generate snapshots of synthetic turbulence using several variants of our method. In order to

quantitatively compare the two, we compute power spectra, increments statistics, transfer functions and we finally assess the quality of the textures we generate by computing scattering transforms statistics (see Allys et al. 2019). We discuss our results in section 5 and develop the prospects of the MuScaTS method in section 6.

## 2. Synthetic models of turbulence: the MuScaTS technique

In the following subsections, we present the ingredients of our framework: a multiscale filtering (2.1), and a generic form (2.2) for partial differential equations which survives filtering under some additional approximations (2.3). This allows us to build a synthetic flow with a sweep from large to small scales (2.4). We then discuss in more details our two main approximations regarding time evolution (2.4.1) and coherence time estimation (2.4.2).

### 2.1. Isotropic filtering

The purpose of this section is to define a continuous series of filters labelled by their respective scales. In particular, this will allow us to reconstruct any zero mean field from its filtered components, and to separate it into a small and a large scale components.

We consider tensor fields (such as the density or the velocity fields, for example) on a space domain $D = \mathbb{R}^d$ of dimension $d$. For a field $f$, we define its filtered version $f_\ell$ at a scale $\ell$ through:

$$f_{\ell} = \varphi_{\ell} * f, \qquad (2.1)$$

where the sign $*$ denotes convolution and $\varphi_\ell$ is a bandpass filter which selects scales on the order of $\ell$.

We define a continuous collection of filters labelled by all possible scales $\ell \in \mathbb{R}^+$ by dilation from a single mother filter $\varphi$ through the formula

$$\varphi_{\ell}(\mathbf{x}) = \varphi(\mathbf{x}/\ell)/\ell^{d}, \qquad (2.2)$$

where the mother filter $\varphi$ selects scales on the order of 1, i.e. its Fourier transform $\hat{\varphi}(\mathbf{k})$ is real positive, decays at small and large wavenumbers and peaks at around $|\mathbf{k}| = 1$, where $|\mathbf{k}|$ denotes the norm of $\mathbf{k}$. Note that since $\hat{\varphi}(\mathbf{0}) = 0$, $\varphi$ has zero mean and so has $f_\ell$. The dilation formula in Fourier space then takes the form

$$\hat{\varphi}_{\ell}(\mathbf{k}) = \hat{\varphi}(\ell\mathbf{k}). \qquad (2.3)$$

We further assume isotropic filtering, i.e. $\hat{\varphi}(\mathbf{k})$ depends only on $|\mathbf{k}|$ through the function $\Phi(s)$ defined for $s \in \mathbb{R}^+$ such that

$$\hat{\varphi}(\mathbf{k}) = \Phi(|\mathbf{k}|). \qquad (2.4)$$

Now, let us define

$$C_{\Phi} \equiv \int_{0}^{+\infty} \frac{ds}{s} \Phi(s), \qquad (2.5)$$

or equivalently, using a short-hand notation that we will use throughout the paper for such integrals,

$$C_{\Phi} \equiv \int_{-\infty}^{+\infty} d \ln s \Phi(s). \qquad (2.6)$$
We will assume that this integral is finite ($0 < C_{\Phi} < +\infty$), which implicitly constrains the shape of our mother filter $\hat{\varphi}$ because $\hat{\varphi}$ should decay fast enough at both $|k| \to 0$ and $|k| \to +\infty$ in order to keep this integral finite. Without loss of generality, we can assume that

$$C_{\Phi} = 1, \qquad (2.7)$$

by redefining $\Phi \to \Phi/C_{\Phi}$. For example, the choice

$$\Phi(s) \equiv s^2 e^{-s^2} / \int_0^{+\infty} x e^{-x^2} dx, \qquad (2.8)$$

results in a normalised $C_{\Phi}$.

In appendix A we show that the normalisation (2.7) allows us to reconstruct simply any field $f$ (with zero mean) from its filtered components $f_{\ell}$, namely we have the following reconstruction formula

$$f = \int_{-\infty}^{+\infty} d \ln \ell f_{\ell}. \qquad (2.9)$$

Finally, we denote

$$f_{<\ell} = \int_{s<\ell} d \ln s f_s \qquad (2.10)$$

and

$$f_{\ge \ell} = \int_{s \ge \ell} d \ln s f_s \qquad (2.11)$$

the components of $f$ at scales respectively smaller and larger than a given $\ell$, so that

$$f = f_{<\ell} + f_{\ge \ell}. \qquad (2.12)$$

## 2.2. Scope of our method: generic evolution equations

We present here a generic form of the partial differential equations for which the approximations inherent to our method can be more easily justified. Consider a $n$-components vector $W$ characterising the state of the gas as a field on our space domain $D$. For example, for incompressible 3D fluid dynamics, we can take $W = (w)$, with $w$ the vorticity vector. Or we could use $W = (w, \rho)$ with $\rho$ the mass density for compressible HD.

We will then consider evolution equations in the form

$$\partial_t W + (v[W].\nabla)W = S[W].W + D[W], \qquad (2.13)$$

where $v$ is an advection velocity vector (i.e. it has the same number of components $d$ as the space dimension), $S$ is a $n \times n$ deformation matrix and $D$ is a term with $n$ components which can (generally but not always) characterise the diffusion or dissipation terms. Here, the brackets $[W]$ indicate that $v, S$ and $D$ are affine functions of the state variables field $W$. In other words, we require the non-linearities to be at most quadratic in the state vector $W$. We keep the freedom of $v$ having a linear dependence on $W$, to accommodate future cases where the advection velocity might be more complicated, but all practical results presented in this paper are for $v[W] = u$, where $u$ is the velocity component in the state vector $W$.

In the following paragraphs, we give a few explicit examples of evolution equations which fit our framework, for the vorticity and divorticity in incompressible 2D hydrodynamics, and for the vorticity in 3D incompressible hydrodynamics. We then list a few more examples without giving proof.

### 2.2.1. Incompressible 2D hydrodynamics

As a first example, consider the evolution equation for 2D incompressible hydrodynamics, which reads:

$$\partial_t w + u \cdot \nabla w = \nu \Delta w, \qquad (2.14)$$

where $w$ is the vorticity component along the normal to the space domain, $u$ is the fluid velocity and $\nu$ is the kinematic viscosity coefficient. It takes the form (2.13) if we simply set $W = (w)$, $v[W] = u$, the fluid velocity as expressed from $w$ by the 2D Biot-Savart formula (without the boundaries' term)

$$u(x) = \int_D d^2 y w(y) \hat{z} \times (x - y) / |x - y|^2 / (2\pi), \qquad (2.15)$$

where $\hat{z}$ is the unit vector normal to the space domain $D$, $S = 0$ and $D[W] = \nu \Delta w$. Note that equation (2.15) can be efficiently computed by using Fourier transforms.

### 2.2.2. Divorticity in 2D hydrodynamics

In 2D incompressible hydrodynamics, it turns out one can also write evolution equations for the curl of the vorticity, a.k.a. the divorticity (Kida 1985; Shivamoggi et al. 2024):

$$B \equiv \nabla \times (w \hat{z}). \qquad (2.16)$$

Taking the curl of the evolution equation of vorticity (2.14), it results:

$$\partial_t B + u \cdot \nabla B = S.B + \nu \Delta B, \qquad (2.17)$$

where $S$ is the velocity gradient matrix, so that $S.B = u \cdot \nabla B$. This equation is in the form (2.13) when we set $W = (B)$, and it is also strikingly similar to the equation for the evolution of vorticity in 3D, see (2.18) below. Note that the velocity field can be recovered from the divorticity by inverting the curl twice, which is a linear operation, and so the matrix $S$ also depends linearly on the field $B$. We will use equation (2.17) as a test bed in 2D on how to handle the respective advection and deformation terms.

### 2.2.3. 3D incompressible hydrodynamics

For 3D incompressible hydrodynamics, the evolution equation for the vorticity $w$ reads

$$\partial_t w + u \cdot \nabla w = w \cdot \nabla u + \nu \Delta w. \qquad (2.18)$$

To show that this equation takes our generic form (2.13), we can set $W = w \equiv \nabla \times u$ the vorticity vector, and then (again, $v[W] = u$ as explained in the introduction of this section 2.2)

$$u = \int_D d^3 y w(y) \times (x - y) / |x - y|^3 / (4\pi) \qquad (2.19)$$

(3D Biot-Savart formula without the boundaries' term). $S$ is the velocity gradient matrix and $D[W] = \nu \Delta w$ is simply the diffusion term.

Note that if we write the incompressible Navier-Stokes equation for the velocity evolution, with the specific pressure $p \equiv P/\rho$ (with $P$ the thermal pressure and $\rho$ the mass density):

$$\partial_t u + u \cdot \nabla u = - \nabla p + \nu \Delta u, \qquad (2.20)$$
<figure><img src="image_2.png" /><figcaption>**Fig. 2.1.** Example of the performance of our method. Panel (a) shows the common initial (Gaussian) vorticity field used in the two other panels, namely a standard 2D incompressible hydrodynamics simulation (panel (b)) and one of our synthesis (panel (c)). Panel (b) shows the vorticity field after the field has been evolved for about one third of a turnover time, using more than 300 standard simulation steps. Panel (c) shows the resulting synthesis after only 4 MuScaTS steps for a total computer processing unit (CPU) cost 17 times smaller.</figcaption></figure>

it does not take the form (2.13). Indeed, for the choice $W = u$, the specific pressure field and its gradient can classically be recovered in the Leray formulation (Majda & Bertozzi 2001, section 1.8) by solving for $p$ in $\Delta p = -\partial_i u_j \partial_j u_i$, but $\nabla p$ then cannot seem to be put into the form $S[u].u$ with $S[u]$ affine in $u$. However, this is not a restriction since we can simply use equations (2.18) with (2.19), but we feel this example clarifies the subtleties in the scope of the generic form (2.13).

### 2.2.4. Other cases

The evolution equations for reduced MHD, incompressible MHD, or self-gravitating isothermal gases can also be massaged into our generic form (2.13). However, we could not make compressible MHD strictly fit our framework, as the Lorentz term combines a product of three primitive variables (and our attempts at changing variables always led to a product of three variables).

In the context of large scale structure formation in cosmology, Euler-Poisson equations in the dust approximation (zero pressure) are easily cast into our generic form (2.13) for the state vector $\mathbf{W} = (u, \rho)$. Indeed, Poisson equation for the gravitational potential is linear in the density. The expansion of the Universe can also be incorporated by switching to properly defined comoving variables. We will use this in section 3.3 to create a link with the multiscale Zel'dovich approximation (Bond & Myers 1996; Monaco et al. 2002; Stein et al. 2019). Nevertheless, the form (2.13) allows for a large variety of fluid evolution equations.

## 2.3. Derivation of the filtered evolution equations

In this section we present a series of approximations leading to a simplified description of the dynamics, namely equation (2.28). From section 2.4 and onward, the latter equation will be the central focus, as, instead of trying to time integrate the more complicated full partial differential equation (2.13), we will mimic the dynamics from (2.28). Thus, we will be able to generate cheap but realistic turbulent fields, as illustrated in figure 2.1.

As a first step, let us recast the evolution equation (2.13) into the form (2.24) below, which spells out how the various scales contribute to the dynamics. To do so, note that because our filtering (2.1) is a convolution product, it com-

mutes with linear functions<fnref n="1" />, so we have in particular

$$ 
\begin{aligned}
\mathbf{v}[\mathbf{W}]_{\ell} &= \mathbf{v}[\mathbf{W}_{\ell}], \\
\mathbf{S}[\mathbf{W}]_{\ell} &= \mathbf{S}[\mathbf{W}_{\ell}], \\
\mathbf{D}[\mathbf{W}]_{\ell} &= \mathbf{D}[\mathbf{W}_{\ell}],
\end{aligned} 
\quad (2.21) $$

and, filtering (2.13) with a given bandpass filter $\varphi_\ell$,

$$ \partial_t \mathbf{W}_\ell = \varphi_\ell * (-\mathbf{v}[\mathbf{W}] \cdot \nabla \mathbf{W} + \mathbf{S}[\mathbf{W}] \cdot \mathbf{W}) + \mathbf{D}[\mathbf{W}_\ell]. \quad (2.22) $$

Then, using the scale decomposition (2.12), we separate the large and small-scale parts of the velocity field as

$$ \mathbf{v} = \mathbf{v}[\mathbf{W}_{<\ell} + \mathbf{W}_{\ge\ell}] = \mathbf{v}[\mathbf{W}_{<\ell}] + \mathbf{v}[\mathbf{W}_{\ge\ell}] = \mathbf{v}_{<\ell} + \mathbf{v}_{\ge\ell}, \quad (2.23) $$

and doing the same for $S$ we rewrite the filtered evolution equation (2.22) as

$$ \begin{aligned} \partial_t \mathbf{W}_\ell \rvert_x = \mathbf{D}[\mathbf{W}_\ell] \rvert_x & - \int_D d^d y \varphi_\ell \rvert_y \mathbf{v}_{\ge\ell} \cdot \nabla \mathbf{W} \rvert_{x-y} \\ & + \int_D d^d y \varphi_\ell \rvert_y \mathbf{S}_{\ge\ell} \cdot \mathbf{W} \rvert_{x-y} + \int_D d^d y \varphi_\ell \rvert_y \mathbf{S}_{<\ell} \cdot \mathbf{W} \rvert_{x-y}, \end{aligned} \quad (2.24) $$

where $\mathcal{D}$ is the space domain. We expressed convolution products with a notation of the form $\dots$ to help the reader keep track of which variable each field depends on. In particular bear in mind that (2.24) is evaluated at position $x$.

The above lengthy but explicit formulation of the dynamics is the viewpoint from which we now proceed to a series of simplifying assumptions, which will result in the approximate evolution equation (2.28).

Our first approximation consists in considering that the factors $\mathbf{v}_{\ge\ell}$ and $\mathbf{S}_{\ge\ell}$ appearing inside two of the integrals in (2.24) are slowly varying in space compared to the filter $\varphi_\ell$, which can be understood as a scale separation approximation. Since these factors presumably remain approximately constant over the local kernel of $\varphi_\ell$, we pull them out of the integrals and write

$$ \begin{aligned} \partial_t \mathbf{W}_\ell \approx & -\mathbf{v}_{\ge\ell} \cdot \nabla \mathbf{W}_\ell + \mathbf{S}_{\ge\ell} \cdot \mathbf{W}_\ell + \mathbf{D}[\mathbf{W}_\ell] \\ & - \int_D d^d y \varphi_\ell \rvert_y \mathbf{v}_{<\ell} \cdot \nabla \mathbf{W} \rvert_{x-y} + \int_D d^d y \varphi_\ell \rvert_y \mathbf{S}_{<\ell} \cdot \mathbf{W} \rvert_{x-y}. \end{aligned} \quad (2.25) $$

In so doing, we have also taken the $\nabla$ sign out of the integrals, an operation allowed by the property of the gradient of a convolution product $\partial_i(f * g) = f * \partial_i g$ for multivariable functions $f$ and $g$. The last two terms on the right hand side are solely due to scales $\ell$ and below, and represent respectively advection and deformation from small scales convolved with the filter at scale $\ell$.

As a second approximation, we consider that these last two terms evolve more rapidly in time than the previous large scale terms, and somewhat chaotically, i.e. scales smaller than $\ell$ introduce some stochastic perturbations which affect the scale $\ell$. Schematically:

$$ \begin{aligned} \partial_t \mathbf{W}_\ell \approx & -\mathbf{v}_{\ge\ell} \cdot \nabla \mathbf{W}_\ell + \mathbf{S}_{\ge\ell} \cdot \mathbf{W}_\ell + \mathbf{D}[\mathbf{W}_\ell] \\ & + \text{"smaller scales stochasticity"}. \end{aligned} \quad (2.26) $$

<footnote n="1">¹ For affine functions, the affine constant becomes filtered, so the notations have to be slightly amended, but the derivation is the same.</footnote>
Nevertheless, we expect the evolution to remain coherent over some time $\tau_\ell$ at a given scale $\ell$. Thereafter, we will either prescribe or estimate this time scale $\tau_\ell$ from the larger scale fields, and we will provide several ways to estimate it in section 2.4.1 below. For now, since $\tau_\ell$ controls the time scale over which two initially adjacent fluid particles diverge from each other, it can be seen as the inverse of the largest positive Lyapunov exponent under the combined effects of the advection flow $v_{\ge \ell}$ and the deformation matrix $S_{\ge \ell}$.

This sensitivity to the initial conditions implies that it is not quite necessary to integrate precisely the flow equations during more than a coherence time. Indeed, if the conditions had been slightly different one coherence time before, the current state of the fluid might now be completely different, to a point these initial conditions might as well be chosen randomly. Therefore, rather than solving a complicated stochastic equation, we will also assume that the effect of the chaos in the small scales is equivalent to randomly reshuffle initial conditions every coherence time $\tau_\ell$.

Thanks to this hypothesis, we do not need to integrate the evolution from the beginning, but only during the length of the last coherent time interval. At each scale $\ell$, we hence assume that the chaotic evolution can be approximated by deterministically integrating to the targeted time $t$ from random conditions at time $t - \tau_\ell$ during a coherence time interval $\tau_\ell$. We denote this change by adding tildes over the vectors subject to this approximation. As a result, we approximate the evolution of the field as:

$$ \partial_t \tilde{W}_\ell = -v_{\ge \ell} \cdot \nabla \tilde{W}_\ell + S_{\ge \ell} \tilde{W}_{\cdot, \ell} + D[\tilde{W}_\ell], \quad (2.27) $$

starting from random conditions $W_\ell^0$ at intermediate time $t_\ell^0 = t - \tau_\ell$ (these intermediate times might hence depend on scale).

The law of $W_\ell^0$ still needs to be specified. The chaos induced by the small scales might well be non trivial, but in a first step we are going to assume that it follows a Gaussian law. Hence we assume that $W_\ell^0$ is a filtered version $W_\ell^0 = \varphi_\ell * W^0$ of a Gaussian field $W^0$, for which the spectrum needs to be prescribed. Furthermore, to reflect the property that neighbouring scales are in fact not independent (we expect the behaviour of the fluid to be continuous across scales), we assume that the fields $W_\ell^0$ are filtered versions of a single random field $W^0$.

We approximate the evolution equation one step further by using the advection velocity field $\tilde{v}_{\ge \ell}$ computed from the synthetic state variable $\tilde{W}_{\ge \ell} \equiv \int_{s \ge \ell} d \ln s \tilde{W}_s$ in equation (2.27), and similarly for $S_{\ge \ell}$, for which we use $\tilde{S}_{\ge \ell}$.

Finally, the line of reasoning detailed above ultimately leads to the following approximate evolution equation for the dynamics at a given scale $\ell$

$$ \partial_t \tilde{W}_\ell \approx -\tilde{v}_{\ge \ell} \cdot \nabla \tilde{W}_\ell + \tilde{S}_{\ge \ell} \cdot \tilde{W}_\ell + D[\tilde{W}_\ell]. \quad (2.28) $$

Interestingly, it is almost in the initial generic form (2.13), except that the advection velocity and deformation matrix now depend on all scales above $\ell$.

In the following, we will introduce yet more approximations, now tailored to solve approximately the approximate evolution equation (2.28) we just derived, i.e. we will present our turbulence synthesis procedure. At this stage,

the reader may legitimately be worried that such an accumulation of (even small) losses of information ought to be too drastic, and can eventually only move us far away from the correct evolution. If so, we invite the reader to recall figure 2.1 (or to have a glimpse at the more detailed figures 4.1 and 4.2) were we showed a non-trivial case were the synthesis and the full simulation do match very well. Therefore, to some extent the whole procedure does keep enough relevant information to be of practical use, i.e. we have meaningful examples where it seems to extract the essence of the dynamics.

## 2.4. Turbulence synthesis

Our last equation (2.28) is in closed form for all scales above $\ell$: the evolution for $\tilde{W}$ at scale $\ell$ can be solved once all scales $s$ above $\ell$ are known. We can hence sweep all scales from the largest down to the smallest and compute $\tilde{W}_\ell$ at all scales $\ell$ (we will discuss how we deal with the largest scale in the numerical implementation below, see section 2.5). The resulting synthesised field can finally be reconstructed as $\tilde{W} = \int d \ln \ell \tilde{W}_\ell$.

Variants of our model will now be defined by the way we approximate the evolution equation (2.28) (detailed in section 2.4.1) and the way we estimate the coherence time $\tau_\ell$ (developed in section 2.4.2).

### 2.4.1. Time evolution approximations

Here, we focus on our approximations to solve the evolution equation (2.28). We assume that the fields $\tilde{v}_{\ge \ell}$ and $\tilde{S}_{\ge \ell}$ evolve slowly compared to the fields at scale $\ell$, as they are based on larger scales. We therefore assume they remain constant during the time interval $\tau_\ell$. For instance, we assume that the fluid parcels move ballistically with the velocity $\tilde{v}_{\ge \ell}$ during a coherence time $\tau_\ell$ to integrate (2.28) in time.

We first focus on a case where $S = D = 0$ and $v \neq 0$. In this case, we only need to randomly draw the initial vorticity field $W_\ell^0$ at time $t_\ell^0$ and remap it at positions $x^0 = x - \tau_\ell \tilde{v}_{\ge \ell}$ where the fluid parcel used to be:

$$ \tilde{W}_{\ell}^{v \neq 0}(x) = W_{\ell}^{0}(x - \tau_{\ell} \tilde{v}_{\ge \ell}). \quad (2.29) $$

Possible shell crossing will be discussed in the next subsection 2.4.2.

We now focus on the case $v = D = 0$ ($S \neq 0$). Since $\tilde{S}_{\ge \ell}$ is constant during $\tau_\ell$, the solution of $\partial_t \tilde{W}_\ell = \tilde{S}_{\ge \ell} \cdot \tilde{W}_\ell$ is obtained as a matrix exponential :

$$ \tilde{W}_{\ell}^{S \neq 0}(x) = \exp(\tau_{\ell} \tilde{S}_{\ge \ell}) \cdot W_{\ell}^{0}(x). \quad (2.30) $$

For the general case, we combine the two previous steps as

$$ \tilde{W}_{\ell}^{D=0}(x) = \exp(\tau_{\ell} \tilde{S}_{\ge \ell}) \cdot W_{\ell}^{0}(x - \tau_{\ell} \tilde{v}_{\ge \ell}), \quad (2.31) $$

and we finally apply the diffusion step as

$$ \tilde{W}_{\ell} \approx \tilde{W}_{\ell}^{D=0} + \tau_{\ell} D[\tilde{W}_{\ell}^{D=0}]. \quad (2.32) $$

We take those three steps in a particular order: advection, deformation, then diffusion. We feel it is justified by
the need to use the quantities $v$ and $S$ at a time as advanced as possible (i.e. as close as possible to the time at which we compute the synthesis flow). Also the computation of the gradients in the deformation matrix requires a synchronous evaluation, so deformation might be better evaluated after advection. We finish by the diffusion step because it somewhat smooths the fronts generated at the convergence points of advection trajectories and where the deformation introduces stretching (a property sought for in adhesion models, for example, see Gurbatov et al. 1989). This also allows to alleviate the small scales generated outside the current filter by the non-linear action of advection and deformation (see discussion in section 5.6, also). However, in the limit where $\tau_\ell$ is very small, time evolution is linear and the order does not matter.

### 2.4.2. Coherence time estimates

So far, we have defined the coherence time in this work as the time during which the evolution of the fluid can be considered deterministic at a given length scale. We will now examine three ways to estimate this time.

1. It is at first tempting to identify it to the correlation time scale. For example, in the direct enstrophy cascade of 2D turbulence, below the injection scale, the correlation time scale is constant, equal to the injection time scale. It would then make sense to adopt $\tau_\ell$ = constant at all scales.

2. The notion of coherence also appeals to the sensitivity to initial conditions. The rate at which fluid trajectories deviate from each other is linked to the symmetric part of the velocity gradient matrix, which characterises the evolution of distances between nearby fluid parcels. We define the stretching time as the inverse of its largest eigenvalue which in incompressible fluids writes as (Haller 2021):

$$ \tau_{\text{stretch}}^{-2}(\mathbf{u}) = |\partial_x u_x \partial_y u_y - \frac{1}{4} (\partial_x u_y + \partial_y u_x)|. \quad (2.33) $$

This eigenvalue can also be seen as the largest instantaneous Lyapunov exponent for fluid trajectories. The above expression for incompressible hydrodynamics (in which $\partial_x u_x = -\partial_y u_y$) can be recast into

$$ 4\tau_{\text{stretch}}^{-2}(\mathbf{u}) = (\partial_x u_x - \partial_y u_y)^2 + (\partial_x u_y + \partial_y u_x)^2 = \tau_{\text{strain}}^{-2}, \quad (2.34) $$

which expresses the stretch rate as a sum of squares of the normal and shear strains. That time scale is hence also naturally related to the deformation of fluid elements. Indeed, in some places the fluid might be rotating in a vorticity clump and its evolution at scales below the size of that vortex could be considered deterministic. On the contrary, fluid parcels in between vortices experience strong deformation through shearing motions. Adopting $\tau_\ell = \tau_{\text{stretch}}(\mathbf{u}_\ell)$ would also ensure that fluid elements retain their integrity (i.e. they are not too deformed and it still makes sense to follow the same fluid particle).

3. Since we resort to a ballistic approximation for the advection of the fluid, our modeled fluid trajectories may cross each other after some finite time. This shell-crossing time can be defined as the inverse of the largest eigenvalue of the velocity gradient (by contrast to its symmetric part above):

$$ \tau_{\text{shell}}^{-2}(\mathbf{u}) = |\partial_x u_x \partial_y u_y - \partial_x u_y \partial_y u_x|. \quad (2.35) $$

Article number, page 6 of 19

Integration of ballistic trajectories for more than this time scale does not make sense. Note that this constraint is more related to our ballistic approximation than to the stochastic behaviour of the fluid, and may be relaxed once we resort to a higher order approximation. However, in the case of 2D incompressible hydrodynamics, this time scale is related to the physically motivated Okubo-Weiss (Okubo 1970; Weiss 1991; Haller 2021; Shivamoggi et al. 2024) criterion $\tau_{\text{shell}}^{-2} = |\tau_{\text{stretch}}^{-2} - w^2|$. That last expression can also be used to show that requesting $\tau_\ell = \tau_{\text{shell}}$ is slightly less stringent than the stretch criterion ($\tau_{\text{stretch}} \le \tau_{\text{shell}}$).

To conclude, because it is not yet clear which of these approaches is best suited, we have tested four choices for the coherence time (from less to more stringent ($\tau_{\text{strain}} < \tau_{\text{stretch}} \le \tau_{\text{shell}}$)):

- Constant $\tau_\ell = t$

- Shell-crossing based $\tau_\ell^{-2} = t^{-2} + \tau_{\text{shell}}(\mathbf{u}_{\ge \ell})^{-2}$

- Stretch based $\tau_\ell^{-2} = t^{-2} + \tau_{\text{stretch}}(\mathbf{u}_{\ge \ell})^{-2}$.

- Strain based $\tau_\ell^{-2} = t^{-2} + \tau_{\text{strain}}(\mathbf{u}_{\ge \ell})^{-2}$.

The minus-two powers in the above expressions take soft minima between the various time scales involved and the integration time $t$. Indeed, when we compare our method to actual simulations, we compute the evolution of the fluid over a finite time $t$ from random initial conditions. In this case, it is not necessary to integrate the evolution of the synthesised flow over times longer than $t$.

As a final note, most of these definitions imply that the coherence time $\tau_\ell(\mathbf{x}, t)$ depends on both scale, position in space and total evolution time $t$ as long as this one is not too long compared to local times of the flow. The initial conditions $\mathbf{W}_\ell^0$ should hence not be understood as starting conditions at a synchronous time. Rather, they are random parameters which reflect the sensitivity to initial conditions.

## 2.5. Numerical implementation

### 2.5.1. Discrete filters

A numerical implementation of the above method requires a space domain with a finite extent and a finite space resolution. As a result, the Fourier domain spans a finite number of wavenumbers, and only a finite number of scales is necessary to recover a given field. We choose to consider a cubic (square in 2D) domain of size $L = 2\pi$ gridded with $N$ pixels and we adopt periodic boundary conditions. The available wavevectors therefore span a regular grid of unit $k_J = 2\pi/L = 1$. We will restrict ourselves to a discrete set of $J + 1$ filters $\{\varphi_j\}_{j=0,...,J}$ logarithmically distributed in scales in such a way that

$$ \hat{\varphi}_{j+1}(k) = \hat{\varphi}_j(k/\lambda), \quad (2.36) $$

with $\lambda \in ]0, 1[$ the ratio between the typical scales of two consecutive filters. This parameter controls the resolution in scales. In the present application, we mostly use $\lambda = 1/2$ but we also tested a finer set of scales with $\lambda = 1/\sqrt{2}$ (hence with twice as many filters), without noticing much change.

We request each filter to be selective in logarithmic scales, in order to be able to use the scale separation arguments at the base of our approximations. The filters $\hat{\varphi}_j$ hence select concentric coronas in Fourier space and are centred around wavenumbers $k_j = k_J \lambda^{j-J}$ with a relative
width of $1/\lambda$. In other words, the typical scale associated with filter $j$ is

$$ \ell_j = L\lambda^{J-j}. \quad (2.37) $$

We point out that these physical scales go in increasing order from the smallest to the largest as $j$ increases.

The discretised version of our normalisation condition (2.7) will simply be to enforce $\sum_{j=0}^J \hat{\varphi}_j = 1$ for any $k$ contained in the computational domain. Indeed, then the reconstruction formula works as

$$ \sum_{j=0}^{J} \hat{f}_j = \hat{f} \cdot \sum_{j=0}^{J} \hat{\varphi}_j = \hat{f}, \quad (2.38) $$

from which we deduce $f = \sum_{j=0}^J f_j$. We present two example filters we used in the present application in the appendix B.

### 2.5.2. Algorithm

The initial Gaussian field $W^0$ is first drawn in real space as a Gaussian white noise. We then convolve it with the desired spectrum by multiplying its Fourier coefficients.

We present below the steps of the algorithm we use to generate turbulent fields $\tilde{W}_j$ at all scales $j$ from the initial Gaussian field $W^0$. We write $\tilde{W}_j = \tilde{W} * \hat{\varphi}_j$ and $W_j^0 = W^0 * \hat{\varphi}_j$.

To initiate the first $j = J$ step of our iterative algorithm, at the largest scale, we simply assume $\tilde{W}_J = W_{j=J}^0$. We then parse the successive scales downward, and we successively apply the following operations at each scale:

1. At stage $j < J$, the fields $\tilde{W}_i$ are known for all $i > j$.

2. We write $\tilde{W}_{\ge j} = \sum_{i>j} \tilde{W}_i + W_j^0$. We compute the advection velocity $v_{\ge j} = v[\tilde{W}_{\ge j}]$ and similarly $S_{\ge j} = S[\tilde{W}_{\ge j}]$.

3. We estimate $\tau_j = \tau_{\ell_j}$ (for which we might need $v_{\ge j}$ if we use the local stretch or shell crossing times, for example).

4. We now replace these values in the construction formulas (2.31) and (2.32). We compute the time-evolved field at scale $\ell_j$ (without diffusion) as

$$ \tilde{W}_j^{D=0}(\mathbf{x}) = \exp(\tau_j \tilde{S}_{\ge j}(\mathbf{x})) \cdot W_j^0(\mathbf{x} - \tau_j \mathbf{v}_{\ge j}) \quad (2.39) $$

where we use a first order linear interpolation scheme to compute $W_j^0$ at positions in between grid points.

5. We finally apply diffusion to the resulting field to get $\tilde{W}_j = \tilde{W}_j^{D=0} + \tau_j D[\tilde{W}_j^{D=0}]$. We now have $\tilde{W}_i$ known for all $i > j - 1$ and we can proceed further back to operation 1 down to the next scale $j - 1$.

The loop stops when we hit the smallest scale labelled by $j = 0$. The final field then results from $\tilde{W} = \sum_{j=0}^J \tilde{W}_j$.

For the interpolation at operation 4, we use RegularGrid Interpolator from the scipy.interpolate Python package. One of the caveats of this interpolation remapping is that it is not divergence and mean preserving. One may need to reproject the divergence free fields and adjust preserved quantities after the interpolation process. In our 2D hydrodynamical application in section 4,

we remove its average from the vorticity field to recover the zero mean. The same holds for the deformation step.

The matrix exponential at operation 4 makes use of the Padé approximant implementation from the EXPOKIT package (Sidje 1998), for which we implemented a vectorised wrapper in order to gain efficiency (essentially to call it from Python in a more efficient way for an array of vectors).

Diffusion at operation 5 is computed in Fourier assuming it proceeds over a homogenised time $\bar{\tau}_j$ averaged over the computational domain. In this simplified case, the effect of incompressible viscous diffusion $\partial_t w = D[w] = \nu \Delta w$ during time $\bar{\tau}_j$ is simply to multiply by $\exp(-\nu k^2 \bar{\tau}_j)$ the Fourier coefficients of the initial field $w$. Note that because $\bar{\tau}_j$ is not necessarily uniform over the computational domain, we resort to using a volume averaged value for $\tau_j$. Otherwise, there is to our knowledge no convenient way to integrate the viscous diffusion term.

### 2.5.3. CPU cost scaling

If we neglect the cost of the matrix exponentiation and the remapping (proportional to the number of pixels), most of the computational cost will amount to Fourier transforms. At each scale index $j > 0$ we need to Fourier transform back and forth because some operations are more efficiently performed in Fourier space (computing the advection velocity, the deformation matrix, or applying the diffusion process) while others are better done in real space (interpolation remap and deformation). To be more precise, the 2D synthesis uses 10 calls to the fast Fourier transform (FFT) for each scale level: 2 for the advection velocity, 4 for the deformation matrix, 2 around the remap step and 2 around the diffusion step. A 3D synthesis would require 3 for the advection velocity, 9 for the deformation matrix, and 6 for each of the remap and diffusion step: a total of 24 calls to a scalar FFT. In general, the method hence requires 10$J$ (2D) or 24$J$ (3D) Fourier transforms to generate one field. Since the CPU cost of a Fourier transform in $d$ dimensions is on the order of $N^d \log_2 N$ operations and $J$ scales like $\log_2 N$, the cost of one field generation scales as

$$ N_{2\text{Dsyn}} = 20N^2(\log_2 N)^2 \quad (2.40) $$

operations in 2D and

$$ N_{3\text{Dsyn}} = 72N^3(\log_2 N)^2 \quad (2.41) $$

in 3D.

The cost of one time step of a pseudo-spectral code is proportional to the order of the time integrator: a Runge-Kutta order 4 such as we use in section 4.1 requires 4 evaluations of the time derivative of the quantities at play. Each of these evaluation costs 5 FFTs in 2D (2 for the velocity, 2 for the diagonal components of the vorticity gradient, 1 to get back to Fourier space), and $3+3+9+3=18$ in 3D where the whole velocity gradient terms are also needed. At a Courant number of 1, the Courant-Friedrich-Lewy condition constrains the time-step to be less than $t_{\text{turnover}}/N$: one requires $N$ time steps to reach one turnover time. The cost of a 2D simulation is hence on the order of

$$ N_{2\text{Dsim}} = 40N^3 \log_2 N \cdot t/t_{\text{turnover}} \quad (2.42) $$

operations for our Runge-Kutta order 4 implementation, and

$$ N_{3\text{Dsim}} = 216N^4 \log_2 N \cdot t/t_{\text{turnover}} \quad (2.43) $$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.2.** CPU cost scaling estimates for 2D (light colors) and 3D (dark colors) hydrodynamics simulations (solid, blue) and MuScaTS syntheses (dashed, red). The scalings are computed according to text, calibrated on a 2D simulation (pseudo-spectral code implemented in Fortran) at $N = 1024$ for one turnover time (right most blue circle). The total CPU cost for simulations at $N = 256$ and $N = 512$ (the other two blue circles) are also indicated. Red squares stand for the CPU cost of a few 2D MuScaTS syntheses.</figcaption></figure>

in 3D.

For a simulation over the length of one turnover time, the CPU time gain factor is hence $2N/\log_2(N)$ in 2D and $3N/\log_2(N)$ in 3D. Both 2D and 3D scalings for simulations and MuScaTS syntheses are displayed on figure 2.2. The scalings are all calibrated on a 2D $N = 1024$ simulation run (rightmost blue circle) using our pseudo-spectral code (see section 4.1) compiled with Fortran on the totoro server at the mesoscale centre mesoPSL. For reference, the CPU hardware on this server is Intel(R) Xeon(R) Gold 6138 CPU, 2.00GHz. Other circles indicate the scaling of simulations at $N = 256$ and $N = 512$: they scale slightly less well than the theory predicts, because we ran our code in parallel on 32 cores, so the communication overhead in the FFT at smaller computational domains (smaller $N$) deteriorates the scaling. The light red squares show the CPU cost on the same hardware for a python 2D implementation of our MuScaTS synthesis.

As a final note, we point out that the current implementation could be further optimised with respect to the multiscale treatment. For a proper choice of filters, we could indeed afford using coarser grids at larger scales than for the smaller scales. This would considerably accelerate the computation at largest scales and lower the number of effective number of scales which enter the CPU cost scaling.

# 3. Comparison to other methods

We investigate how our MuScaTS technique compares to selected existing works.

## 3.1. Comparison to Rosales & Meneveau (2006)

For 3D hydrodynamics, if we neglect the deformation and the diffusion terms, our method is extremely similar to the minimal multiscale Lagrangian map (MMLM) of Rosales &

Article number, page 8 of 19

Meneveau (2006) provided we use $\tau_\ell$ as in their equation (9): $\tau_\ell^{\text{MMLM}} = \ell/u_{\text{r.m.s.}\ell}$ where $u_{\text{r.m.s.}\ell}$ is the prescribed r.m.s. velocity in the initial field at scale $\ell$. The only remaining difference is: we would advect the vorticity field and build the velocity from a Biot-Savart formula, while their Lagrangian map is directly applied to the velocity. They have to reproject their velocity field on $\nabla \cdot \mathbf{u} = 0$ while we get divergence free fields by construction. We hence expect our method should be more realistic, as both vorticity stretching effects and diffusion are retained while in effect, they neglect pressure terms and diffusion.

In a later paper, the same authors Rosales & Meneveau (2008) noticed that they could improve their method by cutting single Lagrangian steps into smaller ones as scales become finer. They considerably improved the realism of the anomalous exponents using this approach, which they call multiscale turnover Lagrangian map (MTLM for short). We note that we could also incorporate this multi-stepping into our MuScaTS method, although it would increase its cost. Our method could perhaps also benefit from the fast algorithms designed by Malara et al. (2016) to improve the efficiency of the MTLM method.

Finally, Subedi et al. (2014) designed an incompressible MHD version of the MMLM method. MuScaTS without diffusion or deformation applied to the choice of variable $\mathbf{W} = (\nabla \times \mathbf{u}, \mathbf{B})$ and using $\tau_\ell = \tau_\ell^{\text{MMLM}}$ should amount to the same model except we advance the vorticity $\nabla \times \mathbf{u}$ rather than the velocity, which directly generates divergence free fields (we would however still have to clean the divergence of the magnetic field). The extension of Subedi et al. (2014) by Lübke et al. (2024) proceeds with advection from an already non-Gaussian field and produced convincing coherent structures. Our own implementation will advance along an already deformed field which contains physically motivated non-Gaussianities: we hence expect that the MHD version of MuScaTS should perform equally well if not better, but this remains to be demonstrated.

## 3.2. Comparison to Chevillard et al. (2011) (CRV)

The link between MuScaTS and CRV is slightly less direct and we need to rework the equations quite a bit to highlight the similarities and differences.

First, let us expose part of our model in a way that will ease the comparison. The work of CRV is in the context of 3D Navier-Stokes incompressible equations, so in our framework we take $\mathbf{W} = (\mathbf{w})$ where $\mathbf{w}$ is the vorticity. We neglect diffusion, or in other words we focus on the inertial range. If we now neglect advection rather than deformation, as in CRV, then only the exponential term remains in equation (2.30) which we rewrite more precisely here as

$$ \tilde{\mathbf{w}}_{\ell}^{S \neq 0}(\mathbf{x}) = \exp(\tau_{\ell} \mathbf{S}[\tilde{\mathbf{w}}_{\ge \ell}]). \mathbf{w}_{\ell}^{0}(\mathbf{x}), \quad (3.1) $$

to make clear the intricate dependence of $\mathbf{S}$ on the larger scales. If to parallel the work of CRV we now approximate the argument of $\mathbf{S}$ by its value ab initio, we get

$$ \tilde{\mathbf{w}}_{\ell}(\mathbf{x}) = \exp(\tau_{\ell} \mathbf{S}[\mathbf{w}_{\ge \ell}^{0}]). \mathbf{w}_{\ell}^{0}(\mathbf{x}) \quad (3.2) $$

or now using equation (2.11), we obtain

$$ \tilde{\mathbf{w}}_{\ell}(\mathbf{x}) = \exp \left( \int_{s \ge \ell} d \ln s \, \tau_{\ell} \mathbf{S}[\mathbf{w}_{s}^{0}] \right) \cdot \mathbf{w}_{\ell}^{0}(\mathbf{x}), \quad (3.3) $$
and, integrating over all scales, we finally get

$$ \tilde{w}^{\text{MuScaTS}} = \int_{-\infty}^{+\infty} d \ln \ell \exp \left( \int_{s \ge \ell} d \ln s \tau_{\ell} \mathbf{S}[\mathbf{w}_s^0] \right) \mathbf{w}_{\ell}^0. \quad (3.4) $$

This formula will be our comparison point to CRV's work below.

Then, let us reformulate part of CRV's model. In their work, they construct a non-Gaussian random velocity field as the result of the vortex stretching of some initial Gaussian vorticity field $\mathbf{w}^0$. Concretely, they consider for $\mathbf{w}^0$ a Gaussian white noise vector (that we will here denote $\eta$) and insert into a Biot-Savart-type formula a proxy for the vorticity field of the form

$$ \mathbf{w}^{\text{CRV}} = \exp(\tau_{\text{CRV}} \mathbf{S}^{\text{CRV}}[\mathbf{w}^0]). \mathbf{w}^0, \quad (3.5) $$

where $\tau_{\text{CRV}}$ is a parameter which controls the amount of stretching, as it relates<fnref n="2" /> to the correlation time-scale of the velocity gradients, and $\mathbf{S}^{\text{CRV}}[\mathbf{w}^0]$ is a matrix field encoding the stretching of $\mathbf{w}^0$ during a typical time $\tau_{\text{CRV}}$. Now, as a general rule, given a vorticity field $\mathbf{w}$ the stretch matrix field reads<fnref n="3" />

$$ \mathbf{S}[\mathbf{w}](\mathbf{x}) \propto \int d^3\mathbf{y} \frac{\mathbf{r} \otimes (\mathbf{r} \times \mathbf{w}(\mathbf{y})) + (\mathbf{r} \times \mathbf{w}(\mathbf{y})) \otimes \mathbf{r}}{r^5}, \quad (3.6) $$

where $\mathbf{r} \equiv \mathbf{x} - \mathbf{y}$ and $\otimes$ denotes the tensor product. In addition, an important feature which is universally observed in both experiments and numerical simulations of turbulence, is that turbulent fields are multi-fractal and have long-range correlations. This formally translates into $\mathbf{S}[\mathbf{w}]$ being log-correlated in space, and a key point in CRV's work which greatly increases the realism of their model, is to implement this feature into their formulae through a simple, effective rule<fnref n="4" />, namely by just changing the exponent in the power-law $r^5$ of the denominator of the general form (3.6). Specifically, by choosing $\eta$ for their $\mathbf{w}^0$ they adopt the rule $r^5 \to r^{7/2}$, i.e. they consider (cf. CRV, equation 12)

$$ \mathbf{S}^{\text{CRV}}[\eta](\mathbf{x}) \propto \int d^3\mathbf{y} \frac{\mathbf{r} \otimes (\mathbf{r} \times \eta(\mathbf{y})) + (\mathbf{r} \times \eta(\mathbf{y})) \otimes \mathbf{r}}{r^{7/2}}, \quad (3.7) $$

and this simple tweak makes their stretch field log-correlated as it should.

Hence, in CRV's work, log-correlation appears through a cunning but somewhat artificial change of exponent, while the stretching time-scale $\tau_{\text{CRV}}$ remains a mere parameter: in particular it is scale-independent. Inspired by their line of thinking, we now suggest an alternative approach in which we retain the physical scaling for the stretch, but we now introduce a scale-dependency of the stretching time-scale to achieve log-correlation. Concretely, starting from expressions (3.5) and (3.6), and imposing log-correlation with a similar scaling argument as CRV's, we now build expression (3.12) below, which is reminiscent of the scale-decomposition (3.4) of our model and may thus be interpreted in terms of a scale-dependent stretching time-scale.

To improve on the physical interpretation of CRV, instead of choosing as initial condition $\mathbf{w}^0$ the simple Gaussian white noise $\eta$ we also adopt a slightly more sophisticated Gaussian field $\mathbf{w}^0$ with the Kolmogorov (1941) power spectrum of the vorticity (i.e. with a spectrum scaling as $E_{w^2}(k) \sim k^{1/3}$). But doing so first raises the question: how is CRV's effective rule then modified? To answer this, we suggest the following rough scaling argument. In Fourier space, for a given corona of radius $k$ and infinitesimal thickness $dk$, by definition of the spectrum $E_{w^2}(k)dk = \hat{w}^2k^2dk$, so our choice of Kolmogorov scaling implies $\hat{w} \sim k^{-5/6}$. Then, since the Fourier transform of a power law $r^{-\alpha}$ is $k^{\alpha-d}$ where $d=3$ is the dimension of space, our scaling in Fourier space corresponds to a scaling $w \sim r^{5/6-3}$ in position space. Hence, compared to CRV's case, we have an additional power $r^{5/6-3}$ in the numerator of the ratio in (3.7), so log-correlation requires that we adjust the exponent in the denominator to compensate it, i.e. we consider

$$ \mathbf{S}^{\text{CRV}}[\mathbf{w}^0](\mathbf{x}) \propto \int d^3\mathbf{y} \frac{\mathbf{r} \otimes (\mathbf{r} \times \mathbf{w}^0(\mathbf{y})) + (\mathbf{r} \times \mathbf{w}^0(\mathbf{y})) \otimes \mathbf{r}}{r^{4/3}}, \quad (3.8) $$

where the denominator should be understood as $r^{7/2+(5/6-3)} = r^{4/3}$. This is our modification to CRV's effective rule.

Now, to draw a parallel with the MuScaTS construction, let us build this stretch matrix $\mathbf{S}^{\text{CRV}}[\mathbf{w}^0]$ from a continuous series of narrow filters, by decomposing it into various scales with (2.9) and using the linearity property (2.21). We thus have

$$ \mathbf{S}^{\text{CRV}}[\mathbf{w}^0] = \int_{-\infty}^{+\infty} d \ln \ell \mathbf{S}^{\text{CRV}}[\mathbf{w}_{\ell}^0], \quad (3.9) $$

but in fact, to be complete, since in CRV's model what drives the dynamics is the full argument of the matrix exponential in (3.5), we also need to take into account the $\tau_{\text{CRV}}$ parameter. Here instead of this mere parameter, we introduce some scale-dependent $\tau_{\ell}$ in the scale-decomposition (3.9) such that

$$ \tau_{\text{CRV}} \mathbf{S}^{\text{CRV}}[\mathbf{w}^0] = \int_{-\infty}^{+\infty} d \ln \ell \tau_{\ell} \mathbf{S}[\mathbf{w}_{\ell}^0], \quad (3.10) $$

where $\mathbf{S}$ is the real, physical stretch (3.6), i.e. with a $r^5$ denominator. Then, to reconcile the $r^{4/3}$ requirement for log-correlation of (3.8) with the real physical $r^5$ dependency of the stretching, we argue (pursuing the hand-wavy reasoning consisting in keeping track and matching length scales) that a power-law dependency

$$ \tau_{\ell} \propto \ell^{5-4/3} = \ell^{11/3}, \quad (3.11) $$

<footnote n="2">² The focus in this section is on the scale dependency of parameters (power-law variation in $\ell$) so throughout this section we omit numerical factors, as in (3.6), but it is worth mentioning that strictly speaking, $\tau_{\text{CRV}}$ does not have the dimension of a time. It is a parameter controlling the intermittency (non-Gaussianity) of the resulting velocity field which appears in CRV's construction from an equation governing the dynamics of stretching, so we choose this $\tau$ notation to help the reader get some intuition, notably since the longer stretching occurs, the larger this parameter is. For more details see CRV's work (Chevillard et al. 2011; Chevillard 2015).</footnote>

<footnote n="3">³ For simplicity and since we focus on the inertial range in this discussion, we omitted the numerical prefactor, the regularization and the principal value present in CRV.</footnote>

<footnote n="4">⁴ In follow-ups of their CRV work, such as Chevillard (2015); Pereira et al. (2018), they do support this with very rigorous but rather lengthy derivations, and we here stick to their simple reasoning yielding the correct result</footnote>
is adequate. Finally, using in (3.5) expression (3.10) and the scale decomposition (2.9) for $w^0$, we rewrite the CRV random field in the following form, which is close to our (3.4),

$$ w^{\text{CRV}} = \int_{-\infty}^{+\infty} d\ln\ell \exp \left( \int_{-\infty}^{+\infty} d\ln s \tau_s S[w_s^0] \right) .w_\ell^0, \quad (3.12) $$

with $\tau_\ell$ given by (3.11) when the vorticity has a Kolmogorov (1941) power-law spectrum. It is in this sense that we mean that the change of power-law to enforce log-correlation in CRV could be interpreted as an integration time varying with scale, with a typical time scale $\tau_\ell \propto \ell^{11/3}$.

This unusual perspective on their theory allows to explicit more of the original physical grounds behind their construction. Thanks to the freedom introduced by the scale dependent stretch parameter, we can both retain the physical expression for the stretch (with $1/r^5$ scaling) and assimilate the generating Gaussian field as the initial vorticity (using Kolmogorov scaling instead of a white noise).

Now, comparing the above reformulations of both formalisms, we reveal three noteworthy differences. First, as can be seen in (3.3), we use a stretch field in the larger scales which has already been advanced, while CRV uses the initial stretch field at all scales. Second, (3.4) confronted to (3.12) shows that in our formulation the effective time scale $\tau_\ell$ applies to all scales above while it is local in scales in CRV. Third, the integration bounds differ: we integrate only the scales above $\ell$ while CRV uses the whole range of scales. In other words, our method suppresses the effects of the smaller scales deformation at a given scale. This might yield less intermittency because the smallest scales which can lead to the largest exponentiated spikes will be suppressed. On the other hand, it may be an improvement for the coherence of the resulting structures compared to Chevillard et al. (2011). Indeed, smoothing avoids strong deformations from very small scales which could result in less correlation between contiguous scales as the strongest peaks at the smallest scales might dominate every scale. In the case of the Zel'dovich Approximation, it was in fact recognised that some smoothing improved the accuracy of the approximation for a given targeted evolution time, which led to the truncated Zel'dovich Approximation (e.g., Coles et al. 1993).

Future investigations of the MuScaTS method in the 3D incompressible HD case will decide which formulation performs best. Finally, note that MuScaTS also retains both advection (neglected in CRV) and diffusion (approximated by a regularisation cut-off at small scales in CRV).

### 3.3. Links to the Zel'dovich multiscale approximation

In cosmology, the concept of multiscale is directly related to hierarchical formation of structures, in particular dark matter halos. In the concordance model of large-scale structure formation, dark matter halos are the hosts of galaxies and clusters of galaxies and their history consists of successive mergers. In this bottom-up scenario, small halos form first and then merge together to form a population of larger halos. The Lagrangian regions occupied by these halos define a smoothing scale, which can itself be associated with a mass and a time scale of halo formation (Press & Schechter

1974). In this picture, mergers are equivalent to connections between Lagrangian regions. As a result, the dynamical process at play can also be described with a multiscale approach accounting for dynamical collapse times of halos of various masses associated with a smoothing level.

The Zel'dovich approximation (ZA) consists in a ballistic integration of initial velocities. Smoothing initial conditions to the ZA was shown to improve the solution, because it avoids the dispersion resulting from shell crossing. This led to the so-called truncated ZA (e.g., Coles et al. 1993) or higher order Lagrangian truncated perturbation theory (e.g. Bernardeau et al. 2002, and references therein), where the level of truncation is the smoothing scale. In the hierarchical formation of halos described above, the truncated ZA (or its higher order equivalent) can be improved further by a multiscale approach: a smoothing length depending on the local position is defined. This multiscale synthesis of dark matter dynamics, has been implemented with success by different authors (see, e.g., Bond & Myers 1996; Monaco et al. 2002; Stein et al. 2019). In these existing multiscale approaches, at every initial position, one finds the largest smoothing scale which does not collapse at the targeted redshift. This identifies both the smoothing length and the evolution time during which this Lagrangian region is evolved.

As mentioned in section 2.2.4, the equations of large scale structure formation including the expansion of the Universe can be cast in our generic form (2.13), and we use a ballistic approximation which in this context amounts to ZA. A multiscale ZA can therefore directly emerge from our MuScaTS framework. In our MuScaTS approach, we suggest to identify the coherence time with the Press & Schechter (1974) collapse time at a given smoothing length, capped by the total time corresponding to the targeted redshift. This results in two subtle differences compared to existing ZA approaches. First, the collapse time at an intermediate scale is computed on a field which has already been advanced at the targeted redshift on the largest scales. Second, every bandpass filtered density field is evolved separately during its own time scale. This may result in different performances compared to existing methods, and it remains to be checked whether they are better or worse. However, the principle of the resulting method is still very close to the original multiscale ZA approaches.

## 4. Validation for 2D incompressible hydrodynamics

In the previous sections, we presented our MuScaTS framework to generate random realisations of a variety of flows. We provided a possible numerical implementation of the procedure. We also connected our synthesis to previously known generating methods in a few specific contexts.

We now turn to validate the method in the simplest possible setup: incompressible 2D hydrodynamics. Our generative procedure deforms an initial Gaussian field under the influence of the flow. We will assess to what extent our synthesis is able to recover the statistics obtained when the flow evolves from an initial Gaussian field with a prescribed spectrum.

To that end, we run accurate simulations of the flow (using a pseudo-spectral code presented below in section 4.1) starting with a large number (here, 30) of independent initial random realisations of Gaussian fields with a
<figure><img src="image_4.png" /><figcaption>**Fig. 4.1.** Vorticity maps for various cases (see text), all based on the *same* original Gaussian field (displayed on panel (a)). (a) initial conditions, (b-c-d) divorticity syntheses with (b) first order deformation (c) advection then deformation (d) deformation only, (e) reference simulation at *t* = 2 (about 1/3rd turnover time scale) (f) basic synthesis for vorticity (g-h) divorticity syntheses with (g) deformation then advection (h) advection only. All syntheses are shown at the same time *t* than for the simulation and a coherence time based on the local strain time (See subsection 2.4.2).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.2.** Same as figure 4.1, but at time *t* = 6, i.e. at about one turnover time scale.</figcaption></figure>

ant of our syntheses generated from 30 *other* independent Gaussian fields. We will thus get quantitative estimates of the performance of our synthesis for each choice of statistical indicators.

### 4.3. *Classical statistical indicators*

After Kolmogorov (1941) discovered his famous power-spectrum law for 3D incompressible hydrodynamics, Kraichnan (1967) was the first to extend it to 2D incompressible hydrodynamics. He predicted that in forced 2D turbulence the velocity power spectrum should scale as
$k^{-5/3}$ above the injection scale (indirect energy cascade) and as $k^{-3}$ below it (direct entrophy cascade). This was later verified both in soap film experiments (Rutgers 1998; Bruneau & Kellay 2005) and in simulations (Boffetta 2007).

Intermittent corrections to these power-law behaviours were shown to be very weak in the indirect cascade (Boffetta et al. 2000). Very weak non-Gaussianities were detected in the velocity increments probability distribution functions (PDFs), whether measured in experiments (Paret & Tabeling 1998) or simulations (Boffetta et al. 2000; Pasquero & Falkovich 2002). The PDF of vorticity increments on the other hand were predicted to display exponential tails by Falkovich & Lebedev (2011). Indeed, significant deviations from Gaussianity were observed for the vorticity increments in simulations by Tan et al. (2014). As a result, we choose to focus on vorticity increments rather than velocity increments since those will provide an easier measurement of intermittency compared to velocity increments.

Analogues of the exact Kármán Howarth Monin relation in 3D exist for both the indirect cascade of energy and the direct cascade of entrophy in 2D. Eyink (1996) showed that in the indirect cascade, the energy transfer function

$$ \mathcal{F}_u = \langle \delta_{\ell} u_{\parallel}^3 \rangle, \quad (4.1) $$

where brackets denote both space and ensemble averaging, has to be positive and proportional to $\ell$. This is linked to slight asymmetries in the PDFs of the longitudinal velocity increments, with a bias towards positive values for $\delta_{\ell}u_{\parallel}$ at large values of $\delta_{\ell}u_{\parallel}$. On the other hand, in the direct cascade, the entrophy transfer function

$$ \mathcal{F}_w = \langle \delta_{\ell} u_{\parallel} (\delta_{\ell} w)^2 \rangle \quad (4.2) $$

was predicted to be negative and proportional to $\ell$ (Paret & Tabeling 1998), see also Bernard (1999). This second result also involves asymmetries of the longitudinal velocity increments, but correlated to the vorticity increment (which is itself symmetric): it is caused by a bias for negative values of $\delta_{\ell}u_{\parallel}$ conditioned on large vorticity increments. In the following we will examine both the energy and entrophy transfer functions $\mathcal{F}_u$ and $\mathcal{F}_w$ as probes of these second order behaviours of the increments PDFs.

The above results were obtained at steady state in driven turbulence. However, here we deal with decaying turbulence, so none of these results may hold. Nevertheless, we will use these classical indicators, power spectra, increments PDFs and transfer functions, as quantitative measurements to assess how our syntheses perform compared to actual simulations.

### 4.3.1. Power spectra

We compare in panel (a) of figure 4.3 the average power spectra over all 30 realisations (reference simulations in blue and synthesis in red, while initial conditions are in green). Both the slope in the inertial range power-law behaviour and the power spectrum shape in the dissipation range seem equally well recovered. On these figures, the shaded area give the error on the mean as estimated from the dispersion over our 30 realisations. It should be stressed that these error estimates are similar between the simulations and the synthesis. In other words, the synthesis is able to reproduce not only the average, but also the variance of

<figure><img src="image_6.png" /><figcaption>Fig. 4.3. Classical statistics averaged over 30 realisations for (a) power spectrum, (b) vorticity increments (for lags of 1, 4 and 32 pixels), (c) energy transfer function $\mathcal{F}_u$ and (d) entrophy transfer function $\mathcal{F}_w$. Green is for the initial conditions (of the syntheses), red for the MuScaTS syntheses (basic run) and blue the simulations, at time $t = 2$. Shaded areas give the 1-sigma error around the mean over the 30 independent realisations (i.e. mean plus or minus standard deviation divided by $\sqrt{30}$).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4.4. Same as figure 4.3, but at time $t = 6$, i.e. at about one turnover time scale.</figcaption></figure>

the power spectrum across the 30 members of our statistical ensemble. After about one turnover time scale, both the power spectrum and its variance are still well accounted for (see figure 4.4).

### 4.3.2. Increments statistics

We compute $\delta_{\ell}X(r) = X(r + \ell) - X(r)$ as the increment of a quantity $X$ over a lag $\ell$ at position $r$ in the map. We
then examine the statistics of all increments which have their lag $\ell$ in a corona of width 1 pixel around a prescribed norm in pixels. For instance, the PDF of vorticity increments at various lags (in pixels) is displayed in panel (b) of figure 4.3. The PDF in the simulation experience large departures from Gaussianity. If we take as metric the logarithm of the PDF wings, at $t=2$ (i.e. a third of a turnover time), a significant fraction of the non-Gaussianity of vorticity increments is reproduced (synthesis takes us more than half way between Gaussian and initial conditions). Also, the relative dispersion of these results is well estimated, too. At one turnover time scale (see figure 4.4), the discrepancy between the PDFs from simulations and synthesis is larger (simulations have slightly increased their departure from Gaussianity while syntheses are still at about the same place).

### 4.3.3. Transfer functions

As mentioned above, in steady state forced 2D turbulence, the energy transfer function $F_u$ is expected to be positive (indirect cascade) for scales above the forcing scale, while the entrophy transfer function $F_w$ should be negative for scales below. Although we have no forcing here, our decaying turbulence simulations display the same signs for the energy transfer functions. Early 3D generative models such as Rosales & Meneveau (2006) or Chevillard et al. (2011) recognised that it is a real challenge to obtain random fields which can account for the non-zero transfer of energy. Here, our synthesis is able to account not only for the correct sign of these transports but also for their shape, as shown in panel (c) and (d) of figure 4.3. Their magnitude however is underestimated in the syntheses. The agreement is actually better at later times (see figure 4.4) because as time proceeds the total energy in the simulation decays while the syntheses become frozen as the strain time is reached everywhere. Finally, most of the relative dispersion of these quantities is reproduced by the synthesis, at all times.

### 4.4. Scattering transform statistics

Scattering transform statistics are a modern and powerful tool characterising textures (Mallat 2012). In the context of cosmological simulations, Allys et al. (2020) have shown that wavelet phase harmonics are able to predict many other non-Gaussian indicators (increment statistics, Minkowsky functionals, bi-spectrum, etc...).

Here, we present an analysis of our results using wavelet scattering transform (WST) coefficients. The first layer WST coefficients of a field $X$ are defined as:

$$S_1(j) = \langle |X * \psi_{j,\theta}| \rangle_{\theta}, \qquad (4.3)$$

where $\langle \rangle_{\theta}$ denotes the combination of ensemble average, spatial average over the computational domain and average over $\theta$. They constitute a statistical probe similar to the power spectrum, with additional sensitivity to sparsity<fnref n="5" /> introduced by the use of the $L^1$-norm instead of the $L^2$-norm.

We do not show the results, as they are very similar to those obtained for the power spectra (see section 4.3.1).

We then define normalised layer 2 coefficients as:

$$S_2(j_1, j_2, \Delta\theta) = \frac{\langle |X * \psi_{j_1,\theta} * \psi_{j_2,\theta+\Delta\theta}| \rangle_{\theta}}{S_1(j_1)S_1(j_2)}. \qquad (4.4)$$

The layer 2 coefficients ($S_2$) estimate the coupling between two angles $\theta_1$ and $\theta_2 = \theta_1 + \Delta\theta$ and two scales $\ell_1 = \ell_0 2^{j_1}$ and $\ell_2 = \ell_0 2^{j_2}$ labelled by $j_1$ and $j_2$, with $\ell_0$ the smallest scale considered (note that by convention scales are labelled by increasing order in WST coefficients). In our case, thanks to isotropy, the coefficients depend only on the difference $\Delta\theta = \theta_1 - \theta_2$ and we have checked that indeed there is very little dispersion of the coefficients at fixed $\Delta\theta$, $j_1$ and $j_2$ when $\theta_1$ varies. The shaded areas on figure 4.5 actually reflect the dispersion resulting from the spatial average over the vorticity fields, the ensemble average over our 30 realisations and from the angular average over $\theta_1$. These shaded area are so thin that they look like solid lines. With our 1024 pixels resolution, we chose 9 dyadic scales along with 16 angles. We plot on figure 4.5 the normalised $S_2$ coefficients at $t=2$ (i.e. about one third of a turnover time).

This powerful tool is much more stringent than the increments and probes much more accurately the defects of our synthesis. The synthesis correctly captures the relative variation of the coefficients with $\Delta\theta$ but fails at reproducing the proper slope of the power-law decay with $\Delta j = (j_2-j_1)$ (although the decay is still somewhat power-law like). At later times (not shown here), the angular variation is also lost by the synthesis. We stress here that although the classical indicators shown in the previous subsection 4.3 were poorly sensitive to the texture differences, the WST coefficients provide a much better quantitative assessment of the discrepancies detected by an eye inspection of figure 4.1. Indeed, WST statistics are known to produce remarkable syntheses (Allys et al. 2019) and constitute a very practical tool to partly assess quantitatively what our eyes are able to probe. First, in this metric, the distance between simulations and syntheses relative to the distance between initial conditions and final result is much larger compared to classical indicators. The synthesis works us a much smaller path towards the actual non-Gaussianities in the eyes of scattering statistics. Second, the discrepancy has a lot more statistical significance: the very small amount of dispersion for these layer 2 coefficients greatly helps disentangling the results between the simulations and our generative models.

We also attempted to find good reduced models for the variation of the layer two coefficients with respect to the relative angle $\Delta\theta$, in the spirit of the reduced WST presented in Allys et al. (2019). We found that the generic shape $S_2(i, j, \Delta\theta) \propto \alpha_{i-j} + \beta_{i-j} \cos(2\Delta\theta) + \gamma_{i-j} \cos(4\Delta\theta)$ provides a good fit to the data, but opted not to show the resulting adjustments to prevent overloading the paper.

## 5. Discussion

We provide in this section a brief account of the parameter space investigation we have conducted, and discuss a few of the shortcomings of our method.

### 5.1. Filter shape

Although all the results of the present paper are computed with a bank of filters based on our cosine filters, we have

<footnote n="5">⁵ The sparsity of a signal refers to its property to have a large fraction of its values being close to 0. It can be measured by comparing its $L^1$ and $L^2$ norms. Such measure can be used to characterize the intermittent nature of turbulent flows (Bruna et al. 2015).</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 4.5.** Mean over 30 realisations of layer 2 WST coefficients $S_2(j_1, j_2, \Delta\theta)$ (see text) for initial conditions (green), the reference simulation (blue) at about a third of a turnover time ($t=2$) and MuScaTS basic syntheses over a time $t=2$ (red). The main length scale runs from its index $j_1=0$ (smallest) to $j_1=J-3=7$. The secondary scale $j_2$ is labelled by $\Delta j = j_2 - j_1$. The relative angle $\Delta\theta$ runs from 0 to $\pi$ from left to right in each broken line of fixed $j_1$ and $\Delta j$. The width of the lines indicates the standard deviation of this mean over our 30 realisations, the computational domain and the absolute angle $\theta_1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 4.6.** Same as figure 4.5 but with a 4 steps bootstrap (see text of subsection 5.4, and see also figure 5.1).</figcaption></figure>

also investigated $B_3$-spline filters (see appendix B). There are very slight differences between the two, hardly significant statistically and we do not show them here. We chose to present the cosine filters because they can cover the Fourier space at a slightly lower cost. In a paper which inspired us for this work, Lewalle (2010) showed how Gaussian based wavelets (Mexican hats) greatly simplify the analytical expressions for the filtered equations. Mexican hats might hence be best suited for theoretical work based on our MuScaTS framework (their weak selective power may hinder the quality of the syntheses, though).

## 5.2. Divorticity

As shown above (see section 2.2.2), the evolution equations for 2D incompressible hydrodynamics for the divorticity also suit our generic formalism (2.13), with both an advection and a deformation term present, very similar to the structure obtained for the 3D vorticity evolution equations.

We find that our method performs very slightly better when we base it on the vorticity evolution equation rather than on the divorticity equation. This might be due to the presence of two terms (advection and deformation) in the divorticity evolution equation. Also, it may be due to the fact that the divorticity power spectrum is even more biased toward the small scales than the vorticity spectrum. Finally, the divorticity remapping introduces spurious non zero divergence in this divergence free field, while the vorticity field does not suffer from this, as it is a scalar in 2D.

We have investigated several variants in the implementation of our method applied to divorticity in 2D. We swapped the order between advection and deformation, without finding significant change (compare panels (c) and (g) on figure 4.1). We tested the full matrix exponential implementation of equation (2.31) and its linearised version in $\tau_\ell$, with a slight preference for the latter implementation (compare panel (b) and (c) on figure 4.1). This might be due to the better consistency between the ballistic approximation for the advection term and a first order approximation in time for the deformation. Classical indicators are almost blind to whether vorticity or divorticity is advanced, with a slight preference for the latter method for the vorticity increments while transfer functions are better reproduced in the former case. Vorticity maps and WST coefficients give a slight preference to the treatment of vorticity rather than divorticity.

The above choices retaining both advection and deformation for the implementation of MuScaTS on the divorticity were not critical. But we found a huge difference when either the advection or the deformation term were omitted, as panels (d) and (h) in figure 4.1 can testify. Even the classical indicators suffer from the loss of either of these terms. This is a good sign for future implementation of our method in the 3D case which should perform even better when compared to existing methods such as Rosales & Meneveau (2006) or Chevillard et al. (2011) and their extensions.
## 5.3. Choosing the coherence time

Most results presented here have used a coherence time based on a smooth minimum between the local strain time (see our last choice in section 2.4.2) and the global target time. However, we have also investigated the usage of a constant global target time, with significantly less accurate results. We also tested a coherence time uniform in space, but dependent on scale as a power-law (which resulted in significantly less accuracy). And finally we tested the other options based on the local shell crossing time or the local stretch time and found the method appeared slightly less accurate in view of the statistical indicators we tested. Perhaps future investigations in compressible or 3D applications will help us better understand which version of the coherence time is most appropriate.

## 5.4. Bootstrapping

Since our method performs better during a smaller fractions of the initial turnover time scale, it is natural to try and start a new MuScaTS step from the resulting snapshot obtained instead of using the initial Gaussian conditions. We have hence split the integration step in two or four smaller successive MuScaTS steps to test the convergence of the method, a process we denote as bootstrapping. The resulting improvement on the vorticity maps is striking, as illustrated on figure 5.1.

Bootstrapping is a poor man's way of probing higher order schemes, hence it augurs well for the performance of future higher order methods. It is also a sign that for steady state turbulence we need better statistics than Gaussian conditions to randomly reshuffle our stochastic terms: including some non-Gaussianity (that of the previous step, when bootstrapping) in these pseudo initial conditions makes an important improvement. This resonates with the findings of Lübke et al. (2024) who found realistic coherent MHD structures when they started the method of Subedi et al. (2014) with a non Gaussian initial field. Finally, the realism of the resulting solution shows that despite the strong approximations we made in deriving our filtered equations, we did not lose too much of the essential physics.

## 5.5. Power spectrum shape

All results shown in this paper are based on an initial velocity power spectrum proportional to $k^{-3}$. However, we have also tested a less steep power spectrum, scaling as $k^{-5/3}$, corresponding to the indirect energy cascade. This last spectrum is less dominated by the large scales, and in particular the vorticity field and more importantly the stretching field is then dominated by the small scales. This may explain why we have found our statistical comparison is significantly less good in this shallower spectrum case (not shown here).

Indeed, in the current application of our method we apply a single sweep through scales from the largest to the smallest, implicitly considering scale causality in a unique direction. In future work, we will investigate the possibility to implement corrections to account for the feedback of the small scales to the largest (see further discussion in subsection 5.6).

<figure><img src="image_10.png" /><figcaption>**Fig. 5.1.** Vorticity maps for one realisation of our basic case at $t = 2$ for (a) the simulation (b) same as panel (d) of figure 4.1 (c) a two-steps ($\Delta t = 1$) bootstrap and (d) a four steps ($\Delta t = 0.5$) bootstrap (see text).</figcaption></figure>

In this paper, we have investigated only decaying turbulence. In principle, it is easy to incorporate a forcing term in our partial differential equations. However, a shortcoming of our method is that we have to introduce a power spectrum for the "initial conditions" $W_l^0$. To estimate this power spectrum, we could resort to theoretical inputs, we could bootstrap our method until we get a steady power spectrum, or there may exist analytical considerations on energy conservation to provide the power spectrum shape from the filtered equations alone.

## 5.6. Small-scale terms

We have currently bundled all small scale terms of our scale filtered equation (2.25) into a "stochastic small scales" term with which we deal with in an extremely expedient manner. A large legacy of literature tackles coarse-graining or Reynolds averaging, which leads to convective terms, turbulent diffusion, small scale dynamo, eddy viscosity which all result from the filtering of the small scales terms. Estimation of these terms requires some sort of closure. Most methods adopt some sort of linear approximation to compute correlations between small scale terms, such as first order smoothing approximation (FOSA), also known as second order correlation approximation (SOCA), as introduced in Krause & Raedler (1980), or mixing length theory approaches (MLT, see Lesaffre et al. 2013; Jermyn et al. 2018, and references therein). Test-field techniques have sometimes been used to validate these approaches in numerical simulations (Brandenburg et al. 2008; Käpylä, P. J. et al. 2020).

Here, we simply have short-circuited these considerations by randomly reshuffling the state variables of the gas over one coherence time-scale. However, once we have a
first guess for the properties of our flow thanks to a first
sweep from large to small scales, we could imagine refin-
ing the scale filtered evolution equations thanks to these
terms which can be readily computed on that first guess.
This would recover some of the causal influence of the small
scales onto the large ones that our method currently com-
pletely neglects.

It seems paradoxical that we can still get correlations
such as energy transfers from small to large scales (as in our
2D example above where we get a positive energy transfer
function). This means that in our 2D HD synthesis, large
scales conspire to arrange velocities in such a way that there
is a net energy transfer towards larger scales. Only, due to
the single direction of our sweep from large to small scales,
that energy transfer is not imprinted on the evolution equa-
tions. In our current investigation, it does not seem to be
critical. However, there may be applications where it could
be important, in particular when the spectrum is shallower,
but also in cases when small-scales dynamos are at play
(i.e. when some small scales instabilities lead to exponen-
tial growth of larger scale averaged quantities).

We finally point here to a slight inconsistency in our for-
malism: the filtered field $\tilde{W}_\ell$ at scale $\ell$ may acquire power
at other scales from its approximate filtered evolution equa-
tion (2.28). This may happen through convergence and di-
vergence of trajectories in the ballistic step, for example,
but any non linear term will generate scales outside the ini-
tial range of the filter at scale $\ell$. At this stage, we are not
sure whether it is a potential drawback or an actual advan-
tage of our method (smaller scales build more efficiently,
and in a natural way). In practice, however, we checked
that the power generated at a scale $\ell$ did not overlap too
much on neighbouring scales.

5.7. Extending the scope of our method

We delve back here on the scope of our method. The generic
form of the evolution equations (2.13) allows already for
quite a few situations which were not dealt with before with
generative methods. We chose this form because it behaves
well under filtering and because it readily expresses advec-
tion and deformation, which clearly are useful concepts in
fluid mechanics. This also allows us to easily link to pre-
vious works which had separately focused on either one of
these two terms. However, it may not be the most general
form which is amenable to our hierarchical filtering. To start
with, we could accept various components of $\mathbf{W}$ to be ad-
vected individually with different speeds (as in the Elsasser
formulation of incompressible MHD where $Z^+$ is advected
by $Z^-$ and conversely) as long as these speeds remain linear
in the state vector $\mathbf{W}$. There may also exist more general
forms which survive filtering and the conservative form is
certainly one which should attract our attention.

In addition, if we now stick to our generic form (2.13), it
is also not always easy to transform common partial differ-
ential equations towards it. For instance in the simple case
of incompressible HD, we could not find a way to write down
the velocity evolution equation (Navier-Stokes) in this form
(there might be one that we could not see, or there might
be a way to still apply the rest of the formalism to it, we
simply have not found it yet).

Finally, even in cases where we could not strictly fit an
evolution equation to our framework, as in the case of com-

pressible isothermal MHD, it may be that hiding isolated
problematic terms into our stochastic reshuffling may not
be too detrimental to the accuracy of our method.

6. Summary and prospects

We have investigated a set of few simple ideas applied to
a large class of partial differential equations: hierarchical
smoothing, first order evolution from random intermediate
initial conditions during a scale and space dependent co-
herence time. With this procedure, we efficiently generate
fields with realistic statistical properties compared to the
results of the full time integration, at least in the simplest
case of 2D incompressible HD.

Previous generative methods applied to 3D HD or MHD
have focused either on advection (Rosales & Meneveau
2006; Subedi et al. 2014; Lübke et al. 2024) or on deforma-
tion (Chevillard et al. 2011; Durrive et al. 2020, 2022). The
present approach retains every physical term in the equa-
tions, including diffusion, and provides a general framework
suitable to a large class of partial differential equations. In
the 2D incompressible HD case, for the evolution of di-
vorticity, we have shown that it is critical to include both
advection and deformation at the same time: this brings
hope that our method when applied to 3D HD or MHD
may significantly outperform those existing techniques.

Classical statistical indicators remain relatively well re-
produced by our approximations during a significant frac-
tion of the turnover time (about one third). A few coherent
structures such as vorticity clumps with spiralling arms are
also reproduced, with the proper winding sense of rotation
according to the vorticity sign. However, more stringent
probes of the texture such as scattering statistics are not
fooled and show discrepancy between simulations and syn-
theses. Yet a limited amount of bootstrap is able to bridge
that gap: higher order integration schemes are likely to yield
syntheses of much better quality at a smaller computational
cost than bootstrap.

With the advent of deep learning, very efficient gener-
ative methods relying on deep learning architectures have
emerged, as for instance diffusion models (see Genuist et al.
2024, and references therein). However, these techniques
suffer from two strong shortcomings: the generative pro-
cedure is often difficult to apprehend, and usually needs
a training set of considerable volume, which in this con-
text means that a lot of full fledged numerical simulations
have to be run prior to training the generative model. Our
method, on the contrary, is directly fed with the physics
contained in the partial differential equations and provides
generated fields at the cost of a handful of simulation steps.
Incidentally, our method could potentially be used to build
large approximate training sets instead of resorting to ac-
tual simulations, or to complement a smaller set of accu-
rate simulations with a large number of MuScaTS generated
snapshots. Other generative methods as maximum entropy
models conditioned by scattering transform statistics, are
also being developed (Bruna & Mallat 2019; Zhang & Mal-
lat 2021; Cheng et al. 2024), which require a very small
training set (only a handful simulations and sometimes even
a single snapshot). However, these methods generate a sig-
nificant extra numerical cost and offer less control of the
distance to the physical solutions compared to the method
we propose.
The original purpose of the synthesis methods as setup by Rosales & Meneveau (2006) was to reduce the time spent by simulations in the initial transient phase from arbitrary initial conditions towards fully developed turbulence. Incidentally, this reduces the carbon footprint of such simulations. Of course, the present model can still be used to the same aim. Here, we propose a range of other future applications. The low generating cost could be used with Monte Carlo Markov Chain techniques to adjust a synthetic field to some observational constraints. For example, one could imagine generating 3D synthetic fields with prescribed projections: to what extent can we recover the underlying 3D structure behind plane-of-sky astrophysical images in such a way ? Similar methods may be used to guess a turbulent field outside regions where it is measured, which could be useful in weather forecasting applications where reliable measurements are sometimes sparsely located. Finally, we hope our simplified models can be used as theoretical basis to understand existing or predict new statistical laws of turbulent flows.

## References

Allys, E., Levrier, F., Zhang, S., et al. 2019, A&A, 629, A115

Allys, E., Marchand, T., Cardoso, J. F., et al. 2020, Phys. Rev. D, 102, 103506

Banerjee, S. & Galtier, S. 2013, Phys. Rev. E, 87, 013019

Banerjee, S. & Galtier, S. 2014, Journal of Fluid Mechanics, 742, 230

Banerjee, S. & Galtier, S. 2016, Phys. Rev. E, 93, 033120

Bernard, D. 1999, Phys. Rev. E, 60, 6184

Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rep., 367, 1

Boffetta, G. 2007, Journal of Fluid Mechanics, 589, 253

Boffetta, G., Celani, A., & Vergassola, M. 2000, Phys. Rev. E, 61, R29

Bond, J. R. & Myers, S. T. 1996, ApJS, 103, 1

Brandenburg, A., Rädler, K. H., & Schrinner, M. 2008, A&A, 482, 739

Bruna, J. & Mallat, S. 2019, Mathematical Statistics and Learning, 1, 257

Bruna, J., Mallat, S., Bacry, E., & Muzy, J.-F. 2015, The Annals of Statistics, 43, 323

Bruneau, C. H. & Kellay, H. 2005, Phys. Rev. E, 71, 046305

Cheng, S., Morel, R., Allys, E., Ménard, B., & Mallat, S. 2024, PNAS nexus, 3, pgae103

Chevillard, L. 2015, Accreditation to supervise research, ENS Lyon

Chevillard, L. & Meneveau, C. 2006, Phys. Rev. Lett., 97, 174501

Chevillard, L., Robert, R., & Vargas, V. 2011, in Journal of Physics Conference Series, Vol. 318, Journal of Physics Conference Series, 042002

Coles, P., Melott, A. L., & Shandarin, S. F. 1993, MNRAS, 260, 765

de Karman, T. & Howarth, L. 1938, Proceedings of the Royal Society of London Series A, 164, 192

Durrive, J.-B., Changmai, M., Keppens, R., et al. 2022, Phys. Rev. E, 106, 025307

Durrive, J.-B., Lesaffre, P., & Ferrière, K. 2020, MNRAS, 496, 3015

Eyink, G. L. 1996, Physica D Nonlinear Phenomena, 91, 97

Falkovich, G. & Lebedev, V. 2011, Phys. Rev. E, 83, 045301

Frisch, U. 1996, Turbulence (Cambridge University Press), 310

Galtier, S. & Banerjee, S. 2011, Phys. Rev. Letters, 107, 134501

Genuist, W., Savin, É., Gatti, F., & Clouteau, D. 2024, working paper or preprint

Goldreich, P. & Sridhar, S. 1995, The Astrophysical Journal, 438, 763

Grauer, R., Krug, J., & Marliani, C. 1994, Physics Letters A, 195, 335

Gurbatov, S. N., Saichev, A. I., & Shandarin, S. F. 1989, MNRAS, 236, 385

Haller, G. 2021, Journal of Fluid Mechanics, 908, A25

Iroshnikov, P. S. 1963, AZh, 40, 742

Jermyn, A. S., Lesaffre, P., Tout, C. A., & Chitre, S. M. 2018, MNRAS, 476, 646

Kida, S. 1985, Journal of the Physical Society of Japan, 54, 2840

Kolmogorov, A. N. 1941, Akademiia Nauk SSSR Doklady, 32, 16

Kolmogorov, A. N. 1962, Journal of Fluid Mechanics, 13, 82

Kraichnan, R. H. 1965, Physics of Fluids, 8, 1385

Kraichnan, R. H. 1967, Physics of Fluids, 10, 1417

Krause, F. & Raedler, K. H. 1980, Mean-field magnetohydrodynamics and dynamo theory (Oxford: Pergamon Press)

Käpylä, P. J., Rheinhardt, M., Brandenburg, A., & Käpylä, M. J. 2020, A&A, 636, A93

Lesaffre, P., Chitre, S. M., Potter, A. T., & Tout, C. A. 2013, MNRAS, 431, 2200

Lewalle, J. 2010, Physica D Nonlinear Phenomena, 239, 1232

Lübke, J., Effenberger, F., Wilbert, M., Fichtner, H., & Grauer, R. 2024, EPL (Europhysics Letters), 146, 43001

Majda, A. J. & Bertozzi, A. L. 2001, Vorticity and Incompressible Flow, Cambridge Texts in Applied Mathematics (Cambridge University Press)

Malara, F., Di Mare, F., Nigro, G., & Sorriso-Valvo, L. 2016, Phys. Rev. E, 94, 053109

Mandelbrot, B. B. & van Ness, J. W. 1968, SIAM Review, 10, 422

Momferratos, G., Lesaffre, P., Falgarone, E., & Pineau des Forêts, G. 2014, Monthly Notices of the Royal Astronomical Society, 443, 86

Monaco, P., Theuns, T., Taffoni, G., et al. 2002, ApJ, 564, 8

Monin, A. S. 1959, Soviet Physics Doklady, 4, 271

Okubo, A. 1970, Deep Sea Research A, 17, 445

Paret, J. & Tabeling, P. 1998, Physics of Fluids, 10, 3126

Pasquero, C. & Falkovich, G. 2002, Phys. Rev. E, 65, 056305

Pereira, R. M., Moriconi, L., & Chevillard, L. 2018, Journal of Fluid Mechanics, 839, 430

Politano, H. & Pouquet, A. 1995, Phys. Rev. E, 52, 636

Press, W. H. & Schechter, P. 1974, ApJ, 187, 425

Richard, T., Lesaffre, P., Falgarone, E., & Lehmann, A. 2022, A&A, 664, A193

Rosales, C. & Meneveau, C. 2006, Physics of Fluids, 18, 075104

Rosales, C. & Meneveau, C. 2008, Phys. Rev. E, 78, 016313

Rutgers, M. A. 1998, Phys. Rev. Lett., 81, 2244

Schekochihin, A. A. 2022, Journal of Plasma Physics, 88, 155880501

She, Z.-S. & Leveque, E. 1994, Phys. Rev. Letters, 78, 336

Shivamoggi, B. K., van Heijst, G. J. F., & Kamp, L. P. J. 2024, Journal of Plasma Physics, 90, 905900214

Sidje, R. B. 1998, ACM Trans. Math. Softw., 24, 130

Starck, J.-L. & Murtagh, F. 2002, Astronomical image and data analysis (Springer)

Stein, G., Alvarez, M. A., & Bond, J. R. 2019, MNRAS, 483, 2236

Subedi, P., Chhiber, R., Tessein, J. A., Wan, M., & Matthaeus, W. H. 2014, ApJ, 796, 97

Tan, H. S., Huang, Y. X., & Meng, J. 2014, Physics of Fluids, 26, 015106

Weiss, J. 1991, Physica D Nonlinear Phenomena, 48, 273

Zhang, S. & Mallat, S. 2021, Applied and Computational Harmonic Analysis, 53, 199

## Acknowledgements

We express gratitude to the referee for having read our paper so carefully and for his constructive comments. The authors acknowledge Interstellar Institute's programs "With Two Eyes" and "I16" and the Paris-Saclay University's Institut Pascal for hosting discussions that nourished the development of the ideas behind this work. PL acknowledges enlightening discussions with Marie Farge who accepted to read through an early version of the paper. PL acknowledges support from the European Research Council, under the European Community's Seventh framework Programme, through the Advanced Grant MIST (FP7/2017-2022, No 742719). Computations were carried out on the totoro machine financed by the same ERC grant and hosted at the meso-center mesoPSL.
## Appendix A: Reconstruction formula

In this appendix, we justify the reconstruction formula (2.9).

Let us take the Fourier transform of the right hand side of (2.9), considering it at a fixed given non zero $\mathbf{k}$. By linearity of the Fourier transform (interchanging the integral of the Fourier definition with the integral over scales $\ell$),

$$ \left( \int_{-\infty}^{+\infty} d\ln \ell f_{\ell} \right) (\mathbf{k}) = \int_{-\infty}^{+\infty} d\ln \ell \hat{f}_{\ell}(\mathbf{k}). \quad (A.1) $$

Moreover, by the convolution theorem (the Fourier transform of the convolution of two functions is the product of their Fourier transforms) relation (2.1) in Fourier space, combined with (2.3) and (2.4), lets us rewrite the integrand on the right hand side of (A.1) as

$$ \hat{f}_{\ell}(\mathbf{k}) = \hat{f}(\mathbf{k})\hat{\varphi}_{\ell}(\mathbf{k}) = \hat{f}(\mathbf{k})\Phi(|\mathbf{k}|\ell). \quad (A.2) $$

Then, as $\hat{f}(\mathbf{k})$ is independent of $\ell$ it can be pulled out of the integral, and using the change of variable $s = |\mathbf{k}|\ell$ (possible since $\mathbf{k}$ is fixed and assumed to be non zero), we get that

$$ \left( \int_{-\infty}^{+\infty} d\ln \ell f_{\ell} \right) (\mathbf{k}) = \hat{f}(\mathbf{k}) \int_{-\infty}^{+\infty} d\ln s \Phi(s) = C_{\Phi} \hat{f}(\mathbf{k}), \quad (A.3) $$

where definition (2.6) has been used in the last step. Given our normalisation (2.7), relation (A.3) demonstrates that both sides of (2.9) have the same Fourier transform for all non zero $\mathbf{k}$. For the special $\mathbf{k} = 0$ case, this statement remains true provided $\hat{f}(0) = 0$ (i.e. $f$ has a zero mean) because all filters have zero means. Finally, by inverse Fourier transform, we obtain the reconstruction formula (2.9).

## Appendix B: Two example filters

In the present application, we tested two sets of filters. The first one is based on a cosine window in the logarithm of the scale. We first define the bump-like function

$$ C_2(s) = \cos^2\left(\frac{\pi}{2}s\right), \quad (B.1) $$

for $s$ in $[-1, 1]$, and zero elsewhere. We then define our list of filters as

$$ \hat{\varphi}_j^C(\mathbf{k}) = C_2(\ln(|\mathbf{k}|/k_j)/\ln[\lambda]), \quad (B.2) $$

where $\lambda$ is the scale ratio between two adjacent scales. Our normalisation condition requires $J = J_C = \lfloor \ln(N/2)/\ln \lambda \rfloor$ (where $\lfloor \cdots \rfloor$ denotes the integral part or floor function) in order to have $\sum_{j=0}^J \hat{\varphi}_j^C(\mathbf{k}) = 1$ for all wavenumbers in the computational domain. Note that the maximum wavenumber represented in our Fourier space is $k_{\max}/k_J = \sqrt{3}N/2$ for a 3D grid, or $k_{\max}/k_J = \sqrt{2}N/2$ for a 2D grid, and since the non-linearities in our method are quadratic, we use the two-third truncation rule in Fourier.

Our second set of filters is based on the box spline of 3rd order (often used in astrophysical applications due to its implementation efficiency, see appendix A of Starck &

<figure><img src="image_11.png" /><figcaption>**Fig. B.1.** Series of cosine (left) and box spline (right) filters for a cubic periodic domain of size $N = 512$ pixels. Their total sum is indicated as well as the wavenumber boundaries for a 3D domain ($k_{\min} = 1$, $k_{\max} = \sqrt{3}N/2$). Colours from red to blue stand for filter indices 0 to $J$. This figure illustrates how many filters are necessary in order for the reconstruction formula to be valid. In practice however, diffusion damps the higher $k$ end of the spectrum, and convergence is reached for a slightly lower number of filters.</figcaption></figure>

Murtagh 2002, for example). Consider the $B_3$-spline cubic function (a box convolved three times with itself)

$$ B_3(s) = \frac{1}{12} (|s-2|^3 - 4|s-1|^3 + 6|s|^3 - 4|s+1|^3 + |s+2|^3). \quad (B.3) $$

$B_3(s)$ is non zero only within the interval $[-2, 2]$. It makes a smooth transition between 2/3 at $s = 0$ and 0 at $s = 2$. We can use this smooth step property to build a smooth logarithmic window:

$$ \hat{\varphi}_j^S(k) = \frac{3}{2} [B_3(|k|/k_{j+1}) - B_3(|k|/k_j)], \quad (B.4) $$

where the 3/2 coefficient is chosen to comply with our normalisation condition. The sum $\sum_{j=0}^J \hat{\varphi}_j^S = \frac{3}{2}B_3(|k|/k_{J+1})$ has limit 1 when $|k|/k_{J+1}$ goes to infinity, i.e. when $J$ is large enough for $k_{J+1}$ to be well above the maximum effective $k$ reached by the wavenumbers ($k_{\max}/k_J = \sqrt{3}N/2$ for a 3D grid, or $k_{\max}/k_J = \sqrt{2}N/2$ for a 2D grid). The normalisation condition is hence only approximately realised in this case.

Figure B.1 displays series of filters for these two choices of mother filter, along with the normalisation condition. The cosine filters are slightly narrower in scales compared to the splines. The spline filters do not match the condition exactly but only to the limit where we take a large number of filters. However, in practice, only the small scales will suffer from lack of power, which is in line with diffusion processes anyway, so it turns out that the use of $J_S = J_C + 1 = \lfloor \ln(N/2)/\ln \lambda \rfloor + 1$ for the spline filters appears quite well converged, with only about 14% of the total power missing compared to the limiting case (estimated from using $J = J_C + 10$), and only 1% missing at $J = J_C + 2$, which we adopted in the present work.