Dear Editor,

Statements that mean and maximal lifespan could increase to 150 years and more, and even without any limit, are not breaking news, as they have been claimed for decades and centuries. For instance, the French mathematician Condorcet [START_REF] De Condorcet | Esquisse d'un tableau des progrès historiques de l'esprit humain (Outlines of an historical view of the progress of the human mind)[END_REF] wrote in 1795 that, in the future, lifespan could continuously increase before reaching, or not, an asymptotic limit. Similarly, the French (formerly Russian) biologist Métalnikov [START_REF] Métalnikov | La lutte contre la mort[END_REF] wrote in 1937 that people could live well beyond 150 years, provided living conditions would be optimal. Robine and Hermann [START_REF] Robine | Maximal human lifespan[END_REF] have reported such claims, including those that endorsed the idea that Patriarchs of the Bible lived for centuries. Nowadays, de Grey [START_REF] De Grey | Foreseeable and more distant rejuvenation therapies[END_REF] has claimed that people born in 2000 could live for 5,000 years, well above the assumed lifespan of Patriarchs. The idea, not to say the hope, that lifespan coud increase without any clear limit, or so far a limit that it is difficult to envisage it, has given rise to many attempts to increase lifespan (for a review see [START_REF] Stambler | History of life-extensionism[END_REF]). Today, only some biogerontologists explicitly accept that lifespan could increase with no or a remote limit, as Guarante and Kenyon [START_REF] Guarente | Genetic pathways that regulate ageing in model organisms[END_REF] who wrote that "we begin to think of ageing as a disease that can be cured". The idea that lifespan could increase at will is let to transhumanists (e.g., [START_REF] Kurzweil | Bridges to life[END_REF]) and to biotechnology companies thriving on the hope to defeat ageing [START_REF] Brooker | Google ventures and the search for immortality[END_REF].

However, there is a long-lasting debate among demographers and biogerontologists on mean and maximal lifespan increases in the future. The 1997 collective book Longevity: to the limits and beyond [START_REF] Robine | Longevity: to the Limits and Beyond[END_REF] confronted the points of view of scholars arguing either that there are practical limits to life expectancy in France [START_REF] Olshansky | Practical limits to life expectancy in France[END_REF] or that the average French baby could live to 100 [START_REF] Vaupel | The average French baby may live 95 or 100 years[END_REF]. In subsequent years, Oeppen and Vaupel [START_REF] Oeppen | Broken limits to life expectancy[END_REF] wrote that country "record life expectancy will reach 100 in about six decades" and Christensen et al. [START_REF] Christensen | Ageing populations: the challenges ahead[END_REF] that the median lifespan of the 2007 Japanese birth cohort could be 107 years, and thus probably longer for women. Recently, Vaupel et al. [START_REF] Vaupel | Demographic perspectives on the rise of longevity[END_REF] wrote that "if the current pace of progress in life expectancy continues, most children born this millenium will celebrate their 100 th birthday", which should probably be understood as "the beginning of the millenium" because, obviously, forecasting what would occur in several centuries from now has not any heuristic value.

Beside forecasts of life expectancy, one can wonder what is the maximal lifespan, i.e., the lifespan of the last survivor (Maximal Reported Age at Death, MRAD), and if the lifespan of the record-woman Jeanne Calment, 122 years and 164 days, will remain the MRAD. Estimating MRAD requires to know mortality rates at extreme ages. It has been proposed that these rates fit to the wellknown Gompertz trajectory, i.e., they increase with age (e.g., [START_REF] Gavrilov | New trend in old-age mortality: gompertzialization of mortality trajectory[END_REF]), or that a plateau is observed above 100 or 105 years of age, with a yearly 50% mortality rate (e.g., [START_REF] Gampe | Mortality of supercentenarians: estimates from the updated IDL[END_REF]). In the first case, MRAD is reached when the death rate is 100%, the last survivor dying in this last time-interval. However, depending on the position and steepness of the slope, which could differ among populations or periods, this limit can be reached at various ages: the inescapable conclusion is that a fixed limit, say 126 years and not a single day more, cannot be practically established even if it theoretically exists.

In the second case, MRAD is depending on the size of the cohort: with a 50% yearly death rate above 100 years of age [START_REF] Gampe | Mortality of supercentenarians: estimates from the updated IDL[END_REF] only one among one million centenarians is expected to live to 120. In other words, there is no theoretical limit to lifespan and the MRAD can slightly increase with the size of the cohort (e.g., 122 years if the number of centenarians is 4 million). Obviously, if the yearly mortality rate at, say, 105 years, is 50%, it follows that MRAD will be higher if mortality rates are constant than if they increase with age.

The issue is that it is difficult to accurately estimate death rates at extreme ages, because the number of subjects is so low that confidence intervals are very high, and hot debates can occur (e.g., [START_REF] Vaupel | Exceptional lifespans[END_REF]) without any firm conclusion. For instance, Gampe [START_REF] Gampe | Mortality of supercentenarians: estimates from the updated IDL[END_REF], studying US supercentenarians (people reaching 110 years of age), wrote that "uncertainty for the annual probabilities of death increases quickly, and is large after age 113, at which point the data are too sparse to allow for an accurate assessment", but nevertheless concluded that "the shape of the survival curve suggests an exponential model, implying a constant hazard". Yet, the Figure 3.2. in [START_REF] Gampe | Mortality of supercentenarians: estimates from the updated IDL[END_REF] could also suggest that death rates slightly increase with age from 110 to 115 years of age, before a decreasing is observed. Because the 95% confidence intervals are so high, covering nearly the 0-100% mortality range at 113 years of age and all the range above this age, one could conclude anything from these results. Similar conclusions could be reached for the European supercentenarians (Figure 3.3 in [START_REF] Gampe | Mortality of supercentenarians: estimates from the updated IDL[END_REF]).

However, MRAD estimates are not really different whether one relies on constant or increasing with age death rates. Robine and Hermann [START_REF] Robine | Maximal human lifespan[END_REF] showed, using the 2013-2070 French forecasts, that in coming decades MRAD would be around 115 years in both cases, the estimate being slightly lower with a Gompertz model. Similar results were observed for the whole planet, MRAD being around 120 years, with a 124 years extreme value. In other words, no matter the model, MRAD will be around 115-120 years in coming decades, which is similar to or slightly higher than current values.

These results also show that increasing the sample size from 65 million to 7.5 billion inhabitants only increases MRAD by no more than 5 years.

Ten years ago, I wrote in a review [START_REF] Bourg | Book review[END_REF] of the book Supercentenarians [19] that "one may regret that this book on supercentenarians did not go beyond the beginning of this debate on maximal longevity". Today, this work has been done, even if the law governing death rates at extreme ages remains intractable, and nearly no authors conclude that maximal longevity will reach 150, 200 or higher values, even if the number of centenarians strongly increases in coming decades. In such conditions, and because it is possible to answer the final question of Vaupel [START_REF] Vaupel | Exceptional lifespans[END_REF], "how long can we live?", by saying that the observed maximal lifespan is 122 years, and maybe slightly longer in next decades [START_REF] Robine | Maximal human lifespan[END_REF], one may argue that time is ripe to relax efforts in searching the best model predicting MRAD, because both models provide the same answer: MRAD is ca 120 years. One can also argue that knowing whether MRAD will be 120, 125, or even 130 years in next decades has no practical importance, because only a very few persons will approach these ages. Thus, there is no concern about the "costs" of these very few supercentenarians. From a theoretical point of view, there is also no actual importance of knowing the best death rate model at the highest ages. As the number of last survivors is obviously very low, estimating mortality rates at these ages will always face the same problems as today: a too low number of subjects. There are other biological traits that appear to have a limit, even if it is not precisely known. Body size has increased in the last decades up to a plateau [START_REF] Marck | Are we reaching the limits of Homo sapiens ?[END_REF], and a very few people reach extreme heights, the world record being 2.72 m [START_REF] De Herder | Acromegaly and gigantism in the medical literature. Case descriptions in the era before and the early years after the initial publication of Pierre Marie (1886)[END_REF] but, as far as I know, there is no passionate debate to know whether there is a limit to human size. Maybe because the valuable question is modulating the growth of children, by increasing or decreasing it, and not to know whether, one day, a new height record will be established.

Therefore, thanks to the work of demographers and biogerontologists, we know how modelling mortality rates at very old ages -two models equally fit -and it is maybe not essential to spend the next ten years in trying to separate these two models. By contrast, knowing whether, before these extreme ages, mortality rates are increasing with age or are constant is useful. Gavrilov et al. [START_REF] Gavrilov | The future of human longevity[END_REF] reported that the death rates of US centenarians have not decreased since the 1940s, in contrast to what is observed for octogenarians and nonagenarians. Does it imply that the death rates of centenarians cannot decrease or is it a transient phenomenon, like was the stagnation of lifespan during the 1970s [START_REF] Gavrilov | Human life span stopped increasing: why?[END_REF]? If death rates are actually constant in centenarians with a yearly 50% value, as hypothesised for supercentenarians [START_REF] Gampe | Mortality of supercentenarians: estimates from the updated IDL[END_REF], no medical progress is expected to occur at these ages. If they increase with age, one could hope to lower the curve of death rates, as done for younger ages, even if its slope cannot be modified. Similarly, selecting the best forecast model of life expectancy at birth in coming decades ( [START_REF] Vaupel | Demographic perspectives on the rise of longevity[END_REF] but see [START_REF] Bourg | New life expectancy forecasts are too optimistic[END_REF]), and thus knowing the number of old and very old people in each country, is of importance to better manage our societies in coming decades.

All these topics are of importance and need the attention of demographers and biogerontologists, but one can surely afford not to use our working time in coming years to separate the two models of mortality rates in supercentenarians, as they roughly provide the same MRAD. As a proverb stands, "It doesn't matter if a cat is black or white, as long as it catches mice": in supercentenarians, MRAD stands for the mice, and the Gompertz and the constant mortality models

Acknowledgements

Thanks are due to Jean-Marie Robine for his helpful comments on a previous draft of this article.

are the black and white cats.


Jean-Baptiste Durrive

Katia Ferrière

Pierre Lesaffre

An intuitive parametric model for 3D compressible hydrodynamical and MHD turbulence

Keywords: turbulence -ISM, magnetic fields -MHD -methods, analytical

), following the same procedure based on the concept of multiplicative chaos. Our model provides a complementary tool to numerical simulations, as it enables us to generate very quickly fairly realistic velocity fields and magnetic fields, the statistics of which are controllable with intuitive parameters. Therefore our model will also provide a useful tool for observers in astrophysics. Finally, maybe even more than the model itself, it is the very procedure that matters the most: our model is modular, in the sense that it is constructed gradually, with intuitive and physically motivated steps, so that it is prone to many further improvements.

INTRODUCTION

Fluids in nature, magnetized or not, are most often in a turbulent state, up to the largest cosmological scales and earliest epochs of the Universe (e.g. Subramanian 2019). Hence, in order to describe natural phenomena, from laboratory to astrophysical environments, it is crucial to be able to model turbulence. But the difficulty to do so precisely is well-known: the solutions of the Navier-Stokes equations and the equations of magnetohydrodynamics (MHD) are still far from reach analytically, and generating realistic data by means of numerical simulations is very costly, both in terms of time and of resources.

In the quest for finding explicit solutions to the incompressible Navier-Stokes equations, Chevillard et al. (2010Chevillard et al. ( , 2011Chevillard et al. ( , 2012Chevillard et al. ( , 2013)); [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]; Pereira et al. (2016Pereira et al. ( , 2018) ) built a powerful (i.e. concise and yet efficient) model for incompressible hydrodynamical turbulence, based on the concept of Gaussian multiplicative chaos. This con-⋆ E-mail:jdurrive@protonmail.com cept was first introduced by the mathematician J.-P. Kahane in Kahane (1985). It is closely related to fractals (e.g. [START_REF] Mandelbrot | Possible refinement of the lognormal hypothesis concerning the distribution of energy dissipation in intermittent turbulence[END_REF], which is why it is a natural tool to model processes comprising self-similarity, such as turbulence. It is suited for many other applications, but our focus here is (magneto)fluid dynamics in view of modeling astrophysical environments. For a recent review on Gaussian multiplicative chaos, see for example Rhodes & Vargas (2014). Mathematically speaking, Gaussian multiplicative chaos consists in considering the exponential of a Gaussian field with logarithmic covariance. Simply said, what lies at the heart of our model is the fact that, while the sum (or more generally a linear operation) of Gaussian processes is a Gaussian process, the product of Gaussian processes is non-Gaussian, i.e. intermittent. Hence the terminology 'multiplicative' in Gaussian multiplicative chaos, and the choice of the nonlinear operation that the exponential is.

Based on the aforementioned series of papers, we proposed in Durrive et al. (2020) three turbulence models for the magnetized case, i.e. for three dimensional MHD turbulence. As in Chevillard et al. (2010), we limited ourselves to the incompressible case. In the present paper, we propose generalizations of this method to include compressibility, for both unmagnetized and magnetized fluids.

Our objective follows that of Durrive et al. (2020): we try to adapt tools developed in the hydrodynamics and MHD communities to the needs of the astrophysics community. In that respect, our priority while designing the present model, has been to try and make it as practical as possible.

Firstly, given its construction, this model should, hopefully, be useful for both numerically-oriented and observationnally-oriented researchers. Indeed, on the one hand, one can use this model to initialize numerical simulations, with quickly generated non-trivial fields. It may also serve as prescriptions for fields at unresolved scales. On the other hand, thanks to the free parameters controlling the properties of the fields, our model may be used to characterize statistically data (stemming from either observations or simulations) or to assess statistical error bars by generating hundreds of realizations. We are currently implementing numerically this model, and in a close future we will make a public release of our code. A crucial feature of this approach is that it is very efficient numerically, because most of the formulae are nothing but convolutions, i.e. simple products in Fourier space. We may generate fields in a few minutes, which, with the same resources and at the same resolution, are obtained in typically several days or weeks by direct numerical simulations (DNS). Obviously DNS results remain more realistic in the details, but with our approach the statistics of the fields are easily controllable. In that sense our approach is complementary to DNS.

Secondly, we made an effort to present the model pedagogically. Hence, while a large part of the derivation is based on mathematical considerations, readers in a hurry or those only looking for a practical tool, do not need to delve into any of these technicalities. Each technical argument is accompanied with a physical meaning, so that the formulae ultimately are very intuitive. Naturally we also added numerous references to the literature, to satisfy as well readers who, on the contrary, are interested in the fundamental aspects of the model. Indeed, practicalities aside, this approach is obviously also interesting per se, from a theoretical point of view, as it takes part in the studies on how the concept of fractals and multiplicative chaos may help modeling turbulence.

Our notations will follow closely those used in Durrive et al. (2020): our work is a sequel of it, and we want to keep the description focused on physical arguments. This is done at the cost of losing some mathematical rigor. The reader caring about the technicalities hidden by these more intuitive notations is invited to have a look at the first appendix of Durrive et al. (2020). Now, in order to make the present paper as self-consistent as possible, let us recall that (i) we use everywhere the shorthand notation r ≡ xy and r ≡ |r|,

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 d, v, ω and m for quantities associated with density fields, velocity fields, vorticity fields and magnetic fields respectively, and (iii) variables with a tilde denote random variables and subscript g is added to random fields with Gaussian statistics.

The paper is organized as follows. First, in section 2, we give a brief summary of the paper, to help readers keep track of the logic behind the various steps of the derivation once they delve into the detailed calculations. Then, in section 3, we present the equations governing the dynamics of the fields of which we want to build a stochastic representation. Section 4 is dedicated to the model in the hydrodynamical case, and section 5 to the magnetized case. We conclude and give a glimpse of our prospects in section 6.

OUTLINE IN SIMPLE WORDS

To improve the presentation, let us summarize, using words only, the adopted method.

Given the observed stochastic nature of turbulence, stochastic calculus1 is a most natural formalism to describe turbulent fields. Starting from the equations of motion, we derive formulae telling us how fluid particles are being deformed over time. We then transform these (approximate) solutions into random fields by setting their initial conditions as Gaussian random fields2 . Randomness is thus implemented through the initial conditions, while non-Gaussian features emerge from the dynamical evolution.

More specifically, the model of Chevillard et al. ( 2010) for an incompressible velocity field is derived as follows. Due to incompressibility, Biot-Savart's law provides a one-to-one correspondence between velocity and vorticity, similar to the one-to-one correspondence between magnetic field (which is divergence-free) and current density. Using this fact, in this model what is constructed ('by hand') is the vorticity field3 , and the velocity field is then deduced by means of Biot-Savart's law. Accordingly the authors start by studying the dynamics of vorticity. They notice that on a short enough timescale τ (namely the correlation timescale of velocity gradients), an efficient approximate solution of the evolution equation for vorticity can be found. Thus they derive an explicit formula for how vortices are being stretched in the fluid. To complete this formula, the initial state needs to be specified: they choose it to be a Gaussian random field, in order to take into account phenomenologically the resultant of all possible fluctuations induced locally by the many surrounding eddies, by a large-scale eddy and/or by some external forcing. Biot-Savart's law then translates this vorticity model into a model for the velocity field. Hence, the velocity field is modeled as the result of the random stretching (dictated by the equation of motion) of an initially Gaussian vorticity field.

Then, the velocity field thus constructed is considered as a random field representative of the stationary state of the fluid. Indeed, τ being a correlation timescale, at each timestep τ the fluid gets deformed and then decorrelates, over and over again, such that the stationary state of the fluid ends up having similar statistics as this velocity field deformed during the time lapse τ . In addition, since the parameter τ corresponds to the time lapse during which stretching occurs, τ controls how much the velocity field deviates from its initial state. The latter being Gaussian, the parameter τ quantifies the intermittency of the velocity field (intermittency being defined as departure from Gaussianity). By simply tuning this intermittency parameter τ , this velocity field model shares many properties of experimental and numerical turbulence (cf. Chevillard et al. 2010, and the following works analyzing these fields).

The models of Durrive et al. (2020) for incompressible MHD turbulence are derived in a similar manner, but starting from the MHD equations. The main new ingredient with magnetic fields is that flux tubes are not only stretched and contracted as vortices are, but they are also being sheared by vorticity.

In the above incompressible cases, the velocity field is fully characterized by its vorticity field, which, by construction, is divergence-free and thus contains only two degrees of freedom. Therefore, in the description, one needs to fix only two degrees of freedom. Indeed, while the deformation itself is obtained from first principles (it is dictated by the equations of the dynamics), the fact that vorticity is modeled as a Gaussian field that is being deformed is a phenomenological input, which fixes these degrees of freedom. Now, in the present work we extend this approach to take compressibility into account. In this case, Biot-Savart's law, generalized to a compressible fluid, expresses the velocity field not only in terms of the vorticity field, but also in terms of the dilatation field. Hence, a third degree of freedom now comes into play, and an additional phenomenological input is required. Guided by astrophysical data, in which density fields tend to be lognormal (e.g. Levrier et al. 2018), we suggest modeling the dilatation field as a fractional Gaussian field (i.e. a Gaussian field of which the statistics are controllable with a few parameters). Having fixed this third degree of freedom, we derive (astrophysically) relevant expressions for compressible hydrodynamical turbulence, as well as MHD turbulence, from a generalized Biot-Savart law.

Finally, it is noteworthy that all the random fields that we build in this model (dilatation, density, velocity, vorticity, deformation, current density, and magnetic fields) are correlated, by construction, because we manage to build them from a single Gaussian white noise vector (the Gaussian 'seed' sg below).

DYNAMICAL EVOLUTION

We aim at constructing an analytic stochastic representation of a density field, a velocity field, a vorticity field and a magnetic field coupled together, as realistically as possible. Mathematically speaking, we are thus aiming at expliciting (approximate) solutions to the system of equations consisting of the continuity equation, momentum equation and induction equation, including viscous and resistive effects. However, a basic feature of the approach of Chevillard et al. (2010) and Durrive et al. (2020) is to start with the aforementioned equations in the limit of vanishing viscosity and resistivity. Indeed, such effects are later taken into account phenomenologically through various regularizations (cf. item (iii) of the turbulization procedure detailed in section 4.1). Therefore, the equations for the dynamics that we are now going to consider for constructing our formulae are the following.

Density field

Denoting the material derivative (i.e. evaluated while moving with the fluid) as

d dt ≡ ∂t + v • ∇, (2)
the conservation of mass of a fluid of density ρ and velocity v can be described by

dρ dt = -ϑρ, ( 3
)
where ϑ is the dilatation field,

ϑ ≡ ∇ • v. (4)
This field ϑ expresses local expansion or contraction, according to whether it is positive or negative respectively. It is equal to zero in the incompressible case of Durrive et al. (2020), and is the new key ingredient of the present compressible case.

Velocity field

The momentum evolution equation for an inviscid fluid of pressure p in a magnetic field B and gravitational field g reads

dv dt = - 1 ρ ∇p + 1 ρ j × B + g. (5)
The terms in the right-hand side correspond to the forces acting on the fluid element, namely the pressure gradient force, the Lorentz force, and the gravitational force respectively. The current density j is given by Ampère's law,

j = 1 µ0 ∇ × B, ( 6
)
where µ0 is the permeability of free space. Equation (5) identifies the physical origin of motion, falling into the category of dynamics. It is at the heart of our model, but in fact a large part of our formalism deals with the fluid deformation and motion, without concern about the causes of these deformation and motion, falling into the category of kinematics. Indeed, our starting point for modeling the velocity field is the Helmholtz decomposition,

v = vc + vr, (7)
where the first term, called the compressional component, satisfies ∇ × vc = 0 and is related to the dilatation field (4) through

vc(x) = 1 4π R 3 ϑr r 3 dV, (8)
while the second term, called the rotational component, satisfies ∇ • vr = 0 and can be written as

vr(x) = 1 4π R 3 ω × r r 3 dV, (9)
(fluid equivalent of the magnetic Biot-Savart law) where the vorticity field is defined as

ω ≡ ∇ × v. (10)
In these expressions we do not consider boundary terms, assuming the fields vanish sufficiently fast at infinity, i.e. far away from the region of interest.

Readers more familiar with electrostatics and magnetostatics than with hydrodynamics will get a feeling for where these expressions come from thanks to the following comparison. The compressional component vc is formally equivalent to the electric field, E, which satisfies Maxwell's equations ∇ × E = 0 and ∇ • E = ρe/ǫ0 (with ρe the electric charge density and ǫ0 the permittivity of free space), and which can be written as

E = 1 4π (ρe/ǫ0) r r 3 dV, ( 11
)
whereas vr is analogous to the magnetic field, B, which satisfies ∇ • B = 0 and ∇ × B = µ0j, and which can be written as the original (magnetic) Biot-Savart law

B = 1 4π (µ0j) × r r 3 dV. ( 12
)

Gradient of the velocity field

A central quantity of our model is the gradient of the velocity field

G ≡ (∇v) T , (13)
where T denotes transposition. Taking the gradient of the Helmholtz decomposition (7), one can write G as [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF])

G = D + Ω, (14)
where the rotation tensor (i.e. the antisymmetric part of G) is a cross product with vorticity, namely

Ω = 1 2 ω×, (15)
and the strain-rate tensor (i.e. the symmetric part of G) may itself be decomposed as

D = Dc + Dr, (16)
where the compressional component reads

Dc = 1 4π R 3 (r 2 I -3rr)ϑ r 5 dV, (17)
and the rotational component

Dr = - 3 8π R 3 r(ω × r) + (ω × r)r r 5 dV. ( 18
)
In these expressions products of vectors of the form xy are tensor products (xy|ij = xiyj ). As in [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF], we do not include a Cauchy principal value in the expressions of ( 17) and (18) (while it is mentioned in Durrive et al. 2020) because in the following this technicality is avoided thanks to the regularizations of the kernels. One should not be misled however by the fact that as such (17) seems to yield a traceless Dc while its trace should equal ϑ. In contrast to D, which has both a compressional and a rotational component (cf. Eq (16)), Ω (cf. Eq (15)) has only a rotational component; its compressional component, Ωc = 1/2ωc×, vanishes trivially because ω ≡ ∇ × vc = 0 (cf. above Eq (8)). In other words, compressibility affects only the strain-rate tensor, D, not the rotation tensor, Ω.

Physically, the matrix D includes vortex stretching, which is an essential kinematic mechanism. This nonlinear effect is responsible for the cascade process in turbulence, by which large-scale eddies break apart into smaller and smaller eddies which rotate faster and faster, i.e. with increasingly stronger enstrophy [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF]. Stretching is most efficient when ω is aligned with the stretching principal axes of D, at a rate given by the eigenvalues of D. Therefore this process is also related to vorticity alignments in turbulent flows. The strain rate that causes stretching can be either a background field induced by other vortices (hence the fact that the expression for D involves a spatial integral), or induced locally by the vortex itself. For a study of the local contribution (induced by a vortical structure in its neighboring vicinity) versus the nonlocal contribution (background strain rate induced in the vicinity of the structure by all the remaining vorticity) in the strain rate tensor, and its implications for vorticity alignment, see for example Hamlington et al. (2008).

The above considerations are kinematic. The dynamics of the velocity gradient, following fluid particles, is obtained by taking the gradient of the equation of motion (5), yielding

dG dt = -G 2 + ∇Π, (19)
where

Π = - 1 ρ ∇p + 1 ρ j × B + g. ( 20
)
The non-linear term G 2 in (19) induces gradient steepening and physically corresponds to self-stretching. This equation is fundamental to analyze turbulent flows. Indeed, the full velocity gradient tensor G characterizes variations of all velocity components, in all directions, which constitutes valuable information, complementary to the usual studies of the statistics and anomalous scaling of longitudinal and transverse velocity increments. A vast literature exibiting its properties is available, especially in the incompressible case (e.g. Zeff et al. 2003;[START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]Wallace 2009;Meneveau 2011) but also compressible (Suman & Girimaji 2011;Parashar et al. 2017). Equation ( 19) is not closed in terms of G because one still needs to prescribe Π. This missing additional constraint on Π is commonly referred to as a closure relation, because it enables us to end up with a closed system of equations. The choice of closure relation is a crucial task, as it impacts the whole model and determines its domain of validity, but it is also a very challenging one. In appendix A we give some examples of closure relations that can be found in the literature specialized in turbulence, which are the relations we took inspiration from for the present work. We postpone to section 4.1 the presentation and justification of the choice of closure for our model, because we need to introduce a few more elements beforehand.

Dilatation field

Taking the divergence of the equation of motion (5), or equivalently taking the trace of ( 19), yields the evolution equation for the dilatation field,

dϑ dt = -tr(G 2 ) -∇ • ∇p ρ + ∇ • j × B ρ + ∇ • g. (21)
In order to better visualize the assumptions underlying our modeling, we write the above equation in full generality, but in fact in this work we will retain only the first term in the right-hand side, i.e. the one stemming from self-stretching.

Let us now expose the reasons why we discard the other terms (for simplicity obviously, but not only), and let us discuss to what extent, a priori, this choice may restrict the domain of validity of our model.

As far as the pressure term is concerned, writing (21) in dimensionless form, it appears that it is inversely proportional to the Mach number squared. Therefore we place ourselves in the high Mach number regime.

Regarding the Lorentz force term, we discard it as a first step, planning to increase the complexity of our model gradually and leaving the task of considering magnetic feedback for future works. That being said, we have more than a simple strategic reason. Firstly it has been observed in hydrodynamics (e.g. Pelz et al. 1985) that, at least in regions of low dissipation, the velocity and vorticity vectors tend to get aligned, statistically speaking (i.e. the probability distribution function of the cosine of the angle between the two fields peaks at ±1, while uncorrelated Gaussian variables would produce a flat distribution). This phenomenon is called Beltrami correlation or Beltramization, in reference to Beltrami vector fields which are parallel to their own curl. Now, similar correlations emerge in MHD turbulence where the fluid has a tendency to produce various correlations between the fields. Among them are Alfvénic correlation, where velocity and magnetic field have a tendency to align, and force-free correlation, where the magnetic field and the electric current density tend to align (Servidio et al. 2008;Miesch et al. 2015;Seta et al. 2020). Hence, the diminishing of nonlinearity in the inertial range by alignment seems common enough to legitimize the fact that we are going to build our model neglecting the feedback from the Lorentz force in (21). Secondly, in the following we will introduce some phenomenological parameters into the dynamical equations of the present section, with the aim of widening their domain of validity. In other words, we will modify our equations a little, or 'post-process' them so to speak, following a procedure detailed in section 4.1. As a consequence, while we neglect the effect of the Lorentz force in the construction at this stage, ultimately our model may still describe to some extent non force-free MHD. This procedure already proved fruitful in the hydrodynamical case (Chevillard et al. 2010), in which the model is first constructed based on the inviscid limit, but in a second step regularization parameters are introduced, which ultimately account for dissipation.

Finally, regarding the gravitational acceleration term, it all depends on the environment one wants to model, because this term as such is not closed. For instance, if we were to model a self-gravitating medium, the relevant closure relation to join to our system of equations would be Poisson's equation ∇ • g = -4πρ. However, we limit the scope of this paper to turbulent interstellar media such as molecular clouds in which self-gravity is subdominant, i.e. not prone to gravitational instability, as is the case in many environments.

Now, since we are considering a situation where equation (21) reduces to

dϑ dt = -tr(G 2 ), (22)
in the following it will matter to have information on the statistical properties of the trace of G 2 . For this purpose, let us relate it to more familiar quantities, whose statistical properties have been previously studied. As detailed in, for example, Ishihara et al. (2009); Meneveau (2011), in smallscale hydrodynamical turbulence, the instantaneous energy dissipation rate defined as (ν is the fluid kinematic viscosity, and summation over repeated indices is implied)

ǫ = 2νDij Dij , (23)
displays large levels of intermittency, being almost lognormal. Furthermore, simulations of [START_REF] Yeung | [END_REF]; Yeung et al. (2006) show that the statistics of the pseudodissipation

ǫ ′ = νGij Gij, (24)
is even closer to log-normality. For this reason, ǫ ′ is often preferred over ǫ in Lagrangian stochastic models (e.g. Sawford 1991;Lamorgese et al. 2007). Meanwhile, the enstrophy,

ξ = ωiωi, (25)
tends to be more intermittent than ǫ (e.g. Zeff et al. 2003).

We mention these three quantities because, with the decomposition ( 14) of G, we have the identities

tr(G 2 ) = ǫ ′ ν -ξ, = 1 2 ǫ ν -ξ , (26)
and given that ǫ, ǫ ′ and ξ are intermittent, these identities suggest that tr(G 2 ) is also an intermittent field.

Vorticity field

Taking the curl of the equation of motion (5) yields the evolution equation for the vorticity field

dω dt = G • ω -ϑω + ∇ρ × ∇p ρ 2 + ∇ × j × B ρ . (27)
Then in view of the decomposition ( 14), together with (15) the first term in the right hand-side simplifies further as

G • ω = D • ω. ( 28
)
This term represents vortex stretching along the vortex axis.

In the incompressible case, stretching along the vortex axis is automatically accompanied by contraction perpendicular to the axis, and by conservation of angular momentum, this perpendicular contraction leads to an increase of the angular velocity. In the compressible case, a new term enters the picture: the term -ϑω, which represents the effect of three dimensional contraction, i.e., contraction along the vortex axis together with contraction perpendicular to the axis. Hence, the sum of the first two terms represents the effect of contraction perpendicular to the vortex axis, which again, by conservation of angular momentum, leads to an increase of the angular velocity. The other terms in the right-hand side correspond to the contribution from baroclinity and from the Lorentz force respectively. We are now going to adopt the same strategy as in Durrive et al. (2020). Indeed, we are going to build our model using simplifying assumptions, but since we will ultimately add several free phenomenological parameters (regularization lengths, Hurst parameters, arbitrary integration regions, and intermittency parameters), the obtained formulae will potentially mimic turbulent (magneto-)fluids in regimes beyond the framework they were inspired from. Sceptic readers at this stage are invited to have a look at the brief review of section 3.4 in Durrive et al. (2020), which stresses that this strategy turns out to be very fruitful in the well-studied incompressible hydrodynamical case of Chevillard et al. (2010). It is these successes that brought us in Durrive et al. (2020) to extend the method of Chevillard et al. (2010) to the magnetized case, and now to the compressible case. Hence as in ( 21), assuming that it will not be as restrictive as it seems, let us neglect the backreaction due to the Lorentz force, and consider a barotropic fluid, i.e. a fluid for which thermal pressure, p, is a function of density ρ only. Then the last two terms of ( 27) vanish and, using mass conservation (3), that equation may be simply written as the Beltrami equation [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF])

d dt ω ρ = D • ω ρ . (29)
Physically, the interpretation of this relation is that the quantity ω/ρ moves as if it were frozen into the fluid, i.e. eddies are transported with the fluid, because this equation has the same form as the equation describing the evolution of the Lagrangian separation vector between two fluid particles (e.g. Rincon 2019).

Magnetic field

The induction equation in the perfectly conducting limit reads

∂tB = ∇ × (v × B). (30)
Expanding the curl of the cross product on the right-hand side, this may be expressed using the velocity gradient and the dilatation field as

dB dt = G • B -ϑB. ( 31
)
With mass conservation (3) this can then be written in a similar way as the evolution equation for the vorticity (29), namely

d dt B ρ = G • B ρ . ( 32
)
This may also be interpreted as a frozen-in condition, but for the magnetic field divided by density this time, which corresponds to Alfvén's theorem (e.g. Rincon 2019). The main difference between ( 29) and ( 32) comes from the fact that Ω • ω = 0, while generally Ω • B = 0, so that (28) does not apply for B. Physically, a vortex cannot rotate its own axis, whereas it can rotate magnetic field lines.

Short-time dynamics

The way we wrote the equations in the previous sections is well-thought-out: The form of the evolution equations for the vorticity field and for the magnetic field, ( 29) and ( 32), most naturally generalize their incompressible counterparts, equations ( 26) and (39) of Durrive et al. (2020), respectively. Indeed, written like this, one simply needs to replace ω by ω/ρ, and B by B/ρ, to go from the incompressible to the compressible case. But this similarity is more than a simple formal convenience: it prompts us to reuse directly the same reasoning as the one used in Durrive et al. (2020) to derive the velocity and magnetic field models. This is what we are going to do, and the first step is to solve the dynamical equations as follows.

Let us consider a short time interval t -t0, where t0 is some arbitrary initial time (only the time difference will matter). By 'short' we mean that it is short enough for the velocity gradient G(t, x) to not vary much during this period, i.e. such that it remains roughly constant, equal to its value at time t0. Now since the above equations, governing the time evolution of the density, the dilatation, the vorticity and the magnetic field, are written in Lagrangian form (with the material derivative), we may integrate them as follows.

The dilatation equation ( 22) is straightforward to integrate since in this framework G is a mere constant. Hence

ϑ(t, x) = ϑ(t0) -(t -t0) tr[G(t0) 2 ]. (33)
Note that in order to improve the presentation we make the x dependencies implicit in the right-hand sides of the expressions throughout this section. Then, plugging this expression for ϑ into the mass conservation (3) yields an equation for ρ which is readily integrated. From this we find that the density field increases (or decreases) locally according to

ρ(t, x) = e -(t-t 0 )ϑ(t 0 )+ (t-t 0 ) 2 2 tr[G(t 0 ) 2 ] ρ(t0). (34)
Similarly, we may integrate equation ( 29) to obtain the vorticity field. Since the tensor D is just the symmetric part of G, which we assume equal to its value at time t0, we may directly set D to its value at t0, thereby finding

ω(t, x) = ρ(t) ρ(t0) e (t-t 0 )D(t 0 ) • ω(t0). (35)
The matrix exponential corresponds to vortex stretching, and the overall factor (the ratio of density at different times, which is absent in the incompressible case), is due to local expansion or contraction. Likewise (32) gives the short-time evolution of the magnetic field

B(t, x) = ρ(t) ρ(t0) e (t-t 0 )G(t 0 ) • B(t0). (36)
The exponential now corresponds to both stretching and also to shear. The overall factor is again a weight due to local expansion or contraction. Note that in this procedure we assume that stretching dominates over advection, as detailed in Durrive et al. (2020). We will come back to this point in section 4.4.

A MODEL FOR COMPRESSIBLE HYDRODYNAMICAL TURBULENCE

In this section we focus on hydrodynamics. Magnetic fields will be discussed afterwards, in section 5.

At the heart of our model is a procedure by which we transform expressions such as those derived in section 3 into random fields 'adapted to mimic turbulence'. Let us illustrate this vague phrase by considering the Helmholtz decomposition (7) as an example. As such, this expression does not exhibit the features we expect from a turbulent field: where is the injection scale, the dissipative scale, or even where is the randomness? Obviously, these would emerge explicitly, somehow, if we were able to completely solve the Navier-Stokes equations, since the latter embody turbulence. But as long as we are unable to do so, we only have the option of mimicking turbulence. Hence, let us try to tweak a little bit our formulae, notably the Helmholtz decomposition, in order to introduce, be it in an ad hoc manner, these important and useful features. This is the very purpose of the procedure we are about to detail: We suggest five slight changes to bring to the dynamical equations from the previous sections, in order to introduce phenomenologically features relevant to stochastic representations of some turbulent states. Now because this procedure consists of several steps, we suggest calling it for short 'turbulization', since it will make our formulae 'look more turbulent'.

First, in section 4.1, we detail what this turbulization procedure consists of, by focusing on a specific example. Then, in section 4.2, we explicit the initial state that we think is most relevant for the astrophysical context, notably of the interstellar medium. Finally, in section 4.3, we show what the full deformed state is (i.e. the result of dynamically evolving the initial state), which constitutes our model for the turbulent state.

turbulization

The ideas underlying turbulization originate from Chevillard et al. (2010). Here we formalize this procedure, as in Durrive et al. (2020), in the sense that we decompose it into five explicit steps. To improve the presentation, we illustrate these steps with the example of the compressible component of velocity, vc, which includes the dilatation field, ϑ. The expressions of vc and ϑ are given by ( 8) and (33), respectively (because the short-time expressions are what matters from now on), which we reproduce here for convenience

vc(x) = 1 4π R 3 ϑr r 3 dV, ϑ(t, x) = ϑ(t0) -(t -t0) tr[G(t0) 2 ]. ( 37
)
Let us make the five following modifications to the expressions in (37): (i) Injection scale parameter: Instead of integrating throughout the whole of three-dimensional space, we limit the integration region to the vicinity of the considered position x, i.e. we make the substitution

R 3 → R, ( 38
)
where the region R enables us to control the typical size of the region of correlation. Indeed, in order to compute the velocity vc at some position x (for simplicity, or if the turbulence is isotropic, say R is a ball of radius L centered at this position x), only the points inside R will matter because points further away than distance L, namely such that r ≡ |x -y| > L, do not enter in the integral (recall that the variable of integration is y). Of course this spatial decorrelation happens smoothly: the ball used to compute vc at a position x ′ close to x, namely such that |x ′ -x| < L, will overlap with the one used to compute vc at x, but the volume of overlap decreases as the separation |x ′ -x| increases.

Hence, this large-scale cut-off parameter L corresponds to some correlation length, and acts, in the turbulence terminology, as the injection scale, controlling the typical size of the largest eddies.

(ii) Roughness/wildness parameter: The expression of vc in (37) is the convolution of ϑ and a power-law kernel. The power of this power law is equal to -3, which stems from the fact that we are considering a three dimensional space. However, fractal geometry demonstrates that considering non-integer dimensions is a powerful tool to describe a wider range of phenomena, in particular non-smooth shapes such as those typically appearing in turbulence. This approach to modeling originates from [START_REF] Mandelbrot | [END_REF]; [START_REF] Mandelbrot | Possible refinement of the lognormal hypothesis concerning the distribution of energy dissipation in intermittent turbulence[END_REF]. In practice, for our purpose, we will make use of this tool in the simplest manner: let us simply generalize the power law of the kernel in vc to an (a priori) arbitrary value h, i.e. we make the substitution

r -3 → r -2h . ( 39
)
As is customary in the literature related to fractals, we call this free parameter the Hurst parameter4 , hence the choice of the letter 'h'. Physically, h quantifies the 'wildness' of the turbulence, in Mandelbrot's wording. Geometrically, this parameter basically controls how rough the field is, how abrupt spatial variations are. Accordingly, the power spectrum of the field is very sensitive to this parameter.

(iii) Dissipation scale parameter: As it happens, the aforementioned kernel is singular, meaning that it diverges as r goes to 0. In order to avoid difficulties due to divergences and to ensure differentiability of the fields we are about to build, we eliminate this singularity by introducing a regularization parameter η, which simply means that we make the substitution

r → (r 2 + η 2 ) 1/2 , ( 40
)
in the denominator of the kernel. Physically, the important advantage of introducing this small-scale cut-off is that we now gain a degree of freedom corresponding to dissipation.

For this reason, we use the notation 5 η by analogy with the Kolmogorov lengthscale, ηK , below which viscosity dominates and the turbulent kinetic energy is dissipated into heat. This scale is usually defined as

ηK = ν 3 ǫ 1 4 , ( 41
)
where ν is the fluid kinematic viscosity and ǫ is the mean rate of dissipation of kinetic energy per unit mass. However, beware that our parameter η is inspired from ηK , but it is not necessarily equal to it. We use it as a free parameter, which broadly takes into account, in a phenomenological manner, all the various dissipative effects not accounted for otherwise.

(iv) Randomization and large-scale features: We choose 'initial' conditions, i.e. prescriptions for the fields ρ(t0), ϑ(t0), ω(t0) and B(t0). Later in the text, we will also introduce the initial current density, j(t0), and we also need to prescribe a value for G(t0), but this will be deduced from the aforementioned fields. This fourth step brings a lot of freedom in the modeling. It may give rise to a wide class of models, depending on the choices made at this level. This versatility is an asset of our approach. We may distinguish two radically different types of choices for the initial state of each field: we may either set it to a random field, in order to introduce fluctuations in the model, or we may choose it to be a smooth field, in order to inject some large-scale features or boundary effects (for example a density profile with a sharp transition between the core and outskirt of a molecular cloud). We will detail all our choices for the initial conditions of the present model in section 4.2. But to finish illustrating turbulization with the example (37), let us mention that in this case this fourth step consists in the substitutions ϑ(t0) → ϑg, G(t0) → Gg.

(42)

We will explicit what ϑg and Gg stand for in section 4.2, and for now we simply point out that the tildes indicate that they are random fields, with Gaussian statistics as indicated by the subscript g.

(v) Intermittency parameter: The expression of vc in (37) is kinematic: it describes how the velocity field varies under a local expansion or contraction, independently of the origin of this expansion or contraction. It is a geometrical relation. On the contrary, the expression of ϑ in (37) is dynamical6 , because it corresponds to the time evolution of ϑ deduced from the equation of motion (5). By construction (cf. section 3.7), this relation holds a priori only during a time lapse of the order of the velocity gradient correlation timescale. However, assuming that after such a period of time the state of the fluid is somehow 'reset' or 'blurred' (by an external forcing or by some complicated non-linear process not modeled here) to a Gaussian state (cf. item (iv) and section 4.2), in the following time step, the fluid particles will once more be deformed according to the equations displayed in section 3.7, during another velocity gradient correlation timescale. With this mental picture in mind, we build a model for a stationary turbulent field from the above dynamical relation, that is to say we replace the time interval t -t0 by a free parameter τ (this substitution is naturally called 'stationarization [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]Durrive et al. 2020) t -t0 → τ.

(43)

Physically, τ corresponds to the time7 during which fluid particles are deformed, i.e. stretched, contracted and sheared. Since the initial state is chosen to be Gaussian (cf. item (iv) and section 4.2), the longer this deformation lasts, the more departure from Gaussianity we expect. Therefore, the parameter τ controls the intermittency of the fluid, i.e. it quantifies how non-Gaussian the field is, and is referred to as the intermittency parameter.

Conclusion: turbulization consists of the five above steps, and applied to the system of equations ( 37), it results in the random fields

vc(x) = 1 4π R ϑr (r 2 + η 2 ) h dV, ϑ(x) = ϑg -τ tr( G 2 g ).
(44) Now that we have illustrated this procedure for the example (37) → (44), we are ready to apply it to all the fields mentioned in section 3, in order to build our models for (magneto-)hydrodynamical turbulence. But before doing so, let us explicit and justify the initial conditions that we are going to choose, such as ϑg and Gg in (44).

'Initial' state (undeformed state)

The initial state represents the state of the fluid before its deformation. Therefore we will also refer to it as the 'undeformed' state. In this section, we propose an initial state that is relevant for the modeling of astrophysical environments.

Inhomogeneity often plays an important role in the turbulence of astrophysical fluids, which is why we choose the initial density field to be a smooth field, i.e. not a random field, by making the substitution

ρ(t0, x) → ρ0(x), (45)
where ρ0(x) is, for instance, the space-averaged density profile of an interstellar molecular cloud or filament, or that of a cluster of galaxies or a filamentary structure of the cosmic web. The substitution (45) provides a way to introduce ordered, large-scale features in the modeling.

Given the stochastic nature of turbulence, it is necessary to introduce randomness at some point as well. To do so, consider two random fields, ϑg and ωg (tildes indicate that they are random). We randomize the initial dilatation field through the substitution

ϑ(t0, x) → ϑg(x), (46)
as well as the initial vorticity field, through ω(t0, x) → ωg(x).

(47)

Let us now expose the reasoning that leads us to choose relevant explicit expressions for these two random fields.

First of all, we decide that these random fields should have Gaussian statistics, and we indicate this property with the subscripts g. We make this decision based on several previous works that show that modeling turbulent fields as dynamically deformed Gaussian fields can be remarkably fruitful. Examples of such works can be found in Appendix A, dedicated to how the question of closure is tackled in the literature.

Then, the simplest option would be to take the initial conditions ϑ(t0, x) and ω(t0, x) as Gaussian white noises. In fact, this is precisely what is done in Chevillard et al. (2010); Durrive et al. (2020), at least for ω(t0), since there is no dilatation field in their incompressible framework. The most evident reason for this is practical: it has the advantage of facilitating analytical calculations of the statistics of the process thus defined (process that we will fully explicit in section 4.3), as performed for instance by Pereira et al. (2016). In addition, this choice is especially natural since Gaussian white noises are structureless (they are 'pristine', in a way), and the whole point of the short-time expansion exposed in section 3.7 is to let the equations of the dynamics determine what deformation, giving rise to the structuring, is physically relevant.

However, interestingly, while Chevillard et al. (2010) ultimately choose ω(t0) as a Gaussian white noise, in fact they first suggest modeling it as a fractional Gaussian field. They do not push this idea further for simplicity, briefly commenting that it is a less straightforward way. Here in contrast, we decide to explore this second possibility because it is actually going to serve our purpose.

Indeed, on the observational side, observations of the interstellar medium show that gas density may often be well represented by a log-normal distribution (e.g. Levrier et al. 2018). At the same time, on the theoretical side, as detailed in Durrive et al. (2020), a convenient manner of modeling an intermittent scalar field with a long-range correlation structure is to take the exponential of a fractional Gaussian field, an idea first introduced by Kahane (1985) (cf. Gaussian multiplicative chaos). Finally, with our approach the density field emerges as (within a corrective function) the exponential of ϑg (cf. equation ( 58) below). For these reasons, we are naturally led to construct ϑg as a fractional Gaussian field (cf. Durrive et al. 2020), namely

ϑg ≡ R d r • sg (r 2 + η 2 d ) h d dV. ( 48
)
This form should now be familiar to the reader: it is reminiscent of the result of turbulization (44) in the example of section 4.1. Accordingly, R d denotes the region of integration, by which the injection (or integral) scale is controlled, η d corresponds to the dissipation scale and h d regulates how rough the field is. The added subscript 'd' stands for 'density' because, as previously mentioned, (48) will serve to control the statistics of the density field through (58) below.

The field sg is a Gaussian white noise vector whose components are independent of one another, zero-averaged, and with unit variance. We will call it the 'seed', because it will ultimately be the only random element in the whole model, and will thus be the very origin of the randomness.

From an observational viewpoint, the choice ( 48) is convenient as it introduces parameters that may readily be fitted to astrophysical data, through (58) below. From a theoretical viewpoint, with this construction we completely anchor this modeling into the multiplicative chaos approach, pursuing in the spirit of the incompressible case (Chevillard et al. 2010;Durrive et al. 2020). An additional advantage is that fractional Gaussian fields are still 'simple enough' to have few structures, such that with this choice we still let the short-time dynamics structure the fields act in a natural way. Indeed, the risk when introducing by hand complex features into the model (through the initial conditions in this case) is to spoil the 'Navier-Stokes driven' features of the fields we build.

It remains to choose ωg. To decide, let us think in terms of degrees of freedom (DoF). We aim at constructing a stochastic representation for a three-dimensional velocity field, therefore we have three DoF to fix. But a velocity field is fully determined by its dilatation field (scalar field, so one DoF) and vorticity field (divergence-free vector field, so two DoF). So far we have fixed one DoF, by prescribing the scalar field ϑg. But to do so we introduced a vector field, sg, that is three DoF. Therefore we need to construct a divergencefree vector field ωg, given that we have at hand a Gaussian white noise vector sg. Seen this way, the most natural choice is to carry on with fractional Gaussian fields, and use the following expression (Robert & Vargas 2008) ωg

≡ Rω r × sg (r 2 + η 2 ω ) hω dV, (49)
which is a divergence-free Gaussian field built with sg, as we were looking for. This expression also elegantly mirrors the expression of ϑg in (48), in the sense that they only differ by the fact that one is derived using a scalar product while the other uses a vector product. Thus we may also infer a physical interpretation for Rω, ηω and hω by analogy, where now the subscript 'ω' means that these are vorticity-related quantities. Finally, expressions ( 48) and ( 49) also have the virtue of being differentiable, without needing to invoke the theory of distributions, contrary to Gaussian white noises.

Having now set the initial dilatation field ϑg and the initial vorticity field ωg, we have in fact also set the initial velocity field, according to the kinematic relation ( 7), together with (8) and ( 9). Now, to be fully consistent, we also need to apply turbulization to this relation. Also, since it will thus be constructed through a linear operation over the Gaussian fields ϑg and ωg, it will also be a Gaussian random field. Let us therefore denote the initial velocity field as vg. Its explicit form is given by

vg = vc,g + vr,g, (50)
with the compressional component

vc,g = 1 4π Rv ϑgr (r 2 + η 2 v ) hv dV, (51)
and the rotational component

vr,g = 1 4π Rv ωg × r (r 2 + η 2 v ) hv dV. ( 52
)
turbulization introduced Rv, ηv and hv which have analogous physical interpretation as the other (R d , η d , h d ) parameters, but related to velocity now, as indicated by the subscript 'v'. Finally, having now set the initial velocity field, we have also set the initial velocity gradient. Indeed, applying turbulization to ( 14) yields an expression for the randomized velocity gradient, namely

Gg ≡ Dg + 1 2 ωg×, (53)
where, following ( 16),

Dg = Dc,g + Dr,g, (54)
with the compressional component

Dc,g = 1 4π Rω (r 2 I -3rr) ϑg (r 2 + η 2 ω ) hω dV. ( 55
)
and the rotational component

Dr,g = - 3 8π Rω r( ωg × r) + ( ωg × r)r (r 2 + η 2 ω ) hω dV. ( 56
)
The resulting fields Dg and Gg are Gaussian, as indicated by the subscript g, because they are defined as a linear operation over a Gaussian white noise. Note that during the turbulization in Dg we take the same Rω, ηω and hω parameters as in (49) because in the following Dg and ωg are both going to be controlling the statistics of the vorticity field (cf.

(59) below).

Deformed state

Now that our initial state is fixed, let us deform (i.e. evolve) the inital fields according to the dynamics described in section 3.7.

The dilatation field evolves according to (33). Now applying turbulization to it, we obtain its randomized version

ϑ = ϑg -τ tr( G 2 g ), ( 57
)
where Gg is evaluated using the expression (53). The density field evolves according to (34), which similarly now results in the random density field

ρ = e -τ ϑg + τ 2 2 tr( G 2 g ) ρ0. (58)
In essence, this expression simply translates mass conservation. The field ρ corresponds to a smooth density profile ρ0(x), introduced in (45), which is perturbed by fluctuations ϑg and Gg, given by ( 48) and ( 53), respectively. The form (58) corresponds to the scalar field built at the beginning of Durrive et al. (2020), cf. their equation ( 9). This latter model was purely formal, inspired from multiplicative chaos only. Here, however, we are led to consider an exponential: it naturally arises from the dynamics (34) itself, such that we can now associate to the Gaussian field entering the description the physical interpretation of a dilatation field.

To construct the field (58), we have replaced G by a Gaussian random variable, Gg. But at the same time the dynamics brought us to introduce the trace of G 2 . Therefore we would like to have an idea of what the statistical properties of the trace of the random variable G 2 g might be.

At this stage, all we can say for sure is that, given that the product of Gaussian processes is in general not a Gaussian process, we at least expect tr( G 2 g ) to be non-Gaussian. This fact is encouraging because, as detailed in section 3.4, we indeed expect the trace of G 2 to be an intermittent field.

An important consequence of the fact that tr( G 2 g ) is intermittent is that our density field (58) is not exactly lognormal. It is rather the product of a lognormal distribution, due to the Gaussian part ϑg, with a non-lognormal distribution, due to tr( G 2 g ). This feature may seem problematic at first sight since in section 4.2 we stated that we needed to end up with a lognormal density field, based on observations. Fortunately, however, here tr( G 2 g ) enters only as a corrective term, since it is a higher order term in the intermittency parameter (it is proportional to τ 2 ). And in fact, it turns out that deviations from lognormality are indeed expected in high Mach number turbulence (Robertson & Goldreich 2018), which is an underlying assumption in our construction. Therefore, this byproduct of our construction is actually an interesting feature.

Likewise, inspired from the short-time dynamics (35) together with (58), we propose to model the turbulent vorticity field as

ω = e -τ ϑg+ τ 2 2 tr( G 2 g ) e τ Dg • ωg. (59)
The exponential of the randomized strain-rate matrix (54) models the random stretching of vortices (in line with the fact that turbulence may basically be defined as randomly stretched vortices, as stated in [START_REF] Wu | Vorticity and Vortex Dynamics[END_REF]), and the exponential containing the dilatation field models the random three-dimensional contraction of vortices (along and perpendicular to their axes). Compressibility also comes into play through ϑg inside the matrix (54). Finally, still in the same spirit, we construct the velocity field from a 'turbulent version' of the Helmholtz decomposition (7), namely

v = vc + vr, (60)
with the compressional component

vc = 1 4π Rv ϑr (r 2 + η 2 v ) hv dV, (61)
and the rotational component

vr = 1 4π Rv ω × r (r 2 + η 2 v ) hv dV, (62)
where the same turbulization was performed, with subscripts v this time. This velocity field v is the deformed (i.e. evolved) version of the undeformed (i.e. initial) velocity field vg given by (50). Indeed, as a simple sanity check, notice that when τ = 0, ϑ given by (57) reduces to ϑg, and ω given by (59) reduces to ωg (the matrix exponential reduces to the identity matrix), and, therefore, v given by ( 60) reduces to vg given by ( 50). As a final remark let us stress that in fact, for pedagogical purposes, throughout the above text we misnamed two quantities. Because we tweaked a little the exact formulae from section 3, strictly speaking ϑ is not the dilatation field corresponding to v in (60) and ω is not exactly its vorticity (a misnomer that was already mentioned in Durrive et al. 2020). Indeed, to be more precise, ϑ and ω are only fields that enter the modified Helmholtz decomposition (60) at the place where the true dilatation and vorticity enter the usual (i.e. unmodified) Helmholtz decomposition (7). Hence ϑ and ω should rather be called something like pseudodilatation and pseudo-vorticity. They only serve as intermediate tools that enable us to model the turbulent velocity field v through (60). The true dilatation field and vorticity field are then given by ∇ • v and ∇ × v, respectively, which can be calculated afterwards.

Advection

The above modeling has an important shortcoming, which was already present, but not tackled, in the incompressible model of Chevillard et al. (2010): advection has been disregarded. Indeed, the relations derived in section 3.7 hold as long as stretching dominates over advection, as detailed in [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF] and in section 3.2 of Durrive et al. (2020). This issue arises because our approach is inspired from works (cf. appendix A) developed in the Lagrangian framework, i.e. following the trajectory of fluid particles, while we work using Eulerian coordinates, i.e. with a fixed frame. Focusing on a regime in which advection is subdominant over stretching is convenient because solving a full advection equation with a variable velocity is not easy. Fortunately, since our approach has the advantage of being intuitive, we may rely on physical arguments to make an educated guess to bypass this limitation at least partially, as follows.

Since during the dynamics of section 3.7 the fluid particles are advected (while being stretched) during a short period of time, let us assume that advection occurs as if the velocity does not vary during this brief time lapse. In other words, we approximate this advection as a simple translation in the direction of the initial velocity field. This is somewhat similar to the Zel'dovich (or ballistic) approximation in cosmology, where the trajectory of particles is approximated as a straight line in the direction of their initial velocity. Of course this analogy has its limitations since for dark matter the dynamics is that of a collisionless gas under gravity, but the fact that this approximation gives rise to the caustics describing the cosmic web is a hint that in our case taking advection into account may help make the shape of the structures in the fields we build more realistic, inducing shock-like structures. Hence, the most natural and convenient way to take this motion into account is to advect the dilatation and vorticity fields by the 'initial' velocity (50), i.e. we propose the expression

v(x) = 1 4π Rv ϑ(y -τ vg)r + ω(y -τ vg) × r (r 2 + η 2 v ) hv dV, ( 63
)
as an improved version of (60).

THREE MODELS FOR COMPRESSIBLE MHD TURBULENCE

Section 4 was dedicated to building a turbulent velocity field. This approach may be extended to describe magnetized turbulence. Under the assumption of incompressibility, this is what Durrive et al. (2020) did. They proposed three magnetic field models, which we are now ready to generalize taking compressibility into account.

Model 1. In many astrophysical environments, density inhomogeneity plays an important role, which is why we introduced the smooth, large-scale density profile ρ0(x) in section 4. Similarly, in general turbulent magnetic fields do not only consist of homogeneous fluctuations, but they also comprise an ordered component. In view of modeling astrophysical fluids, we therefore choose the initial magnetic field entering the time evolution equation ( 36) as a smooth, ordered, field. This means that we make the substitution

B(t0, x) → B0(x), (64)
where the ordered magnetic field B0 provides us with largescale features in the modeling: it can for instance be magnetic spiral arms when modeling discs of galaxies, the Xshaped magnetic field of a galactic halo, the magnetic field compressed in the shell of the Then, applying turbulization to the expression of the magnetic field (36), we construct a first model of turbulent magnetic field, namely

B (1) = e -τ ϑg + τ 2 2 tr( G 2 g ) e τ Gg B0, (65)
where we have used the expression (58) for the random density field. Expression (65) generalizes the incompressible model 1 of Durrive et al. (2020) in two ways: firstly the dilatation field enters explicitly the definition of Gg, and secondly there is now an exponential prefactor due to compressibility, which is identical to the one entering the random vorticity field (59), and which we may interprete similarly as the random three-dimensional contraction of flux tubes along and perpendicular to their axes.

Models 2 and 3. Durrive et al. (2020) introduced another idea, complementary to model 1, in order to model the magnetic field: let us apply the same procedure as above, but starting back from the magnetic Biot-Savart law (12). Of course, since magnetic fields are divergence-free, there is no need to include any dilatation field now. However, this approach introduces the current density field j, for which we need to find a relevant 'randomized' expression j, analogous to ω.

The current density being formally similar to the vorticity of the hydrodynamical case, we may expect a priori that a relevant expression for j will be something similar to (59). Of course, the most appropriate way to pursue would be to derive the full dynamical equation for j, i.e. the generalized Ohm's law, in order to obtain an expression for the matrix that deforms j, in analogy with all that we have done so far. However, we leave this task for future work. Instead, here, we generalize models 2 and 3 of Durrive et al. (2020) by making the same simple reasoning as in their section 4, i.e. we only propose a way to add compressible effects in their framework, without trying to improve their modeling otherwise. Hence, in the line of Durrive et al. (2020), we assume that the matrix field responsible for deforming j is the same as that deforming B, namely G. This assumption is inspired from force-free dynamics in which j is proportional to B. Now, to convince ourselves that this is not as restrictive as it may seem and that the resulting model may potentially be useful beyond this force-free framework, we refer the reader either to Durrive et al. (2020) or to our arguments in sections 3.4 and 3.5 about neglecting the magnetic feedback due to the Lorentz force in ( 21) and ( 27).

On the basis of these considerations, we propose to take inspiration from (65) and to consider as randomized current density field the following expression

j (i) ≡ e -τ ϑg+ τ 2 2 tr( G 2 g ) e τ Gg j (i) 0 , ( 66
)
where i is an index that will be useful right below to distinguish between two models, called models 2 and 3 (hence i = 2 or 3). The dilatation field ϑg is given by ( 48) and the stretching/shearing matrix Gg by (53). By analogy with (65), the two exponential factors in (66) may be interpreted as random three-dimensional contraction and as random stretching and shear, respectively, and the field j (i) 0

corresponds to the initial state of the current density field j i) . We make the following choices for this initial state. As detailed in section (4.2), two extreme cases are of particular interest: taking a random field or on the contrary an ordered field. In the first case, let us consider by analogy with ω, given by ( 49), that the initial current density is a fractional Gaussian field. This will constitute the initial current density for our model 2, namely

j (2) 0 = R j r × sg (r 2 + η 2 j ) h j dV, (67)
where (Rj , ηj , hj ) are the same phenomenological parameters as every time previously, but related to the current density, hence the subscript 'j'. This model is well suited for small-scale turbulence since it does not contain any largescale feature. In the second case, let us consider an ordered, large-scale current density field, j 0 (x). The initial current density for our model 3 is then

j (3) 0 = j 0 , (68)
and if need be, j 0 (x) may be linked to some ordered magnetic field through Ampère's law j 0 = (∇ × B0)/µ0. Finally, applying turbulization to the magnetic Biot-Savart law (12), we naturally build two models (i = 2 or 3) for random magnetic fields, namely

B (i) = µ0 4π Rm j (i) × r (r 2 + η 2 m ) hm dV, (69)
where j (i) is given by (66), and i = 2 if the initial current density j

(2) 0 given by (67) is chosen, and i = 3 if instead j

(3) 0

given by ( 68) is considered. The parameters (Rm, ηm, hm) correspond, in direct analogy with the hydrodynamical case, to the injection scale, dissipation scale and Hurst exponent of the magnetic field. As a final remark, let us highlight an important feature of the above modeling. The link between turbulent and ordered components differs significantly from many other turbulence models in astrophysics. Indeed, our turbulent fields are the result of the deformation of the ordered fields, i.e. the ordered fields are multiplied by some fluctuating function (or matrix), while often turbulence is modeled as a sum, namely B = B ordered + B fluctuating . Our approach thus introduces a natural interaction between scales, without assuming any scale separation. Another important difference is that, as a consequence, our fluctuations are proportional to the amplitude of the ordered component. The benefits and drawbacks of these distinct modelings are worth investigating, but doing so is left for future work.

CONCLUSION AND PROSPECTS

In this work, we have presented an extension to the hydrodynamical turbulence model of Chevillard et al. (2010) and the MHD turbulence models of Durrive et al. (2020) to take into account compression. This new stochastic representation is complementary to the two previous models, as it corresponds to an opposite regime of the dynamics, namely while the former were dealing with the zero Mach number regime (incompressible, highly subsonic), our construction is best suited for the infinite Mach number regime [generalizing to low, intermediate, high or arbitrary Mach numbers is left for future work, but from the construction it is already clear that this will consist in modifying the closure relation, cf. Appendix A, since the Mach number enters the pressure terms of ( 21) and ( 27)]. Even though we have two more fields to consider compared to them due to compressibility, namely the dilatation and density fields, we managed to keep a very strong asset of the method: all our fields (dilatation, density, velocity, vorticity, deformation, current density, and magnetic fields) are correlated, as they should in real turbulence. Indeed, all our fields are derived from a single Gaussian white noise vector, sg, which we call the seed because it is the element that introduces all the randomness in the model. Besides, the fields are correlated in a natural physical way, since their expressions derive from the equations of the dynamics.

Let us mention a few prospects for this work. Several free (physically motivated) parameters emerge from the construction of our model. These parameters are interesting degrees of freedom for practical purposes, as they give more latitude to fit data, but from a theoretical viewpoint the origin of these parameters needs more investigation. We already know at this stage that they are not all independent. For instance, using a one dimensional version of his model, [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF] showed analytically that the intermittency and Hurst parameters are linearly coupled once one imposes the 4/5 law of turbulence. We will develop this important point in a forthcoming paper, as we will explore the parameter space.

A major advantage of the present approach is that it is modular: we construct the fields gradually, adding more and more ingredients in the model, each step being physically meaningful. For simplicity in this paper we made some relatively drastic approximations during the derivation. But because the model is modular, the path to improve it is clearly marked. The construction is based on the physics (the only subtle mathematical tool is multiplicative chaos, but it is not necessary to delve into its details to pursue the modeling), such that many options are already available in the literature to have physically motivated inspiration for new closure relations to input or to change our choices regarding the initial conditions, the large-scale cut-offs and the regularizations. In addition, a vast literature exists (cf. references mentioned in this paper) concerning the kinematics and dynamics of the gradient of velocity. At this stage, we have not taken full advantage of the results derived there. For example we did not push further the analysis of the invariants of the velocity gradients, the so-called (P, Q, R) space, which could guide us towards better choices in the modeling to improve the shape of the structures emerging in the density, velocity and magnetic fields. Nor did we analyze how the properties of our fields vary with the Reynolds number, while work such as Das & Girimaji (2019) could for instance help to do so. Another idea would be to push one order further the expansion we made in the time variable (interpreted as the intermittency parameter) when solving for the vorticity dynamics. Extending our work to take the boundaries of the system into account is important as well, and could for example be done starting with the Biot-Savart operator of a bounded domain (e.g. Enciso et al. 2018). In that sense, improving the present model should be rather straigthforward, and each new choice will result in a new model.

Finally, let us mention that we indeed expect many further improvements in the close future, well beyond our own efforts. The mathematics community is very actively exploring the concept underlying this approach, namely (Gaussian, matrix) multiplicative chaos, notably because it has many fields of applications (e.g. Bröker & Mukherjee 2020;Hager & Neuman 2020;Lacoin 2020;Powell 2020, for most recent work). More specifically, in view of building a realistic stochastic representation of turbulence, ongoing efforts include Chevillard et al. (2019) who manage to build a (onedimensional for now) stochastic field which satisfies the axioms at the core of Kolmogorov's phenomenology, as well as Reneuve & Chevillard (2020) who open the path to developing a non-stationary description.

(ii) which is deformed by the fluid during a short period of time, a period short enough for the velocity gradient to remain constant. Clearly, point (i) is inspired from the EGF (but only the exact analytical part of this model, without requiring the tuning with numerical simulations). This, compared to the RFD, improves the choice of initial condition because thus the initial conditional pressure Hessian is no longer isotropic, which relaxes the strong assumption made in the RFD. In contrast, point (ii) is inspired from the RFD, and is the ingredient that makes the statistics of the field non-Gaussian in a most natural way, i.e. guided by the dynamics.

The above examples were taken from incompressible hydrodynamics studies. As an example of closure constructed in the context of compressible flows, let us mention the so-called Homogenized Euler model by Suman & Girimaji (2009). In this case, a polytropic equation of state is considered, and, most importantly, the velocity gradient G is assumed to be spatially homogeneous, as in Vieillefosse (1982) (hence its name referring to the Restricted Euler closure), and discard all non-local effects. The same authors later introduced an enhanced version of their model, accounting for some critical Mach-number-, Reynolds-number-and Prandtlnumber-dependent physics in Suman & Girimaji (2011). While this model shares the same ambition as our model, namely of tackling the question of compressibility, we preferred nonetheless to take inspiration from the RDGF instead, at this stage at least.

  Local Bubble, or a unidirectional field as the one over the full central molecular zone of the Milky Way (cf. Planck Collaboration et al. (2016), Ferrière & Terral (2014), Alves et al. (2018), and Mangilli et al. (2019) respectively).

© 2020 The Authors

However, let us straightaway reassure any concerned reader that no prior knowledge of this specialized mathematical topic is needed to understand our paper, as detailed in the first appendix ofDurrive et al. (2020).

More precisely, this randomization is only one step out of a five-step procedure that we call 'turbulization'. The other steps introduce an injection scale, a dissipation scale, a parameter quantifying how 'rough' the turbulence is, and a parameter quantifying the intermittency. For details, see section 4.1.

3 Doing so is natural since vortex tubes have long been identified as a key quantity in modeling hydrodynamical turbulence(She et al. 1990;Farge et al. 2001;Apolinário et al. 2020).MNRAS 000, 1-14 (2020)

MNRAS 000, 1-14 (2020)

In fact, in order to lighten our expressions, our definition of the Hurst parameter is shifted by a constant and with a sign reversed compared to the standard definition in, for example, Fractional Brownian Motion.

Compared toDurrive et al. (2020), we now prefer the notation η instead of ǫ, to keep the letter ǫ for dissipation rates.MNRAS 000, 1-14 (2020)

Actually in this particular example we consider a very simple version of the full equation (21), keeping only the quadratic selfamplification of velocity gradients, which is kinematic in nature. But this does not change the fact that the second equation in (37) is the integrated version of the dynamical equation (21).

At least it is so in this picture, useful to construct the model, but the subtlety is that in general tt 0 and τ do not have the same dimensions, because we generalized the value of the Hurst exponent in relation (39).MNRAS 000, 1-14 (2020)

ACKNOWLEDGMENTS

We warmly thank Mathieu Langer and Rodrigo Pereira for fruitful discussions. This research is supported by the Agence Nationale de la Recherche (project BxB: ANR-17-CE31-0022) and the European Research Council (Advanced Grant MIST (FP7/2017-2022, No 742719)).

APPENDIX A: THE CHALLENGE OF CLOSURE

The evolution of velocity gradients along Lagrangian paths (19) is the central equation in our model. Physically, it bears the mechanisms deforming the fluid and giving it its non-trivial structure. In fact, as seen with a first order spatial Taylor expansion, G contains all the information on the small-scale turbulence. Mathematically, this equation is crucial since this is the place where we close our system of equations, i.e. a place where we choose what physics to put into our model. Indeed, equation (19) as such is not closed (not self-consistent): notably we need to specify the pressure term. Now, in order to give the reader a feeling for how this equation may be closed, as well as to illustrate the background material on which we based our own work, we mention in this appendix works taken from the hydrodynamics literature tackling this difficult question.

In the works we detail right below, the authors consider the following analogue of our equation ( 19)

where P is the pressure Hessian tensor

and V is a term corresponding to viscosity, which we do not need to explicit because dissipative effects are taken into account through regularizations, which is the reason why we work in the inviscid limit in section 3. Note, however, that this viscous term is treated in the literature in a very similar way to what we are about to show for the pressure term.

The simplest closure consists in keeping in (A1) only the term -G 2 corresponding to the local self-amplification of the velocity gradient tensor along Lagrangian fluid trajectories, neglecting the non-local pressure contributions and dissipative effects. It was first introduced and studied by Vieillefosse (1982); Vieillefosse (1984); Cantwell (1992) and is called the Restricted Euler approximation. This approximation is theoretically valuable per se because it helps focus on the stretching mechanism, but it is also the starting point for other, more realistic, closures. The problem, however, is that it leads to the development of a singularity in finite time, i.e. a divergence of the velocity gradient. Pressure and dissipation are responsible for opposing the Restricted Euler singularity in real fluids. For instance, to eliminate this singularity, one may add a linear relaxation term as in Martín et al. (1998), though this does not work for arbitrary initial conditions.

In order to deal with the pressure Hessian term, a common approach is to split the latter as

where P ij |G denotes the conditional mean of P ij , and δP ij is a fluctuation about this mean. The former quantity is modeled through statistical approximations (for example the RFD, EGF and RDGF below), and the latter is modeled using Gaussian white noises. With this decomposition, equation (A1) becomes a Langevin equation, where the fluctuating pressure term acts as a stochastic forcing term, accounting for the effects not modeled explicitly, such as perturbations from neighboring eddies, larger-scale eddies or some large-scale forcing. As a first example of this approach, one may have a look at Girimaji & Pope (1990), where the dynamics of the velocity gradient is modeled as a tensor-valued diffusion process, in order to specifically reproduce the log-normal statistics of pseudo-dissipation. But for the direct purpose of the present work, it is the three following examples that matter the most.

At the origin of our model is the Recent Fluid Deformation (RFD) closure, introduced in [START_REF] Chevillard | Habilitation à diriger des recherches[END_REF]; Chevillard et al. (2008); Chevillard & Meneveau (2011). The central idea is to approximate the conditional pressure Hessian P ij |G as:

(i) an initially isotropic tensor, (ii) which is deformed by the fluid during a short period of time, a period short enough for the velocity gradient to remain constant. Point (i) introduces a very simple 'initial' state (i.e. before the deformation), such that the non-trivial statistical properties of P ij |G stem from the deformation of fluid particles by G. Point (ii) clarifies the use of the word 'recent' in the name RFD. The RFD is helpful for understanding the statistics of G (e.g. Moriconi et al. 2014), and to build stochastic models (Chevillard et al. 2010), because it is able to prevent the Restricted Euler singularity and with it one can reproduce many well-known velocity gradient characteristics (geometric properties and anomalous scaling properties).

In parallel, another closure was suggested by Wilczek & Meneveau (2014). These authors derived the exact analytical expression of P ij |G for an isotropic Gaussian velocity field and then, in order to prevent the Restricted Euler singularity due to the self-amplification term -G 2 , they enhanced their model by tuning their parameters with data from direct numerical simulations. Consequently, they gave this closure the name 'Enhanced Gaussian Field' (EGF) closure. While this approach may seem simplistic because of the Gaussian assumption, firstly it demonstrated that a significantly anisotropic conditional pressure Hessian could result from an isotropic Gaussian velocity field (due to the inherent velocity-pressure couplings), and secondly because the resulting pressure Hessian is non-trivial, it turns out that the EGF closure accounts for a certain number of non-trivial properties of the velocity gradient, enough to compete with the RFD.

Soon afterwards, Johnson & Meneveau (2016) weighted the pros and cons of both approaches, and built yet another closure, the so-called Recent Deformation of Gaussian Fields (RDGF) closure, which manages to combine the advantages of both the RFD and the EGF. Now the central idea is to approximate P ij |G as:

(i) initially that of an isotropic Gaussian velocity field,