# Statistical model for filamentary structures of molecular clouds

J.-F. Robitaille, A. Abdeldayem, I. Joncour, E. Moraux, F. Motte, Pierre

Lesaffre, A. Khalil

► To cite this version:

J.-F. Robitaille, A. Abdeldayem, I. Joncour, E. Moraux, F. Motte, et al. Statistical model for filamentary structures of molecular clouds. Astronomy & Astrophysics - A&A, 2020, 641, pp.A138. 10.1051/0004-6361/201937085. [hal-02946117](https://hal.science/hal-02946117v1)

HAL Id: hal-02946117

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

Submitted on 22 Sep 2020

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Statistical model for filamentary structures of molecular clouds
## The modified multiplicative random cascade model and its multifractal nature

J.-F. Robitaille<fnref n="1" />, A. Abdeldayem<fnref n="1" />, I. Joncour<fnref n="1" />, E. Moraux<fnref n="1" />, F. Motte<fnref n="1" />, P. Lesaffre<fnref n="2" />, and A. Khalil<fnref n="3" />

Received 8 November 2019 / Accepted 30 June 2020

## ABSTRACT

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 multifractal 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.

**Key words.** ISM: general – ISM: structure – turbulence – waves – methods: statistical – techniques: image processing

## 1. 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, 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 non-Gaussian segmentation (MnGSeg) technique<fnref n="1" /> that filamentary structures are in fact dominating the scale-free 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 scale-free 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

<footnote n="1">¹ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
e-mail: jean-francois.robitaille@univ-grenoble-alpes.fr</footnote>

<footnote n="2">² Laboratoire de Physique de l'Ecole normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, Paris, France</footnote>

<footnote n="3">³ CompuMAINE Laboratory, Department of Mathematics and Statistics, University of Maine, Orono, ME 04469, USA</footnote>

<footnote n="1">¹ Codes and tutorials are available at https://github.com/jfrob27/pywavan</footnote>
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 Sect. 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: 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 wavelet-based 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 scale-free 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. 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 log-normal density distribution measured in observations and simulations, Elmegreen et al. (2001) proposed exponentiating the fBm following the relation

$$ \rho_{\text{efbm}}(\mathbf{x}) = \exp[\sigma_{\rho} f(\mathbf{x})], \quad (1) $$

where $f(\mathbf{x})$ is the fBm image, $\mathbf{x}$ is the position vector, $\sigma_{\rho}$ is the standard deviation attributed to the fBm, and $\rho$ is the modelled gas density. According to the results of Padoan et al. (1997), $\sigma_{\rho}$ can be related to the Mach number, $M$, as

$$ \sigma_{\rho} = (\ln[1 + 0.5M^2])^{0.5}. \quad (2) $$

Because $\sigma_{\rho} f(\mathbf{x}) = \ln \rho(\mathbf{x})$, the PDF, $P$, can be written as

$$ P(\ln \rho) = \frac{1}{(2\pi\sigma_{\rho})^{1/2}} \exp \left[ -0.5 \left( \frac{\ln \rho - \mu}{\sigma_{\rho}} \right)^2 \right], \quad (3) $$

where $\mu$ is the arithmetic mean value of $f(\mathbf{x})$. The result of the exponentiated fBm (efBm) is shown in Fig. 1b 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.

## 3. 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 \rho$, then the density $\rho$ 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 $\delta$,

$$ \rho(t_n) = \delta_n \delta_{n-1} \dots \delta_1 \delta_0 \rho(t_0). \quad (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 $\epsilon$ is considered first as a nonrandom positive quantity for a cubic volume of side $\ell_0$. This volume is then subdivided into eight equal cubes of side $\ell_1 = \ell_0/2$. For every subdivisions, the dissipation $\epsilon$ 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 $\ell_n = \ell_0 2^{-n}$ becomes $\epsilon_\ell = W_n W_{n-1} \dots W_1 W_0 \epsilon$.

The sequence of products by an independent random variable, $\delta$ 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 Sect. 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. (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 $\delta$ 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
<figure><img src="image_2.png" /><figcaption>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.</figcaption></figure>

(Robitaille et al. 2019),

$$f(\mathbf{x}) = C_{\delta} \sum_{\ell} \sum_{j=0}^{N_{\ell}-1} \ell \tilde{f}(\mathbf{x}, \ell, \theta_j) + \mu_0, \quad (5)$$

where $C_{\delta}$ is a correction factor, $\mu_0$ is the mean value of the original signal, and $\tilde{f}(\mathbf{x}, \ell, \theta_j)$ are the wavelet coefficients calculated from the wavelet transform, defined as

$$\tilde{f}(\mathbf{x}, \ell, \theta) = \mathcal{F}^{-1} \{ \hat{f}(\mathbf{k}) \hat{\psi}_{\ell, \theta}^{*}(\mathbf{k}) \}. \quad (6)$$

In Eq. (6), $\mathcal{F}^{-1}$ denotes the inverse Fourier transform, and $\hat{f}(\mathbf{k})$ and $\hat{\psi}_{\ell, \theta}(\mathbf{k})$ represent the Fourier transform of $f(\mathbf{x})$ and of a daughter version of the mother wavelet, respectively. The latter being dilated to a given scale $\ell$ and rotated by an azimuthal angle $\theta$. 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),

$$\rho_{\text{mc}}(\mathbf{x}) = \exp \left[ \sigma \sum_{\ell, \theta} C_{\ell} \tilde{f}(\mathbf{x}, \ell, \theta) \right] = \prod_{\ell, \theta} \exp \left[ \sigma C_{\ell} \tilde{f}(\mathbf{x}, \ell, \theta) \right], \quad (7)$$

where $C_{\ell}$ is a constant depending on scale $\ell$. 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 squared-amplitude 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 $\mathbf{k}$. In order to meet the random cascade model requirements, all spatial frequencies need to be identically distributed. For this reason, the constant $C_{\ell}$ in Eq. (7) becomes a scale-dependent normalisation factor
<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

$C_ℓ = 1/σ_ℓ$, where $σ_ℓ$ is the standard deviation of wavelet coefficients $\tilde{f}(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 as

$$ 
\begin{aligned}
\hat{\psi}_{\ell,\theta}(\mathbf{k}) &= \ell \cdot e^{-|\ell \cdot \mathbf{k} - \mathbf{k}_0|^2 / 2} \\
&= \ell \cdot e^{-[(\ell \cdot u - |\mathbf{k}_0| \cos \theta)^2 + (\ell \cdot v - |\mathbf{k}_0| \sin \theta)^2]/2}
\end{aligned}
\quad (8) $$

where the vector wavenumber is defined as $\mathbf{k} = (u, v)$ and the constant $|\mathbf{k}_0| = \sqrt{u_0^2 + v_0^2}$ is set to $\pi \sqrt{2/\ln 2} \approx 5.336$ to ensure that the admissibility condition, $\int_{-\infty}^{+\infty} \psi(x)dx = 0$, 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 $M \approx 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.

## 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 star-forming regions. Dynamically, the cumulative effect of the random sequence of $\ln \rho$ 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<fnref n="2" />. Following this new model, Eq. (7) now becomes

$$ \rho_{\text{dmc}}(\mathbf{x}) = \frac{1}{N_{\theta}} \sum_{\theta} \left[ \prod_{l} \exp \left[ \sigma C_{l} \tilde{f}(\mathbf{x}, l, \theta) \right] \right]. \quad (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 $M \approx 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 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

<footnote n="2">² See Robitaille et al. (2019) Sect. 3.1 and Fig. 2 for a representation of the Morlet wavelet sampling in Fourier space.</footnote>
<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 1.** Power laws fitted to the power spectra of Fig. 4.</caption><thead><tr><th></th><th>Power law (γ)</th></tr></thead><tbody><tr><td>fBm</td><td>3.01 ± 0.01</td></tr><tr><td>efBm</td><td>3.12 ± 0.01</td></tr><tr><td>Isotropic random cascade</td><td>1.87 ± 0.01</td></tr><tr><td>Directional random cascade</td><td>1.70 ± 0.01</td></tr></tbody></table>

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.

## 4. 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

$$S_p(\ell) = \langle|\delta C_\ell(\mathbf{r})|^p\rangle \propto \ell^{\zeta_p}, \quad (10)$$

where $\delta C_ℓ$ is the centroid velocity increment, $\delta C_ℓ(\mathbf{r}) = C(\mathbf{r}) - C(\mathbf{r} + \ell)$, $p$ is the statistical moment of the distribution, and the average $\langle\rangle$ is calculated over all possible velocity increments separated by $\ell$. 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 $\zeta_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, $\zeta_p$ becomes quadratic in $p$ such

that $\zeta_p = \frac{1}{2}\mu_l p/3(p/3 - 1)$, where $\mu_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:

$$Z_q(\ell) = \sum_{\mathbf{x}} \sum_{j=0}^{N_\theta-1} |\tilde{f}(\ell, \mathbf{x}, \theta_j)|^q. \quad (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. 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 $\tau(q)$ according to the power-law fit,

$$Z_q(\ell) \propto \ell^{\tau(q)}. \quad (12)$$

In Eq. (12), $\tau(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 $\tau(2)$ is related to the second-order power law $\gamma$ measured by the Fourier power spectrum. The relation between the exponent $H$ and the power spectral index $\gamma$ is given by

$$\gamma = 2 + 2H. \quad (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(\mathbf{x})$. For a monofractal image, the function $\tau(q)$ is a linear function of $q$ of slope $H$. For a multifractal image, the function $\tau(q)$ is non-linear and presents a collection of Hölder exponents, where $h = d\tau(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 \times 1024$ images 64 times and averaged following the quenched averaging method (Hentschel 1994), where

$$e^{\langle\ln Z_q(\ell)\rangle} \sim \ell^{\tau(q)}, \quad \ell \to 0^+. \quad (14)$$

The $Z_q(\ell)$ and $\tau(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, $\tau(q)$ is close to
<figure><img src="image_6.png" /><figcaption>Fig. 5. Partition functions $Z_q(l)$ and the derived $\tau(q)$ functions according to their power-law fits for the quenched averaging of 64 realisations of the four models presented in Fig. 1 and the combined model presented in Sect. 5.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Multifractal spectrum *D*(*h*) derived from the τ(*q*) functions shown in Fig. 5. The dotted red line corresponds to *h* = 0.55, the Hurst exponent imposed on the fBm model.</figcaption></figure>

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 non-linear τ(*q*) function,
as expected for a multifractal geometry.

The multifractal spectrum *D*(*h*), also called the singular-
ity spectrum, is obtained from the Legendre transform of the
partition function scaling exponent τ(*q*):

$$
D(h) = \min_q [qh - \tau(q)]. \tag{15}
$$

The function *D*(*h*) measures the distribution of all the loca-
tions 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),

$$
D(h) = E - \frac{(h - \mu_h)^2}{2\sigma_h}, \quad (16)
$$

where *E* = 2, and *μ*<sub>*h*</sub> and *σ*<sub>*h*</sub> are free parameters. The fitted values
are listed in Table 2. The central value *μ*<sub>*h*</sub> for the fBm corre-
sponds well to the Hurst exponent *H* = 0.55 (*γ* = 2 + 2*H* = 3.1;
the power-law index of the fBm) and has the smallest width
*σ*<sub>*h*</sub> ≈ 0.01 compared to the three other models, which have a
*σ*<sub>*h*</sub> an order of magnitude higher. As expected, for a monofrac-
tal 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 repre-
sents well a non-intermittent turbulent medium with a constant
mean energy dissipation rate on which all statistically averaged
quantities depends. Consequently, because the fluctuating quan-
tities 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 turbu-
lent field with an intermittent energy dissipation rate over the

<table><caption>**Table 2.** Parabola fitted value to *D*(*h*) functions of Fig. 6.</caption>
  <thead>
    <tr>
      <th></th>
      <th>&mu;<sub>h</sub></th>
      <th>&sigma;<sub>h</sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>fBm</td>
      <td>0.57 &plusmn; 0.01</td>
      <td>0.01 &plusmn; 0.01</td>
    </tr>
    <tr>
      <td>efBm</td>
      <td>0.65 &plusmn; 0.01</td>
      <td>0.10 &plusmn; 0.01</td>
    </tr>
    <tr>
      <td>Isotropic random cascade</td>
      <td>0.06 &plusmn; 0.01</td>
      <td>0.13 &plusmn; 0.01</td>
    </tr>
    <tr>
      <td>Directional random cascade</td>
      <td>0.07 &plusmn; 0.01</td>
      <td>0.22 &plusmn; 0.01</td>
    </tr>
  </tbody>
</table>

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 cas-
cade 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 pre-
sented 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 multifrac-
tal singular measures. The authors extrapolated this approach to
continuous functions (i.e. 2D images) by considering Pᵢ 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 mod-
els. Importantly, Khalil et al. (2006) also pointed out that the
calculation of the Z<sub>q</sub>(ℓ) 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 determi-
nation 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 charac-
terisation 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 statis-
tic. 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.

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 small-scale structures due to the multiplicative process, but also due to the normalisation factor $C_t$.

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. (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 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 \ge 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 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 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.

## 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 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:

$$ \rho_{\text{comb}}(\mathbf{x}) = \Gamma(\mathbf{x}) \otimes [A * \rho_{\text{fBm}}(\mathbf{x}) + \rho_{\text{dmc}}(\mathbf{x})] + N(\mathbf{x}), \quad (17) $$

where $\Gamma(\mathbf{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(\mathbf{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
<figure><img src="image_8.png" /><figcaption>**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 power-spectra analysis. Bottom panels: same figures for the Polaris flare region that was previously analysed by Robitaille et al. (2019).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

point of view, the large-scale fractal-dominated medium is con-
nected 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 inter-
mittent behaviour of developed turbulence (Frisch 1995, and
references therein). We demonstrate using the multifractal exten-
sion 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 func-
tion 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 rela-
tion (5), without adding the mean value μ₀ to either of the two
wavelet coefficient subsets, f̃<sub>Gaussian</sub>(x, ℓ, θ<sub>j</sub>) or f̃<sub>coherent</sub>(x, ℓ, θ<sub>j</sub>).
The mean value μ₀ 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 anal-
ysis 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 coher-
ent structures cover a wider range of intensity in the case of the
Polaris region. This difference can be a sign that the multiplica-
tive 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 accre-
tion 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 mod-
els with two different power laws produces no clear break in the
Fourier power spectrum, as for the Polaris flare region. The sym-
bols 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
<table><caption>**Table 3.** Power laws for the power spectra of Fig. 7.</caption><thead><tr><th></th><th>Combined model</th><th>Polaris flare</th></tr></thead><tbody><tr><td>Total Fourier</td><td>2.53 ± 0.01</td><td>2.38 ± 0.02</td></tr><tr><td>Gaussian wavelet</td><td>3.13 ± 0.03</td><td>3.05 ± 0.07</td></tr><tr><td>Coherent wavelet</td><td>2.33 ± 0.05</td><td>2.41 ± 0.04</td></tr></tbody></table>

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 Figs. 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. (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 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.

## 7. 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 hierarchical 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.

Acknowledgements. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie 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) co-funded by CEA and CNES. This research has made use of data from the Herschel Gould Belt survey (HGBS) project (http://gouldbelt-herschel.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, INAF-IFSI Rome and INAF-Arcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC). This research also made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013) and its affiliated package APLpy (Robitaille & Bressert 2012).

## References

Arnéodo, A., Bacry, E., & Muzy, J. F. 1995, Phys. A Stat. Mech. Appl., 213, 232

Arnéodo, A., Bacry, E., & Muzy, J. F. 1998, J. Math. Phys., 39, 4142

Arnéodo, A., Decoster, N., & Roux, S. G. 2000, Eur. Phys. J. B, 15, 567

Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, L6

Arzoumanian, D., André, P., Könyves, V., et al. 2019, A&A, 621, A42

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Bertram, E., Konstandin, L., Shetty, R., Glover, S. C. O., & Klessen, R. S. 2015, MNRAS, 446, 3777

Chappell, D., & Scalo, J. 2001, ApJ, 551, 712

Elia, D., Strafella, F., Dib, S., et al. 2018, MNRAS, 481, 509

Elmegreen, B. G., Kim, S., & Staveley-Smith, L. 2001, ApJ, 548, 749

Falgarone, E., Hily-Blant, P., & Levrier, F. 2004, Ap&SS, 292, 89

Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M.-M. 2010, A&A, 512, A81

Frisch, U. 1995, Turbulence (Cambridge: Cambridge University Press)

Green, D. A. 1993, MNRAS, 262, 327

Hentschel, H. G. E. 1994, Phys. Rev. E, 50, 243

Hily-Blant, P., Falgarone, E., & Pety, J. 2008, A&A, 481, 367

Inutsuka, S.-i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49

Kestener, P., Conlon, P. A., Khalil, A., et al. 2010, ApJ, 717, 995

Khalil, A., Joncas, G., Nekka, F., Kestener, P., & Arnéodo, A. 2006, ApJS, 165, 512

Kirby, J. F. 2005, Comput. Geosci., 31, 846

Kolmogorov, A. 1941, Akad. Nauk SSSR Dok., 30, 301

Kolmogorov, A. N. 1962, J. Fluid Mech., 13, 82

Lis, D. C., Pety, J., Phillips, T. G., & Falgarone, E. 1996, ApJ, 463, 623

Meneveau, C., & Sreenivasan, K. R. 1991, J. Fluid Mech., 224, 429

Miville-Deschênes, M.-A., Lagache, G., Boulanger, F., & Puget, J.-L. 2007, A&A, 469, 595

Muzy, J.-F. 2019, Phys. Rev. E, 99, 042113

Muzy, J. F., Bacry, E., & Arneodo, A. 1993, Phys. Rev. E, 47, 875

Padoan, P., Jones, B. J. T., & Nordlund, Å. P. 1997, ApJ, 474, 730

Panopoulou, G. V., Psaradaki, I., Skalidis, R., Tassis, K., & Andrews, J. J. 2017, MNRAS, 466, 2529

Richardson, L. F. 1922, Weather prediction by numerical process (Cambridge: Cambridge University Press)

Robitaille, T., & Bressert, E. 2012, Astrophysics Source Code Library [record ascl:1208.017]

Robitaille, J.-F., Joncas, G., & Miville-Deschênes, M.-A. 2014, MNRAS, 440, 2726

Robitaille, J. F., Motte, F., Schneider, N., Elia, D., & Bontemps, S. 2019, A&A, 628, A33

Roy, A., André, P., Arzoumanian, D., et al. 2019, A&A, 626, A76

Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998, A&A, 336, 697

Vázquez-Semadeni, E. 1994, ApJ, 423, 681

Vázquez-Semadeni, E., Palau, A., Ballesteros-Paredes, J., Gómez, G. C., & Zamora-Avilés, M. 2019, MNRAS, 490, 3061
## Appendix A: Non-Gaussian segmentation

The non-Gaussian segmentation of the wavelet coefficients is based on an iterative algorithm of denoising. The threshold $\Phi$ splits the non-Gaussian and the Gaussian wavelet coefficients at all scales and directions of the Morlet wavelet transform defined in Eq. (6). The sequence defining $\Phi$ is

$$ \left\{ \begin{array}{l} \Phi_0(\ell, \theta) = \infty \\ \Phi_{n+1}(\ell, \theta) = q\sigma_{\ell,\theta}(\Phi_n(\ell, \theta)), \end{array} \right. \qquad (A.1) $$

where $q$ (that is different from the $q$ in Eq. (11)) is a dimensionless constant that controls how restrictive the definition of non-Gaussianities is. The variance $\sigma_{\ell,\theta}^2$ is defined as

$$ \sigma_{\ell,\theta}^2(\Phi) = \frac{1}{N_{\ell,\theta}(\Phi)} \sum_x \mathbb{L}_{\Phi}(|\tilde{f}_{\ell,\theta}(\mathbf{x})|)|\tilde{f}_{\ell,\theta}(\mathbf{x})|^2, \qquad (A.2) $$

where

$$ \mathbb{L}_{\Phi}(|\tilde{f}_{\ell,\theta}(\mathbf{x})|) = \begin{cases} 1 & \text{if } |\tilde{f}_{\ell,\theta}(\mathbf{x})| < \Phi \\ 0 & \text{else.} \end{cases} \qquad (A.3) $$

and

$$ N_{\ell,\theta}(\Phi) = \sum_{\mathbf{x}} \mathbb{L}_{\Phi}(|\tilde{f}_{\ell,\theta}(\mathbf{x})|). \qquad (A.4) $$

The $q$ parameter varies as a function of $\theta$ and $\ell$. When the algorithm converges to an optimal value for the threshold $\Phi$, 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,

$$ P(\ell) = \frac{\delta\theta}{N_{\theta}N_x} \sum_{\mathbf{x}} \sum_{j=0}^{N_{\theta}-1} |\tilde{f}(\ell, \mathbf{x}, \theta_j)|^2, \qquad (A.5) $$

where $\delta\theta = 2\sqrt{-2\ln 0.75}/|\mathbf{k}_0|$, $N_{\theta} = \Delta\theta/\delta\theta$ is the number of directions $\theta$ needed to sample the Fourier space over the range $\Delta\theta$, and $N_x$ is the number of pixels in the image. The conversion between the spatial scale $\ell$ and the Fourier wavenumber $k$ is made through the relation $k = |\mathbf{k}_0|/\ell$. See Robitaille et al. (2019) for more details.