Jean-Baptiste Durrive

Pierre Lesaffre

Katia Ferrière

Magnetic fields from Multiplicative Chaos

Keywords:

An analytical model for fully developed three-dimensional incompressible turbulence was recently proposed in the hydrodynamics community, based on the concept of multiplicative chaos. It consists of a random field represented by means of a stochastic integral, which, with only a few parameters, shares many properties with experimental and numerical turbulence, including in particular energy transfer through scales (the cascade) and intermittency (non-Gaussianity) which is most conveniently controlled with a single parameter. Here, we propose three models extending this approach to MHD turbulence. Our formulae provide physically motivated 3D models of a turbulent velocity field and magnetic field coupled together. Besides its theoretical value, this work is meant to provide a tool for observers: a dozen of physically meaningful free parameters enter the description, which is useful to characterize astrophysical data.

INTRODUCTION

Magnetic fields are ubiquitous in the Universe, with a stunningly wide range of strengths and coherence lengths, from 10 15 G in magnetars to 10 -6 G in galaxy clusters, or even weaker at even larger scales. Indeed, while historically astrophysical plasmas pertained to solar and stellar physics, magnetism nowadays is also the focus of physicists studying much larger scales. For example, a multitude of highresolution observations indicate that the interstellar medium of our galaxy is a magnetized, multi-phase, highly turbulent and intermittent medium (e.g. Falgarone et al. 1991;Miville-Deschênes et al. 2003;Elmegreen & Scalo 2004;Planck Collaboration et al. 2018). Extra-galactic fields, up to cosmological scales (Mpc size filaments, and, as claimed by an increasing number of groups, inside cosmic voids) are about to be routinely observed as well, thanks for example to the radio-telescopes LOFAR and the SKA. For reviews focused on such fields, see for instance Widrow (2002); Ryu et al. (2012); Subramanian (2019). The bottom line is that, in order to understand the dynamics of astrophysical media, one needs to be able to model MHD turbulence. For astrophysically-oriented reviews on MHD turbulence see e.g. [START_REF] Schekochihin | Turbulence and Magnetic Fields in Astrophysical Plasmas[END_REF]; Brandenburg & Lazarian (2013); Ferrière (2019); Tobias (2019).

Astrophysical systems are fascinating for the same reason that they are so difficult to analyze: what we observe E-mail:jdurrive@irap.omp.eu is the result of extremely rich dynamics, consisting of the intertwining of numerous processes governed by a variety of timescales, operating simultaneously at interdependent lengthscales. For this reason, solving consistently all the equations relevant for the interstellar medium for instance is not reachable analytically, and the relevant tools for this challenge are numerical simulations. Here our approach is to derive analytical expressions which, while capturing as much of the underlying physics as possible, remain simple enough to be discussed intuitively. For example, in our model intermittency is controlled by a single parameter, while in principle it is the result of complex dynamics, namely turbulent dynamo action and random shock compression. Hence, the description introduced here is phenomenological, and may be seen as a complementary tool to numerical simulations. This work was carried out with the intention of providing a tool to describe effectively observations, grasping a fair amount of physics with only a dozen of physically meaningful free parameters. Indeed, our formulae are random processes that contain parameters (and functions) corresponding to: integral scales (at which energy is injected), dissipative scales, the large-scale shape of some mean field, the nature of the process (positively or negatively correlated increments), energy cascades (transfer through scales) and intermittency. Another difficult property to take into account in modeling MHD turbulence, is the coupling of velocity and magnetic fields. Our formulae do account for this in a realistic manner, to a certain extent, as they are built from equations embodying the dynamics. To be able to model such features opens up many applications for astrophysics. For instance, in the interstellar medium the propagation of cosmic rays is very sensitive to the intermittency of the magnetic field in which they move (Shukurov et al. 2017). Nowadays, to explore the correlation between magnetic fields and cosmic rays, one has to run costly numerical simulations (Seta et al. 2018), while with the present model synthetic data may be obtained orders of magnitude faster and controlling the degree of intermittency with a single parameter.

In addition to being a vast subject, astrophysics inherits the intrinsic complexity of the physics itself. Indeed, it is well-known that turbulence per se, even in the incompressible hydrodynamical case, brings serious theoretical challenges, leaving us today with many long-lasting open questions. Exact solutions of the Navier-Stokes equations were found only in quite limited cases, and only a few statistical laws were derived rigorously. Meanwhile, various phenomenologies for MHD turbulence have been suggested in the literature. The study of Alfvénic turbulence has a long history, with pioneering work by Iroshnikov (1963) and Kraichnan (1965) who generalized Kolmogorov's isotropic model. A quantitative picture of the properties of anisotropic cascade emerged with the seminal papers by [START_REF] Frisch | Turbulence. The legacy of A[END_REF]; Ng & Bhattacharjee (1996), at the heart of which lies the concept of critical balance. This picture keeps being refined, with additional ingredients such as dynamical alignment (Boldyrev 2006;Mallet & Schekochihin 2017) and thanks to important on-going efforts to tackle the problem by means of numerical simulations; for reviews see for example Verma (2004); [START_REF] Pouquet | A Review of the Possible Role of Constraints in MHD Turbulence[END_REF]; Rincon (2019); [START_REF] Moffatt | Self-Exciting Fluid Dynamos[END_REF]. Our work relies on results derived by members of the hydrodynamics community, who themselves adapted ideas developed in the mathematics community. The latter physicists also consolidated their formal results by benefiting from close interactions with laboratory experiments. To the best of our knowledge, their approach has not been extended to MHD1 , which is the purpose of the present paper, and has not been mentioned in the astrophysics literature yet either.

More specifically, in a series of papers (essentially Robert & Vargas 2008;Chevillard et al. 2010;[START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]Pereira et al. 2016) a method was proposed to construct analytical expressions of 3D random vectorial fields representing the local structure of homogeneous, isotropic, stationary (fully developed) and incompressible turbulence. The random nature of turbulence is dealt with by means of stochastic calculus, i.e. employing a statistical description of turbulence, modeling fluids in terms of random functions [START_REF] Panchev | Random Functions and Turbulence[END_REF][START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF]. The gist of the method is the following. Starting from a vectorial Gaussian random field, features characteristic of turbulent flows are included one by one, through a series of cunning modifications. Non-Gaussianity, i.e. intermittency, is brought in the model through the (pleasantly named) Gaussian 'multiplicative chaos' introduced by Kahane (1985). Applications of multiplicative chaos go well beyond the theory of turbulence, as it is used not only in mathematics, but also in finance (Duchon et al. 2012) or quantum gravity (Knizhnik et al. 1988;Kupiainen et al. 2018) for example (see Rhodes & Vargas 2014, for a recent review). Then to obtain energy cascades, the key is to extend the scalar theory of Kahane (1985) to matrices. This was first done by Chevillard et al. (2013) in view of applications to turbulence. In fact, one of the important successes of the phenomenology outlined here is that this velocity field is the first stochastic process proposed in the literature that is able to predict a non-vanishing mean energy transfer across scales (Pereira et al. 2016). The long range correlated nature of turbulence is incorporated heuristically using a non linear transformation inspired by the 'recent fluid deformation approximation' proposed by [START_REF] Chevillard | [END_REF]. Finally, note that in the aforementioned literature, both the Lagrangian and Eulerian descriptions are adopted. Here we stick to the Eulerian framework.

The spirit and organization of the paper is the following. In the line of Chevillard et al. (2010), our approach is practical in the sense that we reduce mathematically involved material to intuitive considerations, such that hopefully the reader will feel that the model is built as blocks assembled gradually through a series of intuitive steps. To this end, we begin with the simple case of a scalar field in section 2, which introduces basic ideas. In section 3, the more involved case of a velocity field is treated, adapting previous works to our purpose. This allows us to prepare the ground for our generalization, makes the paper self-consistent, and advocates for digging in this direction. Indeed, at the end of this section, we emphasize how well the results of this procedure agree with direct numerical simulations and experimental data. Next, in section 4, we present our extension to magnetic fields which should appear natural, having carefully detailed the hydrodynamical case. In section 5, we illustrate some properties of the formulae of the paper through numerical examples. Finally, we conclude by giving hints of the prospects of this work.

General comments on the notations: (i) To lighten the paper, we use everywhere the shorthand notation

r ≡ x -y and r ≡ |r|, (1)
as most expressions are convolutions, with r corresponding to the distance from a given position x when integrating with respect to y. (ii) We use subscripts s, v and m for quantities associated with scalar fields, velocity fields and magnetic fields respectively. (iii) The symbol tilde (∼) evoques fluctuations, so we use tildes to indicate random variables (modeling turbulent fields), and we add a subscript g to random fields with Gaussian statistics. For example, adding a noise term in the expression of a smooth non-random field f (x) results in a random field that we denote f (x) if it is intermittent, and fg(x) if it is not.

SCALAR FIELD

Many scalar fields in nature (e.g. density or energy dissipation) are observed to be intermittent, so it is important to find an efficient way of modeling intermittent scalar fields.

To do so, let us start with a formula that we are all familiar with : the gravitational potential φ due to a given distribution of mass ρ is given by

φ(x) = -G R 3 ρ r dV, ( 2
)
where G is Newton's constant. We may reformulate this in terms of the gravitational field g, using Poisson's equation ∇ • g = -4πGρ and the divergence theorem (assuming that the functions vanish sufficiently fast to infinity for surface integrals to be discarded), namely

φ(x) = 1 4π R 3 r • g r 3 dV. ( 3
)
Since g = -∇φ the integral operator in (3) acting on g may be seen as the inverse of the gradient operator, and similarly since ∆φ = 4πGρ the integral operator in (2) acting on ρ may be seen as the inverse of the Laplacian. Expression (3) is an example of a physically relevant scalar field. Now, in order to construct a larger class of such scalar fields, we are going to apply several substitutions to this expression, while keeping its overall structure. First, rather than integrating up to infinite distances, let us consider an arbitrary boundary in the integral, i.e. in (3) make the substitution

R 3 → Rs, (4)
where Rs will now denote the region of integration. Second, let us generalize the power of the power-law kernel

r -3 → r -2hs , (5)
where hs is a free parameter. Third, to avoid the singularity at the origin r = 0 in the integrand brought by this powerlaw behavior, let us make the substitution

r → (r 2 + 2 s ) 1/2 , ( 6
)
where s is another free parameter. This is called a regularization. Mathematically it will spare us difficulties due to divergences and ensure differentiability, while physically the parameter s will efficiently model dissipative effects at small scales, hence this choice of notation traditionally used for small parameters. The final substitution, which we call randomization, consists of replacing the vector field g by a random field

g → η g , (7)
where the components of the vector η g are Gaussian white noises (the subscript g stands for 'Gaussian'), independent of one another, zero-averaged, and with unit variance. This substitution modifies the nature of the integral (3), as it is now a stochastic integral. We give more information on this point and related technical aspects in appendix A.

The above substitutions transform the integral (3) into the following stochastic representation of a scalar random field

sg(x) = Rs η g • r (r 2 + 2 s ) hs dV. (8)
Beware, sg is a general scalar field: it does not bear the physical meaning of a gravitational potential and η g in ( 8) is just some Gaussian white noise vector, not a gravitational acceleration field. Indeed, we started this section mentioning gravity only as an illustration, firstly to build an intuition for where a shape such as (8) may come from, and secondly because this procedure of starting with a well known formula and modifying it with the above substitutions is at the heart of what follows, as we will generalize it to the vector case in section 3.2. As indicated by the notation (tilde and subscript g), expression (8) constitues a Gaussian random process, because it is defined as a linear operation (a convolution) over a Gaussian white noise η g . In fact, this form is basically a spatial version of the fractional Brownian motion (cf. [START_REF] Mandelbrot | [END_REF] in which the exponent of the power-law is called the Hurst exponent, hence our choice of notation h in ( 8) and in all similar integrals throughout this paper. The parameter hs controls the severeness of the singularity and consequently modifies the statistics of the resulting field, as it appears in the exponent of power-law behaving structure functions. More intuitively, while in the classical Brownian motion (also called Wiener process) the increments are independent, in fractional Brownian motion when the Hurst exponent is larger (smaller) than a critical value the increments are positively (negatively) correlated, i.e. there is an increasing (decreasing) pattern in the previous steps (e.g. Rodean 1996, section 13.4).

The fact that ( 8) is a Gaussian field severely restricts its domain of relevance for modeling natural phenomena. Fortunately however, it was shown that from it one may easily construct an intermittent field, by simply applying the exponential function. Indeed, consider the scalar field

s(x) = s0(x) exp [τs sg(x)] , (9)
where τs is a constant and s0(x) is a non-random (i.e. a usual, smooth) scalar field. As illustrated in figures2 B3 and B4, s is intermittent, its degree of intermittency being controlled by τs, hence called the intermittency parameter. Stated in simple words, the field ( 9) is indeed intermittent because the exponential is a non linear function that increases the contrast of the field sg, large values becoming larger and small ones becoming smaller, making large increments more frequent than in the Gaussian case, and resulting in the non-Gaussian wings in the probability density functions of figure B4. This exponential factor then blurs locally, on small scales, the smooth field s0(x) that bears the large-scale behavior.

In a more technical wording, this manner of modeling intermittency (i.e. multifractality) by defining a lognormal process with a long range correlation structure as the exponential of a Gaussian field3 , was introduced in Kahane (1985) and is called 'Gaussian multiplicative chaos'. This terminology refers to Wiener chaos, not to the notion of chaos in dynamical systems, since this theory may be seen as a multiplicative counterpart of the additive Wiener chaos theory (cf Rhodes & Vargas 2014). Such a construction is used for example in hydrodynamics as a stochastic representation of energy dissipation, i.e. a scalar field that shares the same statistical properties as the dissipation field (same expectation value and moments of the filtered fields at a given scale), see Kolmogorov (1962); [START_REF] Mandelbrot | Possible refinement of the lognormal hypothesis concerning the distribution of energy dissipation in intermittent turbulence[END_REF]; [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF] and the Kolmogorov-Obhukov model. Physically this description corresponds to a continuous version of discrete cascades (cf Richardson's picture Richardson 1922) referred to as log-infinitely divisible multifractal processes in the mathematically oriented literature (e.g. Bacry & Muzy 2003;[START_REF] Schmitt | [END_REF].

VELOCITY FIELD

The take home message of the previous section is that applying the simple transformation (9) to a Gaussian scalar field is enough to mimic a real scalar field. Let us now investigate how this may be generalized to a vector field, in order to describe efficiently fluid turbulence. In the same spirit, we will first consider a Gaussian field, and modify it with a physically motivated and intuitive, yet more involved, transformation.

Inspiration from multiplicative chaos:

Introducing matrix Dg

The vectorial version of the Gaussian scalar field ( 8) is obtained by simply substituting a vector product for the scalar product, i.e. symbolically

• → ×. (10)
Indeed, Robert & Vargas (2008) showed that a random vector field that is Gaussian, incompressible, homogeneous, isotropic, and singular can be expressed as the stochastic integral

vg(x) = Rv η g × r (r 2 + 2 v ) hv dV, (11)
where v and hv are constants, and Rv denotes the region of integration. Roughly speaking, the singular power-law shape of the kernel in (11) ensures the scale invariance property (structure functions, i.e. moments of velocity increments, behaving as power-laws) that is needed to satisfy Kolmogorov's 2/3-law. For this reason (11) is a good starting point to model a turbulent velocity field, but it is not enough per se because it lacks two important features. Firstly, it is not intermittent. Secondly, it has a vanishing third-order moment of velocity increments (and actually all odd-order correlators). Therefore it does not reproduce the 4/5-law and is not enough to model energy transfer across scales.

In order to introduce intermittency, Robert & Vargas (2008) suggested disturbing the white noise vector by the measure (they work using measure theory rigorously) given by multiplicative chaos, in the spirit of Kahane (1985). In other words, they replace in (11) the white noise η g by exp[τ dg(y) + c] η g , where τ and c are constants, and dg is a Gaussian scalar field. While this indeed provides an intermittent behavior, it does not yield energy transfer: the skewness of the resulting vector field vanishes at all scales. This is because the velocity increments do not exhibit asymmetry. Soon afterwards, Chevillard et al. (2010) pointed out that if instead of using a scalar field ( dg(y) above) one uses a matrix field ( Dg(y) below), the symmetry is broken and this limitation is bypassed: the resulting vector field contains both intermittency and an energy cascade. Based on these works, let us consider the velocity field

v(x) = Rv [e τω Dg η g ] × r (r 2 + 2 v ) hv dV, (12)
where hv and τω are parameters, whose values will be discussed in section 5, and Dg is the matrix (33) below. The purpose of the next section is to explain the reasons leading to this expression for Dg.

3.2 Inspiration from Biot-Savart's law: Linking Dg to vorticity ω Expression (12) should look familiar to the reader: it can simply be seen as a modified Biot-Savart's law. Indeed, recall that for an incompressible fluid, velocity may be expressed in terms of vorticity (ω ≡ ∇ × v) according to Biot-Savart's law [START_REF] Majda | Vorticity and Incompressible Flow[END_REF])

v(x) = 1 4π R 3 ω × r r 3 dV, ( 13
)
and starting with (13), one obtains (12) with the four following, physically meaningful, substitutions. Firstly, the integration domain is limited Consequently, this allows one to fiddle with the geometry of the field, as illustrated in figure B2 in the scalar field case.

R 3 → Rv. (14)
In section 5, we will take Rv to be a sphere of radius equal to the injection scale, and eddies of that size will be visible (see figures B5 or B9 for example). Secondly, the kernel is generalized:

r -3 → r -2hv , (15)
which allows us to model various types of correlations, as does the Hurst exponent H in fractional Brownian motion. Indeed in fractional Brownian motion this coefficient relates to the persistence of the time series. When H is greater than a half, signals have positive values that persist for a while, before switching to negative values that, in turn, also persist, and so on. On the contrary, for H smaller than a half the signal 'anti-persists', in the sense that it is wildly fluctuating.

Here we are not dealing with time series (fractional Brownian motion) but a spatial version of it (fractional Brownian fields). Therefore the Hurst exponents in our work are not measures of long-term memory but of longe-range correlations. These coefficients enable us to control the smoothness of the fields, by modifying how they are correlated in space.

For instance, Chevillard et al. (2010) choose the power-law entering their rate-of-strain matrix equal to 7/2 which imposes logarithmic long-range correlations. Numerical values are discussed in section 5. Thirdly, the norm is regularized, as was done in (6) in the case of a scalar field. This again mathematically guarantees that ( 12) is differentiable, and physically gives a degree of freedom to model dissipation at small scales (see section 3.4 for more information on dissipation in this model). Finally, in order to complete this analogy, we need the following fourth substitution

ω → e τω Dg η g , (16)
but as such, giving a physical meaning to this is not obvious at all, especially because at this stage we do not even have an expression for Dg. However, exhibiting the substitution ( 16) hints towards a connection between vorticity and matrix Dg. This is surprising a priori: Dg originates from the concept of multiplicative chaos, which was suggested in works unrelated to fluid turbulence. Still, led by this intuition, let us in this section have a closer look at the dynamics of vorticity. As we will see in the next section, interpreting expression (12) as stemming from (13) will pay off, as it will lead us to a relevant explicit expression for matrix Dg, and bring a physical meaning to (12). Momentum conservation for an inviscid fluid of pressure p and constant density ρ0 in a gravitational field g reads

dv dt ≡ ∂tv + v • ∇v = - 1 ρ0 ∇p + g. ( 17
)
For now we are working in the inviscid limit because ultimately the substitution (6) will take dissipative effects into account phenomenologically. The transport equation for vorticity is obtained by taking the curl of ( 17), recalling that ω and v are divergence-free under the present assumptions. Doing so yields

dω dt = G • ω, (18)
where, to improve readability, we denote by G the gradient of velocity:

G ≡ (∇v) T , (19)
with T indicating transposition. Before going further, let us build some intuition of the underlying dynamics, by analyzing the flow locally. Taylor expanding the velocity field at a given point x0 and a given time t, to linear order in |x -x0| we have

v(x) = v(x0) + G(x0)(x -x0). (20)
Next, let us decompose the matrix G into its symmetric and antisymmetric parts as

G = D + Ω, (21)
where

Dij ≡ (∂ivj + ∂jvi)/2, Ωij ≡ (∂ivj -∂jvi)/2. ( 22
)
The symmetric part D is called the deformation (or rate-ofstrain) matrix. Being real and symmetric, it may be diagonalized by an orthogonal matrix, which corresponds physically to stretching and compression along its eigenvectors at a rate given by its eigenvalues. Stretching and compression occur jointly because of incompressibility. Indeed, the condition ∂ivi = 0 in ( 22) implies that D is traceless, and thus the sum of its eigenvalues vanishes. Positive eigenvalues correspond to stretching and negatives eigenvalues to compression. The antisymmetric part Ω is called the rotation matrix. Indeed, using the identity ijk kab = δiaδ jb -δ ib δja, it can be written in terms of the vorticity as

Ωij = ijk ω k /2, (23)
i.e. it is a cross product: In matrix form

Ω = 1 2 ω × . ( 24
)
This reveals that it corresponds to an instantaneous rotation in the direction of -ω with angular velocity |ω|/2. As a result, the expansion (20) becomes

v(x) = v(x0) + D(x0)(x -x0) + 1 2 ω × (x -x0), (25)
from which we see that the flow is locally the sum of a translation, a deformation (stretching/compression) and a rotation (see e.g. [START_REF] Majda | Vorticity and Incompressible Flow[END_REF]. Introducing ( 21), together with (24), into (18) directly gives

dω dt = D • ω. (26)
Physically, this equation describes the effect of vortex stretching: Stretching a vortex along its axis reduces its cross-section due to incompressibility, and this increases its angular momentum due to conservation of angular momentum.

Hence, in order to have a closed evolution equation for ω, it remains to explicit D in terms of ω. Given the definition ( 22) of D, we need to compute the gradient of v. It is tempting to simply differentiate under the integral sign in ( 13), but unfortunately this yields a singularity and for the purpose of the present work we may simply accept the result [START_REF] Majda | Vorticity and Incompressible Flow[END_REF])

D(t, x) = 3 8π P.V. R 3 (r × ω) r + r (r × ω) r 5 dV, ( 27
)
where products like x y are tensor products (x y|ij = xiyj), and 'P.V.' indicates a Cauchy Principal Value integral, defined as

P.V. R 3 f (y) dV ≡ lim →0 |y| f (y) dV. (28)
Equation ( 26) together with ( 27) constitutes the equation satisfied by ω. Rather than trying to solve it exactly, Chevillard et al. ( 2010) expose a cunning method to obtain a very efficient approximate solution, based on the short-time evolution as follows.

First, recall that equation ( 26) is written in Lagrangian form, with the material derivative d/dt = ∂t + v • ∇. Consider a Lagrangian particle at position X at an initial time t0, and at position x(t, X) at a later time t. The vorticity ω(t, X) felt by this particle at time t satisfies (26) where in general D(t) differs from its value at t0. However, after a short enough period of time, say t -t0 = τK where τK is the Kolmogorov timescale which is the characteristic timescale of variation of the correlation of gradients, one may consider that D(t) ≈ D(t0) is a good approximation. Hence, over that period of time, the matrix in equation ( 26) is roughly constant, so we can approximate its solution as ω(t, x(t, X)) ≈ e (t-t 0 )D(t 0 ,x) ω(t0, X).

(29)

This equation shows the two processes at play: the transport of particles from position X to x (see the arguments of ω in the right and left hand sides), and the stretching of the vorticity field by the initial deformation tensor D(t0). Now, it turns out that in the dynamics of turbulence, the components of the matrix τK D(t0) are typically of order unity, while the distance traveled by particles due to advec-

tion τK v(t0) is of order R -1/2 e
where Re is the Reynolds number [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]. Besides, since in astrophysics we are mainly dealing with high Reynolds number fluids, in the following we will consider that the advection of vorticity is negligible compared to its stretching, and finally write ω(t, x) = e (t-t 0 )D(t 0 ,x) ω(t0, x).

(30)

Since this expression stems from an expansion of equation ( 26) near a given time t0, in the literature this procedure is referred to as the 'recent-fluid-deformation' procedure [START_REF] Chevillard | [END_REF]Meneveau 2011).

Putting it all together

Let us now put the final result of section 3.1 together with that of section 3.2: The key feature of expression ( 30) is that this matrix exponential is reminiscent of the matrix exponential appearing in expression ( 12), which completes the analogy between multiplicative chaos and vorticity. Hence, let us finally consider the velocity field

v = 1 4π Rv ω × r (r 2 + 2 v ) hv dV, (31)
where

ω = e τω Dg ωg, (32)
and

Dg = 3 8π Rω (r × ωg) r + r (r × ωg) (r 2 + 2 ω ) hω dV. ( 33
)
In these expressions ωg is a Gaussian white noise vector (just like η g previously, but we change notation now to match its physical interpretation), Rω denotes the region of integration, and hω and ω are constants, with similar roles as hs and s of the scalar case (8). By suggesting this form, a few additional substitutions have been made, which we now detail.

Clearly Dg in (33) stems from D in ( 27), with a few modifications. Firstly, we have introduced a parameter ω in the denominator, i.e. the singular kernel has been regularized, as was done in the scalar field case (6) and in the expression of the velocity field (31). Thanks to this, we may discard the Cauchy Principal Value. Secondly we have generalized the power-law to some free parameter hω, in order to gain an additional degree of freedom in the modeling and to have a larger parameter space for data fitting. The hydrodynamics literature we take the above formulae from is dedicated to the local structure of turbulence and ultimately aims at obtaining exact solutions of the Navier-Stokes equations. In contrast, our own purpose is different, and we do not fix a priori the parameters, e.g. we do not place constraints on the integration regions Rv and Rω. This will give us more degrees of freedom, as illustrated in the scalar field case in figure B2, to fit astrophysical data which include large-scale features. This is somewhat of a poor man's way, as a first approach, to model correlations between large and small scales, in the same spirit as we included a largescale position dependent amplitude s0(x) in the scalar field (9). More refined modeling, which would for instance take into account the advection of vorticity, which is neglected in equation ( 30), is left for future work. Thirdly Dg is a 'randomized' version of matrix D, which is what the tilde notation indicates, as we made the substitution ω → ωg.

(34)

This substitution models the randomness of the stretching occurring in formula (32) (detailed in the next paragraph) due to Dg, since D has the physical meaning of local vortex stretching. Choosing the very same white noise vector ωg in both ( 32) and ( 33) is a thoughtful decision: Chevillard et al. ( 2010) as well as Pereira et al. (2016) showed that if one considers two independent white noises, all odd-order structure functions vanish, which prevents energy transfer. Finally, and most importantly as it is at the heart of the physical interpretation of this model, the velocity field ( 31) is independent of time (it describes a stationary state), which seems in contradiction with the fact that its construction is inspired from dynamics. Indeed, completing the analogy outlined in section 3.2 led us to make this last substitution

t -t0 → τω, (35)
when going from ( 30) to (32). Since this amounts to substituting some parameter for the time variable characterizing the stationary state (more precisely, its intermittency), ( 35) is referred to as stationarization.

With this in mind, the following picture emerges 4 . To begin with, consider a first 'timestep' (of say a Kolmogorov timescale τK ) during which expression (32) may be seen as an evolution equation: an initial vorticity field ωg with Gaussian statistics is stretched according to matrix Dg, resulting in an intermittent vorticity field ω. By incompressibility, through Biot-Savart's law, this translates into the velocity field (31). What the stationarization hypothesis suggests is that this occurs at each timestep: somehow, by some undetermined non-linear process, the non-Gaussianity of the vorticity field is suddenly 'blurred', resetting the underlying vorticity field to a Gaussian field ωg again, such that, at the subsequent timestep, the stretching (32) repeats itself. This repeated succession of stretching/blurring is such that in the end the stationary state of the fully developed turbulent flow has similar statistics as those of a stretched Gaussian vorticity field (it does indeed, cf. the successes of this modeling detailed in section 3.4). The degree of intermittency is then controlled by how long these stretching periods last, which in practice is encoded phenomenologically in the intermittency parameter τω. For example if they are extremely short (τω → 0), there is not enough time for the stretching to occur significantly. The vorticity field then remains Gaussian at each timestep, resulting in a Gaussian stationary velocity field (i.e. (31) tends to (11)). The above picture also helps understand how formula (31) can describe a fully developed state of turbulence, which is the result of highly non-linear dynamics, while it was constructed from a linearization 5 . With this interpretation, taking the very same white noise in ( 32) and ( 33) becomes natural, and not only a technical requirement allowing for energy transfer: In the deformed vorticity field (30) it is the initial value D(t0, x) that comes into play, so it is indeed the 'initial' vorticity ωg that makes sense in Dg.

A word of caution however concerning this interpretation: strictly speaking, ω is not exactly the vorticity field in this model. For example we did not impose ω to be divergence free, while this is required for a vorticity field, and because these are stochastic integrals, one should be a little bit careful with its physical dimension (see appendix A). To clarify: the procedure here is to first suggest (31) as our model for a velocity field v, and only then one defines the vorticity field as the curl of v. All the intermediate steps outlined above should be seen more as ideas useful to construct this model, and build an intuition out of it.

Assets of the adopted method

While part of the procedure presented above is based on rigorous derivations, partly taken from the mathematics literature, some points rely on phenomenological arguments. Yet, by means of both analytical and numerical arguments, previous authors have shown that the model (31) gives a surprising amount of correct predictions when confronted to both direct numerical simulations and experimental data (e.g. Chevillard et al. 2012). What is even more encouraging is that some of its realistic properties are not ingredients introduced by hand but by-products of this construction, which suggests that it really does capture an important part of the underlying physics. In this section we briefly review these successes.

(i) A first important feature is the existence of intermittent corrections. Let us define the velocity increment of lag as

δ v(x) ≡ v(x + ) -v(x). (36)
Longitudinal and transverse increments are the projections of ( 36) respectively along and perpendicular to the direction of . One common way to identify intermittency is to compare the probability density functions (PDFs) of the considered field to those of a Gaussian field. And indeed the PDFs of longitudinal and transverse increments of (31) undergo a continuous deformation as the norm of the lag is decreased, from a Gaussian shape at large lags towards large tails at small lags. These often called 'non-Gaussian wings' are typical signatures of intermittency. Through several examples, we show in section 5 that these wings become increasingly significant as the parameter τω increases, which supports naming it the intermittency parameter. More details may be found in the above-mentioned literature.

(ii) Another classical way to characterize intermittency in isotropic turbulence studies is by analyzing the power-law behavior of structure functions in the inertial range. By definition, the n th order structure function is the n th moment the development of non-Gaussian statistics [START_REF] Chevillard | [END_REF], but considering the short-time dynamics of the transport equation ( 26), we have made its right hand side linear in ω, while in full form it is non-linear in ω.

of velocity increments Sn( ) ≡ (δ v) n where brackets indicate the expectation value [START_REF] Frisch | Turbulence. The legacy of A[END_REF]. In the inertial range, Sn ∝ τn , where the dependence on n of τn quantifies the intermittency: the field is intermittent if and only if τn depends non-linearly on n. For the velocity (31), this non-linearity was confirmed numerically, and also partially analytically: under the simplifying assumption of independence of Dg with the white noise in (31), Pereira et al. (2016) showed that τn = an + bn 2 , where a and b are constants, with b ∝ τ 2 ω , i.e. the quadratic behavior is determined by the intermittency parameter. This again justifies the name given to τω.

(iii) Pereira et al. (2016) further characterized the non-Gaussian behavior of fluctuations by analyzing in details the skewness S ≡ (δ v) 3 /[ (δ v) 2 ] 3/2 and flatness F ≡ (δ v) 4 /[ (δ v) 2 ] 2 of velocity increments [START_REF] Frisch | Turbulence. The legacy of A[END_REF]. The propreties of both longitudinal and transverse increments match the results of direct numerical simulations. In particular, they showed that this velocity field has a higher level of intermittency in the transverse case than in the longitudinal case, which is a feature observed in experiments (see Pereira et al. 2016, and references therein).

(iv) The degree of freedom provided by the Hurst parameter hv allows the velocity field to satisfy the experimental 2/3 law of turbulence.

(v) Another aspect that has been scrutinized in the literature is the statistics of the velocity gradient G, which regulates the dynamics of the turbulent flow. Its eigenvalues are particularly meaningful physically (see section 3.2) and have peculiar statistical properties that (31) reproduces remarquably well. More specifically, these eigenvalues are the solutions of the characteristic polynomial with coefficients given by the two invariants

Q ≡ -1 2 tr(G 2 ) = 1 4 |ω| 2 -1 2 tr(D 2 ), R ≡ -1 3 tr(G 3 ) = -1 4 ω • D • ω -1 3 tr(D 3 ). ( 37
)
Physically, Q is a combination of enstrophy (|ω| 2 ) and dissipation (both per unit viscosity) such that when Q > 0 rotation dominates locally in the flow, and when Q < 0 the region is dissipation-dominated. R is a combination of the local production of enstrophy and dissipation, such that when R > 0 the flow tends to create dissipation, and when R < 0 more enstrophy is produced. Now, it turns out that when representing the joint probability density of R and Q, a characteristic shape appears, known as the 'teardrop shape in the RQ-plane'. See for instance [START_REF] Tsinober | An Informal Introduction to Turbulence[END_REF]; Wallace (2009); Meneveau (2011) for both experimental and theoretical aspects of this feature. In particular, this shape is symmetric with respect to the R = 0 line in the Gaussian case (τω = 0) and tilted anticlockwise otherwise. As detailed notably in [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]; Pereira et al. (2016), this geometry is indeed recovered with (31). This diagnostic tool is yet another piece of evidence of its intermittent behavior, and of the local dissipation occurring in a pretty realistic manner.

(vi) Other essential and non-trivial features are the orientation properties of vorticity with respect to the eigenframe of the rate-of-strain matrix D. In particular, the preferential alignement of vorticity with the eigenvector associated with the intermediate eigenvalue of D does not occur in the Gaussian case, but gets stronger as the intermittency parameter increases. All this is correctly reproduced by (31), cf e.g. [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]; Pereira et al. (2016Pereira et al. ( , 2018)).

(vii) Equation ( 31) is the first stochastic process proposed in the literature to have a non-vanishing mean energy transfer across scales, as observed numerically in Chevillard et al. (2010) and then, by means of a perturbative expansion in the intermittency parameter, Pereira et al. (2016) confirmed analytically the existence of this cascade. As this feature stems from the matrix nature of the multiplicative chaos used, for more details on mathematical aspects we refer the reader to Chevillard et al. (2013) who constructed the theory of matrix multiplicative chaos.

(viii) These formulae are well behaved, notably they are differentiable thanks to the regularization. This aspect is not obvious since they are built from irregular mathematical objects (e.g. sample-paths of Wiener processes are almost nowhere differentiable).

(ix) Last but not least, another success of this approach is that it results in fields that are easy to obtain numerically. Because most equations in this work are convolutions, their implementation is both simple and fast, thanks to Fast Fourier Transform algorithms and because convolutions in Fourier space are just products. To give an order of magnitude, with only 32 CPUs one may generate a 2048 3 resolution cube in only 5 minutes.

As a closing remark, let us add that more achievements related to this description may be found in the abovementioned literature, notably in the Lagrangian framework, studying the time evolution. Here we stick to stationary fields in the Eulerian description. Also, here we aim at constructing expressions that are as handy as possible. For this reason, we chose the simplest regularization, namely (6), but more elaborate regularized norms are used and discussed in the literature. Similarly, instead of our integration regions, it is more common to introduce large-scale cutoff functions in the kernels, which is more convenient to analyze the statistics of the field depending on the properties of these functions.

MAGNETIC FIELD

The above substitutions constitute building blocks to be assembled to construct physically motivated turbulent vector fields. The successes mentioned in the previous section motivate us to extend this approach to magnetized fluids, which are ubiquitous in astrophysics. Per se, this is a very ambitious project, and we will limit ourselves to proposing three straightforward extensions. To do so, we need to adapt the matrix (33) and ensure that the magnetic fields remain divergence-free.

Matrix field of the underlying dynamics

Compared to the hydrodynamical case (17), the momentum equation is now coupled to the induction equation, and contains the additional contribution from the Lorentz force, j × B, where the current density j is given by Ampère's law j = (∇ × B) /µ0 and the constant µ0 is vacuum permeability. In the same spirit as previously, we may omit viscosity and resistivity for now, since dissipative effects will ultimately be taken into account through regularizations similar to (6). We are then left with

∂tv + v • ∇v = - 1 ρ0 ∇p + g + 1 µ0ρ0 (∇ × B) × B, ∂tB = ∇ × (v × B) , (38)
together with the constraint ∇ • B = 0. As a first approach, in this paper we consider incompressible MHD and the case when the back-reaction of the magnetic field, due to the Lorentz force, is negligible. This approximation, called kinematic MHD, is valid as long as magnetic energy is small compared to kinetic energy, i.e. B 2 /(2µ0)

1/2ρ0v 2 . This condition is often not satisfied in astrophysics. However, building our models in this framework is also motivated6 by 3D MHD simulations such as those of Servidio et al. (2008) which show a natural tendency towards force-free behaviors in MHD turbulence. Indeed, as nonlinearity develops, strong alignments spontaneously appear notably between j and B, which is called 'beltramization' in reference to Beltrami vector fields which are fields parallel to their own curl.

In this framework, the induction equation reads

dB dt = G • B. (39)
Introducing ( 21), together with (24), into (39) leads to

dB dt = D • B + 1 2 ω × B. (40)
This equation for the magnetic field is the counterpart of equation ( 26) for the vorticity. The two terms on the right hand side describe flux tube stretching (similar to vortex stretching in (26)) and shear, respectively. By analogy with the hydrodynamical case where we made the substitution (34) to construct Dg, we suggest considering the following 'randomized' velocity gradient

Gg = Dg + 1 2 ωg×, (41)
i.e. in index notation ( Gg)ij = ( Dg)ij + ijk ( ωg) k /2. Expression (41) is the most natural choice to remain consistent with the interpretation outlined in section 3.3: it is the 'initial vorticity' ωg that comes into play in the 'short-time' deformation matrix Gg. Also, with this choice Gg is a Gaussian process, being built from a linear operation on a Gaussian process, and the purpose of the theory of multiplicative chaos is precisely to give a meaning to the exponential of such a field (Rhodes & Vargas 2014).

Three models

Assuming that the matrix field Gg in ( 41) is an accurate representation of the local dynamics in the magnetized turbulent flow, let us now propose three models for MHD turbulence, that are in principle complementary.

Model 1) By analogy with the short time expression of vorticity (30) from the hydrodynamical case, it is natural to consider the short time solution of (39), namely the matrix exponential of Gg. Doing so, the initial value of the magnetic field appears, just like ω(t0, x) in (30). When suggesting the velocity field expression (31), the initial value of vorticity was replaced by a Gaussian field through the substitution ( 16). Here instead, we take advantage of this degree of freedom to introduce some ordered magnetic field B0. Introducing this is important to model astrophysical fluids since the observed magnetic fields in the Universe are usually turbulent with large-scale mean fields (e.g. Ferrière 2007; Jaffe 2019). Hence, our first model is

B (1) = e τm Gg B0, ( 42
)
where τm is a magnetic intermittency parameter. Physically, the magnetic field B

(1)

is the result of the random stretching/compression and shearing (encoded in the matrix exponential, the strength of which is controlled by τm) of a large-scale ordered magnetic field B0(x), due to the turbulent velocity field (31) that it is coupled to. The divergence-free condition will be discussed at the end of this section.

Model 2) By analogy with Biot-Savart's law (13) from the hydrodynamical case, it is natural to consider Biot-Savart's law for magnetic fields

B = µ0 4π R 3 j × r r 3 dV. ( 43
)
To strictly follow the procedure performed with ω, the next step would be to analyze the evolution equation of j, and its short time dynamics would reveal a relevant matrix. However, here we propose a simpler and faster approach. We already found a seemingly relevant generalization of Dg, namely Gg, costlessly obtained from the governing equations (38). And fortunately, the aforementioned beltramization exhibited in simulations of Servidio et al. (2008) prompts us to indeed use Gg: since j tends to be aligned with B, the matrix governing its evolution should in some way resemble the matrix governing the evolution of B, namely (39). Hence, rather than digging further to look for the precise matrix out of the evolution of j, we will directly work with Gg. Admittedly, this is also very convenient. In addition, the simulations of Servidio et al. (2008) also indicate that the alignments are directly related to the self-organization of the magnetofluid giving rise to spatial intermittency, which our formulae do take into account through the intermittency parameters τω and τm. Servidio et al. (2008) furthermore suggest that these dynamically generated correlations suppress nonlinearity, i.e. the equations tend to be quasilinear rather than fully nonlinear. This is reminiscent of the interpretation mentioned in section 3.3 to account for the intriguing fact that our formulae are obtained after a linearization to describe a fully developed turbulent state.

From the above considerations, we propose as a second model

B (2) = µ0 4π Rm j (2) × r (r 2 + 2 m ) hm dV, (44)
where Rm, hm and m have similar physical meanings as Rv, hv and v in the velocity field (31) (for examples see section 5). The field j

(2)

acts like a random current density, which, by analogy with (32), we build as

j (2) = e τm Gg j g , ( 45
)
where j g , the 'undeformed' current density, is a Gaussian white noise. Note that j g must be different from ωg, otherwise, like in the hydrodynamical case, the antisymmetric part of the velocity gradient Gg in (45) would not make any contribution. This magnetic field model would then be redundant with the velocity field model v in (31). On the contrary, one may choose j g to be independent of ωg.

In figure B8 we show what (44) then typically looks like: qualitatively at least, this model is fairly good, especially given its convenient compact expression, and it is worth confronting to data. However, it may have an important theoretical shortcoming: In Pereira et al. (2016) the authors show that considering two independent white noises in the velocity field (31) prevents energy transfer. In our case, we have an additional contribution from shear, so it remains to be seen whether this model has the same limitation or not. A possible improvement may be the following. The two aforementioned cases are two extremes: in one case j g is fully correlated to ωg, while it is completely independent in the other. A balance is to consider a j g 'mildly' correlated to ωg, e.g. a fractional Gaussian vector field such as (11), namely j g = R j ωg × r (r 2 + 2 j ) -h j dV . On the one hand it has the advantage that we control its statistical properties with simple parameters, and this form guarantees it to be divergence-free7 . On the other hand, this introduces yet another set of controlling parameters Rj, hj and j but which are not physically motivated, unlike the other sets of such parameters.

Because ( 44) is formally similar to the velocity field (31), we may give to the magnetic field B

(2)

an analogous physical meaning to that presented in the end of section 3.3, replacing vorticity with current density, and adding the effect of shear (contribution from the antisymmetric part of Gg, due to local rotation). Hence, ( 44) introduces an additional regime compared to its hydrodynamical counterpart (31), namely the regime of shear dominated turbulence, when the local rotation Ω dominates over vortex stretching Dg. Compared to model 1, in this model there is no largescale feature, such as s0 and B0. Therefore, in principle this second model is relevant to describe the local, small-scale structure of a turbulent magnetic field. However, given its phenomenological construction, it may turn out to have a wider range of application depending on what it is used for.

Model 3) Finally, let us propose a model that combines the features of models 1 and 2, by considering again a Biot-Savart law form, but inserting some large-scale feature, namely

B (3) = µ0 4π Rm j (3) × r (r 2 + 2 m ) hm dV, (46)
where

j (3) = e τm Gg j 0 . ( 47
)
At first sight, this third model looks like the second model, but it is fundamentally different. In j

(3)

the undeformed current density j 0 is not a random field as in j

(2)

, but an ordered, large-scale current density field. If need be, it may be linked to some ordered magnetic field through Ampère's law j 0 = (∇ × B0)/µ0. In that sense, ( 46) is somehow 'less' random than (44) since the randomness arises only from the matrix exponential. Thus, technically speaking, this third model is a priori more prone to analytical calculations, since the main mathematical difficulty in ( 44) is that it contains the product of two random distributions (e.g. section 3.4.3 of [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF].

As in model 2, the Biot-Savart form introduces additional degrees of freedom (Rm, hm and m) compared to model 1, while at the same time, as in model 1, it contains a large-scale field (j 0 ). Physically, the magnetic field

B (3)
is the result of the random stretching/compression and shearing (encoded in the matrix exponential, the strength of which is controlled by τm) of a large-scale ordered current density field j 0 , due to the turbulent velocity field (31) that it is coupled to. Pereira et al. (2016) state that the Biot-Savart law form of their hydrodynamical model ensures that their velocity field is divergence-free, as it should since they focus on incompressible dynamics. They also indicate that numerically one has to take a high enough resolution for the size of the grid cells to be small compared to the regularization parameter . Otherwise gradients are not well approximated and the divergence of the field can be, in standard deviation, of the order of the gradient of one of its components, rather than much smaller. In our magnetized case, for models 2 and 3 we also use Biot-Savart forms, and we build model 1 from an evolution operator that conserves the divergencefree condition. Indeed, matrix G stems from the induction equation ( 39) governing the evolution of magnetic fields, and magnetic fields remain divergence-free throughout time. Hence, the very construction of our models is suited to yield divergence-free fields, as it should. However, a word of caution is in order. We generalize to arbitrary exponents and integration regions these Biot-Savart laws, and we randomize the matrix G. At this stage we cannot prove rigorously that the magnetic fields indeed remain exactly divergencefree. That being said, we observe numerically that they are divergence-free up to satisfying levels for the range of parameters that we have tested so far, as long as the resolution is high enough (grid cell size greater than the chosen regularization lengths ), as pointed out in Pereira et al. (2016). In any case, if need be, the simple way out (which is numerically cheap since our formulae are generated in Fourier space) is to apply the projector P which operates on a given

m Lm hm τm j 0 (x)

Table 1. Table of the degrees of freedom (i.e. the free parameters and free functions) appearing in the expressions of the scalar, velocity and magnetic fields introduced in this paper. For simplicity we consider the case of spherical regions R i of radius L i , with i = s, ω, v and m.

vector h as

P(hi) = F -1 δij -kikj/k 2 F(hj) , (48)
where F denotes the Fourier transform. Indeed, for any vector h, P(h) is divergence-free. This projector is frequently used in numerical schemes involving magnetic fields, and it also naturally appears in the rewriting (i.e. without approximation) of incompressible MHD (e.g. Galtier 2016). We summarize in table 1 the many free parameters (and functions) introduced through all the above formulae. These constitute degrees of freedom that may be adjusted to data, in order for example to quantify the degree of intermittency of an observed magnetic field. They also provide a handy tool to model analytically a wide variety of astrophysical environments. Furthermore, this diversity opens a vast parameter space, corresponding to many regimes of the dynamics. For example, since the hydrodynamical and the magnetic regularization parameters v and m correspond to dissipation parameters, it is reasonable to consider that the ratio v / m corresponds to the magnetic Prandtl number Prm (we remain cautious while stating this simply because their precise link with Prm has not been rigorously established yet). Therefore, varying this ratio we may generate turbulent magnetic fields with Prm 1 or Prm 1, which is out of reach of numerical simulations, as it requires a huge dynamical range. In this sense our approach is complementary to numerical simulations.

EXAMPLES

The aim of the present work is to provide a tool for observers or anybody who needs to quickly generate synthetic data with controllable properties (statistics and shape of the structures). We do not claim that our formulae constitute solutions of the MHD equations. Rather, we show a way to build non trivial vector fields with many degrees of freedom that are as physically motivated as possible, in order to be intuitive to handle, based on the MHD equations. It is only in a following paper that we will derive quantitative results on the statistics of the proposed fields, to evaluate how close to, or how far from, real dynamo-generated fields they really are.

The many degrees of freedom summarized in table 1 suggest rich and interesting interplays between the hydrodynamical and magnetic parameters. But an exhaustive anal-ysis of the effects of each of them is also beyond the scope of this paper. Instead, in this section we expose just a couple of selected examples, intending to give the reader a flavor of some key features.

To begin with, let us provide some relevant values for the free parameters, based on Chevillard et al. (2010Chevillard et al. ( , 2011Chevillard et al. ( , 2012)); [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]: (i) in the Gaussian scalar field (8) one may take hs = 3/4 typically, (ii) in the works just quoted, it is shown that hv = 13/12, corresponding to a Hurst exponent H = 1/3 in their notation, makes the velocity field (31) satisfy the experimental 2/3 law of turbulence, (iii) in addition, it is argued that in order for the process to have a finite variance and convenient mathematical properties, one should take 3/4 < hv < 5/4 (corresponding to 0 < H < 1), a constraint that may thus be taken as a prior in data analyses, (iv) in order for the components of Dg to be correlated logarithmically in space, hω = 7/4, according to [START_REF] Mandelbrot | Possible refinement of the lognormal hypothesis concerning the distribution of energy dissipation in intermittent turbulence[END_REF] and Kahane (1985), (v) in the context of laboratory fluids, the value τω = 0.1 is obtained from experiments, after translating the γ 2 = 6.7 × 10 -2 of Pereira et al. ( 2016) into our notations. Indeed, in laboratory and numerical flows, these authors find a universal behavior, independent of the flow geometry and the Reynolds number, for the longitudinal increments, with data well reproduced for that particular value of τω. Finally (vi) regarding the choices of Lv and v , in order to reach a large inertial range, one should maximize the Lv/ v ratio. The small-scale cut-off v is typically taken equal to a few times the resolution length. For instance for a simulation performed in a 1-periodic box with N 3 collocation points, Lv = 1/2 and v = 3/N .

We now show examples of various scalar fields, velocity fields and magnetic fields based on the above formulae. Most figures were produced using Mayavi, an application and library for interactive scientific data visualization and 3D plotting in Python [START_REF] Ramachandran | [END_REF].

Firstly, we need to generate the Gaussian vectorial field ωg, from which all the other fields shown in this paper were generated, making them by construction highly correlated. We display this white noise in figure B1. What is interesting to bear in mind when visualizing this realization, is that obviously it does not have any particular shape or structure, neither in vector form nor in norm. ωg is really only the building block that gives all the other fields their randomness: the structures and non-Gaussianity are produced by the non-linear transformations ( 9), ( 31), ( 42), (44), and (46) applied to it. Thanks to this, if one varies continuously the free parameters, despite the randomness, the fields also vary continuously: we do not generate new realizations for each value of the free parameters.

Secondly, let us illustrate the effect of the regions of integration, in the case of a scalar field, which is easier to visualize. Rs determines the large-scale shape, which may depend on position. In figure B2, we consider a simple ellipsoid, and vary its elongation. The shape of this region roughly determines the trend appearing in the spatial distribution: a spherical region yields rather isotropic 'clouds', a cigare-like region yields filamentary structures and a pancake-like region generates more sheetlike structures. In the rest of this paper, we will stick to spherical regions, namely spheres of radius Li centered at position x:

Ri(x) = y ∈ R 3 /|r| Li , (49)
where i = s, v, m for respectively a scalar field, a velocity field, and a magnetic field. We do so both for simplicity and because in the literature this has been proven to be relevant at least for the local structure of isotropic turbulence. Thirdly, we analyze how the scalar field ( 9) changes as we vary the integral scale Ls together with the intermittency parameter τs. In figure B3, we show 2D slices of 3D realizations of the scalar field, for various values of Ls and τs, in order to show qualitatively how these modify the shape of the structures and the intermittency. Focusing on the modifications along a column, one clearly sees that Ls modifies the global size of the structures formed, i.e. the correlation length. Focusing on the modifications along a row, it appears that basically τs modifies the contrast, and hence controls the intermittency of the field. To complement these intuitive pictures, in figure B4 we provide more quantitative considerations by plotting the PDFs of the increments corresponding to each of the cases shown in figure B3 for three different lags. We normalize the PDFs to unit variance and arbitrarily shift vertically the curves for the sake of clarity (as in Pereira et al. 2018, for example). Focusing on a given row, on the left the PDFs are very close to those of a Gaussian random field, while as we move to the right, the characteristic non-Gaussian wings appear as τs increases. This confirms quantitatively the feeling from figure B3 that τs controls the intermittency of the field. Then considering the modifications along the columns as well, it appears that the 'efficiency' of τs to provide intermittency increases with increasing Ls.

Fourthly, we analyze how the velocity field (31) changes as we vary the integral scale Lv together with the intermittency parameter τω. In figure B5, we show 2D slices of 3D realizations of the velocity field, for various values of Lv and τω, in order to show qualitatively how these modify the shape of the structures and the intermittency. Focusing on the modifications along a column, one clearly sees that Lv modifies the global size of the structures formed, i.e. the correlation length. Focusing on the modifications along a row, it appears that basically τω modifies the contrast of the norm of the vectors, and hence controls the intermittency of the field. To complete these intuitive pictures, in figure B6 we provide more quantitative considerations by plotting the PDFs of longitudinal increments corresponding to each of the cases shown in figure B5 for three different lags. As in the scalar field case, non-Gaussian wings appear as we increase the intermittency parameter τω. These PDFs are very similar to what is usually obtained in direct numerical simulations, but here we may directly control the steepness of the non-Gaussian wings with τω and the fields are generated at very little computational cost. Focusing on a given row, on the left the PDFs are very close to those of a Gaussian random field, while as we move to the right, non-Gaussian wings appear as τω increases. This confirms quantitatively the feeling from figure B5 that τω controls the intermittency. Then considering the modifications along the columns as well, it appears that the 'efficiency' of τω to provide intermittency depends on Lv.

Finally, let us focus on the magnetic fields. To illustrate model 1, we take the following simple divergence-free ordered magnetic field

B0 = x, ( 50
)
where x is the unit vector along the x-direction. This field is then distorted according to (42), as shown in figure B7 for three values of the magnetic intermittency parameter τm. The choice of representation for the arrows is the same as in figure B5. In the left panel τm is very small, such that we almost recover the ordered field (50). Then as τm is increased, the vector field is randomized but it is still clear by eye that its statistics are biased, with a tendency of vectors to be pointing to the right, i.e. along the underlying B0.

To illustrate the second model, we show in figure B8 iso-contours of the norm of the magnetic field (44) in the case when the Gaussian white noise j g in ( 45) is independent of ωg. We do so for four different values of the magnetic Hurst exponent hm, to give a feeling for how this parameter modifies the nature of the field. Just like the Hurst exponent in fractional Brownian motion is a measure of persistence in a time series, it appears visually here that hm enables us to tweak the 'spatial persistence' of the field, i.e. make it smoother or rougher. Note that, although the hydrodynamics literature that we quote (which is very rigorous) identifies a range of validity for the Hurst exponents (notably on hv, cf. discussion above), in order to explore the degrees of freedom provided by our formulae, we also considered values beyond these recommended ranges. For example, in the two bottom panels of figure B8, we use a magnetic Hurst parameter hm smaller than 3/4, which does not belong to the range [3/4, 5/4] suggested for the hydrodynamical parameter hv at the beginning of this section. As can be seen in figure B8, this does result in interestingly looking magnetic fields. In any case, exploring the parameter space freely is even more justified here as we consider the parameters as effective parameters: they are meant to fit astrophysical data, themselves being the result of highly complex dynamics, not captured by our simple initial set of MHD equations, such that ultimately the best fit parameters may not respect the aforementioned priors, which were derived in well-defined, restricted, theoretical frameworks.

For the third model, to exemplify the intertwining of the large-scale (ordered) feature with the smaller-scale turbulence, let us consider the following simple case. Obviously we do not claim that it is realistic, but its simplicity helps the visualization. We leave the modeling with more astrophysically relevant ordered fields for future work. Hence, as an ordered current density field, let us consider a uniform current density wire of radius r0 and norm j0, i.e.

j 0 = j0 ẑ for r r0, 0 for r > r0, (51)
where r 2 ≡ x 2 + y 2 . Pleasingly, in this case the Biot-Savart law (43) can be integrated, so that we know that the corresponding ordered magnetic field is

B0 = µ0j0 2 (-yx + xŷ) 1 for r r0, (r0/r) 2 for r > r0. (52)
In other words, the vectors B0 form concentric circles around the ẑ axis, with a norm increasing linearly with r below radius r0 and decreasing as r -1 beyond r0. In figure B9 we show 3D plots of examples of the turbulent magnetic field (46) using the ordered current density (51), for three different values of the magnetic intermittency parameter τm. This figure is meant to give a visual impression of how the large-scale features of the ordered magnetic field and the turbulence-related scales are intertwined in this model. On the left, τm is very small, such that we recover the ordered magnetic field (52), as expected in the absence of significant distortion. Then as we increase τm, the vectors are more and more randomized. These plots show qualitatively how the statistics of the magnetic field are impacted. The red region in the left panel corresponds to a cylinder of radius r0 (because this is where B0 is maximal, cf ( 52)), while in the middle panel for example we can see that the red region is made of lumps. Their typical size basically (in fact there is a non-trivial dependence on τm) corresponds to the integral scale Lm. Hence we may visualize how the turbulent velocity field distorts the ordered magnetic field B0 and modifies its large-scale properties as τm increases. The right panel shows how this information even gets lost at some point once the distortion becomes really strong and one may discern nothing but eddies.

CONCLUSION & PROSPECTS

As an extension to a recent phenomenology of 3D hydrodynamical turbulence based on the so-called multiplicative chaos, we proposed three models for MHD turbulence. These provide compact, physically motivated formulae for a scalar field, a velocity field and three magnetic fields, which are by construction correlated to one another, because they all derive from the same white noise. Compared to previous works who focused on the local structure of turbulence, we also introduced large-scale features in these fields, to enlarge their possible range of use. Given the successes exhibited in the hydrodynamics literature concerning the velocity field thus constructed, we have hope that our magnetic field formulae will be useful (i) for data analyses, thanks to the free parameters which provide as many fitting degrees of freedom, (ii) for numerical simulations, providing at low computational cost fairly realistic fields to initialize simulation runs, (iii) as building blocks for analytical toy models of astrophysical environments, and (iv) from a deeper theoretical point of view, the phenomenological steps used here may contribute to improve our understanding of MHD turbulence per se. We will soon make publicly available a Python code to generate realizations of the fields presented in this paper. Written with the purpose of being user-friendly, it will also include scripts to conveniently analyze and view the 3D fields using Mayavi.

The present work may be improved and extended in numerous ways. In future works, we will analyze both analytically and numerically the statistics of the magnetic fields introduced here, in the line of Pereira et al. (2016). Then, to better appreciate the scope of these formulae, we will explore more thoroughly the parameter space rather than only a few examples as in section 5. It will also be informative to inspect the various regimes of the dynamics. We will in particular look for conditions under which an inverse energy cascade occurs. Earlier, we stressed that it is important to systematically use the same white noise throughout the formulae, because it was shown to guarantee the existence of a direct cascade in the velocity field. Nevertheless, it remains to assess what happens in the magnetized case. Considering different white noises inside the various formulae, the properties of energy transfers through scales will vary depending on the assumed correlations between the noises. Furthermore, from a practical viewpoint for data analysis, this would provide an additional degree of freedom. Another prospect is to compare the properties of our analytical expressions to direct numerical simulations, notably the correlations between the velocity and magnetic fields, but also regarding the shape of the dissipation regions. Indeed, folded structures (filaments, ribbons, current sheets) naturally form in MHD: In our model considering anisotropic integration regions as well as anisotropic regularization norms (which was taken isotropic here for simplicity) could help generating complex shapes, in addition to the anisotropy introduced by the large-scale ordered fields in the formulae.

Another important task will consist in exhibiting relations between the various free parameters (intermittency, Hurst, and regularization) introduced in our models. In the hydrodynamics literature the intermittency parameter τv is also introduced as a free parameter, and to the best of our knowledge, there is no general, simple, result linking it to the other parameters. So far we have found only results derived in special cases, either under simplifying assumptions (one-dimensional fields, asymptotic computations, low order expansions) or valid only for restricted ranges of the parameter space. For example some relations between these parameters and the Reynolds number may be found in Chevillard & Meneveau (2007); [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]; Pereira et al. (2016). Modeling the transition between viscous and inertial-range is also an important step for this task, which is explored in the multifractal formalism in Meneveau (1996) as well as in Chevillard et al. (2012).

In parallel, it is important to improve our physical understanding of the present phenomenology, and crucial to compare it to others (e.g. Sridhar & Goldreich 1994;[START_REF] Frisch | Turbulence. The legacy of A[END_REF]Politano & Pouquet 1995;[START_REF] Nazarenko | Wave Turbulence[END_REF]Malara et al. 2016;Mallet & Schekochihin 2017). This may possibly require delving back into the references to the mathematics literature that we mentioned in this paper, to acquire a deeper theoretical knowledge about the random measures we handle here. For instance exp(τm Gg) dW is the product of two random distributions, which is a subtle mathematical object per se.

Finally, after surveying the scope and limitations of our description, one may add more ingredients to improve it. The on-going developments in the hydrodynamics community (e.g. Pereira et al. 2018, who study further the dynamics of the velocity gradient) remain a wealth of information. Otherwise, so far we have introduced large-scale features by inputting some ordered magnetic and current density fields, but another important large-scale aspect to take into account are the boundaries of the system. A promising starting point to extend our work in that respect is the Biot-Savart operator of a bounded domain, as studied in Enciso et al. (2018) for example. On the longer term, we will also propose generalizations to compressible turbulence, and go beyond the kinematic regime, to analyse magnetic field dynamos in this formalism.

measure. Here we adopt a less rigorous approach, just like the Dirac δ is commonly handled as a regular function in physics, so we stick to a classical, Riemannian, view of integration. Thirdly, in the literature quoted in the present paper, both notations dW (y) and9 W (dV ) can be found for the infinitesimal element inside stochastic integrals. Because integration is performed over 3D space, the variable y replaces time t of the above one-dimensional example. These two notations are equivalent, but we may give them different interpretations. The first one, dW (y), may be simply interpreted as a field of differential elements with randomly chosen values: to each position y in space we associate three values dWi drawn independently following a Gaussian law. The other notation, W (dV ), may be seen as applying the random process W to the differential element dV at each position y.

In the present paper, as physicists, we choose the notation of (A2), because we consider that overall it makes the method presented in this paper more intuitive and easier to understand. But this is at the cost of introducing an inconvenience in terms of dimensionality. Indeed, for example in Robert & Vargas (2008) the authors, who work in a d-dimensional space and have some lengthscale R in their problem, make their initial Gaussian field dimensionless by introducing a factor R d/2 in their kernel. It is so because the measure dW has dimension R d/2 in R d . Hence in fact, to recover this dimensionality, we should replace in our integrals terms like η g dV by η g √ dV (or the more rigorous dW ), but it is less customary. Actually, many authors (e.g. [START_REF] Rodean | Stochastic Lagrangian Models of Turbulent Diffusion[END_REF]) explicitly link the Wiener process to white noise through dW = η dt, rather than dW = η √ dt as in [START_REF] Chevillard | [END_REF], 2007) for example. In any case, one should bear in mind that these notations have the precise meaning briefly outlined above.

Numerically, in order to build the Gaussian white noise vector field in a box of length B with N collocation points in each direction, one may either consider W (dV ) and generate 3N 3 independent realizations of a Gaussian variable with mean equal to zero and variance equal to dV = ( B /N ) 3 , as in Chevillard et al. (2010); Pereira et al. (2016), or consider dW = η g √ dV , where η g is a normal distribution with mean equal to zero and variance equal to 1, as we choose to proceed in this work. For full details on numerical methods for solving stochastic differential equations, see [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF].

APPENDIX B: FIGURES

In order not to interrupt the reading of the paper with large figures, we gather all the figures in this appendix.

This paper has been typeset from a T E X/L A T E X file prepared by the author. 9) from the same white noise, but varying the integration region Rs. As illustrated in the bottom line, from left to right Rs has the shape of respectively a sphere, a cigare and a pancake. It appears visually that these yield respectively cloud-like, filamentary-like and sheet-like structures, matching qualitatively to the underlying shape of the region. This figure is the only part of the paper where Rs is not spherical. Everywhere else, we use (49) i.e. spheres of radii equal to the injection scales Ls, Lv of Lm. 31) varies with the integral scale Lv and the intermittency parameter τω. Each panel in this array is a 2D slice of 3D velocity fields given by (31), built using the same white noise realization (namely the one shown in figure B1), with the values Lv = (0.08, 0.15, 0.5) (in fraction of the box size) and τω = (0.4, 2.5, 5). Colors indicate the magnitude of vectors, using the light spectrum: red corresponds to the highest velocities and blue to the lowest ones. In 3D all arrows have the same length, so that in these 2D slices the arrows that are small are in fact almost perpendicular to the plane. This way, looking at these 2D representations still gives a feeling of the 3D dynamics. This is the companion figure of figure B6. Figure B8. 3D iso-contours of the norm of the magnetic field given by our second model (44), built using the same white noise realization ωg (namely the one shown in figure B1), in the case when the Gaussian white noise j g is independent of ωg. The left column corresponds to iso-contours in which the norm is equal to the mean minus 1.5 times the root mean square of the field, in the middle column it is equal to the mean, and in the right it is equal to the mean plus 1.5 times the root mean square. In these examples we fix (Lv, τv, Lω, τm, Lm) = (0.5, 2.5, 0.01, 2.5, 0.5), and the Hurst parameters (hv, hω) = (13/12, 7/4), while varying hm: as we go from the top to the bottom row, hm equals 1, 3/4, 1/2 and then 1/4. As we vary hm the field varies smoothly between the configurations shown. 46) using the toy model ( 51) for the ordered current density j 0 , built using the same white noise realization (namely the one shown in figure B1), with from left to right the values τm = (0.01, 0.1, 0.4). The upper left half of each cube is truncated for better visualization, and the choice of representation for the arrows is the same as in figure B5. Lv and Lm are both taken equal to a third of the size of the box.

  Physically, this introduces correlation lengths, since the region Rv controls what points in space (i.e. what y in the integrand, because these are convolutions with r = xy) influence the dynamics at position x. This phenomenologically models the injection scale of the turbulent cascade operating in a given fluid state. If one takes an anisotropic region, several lengthscales are introduced. Points are correlated up to different distances in the various directions.

Figure B1 .

 B1 Figure B1. All the formulae in this paper originate from a single white noise ωg. Therefore to generate the scalar fields, velocity fields and magnetic fields of figures B2, B3, B5, B7, B9 and B8 we used the same realization of a Gaussian random vector field, namely the one shown in this figure. On the left panel, we show a 2D slice of this 3D realization, and its norm on the right.

Figure B2 .

 B2 Figure B2. Three scalar fields s built with formula (9) from the same white noise, but varying the integration region Rs. As illustrated in the bottom line, from left to right Rs has the shape of respectively a sphere, a cigare and a pancake. It appears visually that these yield respectively cloud-like, filamentary-like and sheet-like structures, matching qualitatively to the underlying shape of the region. This figure is the only part of the paper where Rs is not spherical. Everywhere else, we use (49) i.e. spheres of radii equal to the injection scales Ls, Lv of Lm.

Figure B3 .

 B3 Figure B3. This figure shows qualitatively how the scalar field (9) varies with the integral scale Ls and the intermittency parameter τs.Each panel in this array is a 2D slice of 3D scalar fields given by (9), built using the same white noise realization (namely the one shown in figureB1), with the values Ls = (0.08, 0.15, 0.5) (in fraction of the box size) and τs = (0.1, 0.8, 1.6). This is the companion figure of figureB4.

Figure B4 .

 B4 Figure B4. PDFs of the increments δ s(x) ≡ s(x + ) -s(x) of the scalar field (9) for three different lags = 5, 17 and 30 (in fraction of the size of the box), averaged over several white noise realizations. The parameters used in the panels of this array correspond to those used in its companion figure, figure B3, namely Ls = (0.08, 0.15, 0.5) (in fraction of the box size) and τs = (0.1, 0.8, 1.6). This is the companion figure of figure B3.

Figure B5 .

 B5 Figure B5. This figure shows qualitatively how the velocity field (31) varies with the integral scale Lv and the intermittency parameter τω. Each panel in this array is a 2D slice of 3D velocity fields given by (31), built using the same white noise realization (namely the one shown in figureB1), with the values Lv = (0.08, 0.15, 0.5) (in fraction of the box size) and τω = (0.4, 2.5, 5). Colors indicate the magnitude of vectors, using the light spectrum: red corresponds to the highest velocities and blue to the lowest ones. In 3D all arrows have the same length, so that in these 2D slices the arrows that are small are in fact almost perpendicular to the plane. This way, looking at these 2D representations still gives a feeling of the 3D dynamics. This is the companion figure of figure B6.

Figure B6 .

 B6 Figure B6. PDFs of the longitudinal increments of the velocity field (31) for three different lags = 5, 17 and 30 (in fraction of the size of the box), averaged over several white noise realizations. The parameters used in the panels of this array correspond to those used in its companion figure, figure B5, namely Lv = (0.08, 0.15, 0.5) (in fraction of the box size) and τω = (0.4, 2.5, 5). This is the companion figure of figure B5.

Figure B7 .

 B7 Figure B7. These panels show 2D slices of 3D magnetic fields given by our first model (42) using the ordered magnetic field B 0 model (50), built using the same white noise realization (namely the one shown in figureB1), with from left to right the values τm = (0.2, 0.8, 1.7). The choice of representation for the arrows is the same as in figureB5. Lv and Lm are both taken equal to a tenth of the size of the box.

Figure B9 .

 B9 Figure B9. 3D views of magnetic fields given by our third model (46) using the toy model (51) for the ordered current density j 0 , built using the same white noise realization (namely the one shown in figureB1), with from left to right the values τm = (0.01, 0.1, 0.4). The upper left half of each cube is truncated for better visualization, and the choice of representation for the arrows is the same as in figureB5. Lv and Lm are both taken equal to a third of the size of the box.

The closest to our work we found isHater et al. (2011), who also have as a starting point ideas from L.Chevillard and his collaborators, but they generalize to magnetic fields other aspects than we do, namely the time evolution of Lagrangian fluctuations.

MNRAS 000, 1-15 (2020)

In order not to interrupt the reading of the paper with large figures, we gather all the figures at the end, in appendix B.

Which has a logarithmic covariance, i.e. E[ sg(x) sg(y)] ∼ ln(L/|x -y|), with L the integral length scale and E the expectation value.MNRAS 000, 1-15 (2020)

Or at least let us attempt the following intuitive wording, since so far we have not found much discussion about this hypothesis in the quoted literature.

Stretching itself is a non-linear process, and is what accounts for MNRAS 000, 1-15 (2020)

We do not say that dynamo-produced fields will necessarily tend to be force free. What we do is build formulae taking inspiration from force free dynamics. Then the complex intertwining of the free parameters (large-scale cut-offs of the integration regions which may be anisotropic, and the various intermittency, regularization and Hurst parameters) may make these synthetic fields mimic fields which were the result of far more complicated dynamics than the ingredients used to derive our simple expressions. This is the same idea as in the hydrodynamical case, in which the expression for Dg is inspired from the inviscid dynamics, while in the end the result impressively describes flows with dissipation (cf section 3.4), thanks to the phenomenological parameter .MNRAS 000, 1-15 (2020)

Which is not required (cf. the similar discussion on ω at the end of section 3.3), but it is noteworthy, and the impact of considering a divergence-free j or not has to be investigated.MNRAS 000, 1-15 (2020)

This is similar to the fact that in a simple random walk, the typical distance from the origin after n time steps is proportional to √ n.MNRAS 000, 1-15 (2020)

The authors use dy for the volume element, but we prefer dV .MNRAS 000, 1-15 (2020)

ACKNOWLEDGMENTS

We express gratitude to the referee for having read our paper so carefully and for his constructive comments. We thank L. Chevillard for pleasant and fruitful discussions, as well as J.-S. Carrière, R. Pereira, B. Regaldo, T. Richard, and A. Seta for help with the coding part of this work. This research is supported by the Agence Nationale de la Recherche (project BxB: ANR-17-CE31-0022).

APPENDIX A: TECHNICAL ASPECTS

In this appendix, we detail some technical points about stochastic calculus, to clarify and justify our choices of notations.

For illustration, let us consider a basic example of a stochastic differential equation, namely a one dimensional Itô diffusion. In mathematicians' wording, it consists of a time dependent process Xt satisfying dXt = b(Xt) dt + σ(Xt) dWt.

(A1)

In the right hand side, the first term is called the drift term with amplitude b, and corresponds to the deterministic part.

The second term on the contrary is stochastic. It is a noise term with amplitude σ, where Wt is called a Wiener process, and has the following properties. Firstly, while Wt is continuous (almost) everywhere, it is differentiable (almost) nowhere. It is not a regular function, but a distribution. This is why in principle one cannot simply 'divide by dt' the differential relation (A1), but needs to generalize the notion of differentiation, in the framework of distributions. 'Gaussian white noise' corresponds to the (generalized) derivative of the Wiener process. However, it is often convenient to stick to classical notations, and another common way of writing (A1) is

notably in physics (e.g. Vanden Eijnden & Balescu 1996; Bouchet et al. 2018). A historical example in which such a noise term appears is Langevin's equation to model the motion of a particle in a fluid where collisions with the surrounding molecules are modeled as a stochastic force term by means of the white noise η.

Secondly, the noise term in (A1) is a random variable, so it should be thought of in terms of its statistical properties. Now, since Wt is a Gaussian process, it is completely specified by its mean and covariance, respectively E[Wt] = 0 and E[W t Wt] = min(t , t), where E is the expectation value. It can be shown that the covariance can also be specified via E[(W t -Wt) 2 ] = |t -t|. As t moves closer to t, this indicates that dW has a variance equal to the differential element dt. In that sense dW behaves as 8 √ dt. Hence, in the right hand side of (A1), while the deterministic term is an infinitesimal element of order 1, the noise term is an infinitesimal element of order 1/2. Therefore, in Taylor expansions, the term dW 2 is of first order, like dt, and not of second order, like dt 2 (cf. Itô's lemma). For this reason, rules in stochastic calculus, as basic as the chain rule, differ from 'usual' calculus and one has to be cautious with some habits when performing analytical calculations.

In the literature, various notations and terminologies can be found for the above noise term. Firstly, while some call Wt a Wiener process in honor of the mathematician N.Wiener, others refer to it as Brownian motion and write it Bt, in honor of the botanist R.Brown. Secondly, stochastic integrals should in principle be understood in the framework of measure theory, by means of Lebesgue integrals, so mathematically oriented authors call dW a Gaussian white