Statistical model for filamentary structures of molecular clouds -- The modified multiplicative random cascade model and its multifractal nature

We propose a new statistical model that can reproduce the hierarchical nature of the ubiquitous filamentary structures of molecular clouds. This model is based on the multiplicative random cascade, which is designed to replicate the multifractal nature of intermittency in developed turbulence. We present a modified version of the multiplicative process where the spatial fluctuations as a function of scales are produced with the wavelet transforms of a fractional Brownian motion realisation. This simple approach produces naturally a log-normal distribution function and hierarchical coherent structures. Despite the highly contrasted aspect of these coherent structures against a smoother background, their Fourier power spectrum can be fitted by a single power law. As reported in previous works using the multiscale non-Gaussian segmentation (MnGSeg) technique, it is proven that the fit of a single power law reflects the inability of the Fourier power spectrum to detect the progressive non-Gaussian contributions that are at the origin of these structures across the inertial range of the power spectrum. The mutifractal nature of these coherent structures is discussed, and an extension of the MnGSeg technique is proposed to calculate the multifractal spectrum that is associated with them. Using directional wavelets, we show that filamentary structures can easily be produced without changing the general shape of the power spectrum. The cumulative effect of random multiplicative sequences succeeds in producing the general aspect of filamentary structures similar to those associated with star-forming regions. The filamentary structures are formed through the product of a large number of random-phase linear waves at different spatial wavelengths. Dynamically, this effect might be associated with the collection of compressive processes that occur in the interstellar medium.


Introduction
The origin of the scale-free nature of the interstellar medium (ISM) has been debated for several decades. One of the first scale-free signature, the Fourier power spectrum of Galactic HI emission, was attributed to the turbulent nature of the ISM (Green 1993). Turbulence naturally imposes scale-free hierarchical structures to the ISM, where large eddies are extended and subdivided into smaller ones to which they transfer a fraction of their energy up to the dissipation scale (Richardson 1922;Muzy 2019). However, comparisons with statistical models of cloud structure proved that scale-free density fluctuations do not succeed in reproducing the ubiquitous filamentary structures seen in the ISM (Elmegreen et al. 2001;Miville-Deschênes et al. 2007). More recently, some studies on thermal dust emission observed by Herschel raised the question why contrasted filamentary structures, which appear to have a characteristic width (Arzoumanian et al. 2011, do not produce any break or kink in the Fourier power spectrum (Panopoulou et al. 2017). Robitaille et al. (2014Robitaille et al. ( , 2019 showed using the multiscale non-Gaussian segmentation (MnGSeg) technique 1 that filamentary structures are in fact dominating the scale-free Fourier power spectrum. This result demonstrated that dense coherent struc-tures of the ISM also have a hierarchical geometry (most likely multifractal and intimately linked to the formation of stars) and that this geometry is intertwined with the monofractal and diffuse component of molecular clouds. This finding challenges the recent results of Roy et al. (2019), who suggested that only filamentary structures with a sufficiently high area filling factor and/or high column density contrasts have an effect on the scalefree power spectrum of dust-continuum images.
Moreover, Elia et al. (2018) recently emphasised that the monofractal approach to characterising ISM structures underlies a certain degree of degeneracy as the statistical description of the regions is based on a single parameter, the fractal dimension. They proposed to analyse the multifractal properties of ISM structures in Galactic dust-continuum maps using the box-counting approach. The multifractal mathematical framework, which is now applied in many fields, from financial markets to medical imagery analysis, can detect complex local structures and can describe local singularities. Multifractal analysis is also intimately linked to the study of fully developed turbulence (Frisch 1995), where instead of trying to capture the turbulence intermittency from the final state of the dissipation scale, the multifractal model suggests that the overall flow might be described as a superposition of coexisting structures at all scales that are correlated with each other by a cascading intensity (Muzy 2019 The multifractal analysis of ISM structures was also performed by Khalil et al. (2006) on Galactic HI-line emission using a wavelet-based formalism. In contrast with the dust continuum, HI maps revealed a monofractal signature in their study. The monofractal signature of the Galactic HI could be attributed on the one hand to the lower resolution of radio observations, and on the other hand, to the more diffuse nature of the HI emission, which is less inclined to produce local singularities than the dust far-infrared emission.
In this paper, we propose a new statistical model to produce synthetic coherent filamentary structures based on the multiplicative cascade process. This new model reproduces some aspects of the probability distribution function (PDF) of young star-forming regions and the scale-free Fourier power spectrum despite the presence of highly contrasted coherent filamentary structures. Moreover, we propose a multifractal analysis framework based on the wavelet formalism of the MnGSeg technique. As discussed in §4, we strongly suggest using the wavelet approach for the multifractal analysis of continuous images instead of the box-counting formalism developed for discrete sets of points.
The paper is structured as follows: §2 presents the fractal models of the ISM and the general physical properties assigned to them; §3 presents the details of the multiplicative random cascade model and two modified models based on the wavelet decomposition of random fractal images; §4 presents the waveletbased multifractal formalism and an analysis of the synthetic monofractal and multifractal models; §5 discusses the hierarchical nature of the models and compares them with observations; and we finally use this new multifractal framework to produce synthetic ISM structures and lay out prospects for future analyses in §7.

Fractal models
Fractional Brownian motion (fBm) models, as a basic model for spatial intensity fluctuations in the ISM, can reproduce the scalefree nature of structures measured by the Fourier power spectrum. They are produced using the inverse Fourier transform of random phase values in the u-v plane multiplied by a power law for the squared modulus of the complex numbers. For a power law of ∼ −3.0, the distribution function of the fBm model is approximately Gaussian. An fBm model with a power-law index of −3.1 is shown in Fig. 1 (a), and its distribution function is shown in Fig. 2. The fBm model has a zero mean value and a standard deviation of one.
In order to reproduce the typical log-normal density distribution measured in observations and simulations, Elmegreen et al. (2001) proposed exponentiating the fBm following the relation where f (x) is the fBm image, x is the position vector, σ ρ is the standard deviation attributed to the fBm, and ρ is the modelled gas density. According to the results of Padoan et al. (1997), σ ρ can be related to the Mach number, M, as Because σ ρ f (x) = ln ρ(x), the PDF, P, can be written as where µ is the arithmetic mean value of f (x). The result of the exponentiated fBm (efBm) is shown in Fig. 1 (b) for M = 1.5. The distribution function of this model is shown in Fig. 2. As stated by Elmegreen et al. (2001), even if this model succeeds in reproducing the log-normal PDF and the single power law of the ISM, the resulting map lacks the sharp transitions in density that are commonly identified as filaments, shells, and holes in star-forming regions.

Multiplicative random cascade model
As noted by Vazquez-Semadeni (1994), hierarchical structures are naturally expected to arise in flows in which the density has a log-normal distribution. This emergence of a universal distribution can be associated with the central limit theorem, for which the sum of a large number of random variables converges to a normal distribution. If this assumption is attributed to a random sequence of ln ρ, then the density ρ becomes log-normally distributed. Thus, according to this description, after a finite time, density can be considered the product of a large number of independent random fluctuations δ, In this paper, we point out the correspondence between Eq. 4 and the multiplicative random cascade model that is designed to reproduce the multifractal nature of intermittency in developed turbulence (Meneveau & Sreenivasan 1991;Frisch 1995). In this model, the rate of turbulent energy dissipation ε is considered first as a nonrandom positive quantity for a cubic volume of side 0 . This volume is then subdivided into eight equal cubes of side 1 = 0 /2. For every subdivisions, the dissipation ε is multiplied by a series of independent random variables W n that are distributed identically. After n generations, the dissipation value of a cube of side n = 0 2 −n becomes ε = W n W n−1 ... W 1 W 0 ε. The sequence of products by an independent random variable, δ in the case of density fluctuations, has the effect of producing large deviations, that is, an increasingly non-Gaussian statistics towards smaller scales. In the ISM, this effect is measured notably through anomalous scaling of the high-order structure functions of centroid velocity increments (Lis et al. 1996;Hily-Blant et al. 2008;Federrath et al. 2010;Bertram et al. 2015) (see §4). It has been observed in the early 1990s that a probabilistic reformulation of turbulence intermittency analysis in terms of wavelet transforms of the velocity field was a powerful alternative to velocity increments (Muzy et al. 1993;Arnéodo et al. 1995;Frisch 1995). More recently, Robitaille et al. (2014Robitaille et al. ( , 2019) discovered a similar anomalous scaling signature, generally attributed to turbulence intermittency, in Herschel column density maps using PDFs of wavelet coefficients at different spatial scales.

Modified random cascade model
In this paper, we modify the random cascade model by substituting the layered cubes at different scales with the wavelet transforms of the fBm model presented in section 2. Figure 3 gives a schematic representation of the independent random variable δ as multiscale density fluctuations. The wavelet transform of a signal gives the local amount of fluctuations as a function of scales, and, for anisotropic wavelet functions, also as a function of orientations for a 2D signal. Because continuous wavelets are  redundant, a signal can be reconstructed from its wavelet coefficients simply by integrating over the scales and the orientations (Robitaille et al. 2019), where C δ is a correction factor, µ 0 is the mean value of the original signal, andf (x, , θ j ) are the wavelet coefficients calculated from the wavelet transform, defined as In equation 6, F −1 denotes the inverse Fourier transform, and f (k) andψ ,θ (k) represent the Fourier transform of f (x) and of a daughter version of the mother wavelet, respectively. The latter being dilated to a given scale and rotated by an azimuthal angle θ. In order to recover the multiplicative process described  by Vazquez-Semadeni (1994) and modelled by the random cascade, Eq. 1 can be expressed using the wavelet synthesis relation of Eq. 5, where C is a constant depending on scale . Equation 7 demonstrates that when the fractal image is considered as the summation of many spatial fluctuations, the exponentiated fBm model of equation 1 becomes equivalent to the multiplicative process described in equation 4. However, compared to the exponentiated fBm, the layered spatial fluctuations resulting from the fBm wavelet decomposition have to be modified. An fBm image, by construction, follows a squared-amplitude power law for its spatial fluctuations. As shown in Fig. 1 (a), this power law is at the origin of the hierarchal nature of the fractal image, where low spatial frequencies dominate high spatial frequencies k. In order to meet the random cascade model requirements, all spatial frequencies need to be identically distributed. For this reason, the constant C in equation 7 becomes a scale-dependent normalisation factor C = 1/σ , where σ is the standard deviation of wavelet coefficientsf (x, , θ). The construction of a discrete multiplicative cascade, called the 'W-cascades', using an orthogonal wavelet basis was proposed by Arnéodo et al. (1998) to reproduce the multifractal nature of turbulence and of financial time series. More recently, Muzy (2019) propose another grid-free model, closer to the one proposed in the present paper, using continuous wavelet transforms. Our model uses 2D continuous wavelet transforms in order to synthesise the multifractal nature of ISM maps.
To calculate the modified random cascade model, we used the Fan wavelet transform (Kirby 2005;Robitaille et al. 2019). The Fan wavelet is a rotated version of the classical Morlet wavelet and is designed to optimally sample all azimuthal directions in Fourier space. The Fourier transform of the Morlet wavelet is defined aŝ where the vector wavenumber is defined as k = (u, v) and the constant |k 0 | = u 2 0 + v 2 0 is set to π √ 2/ ln 2 ≈ 5.336 to ensure that the admissibility condition, +∞ −∞ ψ(x)dx = 0, is almost met (Kirby 2005). The wavelet transform is performed following equation 6. The result of the modified random cascade model as defined in equation 7 is shown in Fig. 1 (c). This model uses the same fBm model as in Fig. 1 (a). A different standard deviation σ is introduced in Eq. 7 in order to keep the same value of M 1.5 while taking into account the new coefficient C . The distribution function of the modified random cascade model, for the isotropic case, is also shown in Fig. 2. The next section will explore the anisotropic property of the Morlet wavelet in the construction of the modified random cascade model.

Directional modified random cascade model
Finally, we explored the angular dependence of the modified random cascade model. As for the initial fBm and the efBm, the random cascade also fails to reproduce the typical highly contrasted filamentary structures observed in star-forming regions. Dynamically, the cumulative effect of the random sequence of ln ρ could be attributed to the accumulation of shock-like structures in the turbulent ISM (velocity shears, supersonic shocks, etc.). These shock fronts are rarely isotropic, and their preferential direction can depend on the large-scale inhomogeneity of the gas distribution before the passage of the shock front. Moreover, according to the Inutsuka et al. (2015) model, molecular clouds would be preferentially formed in limited regions where the compressional direction of the shock is almost parallel to the local mean magnetic field lines or in regions experiencing a larger number of sequential compressions.
In order to reproduce these anisotropies in the modified random cascade model, we removed the angular dependence in relation 7 and performed the multiplicative process for each angle independently before summing the resulting random cascades over all direction. In our model, the angular directions are discretised in 11 directions by the Morlet wavelet transforms in order to sample the Fourier space optimally 2 . Following this new model, equation 7 now becomes The result of the directional modified random cascade model is shown in Fig. 1 (d). This model succeeds in producing filamentary structures. It uses the same original fBm model as in Fig. 1 (a) and M 1.5. As shown in Fig. 2, the intensity distributions of the isotropic and directional random cascade models and of the efBm model are similar. It is worth noting that increasing the number of angular directions does not change the general aspect of the model, nor the straightness of the filamentary structures. Moreover, because the intensity structures in each model are different, this means that the global intensity distribution alone cannot be used to distinguish the different models. The Fourier power spectra for the four models are shown in Fig. 4. Despite the highly contrasted features and the very different spatial intensity distribution, all power spectra exhibit a single power law. As shown in Fig. 8, the multiplicative operation in our models produces a progressive non-Gaussian contribution  Table 1. towards the small scales at which the Fourier analysis is not sensitive. Compared to the fBm power law, which was chosen to be similar to the power laws measured on Gaussian components of the multiscale segmentation performed by Robitaille et al. (2014Robitaille et al. ( , 2019, the random cascade model power laws have a shallower power law than those measured for the non-Gaussian components. Because of the normalisation factor C , the power law of the random cascade models become independent of the original fBm.

Multifractal analysis
The previous section described that the sequence of products by an independent random variable produces an increasingly non-Gaussian statistics towards smaller scales. In turbulence analysis, this effect is usually measured through structure functions of the centroid of the line-of-sight projected velocity increments, defined as where δC is the centroid velocity increment, δC (r) = C(r) − C(r + ), p is the statistical moment of the distribution, and the average is calculated over all possible velocity increments separated by . When p increases, the velocity increment structure functions give a greater weight to rare events (Bertram et al. 2015). The Kolmogorov (1941) turbulence theory predicts that the constant energy cascade in a turbulent medium yields an exponent ζ p = p/3 for every p. The constant energy cascade of Kolmogorov (1941) for a subsonic non-compressive turbulent medium results in a superposition of random Gaussian fluctuations from large to small spatial scales, which is well represented by stochastic monofractals. This can be completely described statistically with a single power law. The multifractal nature of a medium induced by turbulence intermittency means that the scaling behaviour changes locally in the medium (Frisch 1995). In the case of the log-normal theory of K62 (Kolmogorov 1962), for example, ζ p becomes quadratic in p such that ζ p = 1 2 µ l p/3(p/3 − 1), where µ l quantifies intermittency. For a general description that is not only associated with turbulence, a multifractal scaling can be seen as a collection of interwoven fractal subsets with a range of scaling exponents h, and fractal (Hausdorff) dimension D(h) (Chappell & Scalo 2001).
To estimate the multifractal spectrum D(h), we chose the probabilistic reformulation of turbulence intermittency analysis in terms of wavelet transforms instead of increment structure functions (Muzy et al. 1993;Arnéodo et al. 1995). The partition functions are here defined in terms of the wavelet coefficients: Equation 11 can be considered as the equivalent of the structure function defined in Eq. 10, used here in models of column density maps, where q is analogous to p and gives a greater weight to rare events. This approach has been used to study the multifractal spectrum of HI maps by Khalil et al. (2006) through the formalism of the wavelet transform modulus maxima (WTMM) method (Arnéodo et al. 2000). The main difference with our analysis is that all wavelet coefficients are used to calculate the partitions functions, instead of only a subset of maximum values among the modulus of wavelet coefficients for the WTMM method. We also used the same wavelet as was used in the current work and in the MnGSeg technique (Robitaille et al. 2019), the complex Morlet wavelet defined in Eq. 8, instead of the derivative of a Gaussian, the DoG wavelet. From these partition functions, we can define the scaling exponent τ(q) according to the power-law fit, In equation 12, τ(1) corresponds to the Hurst exponent H that is related to the fractal dimension d f = E − H, where E is the Euclidean dimension of the image (Stutzki et al. 1998;Khalil et al. 2006). The index τ(2) is related to the second-order power law γ measured by the Fourier power spectrum. The relation between the exponent H and the power spectral index γ is given by Properly speaking, the Hurst exponent H is generally used to characterise a monofractal image or surface, as is the case for fBms. In the case of a multifractal image, where the scaling exponent changes from point to point, this exponent becomes a local quantity and is generally referred to as the Hölder exponent h(x). For a monofractal image, the function τ(q) is a linear function of q of slope H. For a multifractal image, the function τ(q) is nonlinear and presents a collection of Hölder exponents, where h = dτ(q)/dq. As noted by Khalil et al. (2006), two realisations of a given stochastic process are not a priori expected to have the same set of Hölder exponents. For this reason, multifractal analyses should be performed on several realisations of the same process. For this analysis, each model is realised on 1024 × 1024 images 64 times and averaged following the quenched averaging method (Hentschel 1994), where The Z q ( ) and τ(q) functions for the four models presented in this paper and the combined model presented in §5 are plotted in Fig. 5. As expected for the first fBm model, τ(q) is close to a linear function of q of slope H 0.55, which corresponds to γ 3.1 (see Eq. 13). On the other hand, the efBm model and the two random cascade models have a nonlinear τ(q) function, as expected for a multifractal geometry. The multifractal spectrum D(h), also called the singularity spectrum, is obtained from the Legendre transform of the partition function scaling exponent τ(q): The function D(h) measures the distribution of all the locations with a given scaling exponent h(x), that is, the Hausdorff dimension of that set of points. In other words, D(h) shows the manner in which regions of different scaling exponents fill space (Chappell & Scalo 2001). The resulting multifractal spectra for the models are shown in Fig. 6. The spectra were fitted with a parabola according to the so-called log-normal model (Kestener et al. 2010;Frisch 1995), where E = 2, and µ h and σ h are free parameters. The fitted values are listed in Table 2. The central value µ h for the fBm corresponds well to the Hurst exponent H = 0.55 (γ = 2 + 2H = 3.1; the power-law index of the fBm) and has the smallest width σ h 0.01 compared to the three other models, which have a σ h an order of magnitude higher. As expected, for a monofractal model, the distribution almost collapses to one point, but the multiplicative models of multifractal nature display a wide range of scaling exponents with a range of fractal dimensions. An independent multifractal analysis of the models has also been performed using the WTMM method (Arnéodo et al. 2000;Khalil et al. 2006). This second analysis gave the same results within the numerical uncertainties.
According to the multifractal analysis, the fBm model represents well a non-intermittent turbulent medium with a constant mean energy dissipation rate on which all statistically averaged quantities depends. Consequently, because the fluctuating quantities modelled by an fBm model are statistically well described by a single scaling exponent, H or γ, its scaling index can also be estimated with its Fourier power spectrum. A turbulent field with an intermittent energy dissipation rate over the field affects other averaged quantities, such as the velocity or the density fluctuations, and they cannot be measured with a single Fourier power-spectrum index. When the multifractal analysis is applied on the efBm, the isotropic and directional random cascade models demonstrate that these models can reproduce the intermittent behaviour of fully developed turbulence. However, the directional random cascade model is the only model of those presented here that is able to visually reproduce the filamentary aspect of the ISM.
The multifractal analysis we performed on our synthetic models has many differences compared to the recent study presented by Elia et al. (2018), and we caution against comparing the two approaches. As pointed out by Khalil et al. (2006), the box-counting approach used by Elia et al. (2018) and Chappell & Scalo (2001) is usually applied for the analysis of multifractal singular measures. The authors extrapolated this approach to continuous functions (i.e. 2D images) by considering P i as the sum of the pixel brightness for all pixels contained in the ith box. One consequence of this approach is the trivial estimate of the scaling exponent h = 2.0 for all fields instead of its relation with the Fourier spectrum power-law index, as demonstrated in our study. Alternatively, Arnéodo et al. (1995) suggested the use of wavelet functions for the multifractal analysis of continuous functions as a natural oscillatory variant and generalisation of the box function.
Furthermore, Khalil et al. (2006) showed that a multifractal analysis must be performed on several realisations of the same process with a well-defined averaging protocol. Our analysis was therefore realised on 64 images of similar synthetic models. Importantly, Khalil et al. (2006) also pointed out that the calculation of the Z q ( ) function is increasingly less accurately estimated at large scales, especially for high values of |q|. For this reason, great care should be exercised when considering the range of q values used in a multifractal analysis. The determination of this range of statistical order moments should not be arbitrary, but rather dictated by the availability of statistics in the analysed dataset to ensure statistical convergence. The range of q values can be robustly and objectively determined (see Khalil et al. (2006) Eq. 32 and Figures 8 and 25). The multifractal characterisation of a medium or a surface is a very sensitive analysis, and we advise care in the choice of the approach according to the type of data and the size of the sample exhibiting a similar statistic. Unfortunately, several multifractal studies use q values are high as |q| = 20 without justification (Elia et al. 2018;Chappell & Scalo 2001), which can lead to numerical instabilities that may very well create artificial multifractal signatures.

Hierarchical nature of filaments
The Fourier power spectra for the four models presented in sections 2 and 3 are shown in Fig. 4. All models follow a single scaling power law. The fBm and efBm models have a similar power law, which is the power law that was initially imposed on the complex noise in the Fourier space for the fBm. The two random cascade models have shallower power laws of ∼ 1.8, indicating an enhancement of small-scale structures due to the multiplicative process, but also due to the normalisation factor C .
We recall that the models presented here are purely statistical and do not attempt to reproduce all the complex non-linear relations of physical processes that shape the ISM. However, without perfectly reproducing observational features, such as the curved shape of some filamentary structures, these models succeed in reproducing fundamental statistical properties of star-forming regions, such as the connection between the diffuse background density fluctuations, filamentary structures, and the multifractal nature of these regions. Robitaille et al. (2014Robitaille et al. ( , 2019   tures of the ISM also possess a hierarchical nature. They also demonstrated for two different regions of the ISM that the power spectrum of filamentary structures dominates intermediate and small scales without producing any break in the total Fourier power spectrum of the region. Because these contrasted features are the result of progressive non-Gaussian contributions towards the small scales, no kink or break is visible in the Fourier power spectrum. These progressive non-Gaussian contributions may or may not change the total power law measured by the Fourier spectrum because second-order moments alone cannot fully describe the statistical properties of non-Gaussian fields. A similar progressive non-Gaussian contributions is responsible for the anomalous scaling of structure functions of centroid velocity increments for orders p 4. This statistical description of the filamentary dense coherent structures observed notably in the dust continuum images differs from the model suggested by Roy et al. (2019), where the synthetic filaments were modelled with a Gaussian or Plummer profile. This model adds a contribution in terms of power only to the spatial scale corresponding to the width of the filament profile. In this case, the power law measured by the Fourier power spectrum is only produced by the fluctuating background. The filament model proposed here contributes in terms of power to a broad range of spatial scales and has a power law that is only limited in the case of real observations by the telescope transfer function. Moreover, the multiplicative operation of the model imposes a spatial coherence through spatial scales to the filamentary structures. The directional modified random cascade model is the first statistical model that can produce filamentary structures that follow a global log-normal density distribution and has a hierarchical nature that is measured with its Fourier power spectrum.
The coherent structures are formed through the product of a larger number of random phase linear waves at different spatial wavelengths. Dynamically, this effect might be associated with the collection of compressive processes that occur in the ISM. The final amplitude of the resulting filamentary structure depends on the phase coherence of the linear waves that model the shock-like wave fronts. The initial fBm image we used to construct the directional multiplicative random cascade is built with a random phase, which intends to be representative of the random mixing of the non-intermittent turbulence. As shown in Figures 1 (d), the hierarchy of filamentary structures is formed through the random phase coherence of linear waves at different spatial wavelengths. We showed with the modified random cascade model that the product of initially random fluctuations is sufficient to reproduce some of the complex statistical properties of molecular clouds.
The fBm and the efBm models have a similar power law in Fig. 4. This is due to the sensitivity limit of the Fourier power spectrum to non-Gaussian distributions as a function of scales. Robitaille et al. (2019) proposed performing the segmentation of non-Gaussianities as a function of scales using complex wavelet transforms in order to detect the contribution of these structures and extract them. Previously, Robitaille et al. (2014) showed that the efBm model failed to reproduce the shallower power law of non-Gaussianities as it is measured in observations. As shown in Fig. 4, through the multiplicative process and the normalisation as a function of scales, the modified random cascade models succeed in producing a shallower power law than the Gaussian fBm by introducing increasingly more non-Gaussianities towards the small scales. Moreover, the distribution function and the Fourier power spectrum, when averaged over angles, are blind to the directional aspect of the modified random cascade model, which makes these measurements highly degenerate.

Comparison with observations
This study demonstrates that common statistical tools, such as the Fourier power spectrum and PDFs, present a high level of degeneracy. These techniques fail to capture the level of complexity of many hierarchical structures. These hierarchical properties observed in many star-forming regions can be represented by a collection of interwoven scaling exponents with different fractal dimensions. Unfortunately, because the multifractal analysis is high sensitive, the direct application of such methods on a limited number of observations that are affected by very similar physical processes is hazardous. However, knowing that these regions possess multifractal properties opens exciting new opportunities with regard to the development of new approaches to analysing their hierarchical structures and to understand how it affects the star formation activity of the regions. The development of such approaches is beyond the scope of this paper, but this section attempts to compare some fundamental properties observed in ISM structures with our models.
The top panels of Fig. 7 show the results of the addition of two models, the initial Gaussian fBm and the directional modified random cascade produced with the same fBm. The fBm power law was kept at 3.1, but the directional modified random cascade power law was artificially changed to 2.5 to match observations (Robitaille et al. 2019). This combined model is defined by the following equation: where Γ(x) modelled the telescope transfer function with a Gaussian kernel with a full width at half maximum (FWHM) of 3.0 pixels. The crossed circle is the convolution operation. A is a constant, fixed to 5.0 in our model, to adjust the fBm power contribution compared to the multiplicative cascade. The standard   deviation and the mean value of the combined model are also adjusted in order to compare the model with the Polaris flare region (lower panels of Fig. 7). The noise level is modelled by N(x) and also based on the noise level of the Polaris region. This combined model is an attempt to statistically reproduce the dual nature of molecular clouds described by Falgarone et al. (2004). In this description, star-forming clouds are seen as the combination of a diffuse component that dominates the large scales and has a fractal geometry, and a coherent component, which is well described by a network of filaments and dense cores. This description suggests a connection between the highly dynamical diffuse medium and dense molecular clouds. From a turbulence point of view, the large-scale fractal-dominated medium is connected to the non-intermittent and monofractal turbulence of Kolmogorov (1941), and the hierarchical coherent component is attributed to a multifractal subset of the field due to the intermittent behaviour of developed turbulence (Frisch 1995, and references therein). We demonstrate using the multifractal extension of the MnGSeg technique (Robitaille et al. 2019) that the multiplicative models, including the efBm, are multifractal and thus have statistical properties similar to intermittent velocity fields of fully developed turbulence.
The top middle panel of Fig. 7 shows the distribution function of the combined model. The distribution function follows a log-normal distribution as for the efBm model. The segmented PDFs, according to the separation made by MnGSeg, correspond to the coloured distributions shown in Fig. 7. The MnGSeg technique was applied on the wavelet coefficients in order to segment the non-Gaussianities as a function of the scale and the orientation of the Morlet wavelet transforms (see Appendix A for more details). The Gaussian and the non-Gaussian (also called coherent) components of the map are reconstructed using relation 5, without adding the mean value µ 0 to either of the two wavelet coefficient subsets,f Gaussian (x, , θ j ) orf coherent (x, , θ j ). The mean value µ 0 was also subtracted from the total image PDF. Because the segmentation algorithm is applied in the wavelet space at multiple scales and for different orientations, the position of both components for the reconstructed maps is not exclusive. This property of the segmented map explains why their PDFs are overlapping. It is clear from the PDF analysis of the combined model (top middle panel of Fig. 7) and for the Polaris flare molecular complex (bottom middle panel of Fig. 7) that the exponential part of the PDFs is associated with the coherent structures in the maps. However, the coherent structures cover a wider range of intensity in the case of the Polaris region. This difference can be a sign that the multiplicative cascade alone cannot explain the highest density structures even for a low stellar activity region such as Polaris. This high density compared to the multiplicative cascade model might be explained by the effect of a global collapse and/or accretion mechanisms, which in addition to multiplicative processes have the ability to concentrate more the mass of the coherent structures (Vázquez-Semadeni et al. 2019).
The right top panel of Fig. 7 presents the power spectrum analysis of the combined model. The Fourier analysis gives the total power spectrum of the model. The addition of the two models with two different power laws produces no clear break in the Fourier power spectrum, as for the Polaris flare region. The symbols are the averaged squared absolute value of complex wavelet coefficients at a certain spatial scale (see Appendix A). The segmentation shows that the Fourier power spectrum is dominated by the filaments structures over a wide range of scales. The fitted power laws for the combined model and Polaris flare are shown in Table 3. Robitaille et al. (2019) showed that the small-scale flattening in the Gaussian wavelet power spectrum was caused by the cosmic infrared background noise, which here was not modelled by the combined model.
As shown in the turbulence review by Frisch (1995), the models presented in Figures 1 and 7 satisfy the statistical spatial distribution of turbulent velocity fields and energy dissipation associated with non-intermittent monofractal media and/or intermittent multifractal media well. Robitaille et al. (2014Robitaille et al. ( , 2019 showed that these statistical properties, like the single power law of the Fourier power spectrum that is dominated by the contrasted coherent structures in field, are also characteristic of the gas column-density distribution observed through the infrared thermal dust emission. The direct dependences of the velocity and density fields in turbulence are difficult to determine. It is likely, however, that the multifractal nature of one may affect the statistical properties of the other. There is also a possibility that because of the non-linear nature and the hierarchical structures involved in this process, the mechanism of hierarchical collapse of molecular clouds such as described by Vázquez-Semadeni et al. (2019) has the capacity of enhancing the multifractal nature of star-forming regions.

Conclusion
We presented the first statistical model that can reproduce the fundamental gas distribution properties that have been seen in observations and numerical simulation of the ISM. It reproduces the ubiquitous log-normal distribution, the scale-free power spectrum, some aspects of the filamentary structures, and the multifractal geometry that is observed and measured in modern studies. The modified multiplicative random cascade model is easy to produce. The multifractal nature of star-forming regions has great implications on the hierarchal properties of their density distribution. The statistical modelling proposed in this paper indicates that dense coherent structures that are isolated through multiscale non-Gaussian segmentation exhibit a continuum of scaling exponents of different fractal dimensions. Multiplicative cascades can explain the origin of these structures. However, comparisons with observations seem to indicate that additional compressive effects, such as the gravitational collapse of the region, are necessary to explain the high density of observed star-forming regions compared to their diffuse monofractal background. The multifractal nature of these particular regions as a function of their derived properties or degree of evolution will be investigated in future works.
The q parameter varies as a function of θ and . When the algorithm converges to an optimal value for the threshold Φ, the skewness, the third moment of the distribution, of the Gaussian wavelet coefficient distribution is calculated. If the skewness is higher than a certain value, the parameter q is diminished by 0.1. This operation is repeated until convergence of the parameter q. The segmented wavelet power spectra can be calculated by averaging the square absolute value of complex wavelet coefficients as a function of spatial scales, where δθ = 2 √ −2 ln 0.75/|k 0 |, N θ = ∆θ/δθ is the number of directions θ needed to sample the Fourier space over the range ∆θ, and N x is the number of pixels in the image. The conversion between the spatial scale and the Fourier wavenumber k is made through the relation k = |k 0 |/l. See Robitaille et al. (2019) for more details.