1 Introduction
Stars are born within the cold (~10 K) and dense (>10⁴ cm⁻³) interstellar medium (ISM), a phase representing a small fraction of galactic mass. This rarity of star-forming gas underscores the importance of understanding its origins within the broader ISM physics, from galactic scales down to the scales of star formation itself, encompassing processes like heating, cooling, magnetic fields, cosmic rays, and gravity. Deciphering the emergence of this dense gas from the predominantly diffuse and warm ISM is crucial to fully understanding the initial conditions of star formation.
Historically, models from the 1980s and 1990s depicted the star-forming ISM as governed by a quasi-static evolution of clouds and clumps, potentially collapsing into dense cores. However, the significance of supersonic motions and dynamic processes has since gained recognition, leading to the gravo-turbulent paradigm. This paradigm posits that pockets of quiescent dense gas can arise from a turbulent sea of low-density gas. More recently, thanks to the Herschel Space Telescope, observations have revealed that dense star-forming gas predominantly resides in filamentary structures. This discovery complicates the traditional view of isotropic, thermal pressure-driven turbulence. In parallel, numerical simulations have also shown the formation of coherent structures like filaments and hubs, considering factors such as gravity, turbulence, magnetic fields, radiation, and thermodynamics. These simulations are essential for evaluating the statistics of intensive variables exhibiting multifractal behavior.
The turbulent nature of the ISM is firmly established by high Reynolds numbers. The need for a more precise method to quantify turbulence properties has led to fractal approaches. Initial studies using monofractal descriptors to describe the ISM’s self-similarity were superseded by the realization that observations deviate from pure monofractal characteristics, aligning with advanced turbulence descriptions. Characteristic scales, typically from subparsec to a few parsecs, were identified, challenging the notion of pure self-similarity. The transition from incoherent to coherent structures is often linked to dense structure formation and turbulence dissipation. If coherent structures display scale invariance, they exhibit multifractal properties, evident in power-law statistics for spatial and temporal correlation functions. This points to critical manifolds, particularly in turbulent dynamical systems. To understand ISM turbulence mechanisms and differentiate between magnetic pressure and gravity effects on structure evolution, advanced multiscale and multifractal analysis techniques are essential.
Computational tools have been developed to extract key multifractality signatures from data and connect them to multiplicative processes and energy cascades in turbulence. Multifractality in physical systems is closely tied to singular measures defined by observables. For instance, in 3D turbulence, the kinetic energy dissipation rate per unit volume, ε˙(x), defines the total energy dissipation E as a measure with singular behavior around points x, characterized by singularity exponents h(x). The energy dissipation E within a ball B(x,r) centered at x with radius r scales as *r**h*(x) as r approaches zero. The spatial intermittence of energy transfer support is directly linked to critical manifolds defined by singularity exponent distributions, implying complex energy partitioning across scales and power-law behavior in physical variable statistics.
Canonical formalisms are commonly used to apply multifractal analysis. These methods rely on statistical averages from large ensembles of system realizations to derive singularity spectra and values. However, astronomy often lacks multi-epoch data, limiting canonical approaches to MHD simulations. Single-realization observations, typical in astronomy, cannot be analyzed using canonical methods.
This study overcomes canonical approach limitations by employing a “microcanonical” multifractal signal analysis. This method evaluates individual microstates within a single realization. We apply this to study turbulence in the Musca filament, observed by Herschel. Musca is a prime example of an isolated filamentary structure within the Chamaeleon-Musca cloud complex, unaffected by stellar feedback. Herschel‘s sensitive space-based observations provide unprecedented spatial and flux dynamic ranges, crucial for statistically studying turbulent flows potentially leading to dense gas and star formation. Herschel‘s photometric surveys cover 70–500 μm with 6″ to 36″ angular resolution, spanning several square degrees per region.
This paper proceeds as follows: Section 2 details the Herschel data used. Section 3 introduces multifractal formalism, its relation to turbulence, intermittency, multiplicative cascades, and both canonical and microcanonical approaches. Section 4 investigates the scale invariance of observation maps. Section 5 introduces sparse filtering for revealing filaments obscured by noise and the 2D structure function approach. Section 6 describes the MHD simulations used. Section 7 presents results from multifractal analysis on Herschel Musca observations and MHD simulations. Section 8 discusses these results, and Section 9 concludes the study.
2 The Musca Cloud: Herschel Flux Density Map
Musca, a prominent 6 pc long filament with a high aspect ratio, is situated at a distance of 140–150 pc. It exhibits a low average column density N, with N(max) ~ 4–8 × 10²¹ cm⁻², and harbors only one protostellar source at its northern end. Herschel observations reveal striations orthogonal to the filament, interpreted as mass inflow along magnetic field lines. Evidence of continuous mass accretion toward Musca is supported by observed low-velocity filament accretion shocks. Pointed observations of ¹³CO and C¹⁸O suggest the Musca filament crest is a single velocity-coherent structure. While some proposed Musca as an edge-on sheet, detailed spectroscopic observations from APEX and SOFIA indicate that the crest is a dense (nH₂ ~ 10⁴ cm⁻³) cylindrical structure, not a low-density sheet. Mass input is attributed to large-scale flows, without clear evidence of striations as inflow channels. Striations might be a result of gas compression by magnetosonic waves, consistent with recent theoretical proposals.
Fig. 1. Musca flux density map from Herschel at 250 μm. The filamentary structure with a high aspect ratio is obvious. Perpendicular to the main ridge of emission are fainter hair-like structures (called striations in Cox et al. 2016) that are mostly attached to the main filament.
The Herschel SPIRE 250 μm data of Musca, shown in Fig. 1, are from the Herschel Gould Belt Survey and published in Cox et al. (2016). This study uses the same dataset, obtained in parallel mode at 10 Hz sampling and 60″ s⁻¹ speed, reduced using modified HIPE version 10 pipeline scripts. Level 1 contexts for orthogonal scan directions were combined using the naive map maker in the de-striper module. Conversion to surface brightness (MJy sr⁻¹ from Jy beam⁻¹) used beam-areas from Neptune measurements (March 2013). ZeroPointCorrection task calibrated maps using Planck satellite data. The angular resolution at 250 μm is 18″. We utilize the 250 μm flux map, not the column density map from prior studies, due to its higher angular resolution. The 250 μm map serves as a good column density proxy in regions minimally affected by stellar feedback, tracing cold dust mixed with cold molecular gas.
3 Canonical and Microcanonical Approaches to Multifractality
Numerous galactic studies using continuum data or atomic/molecular emission lines have demonstrated that the Fourier power spectrum of observed line intensity often follows a power law, at least within certain size ranges. This is commonly interpreted as a consequence of the scale-free, turbulent nature of the ISM. This observed scale invariance suggests an energy cascade, where energy injected at larger scales transfers to smaller scales, naturally explaining the ISM’s complex structure at all scales. However, the Fourier power spectrum, with its single descriptive parameter (spectral slope), is too simplistic to capture the full range of turbulent and scale-free phenomena in the ISM, particularly filamentary structures. Monofractal models, relying on a single fractal dimension parameter characterizing energy transfer geometries, represent the simplest models to describe energy partitioning across turbulent medium scales.
For astronomical modeling, monofractal signals can be generated using stochastic processes like fractional Brownian motion (fBm), especially 2D fBms. These fBms are defined by their Hurst exponent, determining their monofractal properties. Khalil et al. (2006) found a monofractal signature in Canadian Galactic Plane Survey atomic hydrogen data, alongside anisotropic features. However, the complex ISM organization necessitates the multifractal formalism, as crucial in magnetohydrodynamic turbulence as in hydrodynamic turbulence.
Multifractal formalism finds applications in various astronomical contexts, including heliophysical turbulence, the ISM, the Universe’s large-scale structure, galaxy mergers, and gravitational wave detection. This formalism exists in two main forms: canonical and microcanonical. The canonical approach, more prevalent, includes methods like Wavelet Transform Modulus Maxima (WTMM) and cumulant approaches. Microcanonical approaches evolved from less precise methods like counting box and histogram methods to the more efficient geometrico-statistical formulation, which we adopt in this study. The following sections review key multifractal approach aspects, without exhaustive detail.
Hydrodynamic or magnetohydrodynamic turbulence exhibits complex energy partitioning across scales. This complexity is attributed to intermittency, stemming from the intricate geometrical organization of energy support at small scales. Observational evidence supports this, with detected strong velocity shears at subparsec scales. Multifractal models were introduced to describe intermittency, linking physical variable fluctuations (e.g., velocity projection) to the geometry of abrupt change points. Consider a fluid in 2D or 3D medium, with velocity vector v(x) at point x. Measurements can be taken over time at a fixed x or spatially at a fixed time. Experimentally, measuring turbulent velocity field projection v in a fixed direction u is often simpler. This allows evaluation of statistical moments:
〈|(v(x+ru)−v(x))⋅u|q〉,
of the turbulent velocity field v. Averages in Eq. (1) can be temporal or spatial, depending on the setup. Experimentally, moments in Eq. (1) scale as *r*ξ(q) within an inertial range [*r*₁, r₂]:
〈|(v(x+ru)−v(x))⋅u|q〉~rξ(q).
Linear ξ(q) indicates monofractality, while experimentally observed nonlinear ξ(q) signifies multifractal behavior, attributed to intermittency of geometric origin. Let Fh denote the set of points x where velocity increments in direction u scale as *r**h*:
Fh = {x s.t. |(v(x+ru)−v(x))⋅u|~rh when r→0}.
These sets Fh, particularly for small or negative h, have complex geometry related to the nonlinearity of ξ(q). Quantifying Fh geometry involves Hausdorff fractal dimension, D(h):
D(h) = dimH(Fh).
The link between geometry and statistics is given by:
ξ(q) = infh>0(qh+d−D(h)),
where d is embedding space dimension (d = 2 or 3). ξ(q) is the Legendre transform of the singularity spectrum, h↦D(h). In monofractality, ξ(q) is linear, and the spectrum reduces to a single point. Multifractality yields a distribution of fractal dimensions with h. Generally, h = dξdq(q).
This concept extends to any signal s(x) or random process (Xt)t≥0 displaying multiscale behavior (nonlinear moment scaling) or more general multifractal processes. It also applies to measures. The singularity spectrum, therefore, encapsulates substantial turbulence statistics, but its computation presents a fundamental challenge.
Velocity field intermittency, driving scaling and energy partitioning, is described by random multiplicative cascades, proposed by the Russian school. Denoting (δv)r as longitudinal increments |(v(x+ru)−v(x))⋅u|, the cascade can be modeled through the probability distribution for (δv)r:
P(δv)r1 = ∫ℝGr1r2 (u)e−uPe−u(δv)r2 du,
expressing probability distribution scale change via a kernel measure Gr1r2 (r₁ < r₂). Dirac measure kernels would imply scale-invariant distributions. Equation (7) holds if the kernel G satisfies the convolution law (for monotonic scale sequences r₁ < r₂ < … < rn):
Gr1rn = Gr1r2*⋯*Grn−1rn,
implying indefinitely divisible laws P(δv)r, physically interpretable as multiplicative cascades.
3.1 Log-normal and log-Poisson Processes
Exact singularity spectrum computations are limited to a few cases. We present two examples, relevant for physical interpretation and this study: log-normal and log-Poisson processes.
The energy cascading model for turbulence intermittency implies a multiplicative form for energy dissipation rate ε˙(x). Recursive division of a spatial volume and multiplication of dissipation rate by random variables leads to a multiplicative form at each scale. The logarithm becomes a sum of random variables, approximating energy dissipation rate as a log-normal process under central limit theorem assumptions. In this model, space intermittency arises from numerous regions with similar dissipation rates. Kolmogorov introduced the log-normal model in 1962.
However, intermittency models vary based on energy dissipation rate assumptions. Gledzer et al. (1996) introduced an intermittency model involving rare, localized regions of both high and low energy dissipation, similar to She & Leveque (1994), resulting in log-Poisson statistics. This model features dissipatively active and passive localized regions, allowing “dissipation holes.” Log-Poisson processes are particularly relevant here due to their potential to describe filament distributions, as demonstrated in this work.
Singularity spectrum h↦D(h) has closed-form solutions for log-normal and log-Poisson processes. Log-normal processes in ℝᵈ have a parabolic singularity spectrum:
D(h) = d−12(h−hmσh)2,
where ξ(q) = hmq−12σhq2, hm is the mean singularity, and *σ*h is singularity dispersion. Log-Poisson processes in ℝᵈ, translationally invariant and indefinitely divisible, have a singularity spectrum:
D(h) = Dmin+(d−Dmin)ω(h)(1−logω(h)),
where ω(h) = h−hmin(d−Dmin)(−logβ), hmin is minimum singularity, Dmin is the dimension of Fhmin, and β is the dissipation parameter (0 < β < 1), with β = 1+hmind−Dmin. Accordingly, ξ(q) = hminq + (d − Dmin)(1 − *β*q).
Typically, singularity spectrum D(h), Legendre transform of ξ(q), is concave (Fig. 2, left). Figure 2 also displays log-normal (middle) and log-Poisson (right) spectra. Log-normal spectra, being parabolic, are symmetric, while log-Poisson spectra are asymmetric. Asymmetry indicates distinct dynamic properties compared to log-normal models. This study focuses on log-normal and log-Poisson processes due to their geometrical models. Other models exist, like log-α-stable, also with asymmetric spectra.
Fig. 2. Singularity spectra of log-normal and log -Poisson stochastic processes. Left: general form of the singularity spectrum D(h), Legendre transform of ξ(q). See text for a definition of D(h), ξ(q), and Fh. Middle and right: singularity spectrum of a log-normal process with d = 2, hm = 0, and *σ**h** = 0.5 (middle panel) and of a log-Poisson process with d = 2, Dmin = 1.07, and h**min = −0.5 (right panel*).
3.2 Existing Computational Approaches
This section summarizes prior astronomical multifractal approach applications, building upon introductions in Elia et al. (2018), Khalil et al. (2006), McAteer et al. (2007), Salem et al. (2009), Kestener et al. (2010), Macek et al. (2014), Robitaille et al. (2020), and more fundamental presentations in Arneodo et al. (1995) and Turiel et al. (2008) on signal processing.
Equation (3) defines scaling exponents h from the true turbulent velocity field v(x), inaccessible in astronomy and geophysics. Computations rely on observational maps, original signals denoted s. Multifractal formalism in signal processing involves defining and computing singularity exponents h(x) and spectra from s, using general positive measures μ instead of v(x). These can be probability measures from s or additive quantities in s domain.
For a positive measure μ and ball B(x,r), the singularity exponent h(x) is:
h(x) = limr→0logμ(B(x,r))log(r),
encoding limiting scaling information at each point x. It is a pointwise generalization of global scaling parameters like power spectrum slope or Δ-variance. Computing h(x) at each x allows generalizing Eq. (3) and defining sets:
Fh = { x | h(x) = h}.
As in Eq. (4), Fh complexity is measured by fractal dimension D(h). Covering μ support with size r balls, the histogram *N**h**(r), number of balls scaling as r**h***, behaves as:
Nh(r)~r−D(h),
(see Arneodo et al. 1995). Singularity spectrum h ↦ D(h) represents the limiting behavior of *N**h**(r) histograms as r*** → 0, a fundamental tool for analyzing complex, turbulent signals. The main challenge is computing this spectrum from astronomical maps.
Direct methods like box-counting (estimating logμ(B(x,r)) slope vs. log r) are inefficient and error-prone (Arneodo et al. 1995; Chappell & Scalo 2001). Elia et al. (2018) use generalized fractal dimensions *Dq. Covering signal domain with boxes B(xi,r), μi(r*) is the signal proportion inside B(xi,r), treating signal as probability distribution. Partition function is defined:
Z(q,r) = ∑i = 1N(r)μi(r)q.
As r → 0, Z(q,r)~rτ(q), with generalized dimensions Dq = τ(q)q−1. D₀ is box dimension of μ support, D₁ is information dimension, and *Dq encodes correlation integral scaling for q ≥ 2. Arneodo et al. (1995) showed the relation between τ(q) = (q − 1)Dq and singularity spectrum D*(h):
τ(q) = minh(qh−D(h)),
a Legendre transform. Multifractal formalism is thus thermodynamically interpretable: q as inverse Boltzmann temperature, studying self-similarity phases of μ. Limit r → 0 is thermodynamic limit at infinite volume (V = log(1/r)), and hi = logμilogr = −logμilog(1r) corresponds to microstate i‘s energy per unit volume *E*i. Partition function becomes:
Z(q,r) = ∑ie−qEi.
Singularity spectrum h↦D(h) is entropy per unit volume, so spectrum computation is equivalent to entropy per internal energy calculation in a large multibody system. This is a microcanonical formulation, valid in the thermodynamic limit, implemented numerically via box-counting or histogram methods, used in Chhabra et al. (1989), Chappell & Scalo (2001), Elia et al. (2018), and Movahed et al. (2006). However, it yields poor results due to finite size effects (Arneodo et al. 1995). Crucially, result quality verification requires testing on rescaled signal versions. Singularity spectra should be scale-invariant, which box-counting method fails, as shown in Section 4.
Canonical formulations address microcanonical approach issues by evaluating h and D(h) as averages across realizations, i.e., in a canonical ensemble, requiring multiple system realizations. Advanced canonical numerical implementations include WTMM, cumulant analysis (used here for log-correlations), wavelet leaders, and multifractal detrended fluctuation analysis (MDFA).
For example, WTMM on observational map signal s(x) uses analyzing wavelet ψ with vanishing moments to filter long-range correlations. Wavelet-projected signal Tψ(x,r) is evaluated at position x and scale r. A qth order partition function Z(q, r) is built from Tψ such that Z(q, r) ~*r*τ(q) as *r* → 0. Singularity spectrum is D(h) = minq(qh−τ(q)). This relies on accurate τ(q) evaluation via log-log regression with large ensemble data.
ISM studies often lack large observational map ensembles, prompting Elia et al. (2018) to favor microcanonical approaches. However, box-counting/histogram methods poorly estimate singularity spectra. This study develops a more reliable microcanonical implementation, introduced next.
Fig. 3. Increasing sequence of lattice nets Ω₀ ⊂ Ω₁ ⊂⋯ ⊂ Ω ⊂⋯ inside the unit square.
3.3 Singularity Spectrum from Microstates
A key goal of microcanonical multifractal formalism is precise evaluation of h(x) for a measure μ and singularity spectrum h↦D(h) derivation to understand complex, turbulent signal statistics. This section and Appendix A detail singularity spectrum/exponent determination in a microcanonical formulation based on predictability theory in complex systems, overcoming box-counting/histogram method limitations.
Astronomical maps and turbulence data exhibit scaling only within an inertial scale range. Measures must be defined from finite, discrete map values. We follow the multifractal ansatz in disordered systems physics for rigorous discrete signal presentation (Fyodorov 2010; Fyodorov et al. 2012). First, define a measure from observational map data. Prior works often defined measure from signal s itself. Turiel et al. (2008) showed better accuracy from discrete gradient information. Consider an image domain inertial range [r₁, r₂ ], identified with unit square [0,1]². Define refined lattice point collection. For each n∈ℕ, consider Ω set of points with coordinates as integer multiples of 2⁻ⁿ within the unit square, forming an increasing sequence Ω₀ ⊂ Ω₁ ⊂⋯ ⊂ Ω ⊂⋯ (Fig. 3).
For each x ∈ Ω, define a finite set of neighboring points Vn(x) (Fig. 4).
For x ∈ Ω, signal data s(x) and discrete gradient norms ∇d(x,y) = |s(x)−s(y)|‖x−y‖ are available for y∈Vn(x) vicinity points around x at scale 2⁻ⁿ. Define a measure *μn at scale 2⁻ⁿ from discrete gradient information. Appendix A provides exact definitions. μn is a gradient measure associated with observational map data at scale 2⁻ⁿ, meaning μ*n-measure of a ball is the sum of discrete gradient norms within.
Multifractal formalism application to signal s and *μn is valid only if μn* satisfies scaling hypothesis in Eq. (A.9). Canonical multifractal implementations require scaling verification, generally via log-regression on partition functions. Existing microcanonical implementations show scaling only for fractional Brownian motion simulations. Section 5.3 uses 2D structure function methodology to check data scaling.
Key idea in microcanonical singularity exponent computation based on predictability is relating h(x) to signal reconstruction. Singularity exponent values are bounded; lowest value *h*∞ = min{*h*(x)} corresponds to sharpest signal transitions (Eq. (12)). Points x with h(x) = *h*∞ form a fractal subset F∞ encoding strongest transitions. In 2D signals, F∞ encodes strongest edges, while other Fh sets represent smoother transitions as *h* increases. Assuming F∞ coincides with most unpredictable points, h(x) computation at scale 2⁻ⁿ involves only immediate neighboring points around x. Appendix B rigorously formulates this, leading to h(x) computation using correlation measure H, evaluated locally around x, leading to Eqs. (B.5) and (B.7), our fundamental methodology for microcanonical singularity exponent computation.
Fig. 4. Discrete neighborhood sets Vn(x) at scale 2⁻ⁿ around point x.
4 Experimental Comparison of Microcanonical Methodologies: Scale Invariance
This section tests scaling hypothesis validity (Eq. (A.9)) by comparing singularity spectra computed in microcanonical multifractal formalism using two methods: enhanced counting-box method (gradient modulus wavelet projection, Turiel et al. 2006) and the method from Section 3.3. Experiment performed on Musca Herschel 250 μm map.
Scale invariance in Eq. (A.9) implies that singularity spectra from scaled versions of an observational map should coincide. Experiment: take Musca map, generate two downscaled versions, compute singularity spectra with both methods on downscaled versions, check for spectrum coincidence.
4.1 Gradient Modulus Wavelet Projection Method
Algorithm detailed in Turiel et al. (2006). Enhanced counting-box method for computing singularity exponents h(x) and spectrum h↦D(h) via log-regression, not directly on scaling hypothesis (Eq. (A.9)), but on wavelet projections for better accuracy. Allowed because wavelet projections of measure μ scale with same exponents h(x) if analyzing wavelet has n vanishing moments with n > h(x) (Venugopal et al. 2006a; Turiel et al. 2008): if r > 0 is scale, μ measure on ℝ², ψ real wavelet, and *λ**r** measure 1rψ(−xr) dx, wavelet projection of μ at scale r is Tψ(μ,r), convolution of μ and λ**r***:
Tψ(μ,r) = λr*μ.
If μ has density, Tψ(μ,r) density is the continuous wavelet transform of original density with wavelet ψ. Measure *μ*n is from Eq. (A.1) with Vn(x) (Fig. 4, left part, two neighbors). Analyzing wavelet is β-Lorentzian:
Lβ(x) = 1(1+‖x‖2)β,
with β = 3, support scaled to Nyquist frequency. Log-regression over 30 scales.
Scaling check: downscale signal on two consecutive scales below original resolution using discrete wavelet transform (reverse bi-orthogonal projection order 4.4). For each downscaled projection, compute singularity exponents via log-regression with Lorentzian wavelet to get scalar fields h(x) at two resolutions. Compute h ↦ D(h) using Eq. (A.11) from each scalar field. If scale invariance and correct exponents, mappings should coincide. Figure 5 shows log-regression result. For informative part (h ≤ 0), mappings do not coincide, indicating poor log-regression method results and lack of scaling.
Fig. 5. Checking scale invariance: gradient modulus wavelet projection method, log-regression. Top row: visualization of the singularity exponents computed through log-regression using a Lorentz wavelet of order 3 on the Musca 250 μm Herschel flux map presented in Sect. 2. The singularity exponents are computed with the gradient modulus wavelet projection method algorithm. The two upper images correspond to two consecutive downscaled wavelet projections of the original signal using a discrete wavelet projection. Bottom: resulting mappings h ↦ D(h) using the singularity exponents computed on the two consecutive scales (red: one-half scale of the original Musca 250 μm Herschel flux map; green: one-fourth scale of the original Musca 250 μm Herschel flux map). Inside the domain corresponding to h ≤ 0 the two mappings do not coincide. Consequently, the log-regression method for computing the singularity exponents provide poor result in this case.
4.2 Local Correlation Measure of Section 3.3
Figure 6 repeats the Musca 250 μm map experiment with singularity exponents computed using Eqs. (B.7) and (B.5). Resulting h↦D(h) maps show scaling behavior more satisfactorily, coinciding closely for two consecutive scales over singular values. This agreement is stronger than log-regression exponents in Fig. 5. Local correlation measure algorithm exhibits scaling better. Resulting graphs differ from log-regression in Fig. 5, being less parabolic. Log-regression method (gradient modulus wavelet projection) tends to produce log-normal type singularity spectra in this case.
Fig. 6. Checking scale invariance: local correlation measure. Top left and right: visualization of the singularity exponents computed with the local correlation measure described in Eqs. (B.7) and (B.5), respectively, on the Musca 250 μm Herschel flux map presented in Sect. 2. The two images correspond to two consecutive downscaled wavelet projections of the original signal using a discrete wavelet projection. Bottom: resulting mappings h ↦ D(h) using the singularity exponents computed on the two consecutive scales (red: one-half scale of the original Musca Herschel map; green: one-fourth scale of the original Musca Herschel map). The two mappings coincide much better than with the exponents computed from log -regression and displayed in Fig. 5. Consequently, the local correlation measure algorithm exhibits the scaling in a much betterway.
5 Application to the Musca Observation Map, Sparse Filtering, and Determination of the Inertial Range
5.1 Map of the Singularity Exponents
In Section 3, multifractal formalism was introduced via velocity field intermittency analysis, standard in turbulence literature. This study analyzes thermal continuum emission using discrete gradient measures from Section 3.3. Figure 7 displays singularity exponents h(x) on the Herschel 250 μm map, computed with local correlation measure (Eqs. (B.7) and (B.5)), on a filtered map version (Section 5.2) to reduce background noise. Noise reduction improves visualization and singularity spectrum precision. Wedge in Figure 7 makes negative exponents brighter. Singularity exponents trace signal transitions, particularly along filamentary structures, well-rendered by negative h(x) values. Thin dark region along main filament crest indicates constant flux along a cylindrical structure generatrix, also seen in MHD simulations. Musca crest shows small elongated structures perpendicular to the crest, adding complexity. Small circular features likely unremoved background galaxies.
Fig. 7. Map of the singularity exponents computed on filtered Musca 250 μm data. The image is a magnification of the singularity map over the central part of the Musca observation map to better show how the singularity exponents reflect the complex distribution of filamentary coherent structures.
5.2 Sparse Filtering of an Observational Map
Musca Herschel 250 μm flux map contains point-like sources (mostly galaxies), cosmic infrared background (CIB), and cosmic microwave background (CMB), isotropic low-amplitude values close to small-scale filamentary structures. Noise reduction is crucial for multifractal analysis. Robitaille et al. (2019) used filtering to separate monofractal noise from multifractal coherent structures via wavelet coefficient PDF thresholding, identifying monofractal noise and turbulent ISM components. However, this method is not used here as small-scale wavelet decomposition thresholding blurs filaments after reconstruction, potentially impacting singularity spectra. Edge-aware noise filtering (Badri 2015) is preferable, preserving coherent structures while reducing noise.
Compute filtered image sf(x) close to original data s with sparse gradients. Solve optimization problem:
argminsf ‖s−sf‖+λ‖∇sf‖.
Candès et al. (2008) show L¹ norm better promotes sparsity than L² norm while remaining convex, facilitating numerical implementation and global minimum existence. Use L¹ norm ‖∇s‖₁ = ∫ℝ² (|∂s∂x₁|+|∂s∂x₂|)dx₁ dx₂ and solve:
argminsf ‖s−sf‖1+λ‖∇sf‖1,
with λ > 0 tuning parameter. First term ∥s − sf∥₁ ensures filtered image similarity to original (data fitting), second term λ ∥∇sf∥₁ minimizes gradient sparsity, reducing noise while preserving gradient information. L¹ minimization solved using half-quadratic splitting (Geman & Chengda Yang 1995; Schmidt & Roth 2014). Figure 8 shows sparse filtering on flux density data as standard noise reduction. Figure 9 shows enhanced filamentary structures in singularity exponents after filtering compared to unfiltered version. Filtering reduces CIB and CMB, exposing filaments obscured by Gaussian noise. Figure 10 shows noise suppressed by filtering (s − sf). Figure 11 compares singularity spectra with and without filtering. Noise overestimates low-dimensional filaments and Gaussian noise makes unfiltered spectrum more parabolic. Noise reduction is necessary preprocessing for accurate singularity spectrum computation.
Edge-aware filtering (Eq. (22)) does not introduce spurious gradients. Minimization algorithm implementation sets low-level gradients to 0, keeps gradients above a threshold, and generates new image with weighted sum of original and previous step gradient images, not introducing new gradients.
Section 4 and Appendix B checked gradient measure scaling by computing singularity spectra at consecutive wavelet projection scales, observing scaling when exponents are computed using Appendix B method on filtered Musca map (λ = 0.7). Figure 12 shows same test on non-filtered Musca map. Scaled versions generated using reverse bi-orthogonal discrete wavelet transform of order 4.4. Noise alters measure scaling, advocating for filtered maps.
Figure 13 shows q↦ ζ(q) function (structure function method, Eq. (25)) for unfiltered and filtered Musca maps. Distinct graphs indicate noise effects on scaling properties.
Singularity spectrum is estimated from singularity exponent histogram (Eq. (A.11)). Error bars can be estimated. Discretizing histogram with bins, probability *pα of exponent in bin Bα is Nα∕N. For large N, Nα is Gaussian, assigning error bars with 99% confidence interval: ΔD(h) = 3logr1⋅1Nh. Figure 14 shows singularity spectrum with error bars for filtered Musca 250 μm map (λ* = 0.7).
Fig. 8. Visualization of the filtering (computed with λ = 0.7) zoomed in on the central part of the Musca 250 μm map. The level was chosen to emphasize low values. The left panel shows the filtered data and the right the unfiltered data. The positive effect of the filtering on the flux density data becomes obvious. The visualization performed on the singularity exponents (see Fig. 9) demonstrates the importance of filtering on filamentary structures, which are oflow fractal dimensions.
Fig. 9. Illustration of background noise reduction and enhancement of the filamentary structures in the Musca 250 μm map by edge-aware nonlinear filtering while promoting a signal’s gradient sparsity with L¹ norm. The images show the singularity exponents of a region in Musca that is rich in filamentary structures. The right panel displays the singularity exponents derived from the unfiltered original map, while the left panel shows the singularity exponents after filtering the map with λ = 0.7 (Eq. (1)). The same grayscale wedge is used for both images.
Fig. 10. Visualization of noise reduction by the edge-aware noise filtering. The data displayed are s − sf in the same area as shown in Fig. 9.
Fig. 11. Singularity spectrum of the Musca 250 μm Herschel map, computed with and without edge-aware filtering. The orange curve shows the singularity spectrum computed on raw unfiltered data, while the blue curve displays the singularity spectrum computed on edge-aware filtered data using the filtering introduced in Sect. 5.2 with parameter λ = 0.7. The horizontal axis gives the h-values of the scaling exponents, and the vertical axis gives the fractal dimension D(h) = dim(Fh) (see Sect. 3.3). First, it becomes obvious that the part of the spectrum corresponding to h ≤ 0, i.e., the strongest transitions, is not accurately computed on the unfiltered data; the presence of the background noise leads to an overestimation of the multifractal spectrum. Second, the slope of the right part of the graph is poorly evaluated in the presence of noise. As a result, the singularity spectrum computed on the filtered observation map is less symmetrical with regards to the vertical axis h = 0.
Fig. 12. Mappings h↦D(h) computed for two consecutive scales of the unfiltered Musca 250 μm Herschel observational map using a reverse bi-orthogonal discrete wavelet transform of order 4.4. The graphs show that the presenceof the noise alters the scaling of the measure.
Fig. 13. Plot of ζ(q) for the Musca 250 μm Herschel observational map with filtering (λ = 0.7, curve in red)and without filtering (curve in green); the mapping q↦ζ(q) is defined by the structure function method (Eq. (25)). This graph shows that the background noise does affect the statistics of scaling.
Fig. 14. Singularity spectrum of the Musca Herschel 250 μm map with error bars as defined by Eq. (24); the observational map is edge-aware filtered with a tuning parameter λ = 0.7.
5.3 Inertial Range and the 2D Structure Function Method
Canonical multifractality approaches determine scaling laws by statistically analyzing wavelet projection moments over scales and realizations (Venugopal et al. 2006a). Verifying scaling laws via grand ensembles is infeasible with a single Musca image.
Use spatial approach based on 2D structure functions (Renosh et al. 2015) to check for inertial range. Scaling laws exist for spatial distances when:
〈|s(x1)−s(x2)|q〉~‖x1−x2‖ζ(q).
Spatial moments ⟨|s(x₁) − s(x₂)|q⟩ are log-log plotted against distances ∥x₁ −x₂∥ for large (x₁, x₂) ensembles. Experiments on original and filtered data, and MHD simulations (Section 6), used 5 × 10⁸ (x₁, x₂) pairs randomly chosen in Musca map. Figure 15 shows structure function method results for filtered Musca 250 μm map (λ = 0.7). Top left image: log-log plot with moment values q colored, linear regression fits (dotted lines) in inertial range [13, 160] pixels ([0.053, 0.65] pc). Middle image: q↦ζ(q) map. Right images: linear regression quality (sse, Pearson correlation r).
Inertial range is likely affected by beam size (~ 3 pixels, 0.012 pc) for lower boundary (0.05 pc) and map sizes (1786 x 2135 pixels, 7.144 x 8.54 pc) for upper boundary (0.65 pc).
Fig. 15. Result of the 2D structure function method for the determination of an inertial range of spatial scales. Computations are done on the Musca 250 μm map filtered data with λ = 0.7 (Eq. (22)). A total of 5 × 10⁸ couples of points (x₁, x₂) are chosen uniformly in the spatial domain; moments q are between 0.1 and 5, with increment steps of 0.1 (for clarity, only some of the moments are shown). The left panel shows the graph of a log −log plot (Eq. (25)), with colors indicating some moment values, q. The image also shows the graphs of linear regression fits (dotted lines), performed in an inertial range covering the interval [13, 160] pixels, corresponding to the range of distance [0.053, 0.65] pc. The middle panel displays the resulting map q↦ζ(q). The two panels on the right show the quality of the linear regression: the sse and the r correlation coefficient. Hence, the quality of the fit decreases with the moment order q.
6 MHD Simulation Data
Simulations used are from Dib et al. (2007, 2008). Ideal MHD equations solved on uniform 3D cubic grid using TVD scheme (second-order upwind, Kim et al. 1999). Periodic boundary conditions in all directions. Poisson equation solved for self-gravity using Fourier algorithm. Second-order time accuracy via momentum density update due to gravity (Truelove et al. 1998).
Turbulence driven continuously in simulation box, kinetic energy input rate adjusted to maintain constant rms sonic Mach number Ms = 10, as in Stone et al. (1998). Kinetic energy injected at large scales, wavenumber range k = 1–2. Physical units: 4 pc box size, average number density 500 cm⁻³. Column density ~ 5 × 10²¹ cm⁻², similar to molecular clouds except massive star-forming regions. Temperature 11.4 K, sound speed 0.2 km s⁻¹, initial rms velocity 2 km s⁻¹ (Ms = 10). Four simulations vary by initial magnetic field strength, from magnetically subcritical to non-magnetic cloud. Initial magnetic field for subcritical, moderately supercritical, and strongly supercritical models are B₀ = 45.8, 14.5, 4.6 μG, corresponding to β plasma and mass-to-magnetic flux values of β = 0.01, 0.1, 1, and *μ*box = 0.9, 2.8, 8.8. Simulations start with uniform density, evolve for half sound crossing time (ts = 20 Myr), equivalent to five turbulent crossing times (tc = 2 Myr), before gravity is turned on. Common practice for turbulent cascade development. Figure 17 (left) shows integrated column density maps for hydrodynamical case (no magnetic field) at time t = 0.422 tc Myr after gravity onset.
Figure 16 shows inertial range determination using 2D structure function method (Section 5.3) on MHD simulation output. Differences from Musca 250 μm map results are considerable, especially in linear regression quality. Linear regression range [0.18, 1.5] pc. Other MHD outputs show similar plots. MHD simulation output scaling properties differ from observations.
Figure 17 (middle and right) shows singularity exponents and spectra (Section 3.3). Error bars show more uncertainty than Musca map (Fig. 14) due to smaller sample size (256 x 256 image data). Figure 18 shows singularity spectrum (orange) and fitted log-normal spectrum (blue) D(h) = 2−12(h−hmσh)2 (Eq. (9)) for t = 0.62 tc Myr output. Quasi-Newton minimization algorithm used for fitting, yielding hm = −0.0223, *σ*h = 0.2037, so D(h) = 2−(h+0.0223)20.0829. Fit is good, indicating log-normal process approximation. Good fit observed for all simulation outputs. Musca 250 μm map singularity spectrum (Figs. 14 and 19) does not fit log-normal as well. Discussed further in next sections.
Fig. 16. Result of the 2D structure function method, applied over one of the MHD simulation outputs presented in Sect. 6, for the determination of an inertial range of spatial scales. A total of 5 × 10⁸ couples of points (x₁, x₂) are chosen uniformly in the spatial domain; the moments q have values between 0.1 and 5, with increment steps of 0.1 (for clarity, only a few of the moments are shown). The first panel shows the graphs of the log − log plot (Eq. (25)), with colors indicating some moment values, q). As in Fig. 15, the graphs show linear regression fits (dotted lines) performed in an inertial range covering the interval of distance [0.18, 1.5] pc. The second panel displays the map q↦ζ(q). The third and fourth panels show the quality of linear regression: sse and the r correlation coefficient.
Fig. 17. Multifractal analysis of MHD data. Left panel: one snapshot of the column density output of the MHD simulation data used in this work, integrated along one axis of the cube. The output shown here corresponds to the hydrodynamical case with no magnetic field at time t = 0.422 tc. Middle panel: map of the singularity exponents. Right panel: singularity spectrum with its error bars computed using Eq. (24).
Fig. 18. Log-normal fit of the singularity spectrum for the data corresponding to time t = 0.62tc Myr after gravityis turned on. Shown is the singularity spectrum computed in the microcanonical framework (orange) and a fitted quadratic log-normal singularity spectrum (blue). The fit is very good, which indicates a good approximation by a log -normal process.
7 Multifractal Analysis of Data
7.1 Detection of a Multiplicative Cascade
Multiplicative cascade existence is crucial in turbulence phenomenological descriptions (Section 3). Multiplicative cascade models generally have log-correlations C(r,Δx) = −σ²log(Δx), where σ² = Var(log|W|) and W is the cascade-generating random variable (Appendix D). For log-normal cascades, C(r,Δx) = −c₂log(Δx), with c₂ as the second cumulant (Arneodo et al. 1998b). Appendix D summarizes relevant tools. Investigate multiplicative cascade existence in Musca map using cumulant approach (Appendix C). Cumulant analysis on 1D signals from 2D map columns (Fig. 20, left). Figure 20 example shows 1D signal from Musca map column crossing low- and high-flux regimes. Cumulant analysis uses Mexican hat wavelet (Eq. (26)) with deviation σ = 3.2. Wavelet: ψ(x) = 2π¼3σ(x²σ²−1)e−x²2σ² (second Gaussian derivative). Continuous wavelet projection (CWT) of 1D signal f(x) at scale r > 0 is Tψ(f)(x,r) = 1r∫ℝf (u)ψ(x−ur) du. L¹ normalization used instead of L²; energy conservation not necessary for correlations, L¹ normalization better adapts to signal variations (Venugopal et al. 2006a). Use Herschel 250 μm Musca dataset filtered with λ = 0.7 (Eq. (22), Section 5.2). Figure 20 shows Musca column and 1D signal f(x) with CWT projections Tψ(f)(x,r) at scales r₁ and r₂ (r₂ > r₁).
Cumulant analysis on each 1D signal (columns) by computing log-correlations C(r,Δx) (Eq. (D.1)) and plotting vs. log(Δx). Look for long-range correlation fingerprints and multiplicative cascade evidence (C(r,Δx)~logΔx, Eq. (D.2)). Use 500 scales from r = 1 to 500, 400 Δx values from 1 to 400. Slope c₂ computed via linear regression (Eq. (C.3)) for scales 10-200. Long-range correlations found for all 1D columns, multiplicative cascade exists in Musca. Figure 21 (left) example graph (Col. 970): log-correlations for different scales decrease with y = −c₂log(Δx) line, indicating multiplicative cascade. However, some columns (rectangular area, Fig. 22, left) show different behavior (Fig. 21, right). Multiplicative cascade still present, but an area perturbs statistics, suggesting different processes in Musca cloud.
Fig. 20. CWT analysis performed on a 1D signal extracted from the Musca 250 μm map. Left panel: Musca 250 μm map in which the white column defines a 1D signal. Right panel: original 1D signal (centered and rescaled) and CWT projections of the selected 1D signal with the Mexican hat wavelet at scales 5 (red) and 60 (blue).
Fig. 21. Plots of the log-correlations C(r,Δx) vs. log (Δx) for two 1D signals extracted from the Musca 250 μm map according to the scheme shown in Fig. 20. The analyzing wavelet is given by Eq. (27) with σ = 3.2. We use 500 scales from r = 1 to r = 500 and 400 values for Δx between 1 and 400. The chosen values for the scales are enough to highlight more than two nodes in a cascade, if such a cascade exists. Each image shows the graphs of the log-correlations for a few scales (curves from green to purple, corresponding to increasing scales) as well as the line y = − c₂ log(Δx) with the slope c₂ computed through linear regression, according to Eq. (C.3). The linear regression used to compute c₂ is performed for scales between 10 and 200. The left panel shows the log-correlations corresponding to Col. 970 in the Musca observational map. The graphs of C(r,Δx) for different scales follow the slope given by the line y = −c₂log(Δx). The right panel, corresponding to Col. 947, displays a different behavior: long-range correlations are observed, but the behavior seen in a log-normal multiplicative cascade is not present.
Fig. 22. Definition of areas. In the left panel, the statistics of the log-correlations are differentfrom those of the other parts of Musca in the rectangular area shown in the picture. In the right panel, a definition of six ROIs inside the Musca ISM are shown.
7.2 Musca Spectrum: Multifractality, Log-normality
Final singularity spectrum for Musca region in Figure 14. Unlike previous ISM spectra, this spectrum isn’t symmetric, showing a steep descent on the h ≤ 0 side and a rounder Gaussian shape on the h ≥ 0 side. Asymmetry enhanced by noise filtering (Fig. 11). Despite asymmetry, Figure 23 (left) fits a parabolic curve representing log-normal behavior, finding hm = 0.0206, *σ**h** = 0.4564. Spectrum for h ≤ 0 deviates from parabolic fit, suggesting non-log-normal process in Musca turbulence, resembling log-Poisson spectrum (Fig. 23, right). Log-Poisson spectra are asymmetric with steep descent on negative h*** side (Fig. 2). Fit results can be compared to unfiltered case (Fig. 19, bottom). Noise filtering enhances filamentary structures, suggesting asymmetry and log-Poisson behavior relate to 1D filaments, as expected for most singular turbulent cascade structures (She et al. 1990, Section 3.1). Propose ISM in Musca region shows turbulence better represented by dynamical intermittency (localized dissipation) than space intermittency (widespread dissipation). Filamentary structures in all Musca ROIs and ROI1 singularity exponent maps could affect log-normality (Fig. 26). According to She et al. (1990), most intermittent incompressible hydrodynamic turbulence structures are filamentary. In compressible, magnetized turbulence, filaments might be more stable.
Propose ISM with Musca filament shows turbulence better represented by dynamical intermittency (localized enhanced dissipation) than space intermittency. From Bonne et al. (2020a) and filamentary nature of intermittent regions (She et al. 1990; She & Leveque 1994), speculate localized dissipation links to Musca filament and striations, possibly due to magnetic field guiding and focusing matter accumulation, enhancing dissipation in accretion events. Local shocks providing strong turbulence dissipation may explain log-Poisson statistics of large singularities (large h values, large local gradients, white in Fig. 7). Most non-symmetric (potentially log-Poisson) spectra found in Musca filament ROIs, reinforcing dynamical intermittency association with filament formation.
Intriguingly, microcanonical analysis of simulations did not reveal log-Poisson behavior, always yielding log-normal spectra. Could reflect continuous injection of well-behaved large-scale turbulence in simulations. Self-gravity inclusion doesn’t affect spectra. Simulations may lack elements to reproduce observed data, but statistical tool is sensitive enough to differentiate simulations and real datasets.
Fig. 23. log-normal and log -Poisson fitting of Musca singularity spectrum. Left: singularity spectrum of the filtered (λ = 0.7) Musca 250 μm observational map (in orange) fitted against a log-normal parabolic spectrum (blue) defined by Eq. (9). The fit is performed using the quasi-Newton minimization algorithm implementation available in Matlab. Right: singularity spectrum of the filtered (λ = 0.7) Musca 250 μm observational map (in orange) fitted against a log-Poisson singularity spectrum (blue) defined by Eq. (10). The minimization algorithm is the same. Both images display the values found by the minimization algorithm and the error fit.
7.3 Distinct Statistical Properties Observed in the Data
Compare singularity spectra computed over different regions of interest (ROIs) to reveal potential turbulent behavior variations within cloud. Figure 22 shows six ROIs in Musca, selected based on physical properties. ROI1: northern filament end with protostar. ROI2: eastern cloud with striations. ROI3 & ROI4: highest density crest regions; ROI4 southern, fragmented; ROI3 northern, homogeneous crest. ROI5: southern filament end, less organized. ROI6: western embedding cloud, less prominent striations. For each ROI, compute singularity spectrum and log-normal fit parameters hm and *σh** of Dlogn-fit(h) = 2−12(h−hmσh)2. Fit quality estimated by L² error ‖D(h)−Dlogn-fit(h)‖₂². Best-defined curves: ROI2 and ROI5 (most pixels). Figure 24 and Table 1 show results. Graphs and values confirm processes in each ROI can deviate from log-normal. Two curve classes: close to log-normal (ROI1, ROI4, errors 0.028 and 0.19) and deviating significantly. ROI1 and ROI4 have largest fitted log-normal singularity dispersions σh (0.62 and 0.44) and most symmetric spectra. Other ROIs show steeper slope at negative h*** values, ROI2 worst log-normal fit, related to striations in ROI2, rectangular area in Fig. 22, and log-correlation deviations in Section 7.1. For ROIs except ROI1 and ROI4, spectrum closer to log-Poisson than log-normal, also true for full Musca map spectrum (Figs. 14, 19, 23).
Fig. 24. Singularity spectrum and log-normal fit for each ROI. Shown for each region of interest in the Musca cloud delimited in Fig. 22 are the singularity spectrum h↦D(h) of an ROI (in orange), a log-normal fitted approximation singularity spectrum (in blue; Eq. (9)), the obtained numerical values hm and *σ**h**, and the error of the fit. The various values of hm and σ**h*** for each ROI are also given in Table 1.
7.4 Application on Data from Simulations
Figure 18 shows typical log-normal fit result for MHD simulation data projection (Section 6). Fit is very good in all directions and for all simulation outputs, indicating MHD outputs are log-normal processes. Track log-normal spectrum width variation (coefficient *σh** of Eq. (9)) over time. Figure 25 shows dispersion σh vs. time (Myr) for β = 0.1. Dispersion widens as gravitational effects increase (more bound, collapsing structures), as in Elia et al. (2018), where gravity broadens spectrum, also true for compressive turbulence driving, generating larger density contrasts and mimicking gravity. Log-normal dispersion σ**h*** can indicate gravitational effects.
Fig. 25. Variation in time of the fitted *σ*h of the log-normal singularity spectrum (Eq. (9)) for the simulation outputs described in Sect. 6. The value of the magnetic pressure β = 0.1.
8 Discussion of Results: Musca Turbulent Dynamics and Implication on Star Formation Processes
Bonne et al. (2020a) and others suggest several physical processes explain Musca filament formation and ISM behavior: large-scale turbulence injection, small-scale self-gravity, magnetic field pressures/tensions, and local shocks. Long-term goal: seek statistical signatures of these processes using multifractal analysis and Herschel maps unprecedented depth. Facing challenges in deriving clean statistical probes and real data complexity, this paper is a first step toward this goal, providing hints requiring further confirmation via simulations and star-forming region studies.
Multiplicative cascade with significant inertial range (0.05–0.65 pc) found, supporting turbulence interpretation of Herschel dust emission in low column density Musca cloud. Turbulent motions dominant, originating from at least parsec scales.
New microcanonical approach, not requiring multiple realizations, yielded precise singularity spectrum to characterize turbulence. Scale invariance verified, crucial for deriving accurate singularity spectra and quantifying differences between regions and observations/simulations to understand cloud and filament formation physics.
Singularity spectrum confirms Musca ISM and surroundings have multifractal structure, implying intermittency and non-Gaussian behavior, critical for densest, star-forming structure formation. Unlike previous log-normal spectra, Musca spectrum deviates significantly, showing asymmetry. Deviation enhanced by noise filtering (Fig. 11). Spectrum asymmetry suggests non-log-normal process, potentially log-Poisson behavior, linked to dynamical intermittency (localized dissipation). Log-Poisson behavior expected for energy cascade models with rare localized regions of high/low dissipation (dynamical intermittency), while log-normal implies space intermittency (widespread dissipation). Geometric models for log-Poisson spectra are relevant in astronomy, requiring further analysis to link observed statistics with specific processes. Filamentary structures in Musca ROIs and ROI1 singularity exponent maps could influence spectrum log-normality (Fig. 26). Most intermittent incompressible hydrodynamic turbulence structures are filamentary for stability (She et al. 1990). In compressible, magnetized turbulence, filaments likely more stable.
Propose Musca region ISM turbulence better represented by dynamical intermittency with localized enhanced dissipation than space intermittency. From Bonne et al. (2020a) and filamentary intermittent regions (She et al. 1990; She & Leveque 1994), speculate localized dissipation links to Musca filament and striations, possibly due to magnetic field guiding and focusing matter accumulation and enhanced dissipation in accretion events. Local (accretion) shocks, providing strong local turbulence dissipation, may explain log-Poisson statistics of large singularities (large h, large local gradients). Most non-symmetric (possibly log-Poisson) spectra found in Musca filament ROIs, reinforcing dynamical intermittency-filament formation link.
Microcanonical analysis of simulations did not yield log-Poisson behavior, always log-normal spectra. Could reflect continuous injection of well-behaved large-scale turbulence in simulations. Self-gravity inclusion doesn’t affect spectra. Simulations may lack elements to reproduce data, but statistical tool can differentiate simulations and datasets.
Fig. 26. Distribution of the singularity exponents in ROI1 (left) and in the southern region of Musca ROI5 (right). Both regions show complex linear structures, but of different sizes and organizations.
9 Conclusion
This work used a new microcanonical multifractal formalism, based on predictability in complex systems, to analyze ISM turbulence structure in a Musca Herschel map. Effective multifractal formalism use requires checks on observational maps: inertial range determination, scale invariance (spectrum coincidence at different scales), and noise reduction preserving weak coherent structures.
Self-contained study detailed experiments. Scale laws and inertial range verified using 2D structure function methodology. Microcanonical formulation based on direct singularity exponent computation. Singularity spectra derived, scale invariance checked. Background noise log-normalizes signals, making spectra parabolic. L¹ denoising algorithm (edge-aware filtering) preserves low-dimensional features. Cumulant theory used to determine multiplicative cascade presence using log-correlations.
Results show multiplicative cascade with significant inertial range (0.05–0.65 pc), indicating dominant turbulence from scales larger than a parsec. Intermittency precisely studied. Data and analysis challenge log-normality of turbulence singularity spectrum, major result. Suggests Musca filament turbulence exhibits dynamical intermittency with localized dissipation more than space intermittency. Hints of different turbulent behavior and missing simulation physics indicate this study is a first step toward statistical signatures of ISM turbulent physical processes.
Acknowledgements
This work is supported by the GENESIS project (GENeration and Evolution of Structure in the ISm), via the french ANR and the german DFG through grant numbers ANR-16-CE92-0035-01 and DFG1591/2-1. N. Schneider acknowledges support by the German Deutsche Forschungsgemeinschaft, DFG project number SFB 956. L. Bonne acknowledges support by the Région Nouvelle-Aquitaine, France.
Appendix A Definition of the Measure
Appendix details measure definition, using notations from Section 3.3. Consider Vn(x) (Fig. 4, left), each with two points, and discrete measure *μ*n on unit square:
μn = (12n)2∑x∈Ωn∑y∈Vn(x)∇d(x,y)δx.
For each n, *μn is sum of local measures μn = (12n)2∑x∈Ωnμnx with μnx = ∑y∈Vn(x)∇d(x,y)δx. Examine μn behavior as n → +∞. If s* is differentiable:
|s(x)−s(y)| = |〈∇s(x) | (x−y)〉+‖x−y‖ε(x−y)|,
so,
|s(x)−s(y)|‖x−y‖ = |〈∇s(x) | (x−y)‖x−y‖〉+ε(x−y)|.
If f is continuous bounded function, ∫f dμnx→f(x)(|∇s(x)₁|+|∇s(x)₂|) as n → +∞, with ∇s(x) = (∇s(x)₁,∇s(x)₂), because (x−y)‖x−y‖ is unitary vector and neighbor points in Fig. 4. Hence, ∫f dμnx→f(x)‖∇s(x)‖L₁. To get f(x)‖∇s(x)‖L₂, use:
μn = (12n)2∑x∈Ωn∑y∈Vn(x)∇d(x,y)2δx.
Using second set of V(x) (Fig. 4, right), use:
μn = (12n)2∑x∈Ωn12(∑y∈Vn(x)∇d(x,y))δx,
to converge to f(x)‖∇s(x)‖L₁ for ∫f dμnx. If f is bounded continuous, as n → +∞, ∫f dμn→∫f ‖∇s‖ dλ with λ Lebesgue measure on unit square. Measures *μ*n vaguely converge to gradient norm measure density when s is differentiable.
When s is any signal, x ∈ Ω = ⋃n≥₀Ωn, x ∈ Ω for some n, and x ∈ Ω for all k ≥ n. Ball B(x,2¹⁻ⁿ) with n ≥ 1, measure μ(B(x,2¹⁻ⁿ)) = limk→+∞μk(B(x,2¹⁻ⁿ)) = limk→+∞(12k)²(∑z∈Ωk∩B(x,2¹⁻ⁿ)∑yk∈Vk(z), k>n∇d(yk,z)). As n → +∞, last term is sum of discrete differences in B(x,2¹⁻ⁿ). Scaling hypothesis:
μ(B(x,21−n))~A(x)(12n)h(x) as n→+∞,
expressing measure scaling as n → +∞, defining singularity exponent scalar field on dense Ω ⊂ unit square, with continuous extension over unit square if limx→y,x∈Ωh(x) exists.
To count sites with same scaling, define exponent density at resolution 12n:
ρn(h) = ∑x∈Ωnδ(logμ(B(x,21−n))log2n−h).
Histogram of exponents. As n → +∞, ρn(h)~cnr0−D(h), with *cn > 0, h↦D(h) singularity spectrum of gradient norm density measure (d μ = ∥∇s∥ dx**) for differentiable s*.
Sets Fh = { x | h(x) = h}. For complex signals s like turbulence, Fh sets are dense. If some Fh not dense, open ball exists with no points having exponent h, implying forbidden gradient norm in that ball for differentiable case.
Appendix B Singularity Exponents and Predictability
Appendix details singularity exponent h(x) evaluation method yielding better singularity spectra. Based on predictability in complex systems. Predictability time in dynamical systems: Tp = 1λdlog(Δδ), with *λd leading Lyapunov exponent. Aurell et al. (1997) extended this to turbulence, relating to singularity spectrum h ↦ D(h***). Singularity exponents, via spectrum, relate to predictability. Pont et al. (2006, 2013), Turiel et al. (2006, 2008) proposed relating exponent computation to predictability, information content, and reconstructability. Unpredictable points encode system information, forming a set for complete system reconstruction.
Subset F reconstructs signal s if ∇s(x) = G(∇Fs(x)), with G reconstruction functional, ∇F gradient operator restricted to F. Turiel et al. (2008) showed, assuming deterministic, linear, translationally invariant, isotropic G, reconstruction kernel g exists, leading to Fourier space reconstruction:
s^(k) = g^(k)⋅∇Fs^(k).
If F satisfies Eq. (B.3), then div(∇Fcs) = 0, with Fc complementary set of F. Locality of gradient/divergence means point x∈F decision is local. Assumption:
Most unpredictable point set F∞ gives perfect reconstruction (Eq. (B.3)); x∈F∞ decision is local, F∞ is lowest singularity exponent set (F∞ = Fh∞ = {x | h(x) = h∞}).
Under assumption, exponent h(x) computation at lowest scale r₀ = 2⁻ⁿ, instead of log-regression. For differentiable s with gradient ∇s, at scale r₀ = 2⁻ⁿ, ‖∇s‖(x)~r02 μ(B(x,2¹⁻ⁿ))~A(x)(12n)h(x)+2 = A(x)r0h(x)+2. Translational invariance gives average 〈μ(B(x,2¹⁻ⁿ))〉 = C. Then,
h(x)+2 = log(‖∇s‖(x)/〈‖∇s‖〉)logr0−1logr0log(A(x)C).
Quantity h˜(x) = log(‖∇s‖(x)/〈‖∇s‖〉)logr0 is exponent estimate at x if r₀ is small enough for log(A(x)C) variations to be negligible. Then h(x) = h˜(x)−2.
For arbitrary signal s and lowest scale r₀ = 2⁻ⁿ, exponent h(x) is similarly derived from h˜(x) = log(H(μn,x,r0)/〈H(μn,x,r0)〉)logr0, with H function of measure *μn, using local information around x at scale r₀. 〈H(μn,x,r0)〉 is spatial average of H over map, lessening μn* fluctuations at scale r₀. H measures local unpredictability, subtracting signal value at x from neighbor-inferred value, after detrending. H is local correlation measure. Computation:
Consider point x for h(x) and local neighborhood W(x) at scale r₀ = 2⁻ⁿ (Fig. B.1).
Vector (s(x), s(x₁), s(x₂), s(x₃), s(x₄)) is detrended using d = 13(s(x)+s(x1)+s(x2)+s(x3)+s(x4)), yielding detrended vector (p₀(x), p₁(x), p₂(x), p₃(x), p₄(x)). Define εx(x) = p₂(x)−p₁(x), εy(x) = p₄(x)−p₃(x), related to Wx(x), Wy(x) (Fig. B.1). Local correlation measure:
H(μn,x,r0) = ((εx(x)²+εy(x)²)(A(x,r0)μn(W(x))))½,
with A(x,r0) = |εx(x)(∑y∈W(x)εx(y))+εy(x)(∑y∈W(x)εy(y))|. Value h(x) computed using Eq. (B.5) as h(x) = h˜(x)−2. Exponents define sets Fh = {x | h(x) = h}, important in turbulent system description (Frisch 1995). Turiel et al. (2008) derived log-Poisson process argument generalizing She & Leveque (1994), describing multiplicative cascade from geometrical Fh organization. General justification for physical processes is still lacking.
Fig. B.1. Neighborhood system W(x) considered at a point x for the computation of local correlation measure H(μ,x,r0). It consists of a cross made of the horizontal and vertical sets Wx(x), Wy(x) with other neighboring points, including x.
Algorithm for Singularity Spectrum Computation
Equation (A.11) relates singularity spectrum h ↦ D(h) to logarithm of histogram *ρn(h), but practical computation requires determining cn* in Eq. (A.11). Physical observation: there’s a singularity exponent value h₀ such that Fh₀ has dimension of observational map (D(h₀) = 2). From (A.11), ρ₀ = ρn(h₀) = cnr0−2. D(h) ≤ 2, so h₀ is most probable exponent, maximum of *ρn (h***) (or ρn(h˜)). From Eq. (A.11), ρn(h) = ρ₀r02−D(h), so D(h) = log(ρn(h)/ρ₀)−log(r₀)+2. Algorithm 1 computes singularity spectrum.
Algorithm 1:
- Compute singularity exponents h(x) for each point x in observational map using Eqs. (B.5) and (B.7).
- Compute histogram *ρn(h) of exponents h(x***).
- Find maximum ρ₀ of *ρn(h***).
- Compute singularity spectrum using D(h) = log(ρn(h)/ρ₀)−log(r₀)+2, with r₀ = 2⁻ⁿ lowest signal resolution scale.
Appendix C Cumulant Analysis Method
Appendix details canonical multifractality description using cumulant approach (Delour et al. 2001). Used to determine multiplicative cascade existence (Arneodo et al. 1998a, 1999b) and long-range correlations (Venugopal et al. 2006a). For probability measure μ on line, characteristic function fμ(z) = Eμ(eiz) = ∑n = 0∞Mn(iz)nn!, with moments *Mn = ∫ xndμ(x). Cumulant generating function gμ(z) = logfμ(z) = ∑n = 0∞Cn(iz)nn!, with cumulants Cn. First cumulants: C₁=M₁, C₂=M₂−M₁², C₃=M₃−3M₂M₁+2M₁³. Recurrence: Cn = Mn−∑l = 1n−1(n−1l−1 )ClMn−l. Relation between Cn and multifractal spectrum when μ is multifractal measure. Let s be signal, ψ wavelet, Tψ(x,r) wavelet projection, Z(q*, r) partition function. Measure μ density log|Tψ(x,r)|. Let *Cn (r) be its cumulants. Slope of Cn(r*) vs. log(r) gives, as r → 0, Cn~(−1)n+1cnlog(r), with τ(q) = −c₀ + c₁q − c₂q²/2! + c₃q³/3! + ⋯. For 1D log-normal process, *cn = 0 for n > 2, D(h) = c₀−(h−c₁)2c₂, with c₀ = 1, c₁ = h***m, c₂ = σh².
Appendix D Two-Point Magnitude Statistical Analysis and Multiplicative Cascade
Two-point correlation function for scale r and spatial interval Δx: C(r,Δx) = 〈 (log|Tψ(x,r)|−〈log|Tψ(x,r)|〉 )⋅ (log|Tψ(x+Δx,r)|−〈log|Tψ(x,r)|〉)〉. If C(r,Δx) is logarithmic in Δx and scale-independent for Δx > r (C(r,Δx)~logΔx), long-range dependence exists (Arneodo et al. 1998a,b; Venugopal et al. 2006a). For log-normal multiplicative cascade, C(r,Δx)~−c₂logΔx. Self-similar processes can exist without multiplicative cascade (Arneodo et al. 1999a,b; Venugopal et al. 2006a), often neglected in astronomical multifractal analysis, crucial for proving multiplicative cascade existence in ISM.
