Issue 
A&A
Volume 641, September 2020



Article Number  A138  
Number of page(s)  11  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201937085  
Published online  22 September 2020 
Statistical model for filamentary structures of molecular clouds
The modified multiplicative random cascade model and its multifractal nature
^{1}
Univ. Grenoble Alpes, CNRS, IPAG,
38000
Grenoble, France
email: jeanfrancois.robitaille@univgrenoblealpes.fr
^{2}
Laboratoire de Physique de l’Ecole normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Diderot,
Sorbonne Paris Cité,
Paris, France
^{3}
CompuMAINE Laboratory, Department of Mathematics and Statistics, University of Maine,
Orono,
ME
04469, USA
Received:
8
November
2019
Accepted:
30
June
2020
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 lognormal 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 nonGaussian 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 nonGaussian 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 starforming regions. The filamentary structures are formed through the product of a large number of randomphase linear waves at different spatial wavelengths. Dynamically, this effect might be associated with the collection of compressive processes that occur in the interstellar medium.
Key words: ISM: general / ISM: structure / turbulence / waves / methods: statistical / techniques: image processing
© J.F. Robitaille et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The origin of the scalefree nature of the interstellar medium (ISM) has been debated for several decades. One of the first scalefree signature, the Fourier power spectrum of Galactic HI emission, was attributed to the turbulent nature of the ISM (Green 1993). Turbulence naturally imposes scalefree 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 scalefree density fluctuations do not succeed in reproducing the ubiquitous filamentary structures seen in the ISM (Elmegreen et al. 2001; MivilleDeschê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, 2019), do not produce any break or kink in the Fourier power spectrum (Panopoulou et al. 2017). Robitaille et al. (2014, 2019) showed using the multiscale nonGaussian segmentation (MnGSeg) technique^{1} that filamentary structures are in fact dominating the scalefree Fourier power spectrum. This result demonstrated that dense coherent structures 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 dustcontinuum 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 dustcontinuum maps using the boxcounting 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 HIline emission using a waveletbased 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 farinfrared 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 starformingregions and the scalefree 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 Sect. 4, we strongly suggest using the wavelet approach for the multifractal analysis of continuous images instead of the boxcounting formalism developed for discrete sets of points.
The paper is structured as follows: Sect. 2 presents the fractal models of the ISM and the general physical properties assigned to them; Sect. 3 presents the details of the multiplicative random cascade model and two modified models based on the wavelet decomposition of random fractal images; Sect. 4 presents the waveletbased multifractal formalism and an analysis of the synthetic monofractal and multifractal models; Sect. 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 Sect. 7.
2 Fractal models
Fractional Brownian motion (fBm) models, as a basic model for spatial intensity fluctuations in the ISM, can reproduce the scalefree nature of structuresmeasured by the Fourier power spectrum. They are produced using the inverse Fourier transform of random phase values in the uv 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 powerlaw index of − 3.1 is shown in Fig. 1a, 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 lognormal density distribution measured in observations and simulations, Elmegreen et al. (2001) proposed exponentiating the fBm following the relation (1)
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, , as (2)
Because σ_{ρ}f(x) = lnρ(x), the PDF, P, can be written as (3)
where μ is the arithmetic mean value of f(x). The result of the exponentiated fBm (efBm) is shown in Fig. 1b for . 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 lognormal 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 starforming regions.
3 Multiplicative random cascade model
As noted by VazquezSemadeni (1994), hierarchical structures are naturally expected to arise in flows in which the density has a lognormal 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 lognormally distributed. Thus, according to this description, after a finite time, density can be considered the product of a large number of independent random fluctuations δ, (4)
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 nonGaussian statistics towards smaller scales. In the ISM, this effect is measured notably through anomalous scaling of the highorder structure functions of centroid velocity increments (Lis et al. 1996; HilyBlant et al. 2008; Federrath et al. 2010; Bertram et al. 2015) (see Sect. 4). It has been observed in the early 1990s that a probabilisticreformulation 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. (2014, 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.
3.1 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 Sect. 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), (5)
where C_{δ} is a correction factor, μ_{0} is the mean value of the original signal, and are the wavelet coefficients calculated from the wavelet transform, defined as (6)
In Eq. (6), denotes the inverse Fourier transform, and and 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 VazquezSemadeni (1994) and modelled by the random cascade, Eq. (1) can be expressed using the wavelet synthesis relation of Eq. (5), (7)
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 Eq. (1) becomes equivalent to the multiplicative process described in Eq. (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 squaredamplitude power law for its spatial fluctuations. As shown in Fig. 1a, 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 Eq. (7) becomes a scaledependent normalisation factor C_{ℓ} = 1∕σ_{ℓ}, where σ_{ℓ} is the standard deviation of wavelet coefficients . The construction of a discrete multiplicative cascade, called the “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 gridfree 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 as (8)
where the vector wavenumber is defined as k = (u, v) and the constant is set to to ensure that the admissibility condition, , is almost met (Kirby 2005). The wavelet transform is performed following Eq. (6). The result of the modified random cascade model as defined in Eq. (7) is shown in Fig. 1c. This model uses the same fBm model as in Fig. 1a. A different standard deviation σ is introduced in Eq. (7) in order to keep the same value of 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.
Fig. 1 Spatial intensity distribution of the fractal models described in Sect. 2. Panel a: fBm model and panel b: exponentiated fBm model, and the random cascade models described in Sect. 3. Panel c: isotropic and panel d: directional random cascade models. 
Fig. 2 Distribution functions for the four models presented in Fig. 1. Distribution functions for the efBm, the isotropic, and the directional random cascade are normalised by the mean pixel value of the image. 
Fig. 3 Schematic representation of the modified random cascade model, where each level is the normalised spatial fluctuations of an fBm as a function of scales obtained from wavelet transforms. 
3.2 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 starforming regions. Dynamically, the cumulative effect of the random sequence of ln ρ could be attributed to the accumulation of shocklike structures in the turbulent ISM (velocity shears, supersonic shocks, etc.). These shock fronts are rarely isotropic, and their preferential direction can depend on the largescale 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, Eq. (7) now becomes (9)
The result of the directional modified random cascade model is shown in Fig. 1d. This model succeeds in producing filamentary structures. It uses the same original fBm model as in Fig. 1a and . 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 nonGaussian contribution 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. (2014, 2019), the random cascade model power laws have a shallower power law than those measured for the nonGaussian components. Because of the normalisation factor C_{ℓ}, the power law of the random cascade models become independent of the original fBm.
Fig. 4 Fourier power spectrum for the four models presented in Fig. 1. The dashed lines show the fitted power laws. The fitted power indices are listed in Table 1. 
4 Multifractal analysis
The previous section described that the sequence of products by an independent random variable produces an increasingly nonGaussian statistics towards smaller scales. In turbulence analysis, this effect is usually measured through structure functions of the centroid of the lineofsight projected velocity increments, defined as (10)
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 noncompressive 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 lognormal theory of K62 (Kolmogorov 1962), for example, ζ_{p} becomes quadratic in p such that 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: (11)
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. Wealso 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 tothe powerlaw fit, (12)
In Eq. (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 secondorder power law γ measured by the Fourier power spectrum. The relation between the exponent H and the power spectral index γ is given by (13)
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 (14)
The Z_{q}(ℓ) and τ(q) functions for the four models presented in this paper and the combined model presented in Sect. 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): (15)
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 socalled lognormal model (Kestener et al. 2010; Frisch 1995), (16)
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 powerlaw 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 nonintermittent 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 powerspectrum 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 boxcounting 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 powerlaw 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 welldefined 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 Figs. 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.
Fig. 5 Partition functions Z_{q}(ℓ) and the derived τ(q) functions according to their powerlaw fits for the quenched averaging of 64 realisations of the four models presented in Fig. 1 and the combined model presented in Sect. 5. 
Fig. 6 Multifractal spectrum D(h) derived fromthe τ(q) functions shown in Fig. 5. The dotted red line corresponds to h = 0.55, the Hurst exponent imposed on the fBm model. 
5 Hierarchical nature of filaments
The Fourier power spectra for the four models presented in Sects. 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 smallscale 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 nonlinear 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 starforming regions, such as the connection between the diffuse background density fluctuations, filamentary structures, and the multifractal nature of these regions.
Robitaille et al. (2014, 2019) showed using the MnGSeg technique that the dense coherent and often filamentary structures 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 nonGaussian contributions towards the small scales, no kink or break is visible in the Fourier power spectrum. These progressive nonGaussian contributions may or may not change the total power law measured by the Fourier spectrum because secondorder moments alone cannot fully describe the statistical properties of nonGaussian fields. A similar progressive nonGaussian 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 lognormal 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 shocklike 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 nonintermittent turbulence. As shown in Fig. 1d, 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 nonGaussian distributions as a function of scales. Robitaille et al. (2019) proposed performing the segmentation of nonGaussianities 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 nonGaussianities 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 nonGaussianities 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.
6 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 starforming 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: (17)
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, starforming 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 largescale fractaldominated medium is connected to the nonintermittent 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 lognormal 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 nonGaussianities as a function of the scale and the orientation of the Morlet wavelet transforms (see Appendix A for more details). The Gaussian and the nonGaussian (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, or . 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ázquezSemadeni 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 smallscale 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 Figs. 1 and 7 satisfy the statistical spatial distribution of turbulent velocity fields and energy dissipation associated with nonintermittent monofractal media and/or intermittent multifractal media well. Robitaille et al. (2014, 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 columndensity 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 nonlinear nature and the hierarchical structures involved in this process, the mechanism of hierarchical collapse of molecular clouds such as described by VázquezSemadeni et al. (2019) has the capacity of enhancing the multifractal nature of starformingregions.
Fig. 7 Top left panel: combined model, which is the addition of the fBm model in Fig. 1a and the directional modified random cascade model in Fig. 1d following Eq. (17). Top middle panel: its segmented PDFs. Top right panel: segmented powerspectra analysis. Bottom panels: same figures for the Polaris flare region that was previously analysed by Robitaille et al. (2019). 
Fig. 8 PDFs as a function of scales for the combined model defined by Eq. (17) that is shown in the top left panel of Fig. 7. 
7 Conclusion
We presented the first statistical model that can reproduce the fundamental gas distribution properties that have been seen in observationsand numerical simulation of the ISM. It reproduces the ubiquitous lognormal distribution, the scalefree power spectrum, some aspects of the filamentary structures, and the multifractal geometry that is observed and measured in modern studies. Themodified multiplicative random cascade model is easy to produce. The multifractal nature of starforming 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 nonGaussian 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 starforming 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.
Acknowledgements
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie Grant Agreement No. 750920 and from the programme StarFormMapper under grant agreement No. 687528. This work was also supported by the Programme National de Physique Stellaire and Physique et Chimie du Milieu Interstellaire (PNPS and PCMI) of CNRS/INSU (with INC/INP/IN2P3) cofunded by CEA and CNES. This research has made use of data from the Herschel Gould Belt survey (HGBS) project (http://gouldbeltherschel.cea.fr). The HGBS is a Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAFIFSI Rome and INAFArcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC). This research also made use of Astropy, a communitydeveloped core Python package for Astronomy (Astropy Collaboration 2013) and its affiliated package APLpy (Robitaille & Bressert 2012).
Appendix A NonGaussian segmentation
The nonGaussian segmentation of the wavelet coefficients is based on an iterative algorithm of denoising. The threshold Φ splits the nonGaussian and the Gaussian wavelet coefficients at all scales and directions of the Morlet wavelet transform defined in Eq. (6). The sequence defining Φ is (A.1)
where q (that is different from the q in Eq. (11)) is a dimensionless constant that controls how restrictive the definition of nonGaussianities is. The variance is defined as (A.2)
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, (A.5)
where δθ = 2 , 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.
References
 Arnéodo, A., Bacry, E., & Muzy, J. F. 1995, Phys. A Stat. Mech. Appl., 213, 232 [NASA ADS] [CrossRef] [Google Scholar]
 Arnéodo, A., Bacry, E., & Muzy, J. F. 1998, J. Math. Phys., 39, 4142 [CrossRef] [MathSciNet] [Google Scholar]
 Arnéodo, A., Decoster, N., & Roux, S. G. 2000, Eur. Phys. J. B, 15, 567 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, L6 [Google Scholar]
 Arzoumanian, D., André, P., Könyves, V., et al. 2019, A&A, 621, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bertram, E., Konstandin, L., Shetty, R., Glover, S. C. O., & Klessen, R. S. 2015, MNRAS, 446, 3777 [NASA ADS] [CrossRef] [Google Scholar]
 Chappell, D., & Scalo, J. 2001, ApJ, 551, 712 [NASA ADS] [CrossRef] [Google Scholar]
 Elia, D., Strafella, F., Dib, S., et al. 2018, MNRAS, 481, 509 [NASA ADS] [CrossRef] [Google Scholar]
 Elmegreen, B. G., Kim, S., & StaveleySmith, L. 2001, ApJ, 548, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Falgarone, E., HilyBlant, P., & Levrier, F. 2004, Ap&SS, 292, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., RomanDuval, J., Klessen, R. S., Schmidt, W., & Mac Low, M.M. 2010, A&A, 512, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frisch, U. 1995, Turbulence (Cambridge: Cambridge University Press) [Google Scholar]
 Green, D. A. 1993, MNRAS, 262, 327 [NASA ADS] [CrossRef] [Google Scholar]
 Hentschel, H. G. E. 1994, Phys. Rev. E, 50, 243 [CrossRef] [Google Scholar]
 HilyBlant, P., Falgarone, E., & Pety, J. 2008, A&A, 481, 367 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Inutsuka, S.i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kestener, P., Conlon, P. A., Khalil, A., et al. 2010, ApJ, 717, 995 [NASA ADS] [CrossRef] [Google Scholar]
 Khalil, A., Joncas, G., Nekka, F., Kestener, P., & Arnéodo, A. 2006, ApJS, 165, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Kirby, J. F. 2005, Comput. Geosci., 31, 846 [NASA ADS] [CrossRef] [Google Scholar]
 Kolmogorov, A. 1941, Akad. Nauk SSSR Dok., 30, 301 [Google Scholar]
 Kolmogorov, A. N. 1962, J. Fluid Mech., 13, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Lis, D. C., Pety, J., Phillips, T. G., & Falgarone, E. 1996, ApJ, 463, 623 [NASA ADS] [CrossRef] [Google Scholar]
 Meneveau, C., & Sreenivasan, K. R. 1991, J. Fluid Mech., 224, 429 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M.A., Lagache, G., Boulanger, F., & Puget, J.L. 2007, A&A, 469, 595 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Muzy, J.F. 2019, Phys. Rev. E, 99, 042113 [CrossRef] [Google Scholar]
 Muzy, J. F., Bacry, E., & Arneodo, A. 1993, Phys. Rev. E, 47, 875 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Jones, B. J. T., & Nordlund, Å. P. 1997, ApJ, 474, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Panopoulou, G. V., Psaradaki, I., Skalidis, R., Tassis, K., & Andrews, J. J. 2017, MNRAS, 466, 2529 [NASA ADS] [CrossRef] [Google Scholar]
 Richardson, L. F. 1922, Weather prediction by numerical process (Cambridge: Cambridge University Press) [Google Scholar]
 Robitaille, T., & Bressert, E. 2012, Astrophysics Source Code Library [record ascl:1208.017] [Google Scholar]
 Robitaille, J.F., Joncas, G., & MivilleDeschênes, M.A. 2014, MNRAS, 440, 2726 [NASA ADS] [CrossRef] [Google Scholar]
 Robitaille, J. F., Motte, F., Schneider, N., Elia, D., & Bontemps, S. 2019, A&A, 628, A33 [CrossRef] [EDP Sciences] [Google Scholar]
 Roy, A., André, P., Arzoumanian, D., et al. 2019, A&A, 626, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998, A&A, 336, 697 [NASA ADS] [Google Scholar]
 VazquezSemadeni, E. 1994, ApJ, 423, 681 [NASA ADS] [CrossRef] [Google Scholar]
 VázquezSemadeni, E., Palau, A., BallesterosParedes, J., Gómez, G. C., & ZamoraAvilés, M. 2019, MNRAS, 490, 3061 [NASA ADS] [CrossRef] [Google Scholar]
Codes and tutorials are available at https://github.com/jfrob27/pywavan
See Robitaille et al. (2019) Sect. 3.1 and Fig. 2 for a representation of the Morlet wavelet sampling in Fourier space.
All Tables
All Figures
Fig. 1 Spatial intensity distribution of the fractal models described in Sect. 2. Panel a: fBm model and panel b: exponentiated fBm model, and the random cascade models described in Sect. 3. Panel c: isotropic and panel d: directional random cascade models. 

In the text 
Fig. 2 Distribution functions for the four models presented in Fig. 1. Distribution functions for the efBm, the isotropic, and the directional random cascade are normalised by the mean pixel value of the image. 

In the text 
Fig. 3 Schematic representation of the modified random cascade model, where each level is the normalised spatial fluctuations of an fBm as a function of scales obtained from wavelet transforms. 

In the text 
Fig. 4 Fourier power spectrum for the four models presented in Fig. 1. The dashed lines show the fitted power laws. The fitted power indices are listed in Table 1. 

In the text 
Fig. 5 Partition functions Z_{q}(ℓ) and the derived τ(q) functions according to their powerlaw fits for the quenched averaging of 64 realisations of the four models presented in Fig. 1 and the combined model presented in Sect. 5. 

In the text 
Fig. 6 Multifractal spectrum D(h) derived fromthe τ(q) functions shown in Fig. 5. The dotted red line corresponds to h = 0.55, the Hurst exponent imposed on the fBm model. 

In the text 
Fig. 7 Top left panel: combined model, which is the addition of the fBm model in Fig. 1a and the directional modified random cascade model in Fig. 1d following Eq. (17). Top middle panel: its segmented PDFs. Top right panel: segmented powerspectra analysis. Bottom panels: same figures for the Polaris flare region that was previously analysed by Robitaille et al. (2019). 

In the text 
Fig. 8 PDFs as a function of scales for the combined model defined by Eq. (17) that is shown in the top left panel of Fig. 7. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.