11institutetext:  Laboratoire de Physique de l’École normale supérieure, ENS,
Université PSL, CNRS, Sorbonne Université, Université Paris Cité, Paris,
France 22institutetext:  Observatoire de Paris, Université PSL, Sorbonne
Université, LERMA, CNRS UMR 8112, 75014 Paris, France

# Molecular clouds: do they deserve a non-Gaussian description?

Pablo Richard Molecular clouds: do they deserve a non-Gaussian
description?Molecular clouds: do they deserve a non-Gaussian description?
Erwan Allys Molecular clouds: do they deserve a non-Gaussian
description?Molecular clouds: do they deserve a non-Gaussian description?
François Levrier Molecular clouds: do they deserve a non-Gaussian
description?Molecular clouds: do they deserve a non-Gaussian description?
Antoine Gusdorf Molecular clouds: do they deserve a non-Gaussian
description?Molecular clouds: do they deserve a non-Gaussian
description?Molecular clouds: do they deserve a non-Gaussian
description?Molecular clouds: do they deserve a non-Gaussian description?
Constant Auclair Molecular clouds: do they deserve a non-Gaussian
description?Molecular clouds: do they deserve a non-Gaussian description?

(Submitted to A&A on July 13, 2024.)

Molecular clouds are astrophysical objects whose complex non-linear dynamics
is reflected in their complex morphological features. Many studies,
investigating the bridge between higher-order statistics and physical
properties, have shown the interest provided by non-Gaussian morphological
features to grasp physical information. Yet, as this bridge is usually
characterized in the supervised world of simulations, transferring it onto
observations can be hazardous, especially when the discrepancy between
simulations and observations remains unknown. In this paper, we aim at
evaluating, directly from the observation data, the discriminating ability of
a set of statistics. To do so, we develop a test that allows to compare, for a
given unlabeled dataset, the informative power of two sets of summary
statistics, relying on a notion of statistical compatibility. Contrary to
supervised approaches, this test does not require the knowledge of any class
or parameter corresponding to the data, but focuses instead on comparing the
degeneracy levels of these descriptors, relying on a notion of statistical
compatibility. We apply this test to column density maps of 14 nearby
molecular clouds observed by Herschel, and iteratively compare different sets
of usual summary statistics. We show that a standard Gaussian description of
these clouds is highly degenerate but can be substantially improved when being
estimated on the logarithm of the maps. This illustrates that low-order
statistics, properly used, remain a very powerful tool. We then further show
that such descriptions still exhibit a small quantity of degeneracies, some of
which are lifted by the higher order statistics provided by reduced wavelet
scattering transforms. This property of observations quantitatively differs
from state-of-the-art simulations of dense molecular cloud collapse and is not
reproduced by log fractional Brownian motion models. Finally we show how the
summary statistics identified can be cooperatively used to build a
morphological distance, which is evaluated visually, and gives very
satisfactory results.

###### Key Words.:

ISM: clouds – ISM: structure – Submillimeter: ISM – Methods: statistical

##  1 Introduction

Molecular clouds (MCs) play a key role in star formation. Their highly non-
linear dynamics tends to couple spatial scales over a wide range, to create
dense filamentary structures, in which clumps form that will eventually lead
to prestellar cores (McKee & Ostriker, 2007). However, the precise role of
each physical ingredient (such as turbulence, magnetic fields, and gravity, to
name but a few) in this dynamics still remains to be fully understood.

Figure 1: Column density map of a region in Perseus (top row) and sample of a
log-Fractional Brownian Motion field (bottom row). Each image has 256 $\times$
512 pixels but is tiled into 2 complementary patches of size 256 $\times$ 256
on which the statistics shown in Fig. 2 are computed. While these one-point
and two-point statistics (some being non-Gaussian) are clearly compatible,
these two images have manifestly different morphologies.

Figure 2: Statistical properties of the maps from Fig. 1: PDFs (top), power
spectra (center), and power spectra of the logarithms of the maps (bottom).
The dashed vertical line in the power spectra plots is located at the
$18.2\arcsec$ resolution of the column density map from the observational
dataset. Power spectra are apodized as explained in Sec. 4. We see that the
logFBM process is compatible with this observed portion of MC for all these
statistics (some of theme being non-Gaussian), while having a manifestly
different morphology, as revealed in Fig. 1. All these statistics are thus
degenerate for such a comparison.

A key step towards such an understanding certainly resides in our ability to
decipher, in the morphology of these structured clouds, some signature of the
physical processes at play. For instance, the power-law shape of the power
spectrum in emission maps may trace properties of the turbulence (Miville-
Deschênes et al. (2007); Federrath & Klessen (2013)), or the shape of the
probability distribution function (PDF) of column density maps (Vázquez-
Semadeni et al., 1997; Hennebelle & Chabrier, 2008; Kainulainen et al., 2009;
Federrath & Klessen, 2013; Schneider et al., 2022; Appel et al., 2022) may
trace the impact of gravitational collapse or stellar feedback as it
transitions from a log-normal shape to the development of a heavy tail.

However, the strong non-linear interplay between the physical processes at
work leads to the emergence of non-Gaussian features in interstellar
structures, such as filaments. This has motivated numerous studies to capture
physical information beyond one-point and two-point statistics. Many
approaches have been investigated including, for instance, diagnostics of the
phase coherence of the Fourier modes (Levrier et al., 2006; Burkhart &
Lazarian, 2016), bispectrum (Burkhart et al., 2009), structure functions
(Heyer & Brunt, 2004), dendograms (Goodman et al., 2009), multiscale
segmentation (Robitaille et al., 2019), scattering transforms (Allys et al.,
2019; Regaldo-Saint Blancard et al., 2020; Saydjari et al., 2021), and neural
networks (Peek & Burkhart, 2019; Zavagno et al., 2023).

If higher-order statistics represent an appealing way of refining the
geometrical description of complex ISM structures, they are nevertheless
subject to the following caveats:

  1. 1.

They should not screen out the potential of low-order statistics, but should
come as complementary diagnostics (Burkhart & Lazarian, 2016). Low-order
descriptors should not be underestimated: it is not because the processes
under study are highly non-Gaussian that low-order statistics are only
marginally informative. Yet, in numerous studies, such easy-access information
is not considered, or even intentionally discarded111The mean and variance of
the data may be normalized, the power spectra flattened, and histograms (non-
linearly) equalized., to emphasize the contribution of high-order statistics.

  2. 2.

Refining the statistical description with higher-order terms complicates the
connection between the statistical descriptions and the physical processes and
associated parameters. Such a connection is thus usually learned in the
supervised world of simulations, but transferring it to observations can be
hazardous, especially when their distance to simulations remains unknown
(Falgarone et al., 2004; Peek & Burkhart, 2019).

Figure 3: Footprints of the Herschel Gould Belt Survey (HGBS) fields used in
this study, overlaid on the total thermal dust intensity at
$353\,\mathrm{GHz}$ from the GNILC (Remazeilles et al., 2011) variable
resolution data of Planck (Planck Collaboration, 2020).

In this paper, we aim at identifying what sets of summary statistics are
relevant to characterize the diversity of observed molecular clouds, without
relying on simulations nor prior knowledge. In such an unsupervised framework,
various observations cannot be grouped a priori in a similar class, so one
often has to work in a very low data regime, which restricts dramatically the
range of tools that can be used. This unsupervised and prior-knowledge free
framework excludes features that require significant tuning to extract
information, such as neural networks. This priority shift, from theoretical
informativeness to actual information retrievable in this framework, might
reshuffle the cards among the various sets of summary statistics to consider.
This leads us to wonder to what extent non-Gaussian statistics can actually
bring meaningful contributions in a fully observation-driven and unsupervised
framework.

To answer this question, we rely on a dataset of column density maps
constructed from a survey of nearby molecular clouds. Our approach is to
identify the amount of degeneracy that a set of summary statistics can have on
this dataset, i.e., finding a pair of maps in this dataset that are compatible
according to these statistics, but actually have a different morphology. A
simplified illustration of such a degeneracy is given in Figs. 1 and 2, where
we compare an observed column density map with a realization of a field whose
logarithm is sampled from a Fractional Brownian Motion, that is a specific
type of Gaussian process with power-law type constrained power spectrum. In
this example, the PDFs of these maps, their power spectra, and the power
spectra of their logarithms, all seem to be compatible, even though the maps
clearly have different structures. This illustration underlines the limitation
of those statistics, but only for this specific comparison, where we compare
one observation to a specific known model. However, in an unsupervised
framework, when working solely with unlabeled observations, we will exhibit
these limitations by introducing additional statistics that could be used to
lift these degeneracies.

In this work, we use such a degeneracy diagnostic on the observed molecular
clouds dataset to progressively build a robust set of summary statistics. We
start with basic statistics, such as the power spectrum and one-point
statistics. We evidence strong limitations that we investigate and
substantially mitigate by taking the logarithm of the maps beforehand. We
evidence further but moderate limitations of this improved set, leveraging
higher-order statistics through the reduced wavelet scattering transform
(RWST). This study, initially carried out only on the observational dataset,
is then extended to logFBM data, as well as to a set of magnetohydrodynamics
(MHD) simulations intended to reproduce observations of dense star-forming
molecular clouds. We conclude by discussing the relevance and limitations of
the final set of statistics we have constructed, and show how it can be used
to assess the distance, in a statistical sense, between different maps of the
interstellar medium.

The structure of the paper is as follows:

  * •

In section 2, we present the set of MCs considered in this work and the
corresponding column density maps, as well as the numerical simulations and
the logFBM models.

  * •

We present in section 3 the diagnostic of statistical compatibility that we
will rely on, and our general methodology.

  * •

We present in section 4 the sets of statistics that will be confronted.

  * •

We apply this methodology in section 5 to confront these sets of statistics on
observations, and from these results we design an informative and low-
dimensional set of summary statistics $\phi_{\text{final}}$.

  * •

We finally define, in section 6, a distance based on $\phi_{\text{final}}$
that allows us to compare datasets, such as observations and simulations.

##  2 Data

In this section, we present the main dataset that we use for this study: an
ensemble of $\sim$550 molecular hydrogen ($\mathrm{H}_{2}$) column density
maps from a survey of nearby molecular clouds (MCs). Three other datasets will
be used in this paper:

  * •

a set of $\sim$230 total gas column density maps built from
magnetohydrodynamical (MHD) numerical simulations of dense molecular clouds,
classified into three subsets with varying values of the magnetic field, from
null (_hydro_), to medium (_MHD_), and high (_MHD high B_). These simulations
are state-of-the art attempts to reproduce the physics at play in the early
stages of star formation, including self-gravity and decaying turbulence, and
are therefore well adapted to compare to our observations of the HGBS clouds.

  * •

A set of $\sim$500 synthetic maps sampled from logFBM models, whose parameters
reproduce the diversity of the observational dataset. This type of purely
synthetic fields has already been extensively used to model the interstellar
medium (see, e.g., Elmegreen, 2002; Brunt & Heyer, 2002; Miville-Deschênes et
al., 2007; Levrier et al., 2018).

  * •

A set of $\sim$1000 images from a large collection of everyday textures, the
Describable Texture Dataset (DTD, Cimpoi et al. (2014)). We use these to
emphasize the specificity of ISM fields in the context of image texture
analysis.

More details about these sets are given in appendix A. All these maps have a
size of 512 $\times$ 512 pixels.

###  2.1 Observations: column density maps from the HGBS

We focus our study on a set of MCs targeted by the Gould Belt Survey (HGBS)
(André et al., 2010) with the Herschel Space Observatory (Pilbratt et al.,
2010), whose footprints on the sky are shown in Fig. 3. The HGBS combines two
main criteria we require for this work:

  * •

it sampled numerous MCs with a diversity of physical and environmental
conditions, from diffuse and quiescent regions with no sign of star-formation
activity such as the Polaris Flare (Heithausen & Thaddeus, 1990; André et al.,
2010; Miville-Deschênes et al., 2010) to very dense and active ones such as
Orion B (Schneider et al., 2013) or the Aquila Rift cloud (Könyves et al.,
2015). Examples of molecular column density maps from the HGBS are shown in
Fig. 4. We however emphasize that this survey is limited to local clouds
(distances $d\leq$ 500 pc) and does not cover the whole range of conditions
expected in Galactic molecular clouds222A further study could target more
distant, more massive star forming regions with interferometric observations.
We postpone such a study to a future paper..

  * •

It imaged a broad range of scales, from the full cloud size ($\sim 10\,\,$pc,
corresponding to a few degrees for these nearby clouds) down to the $\sim
0.1\,\,$pc scales of filaments (Arzoumanian et al., 2011; André et al., 2014),
which are spatially resolved for the nearest clouds333The typical resolutions
range from $10\,\arcsec$ at 70 $\mu$m to $36\,\arcsec$ at 500 $\mu$m., thus
covering more than 2 orders of magnitude in spatial scales. This allows us to
perform an in-depth morphological analysis, based on a local description of
multi-scale interactions.

We consider the high resolution (18 ″) column density maps of 14 regions,
produced with the procedure described in Appendix A of Palmeirim et al. (2013)
and publicly available from the Herschel Gould Belt Survey
Archive444http://www.herschel.fr/cea/gouldbelt/en/Phocea/Vie_des_labos/Ast/
ast_visu.php?id_ast=66. This 18 $\arcsec$ angular resolution corresponds to a
$12\,$mpc spatial resolution for the nearest clouds, such as Ophiuchus and
Taurus ($d\sim 140\,$pc), and up to $40\,$mpc for the most distant clouds,
such as Orion B ($d\sim 450\,$pc). The main properties of these clouds are
given in Table 1.

Figure 4: Column density maps of three regions from the HGBS. The spatial
resolution of the maps is $18\,\arcsec$, sampled with a pixel size of 3
$\arcsec$. A white dot with diameter $5\times 18\,\arcsec$ is shown at the
bottom right corner of each map. Estimates of spatial scales, also shown on
the maps, are based on reported distances (see Table 1). The Polaris Flare is
an example of a diffuse and quiescent cloud, while Aquila and Orion B are
dense and very active star-forming regions.

###  2.2 Subsampling and tiling

The presence of physical processes such as gravity in the dynamics of MCs
implies spatial variations of their statistical properties, which prevents us
from modeling them as stationary stochastic processes. For instance, they
usually have strong local overdensities, while being surrounded by diffuse
borders. This non-homogeneity makes it difficult to compare such objects as a
whole, as well as to properly estimate their statistical properties, when seen
as realizations of a random process, as these properties will for instance
depend heavily on the identification of the cloud boundary, or on the precise
definition of what a cloud is. Furthermore, variance estimates for random
processes generally rely on a homogeneity assumption to estimate the intrinsic
variability of a process from its spatial variations. More broadly, the
estimation of the statistical moments and properties of a process is usually
based on spatial averaging, which assumes the homogeneity of the sample.

To avoid this problem, we choose to restrict the comparison to local patches
of these MCs, assuming statistical homogeneity within each individual patch.
This requires the identification of a characteristic stationarity length over
which it can be assumed that the statistical properties of the cloud do not
vary significantly. Different patches extracted from the same cloud may have
different statistical properties, which will be representative of the non-
homogeneity of the cloud as a whole. This encourages us to use small patches.
However, an additional difficulty is that these patches must be large enough
to make statistical estimates. In particular, variance estimates require
sampling beyond the correlation length of the process under study.

In this paper, we have chosen to cut patches with angular sizes ranging from
$0.85$° to $2.55$°, that correspond respectively to $3\,$pc and $9\,$pc for a
cloud at a typical distance of $200\,$pc. We believe this represents a good
compromise with regard to the trade-off expressed above. However, we are aware
that there may not be an ideal solution, especially as we have no precise
estimate of neither the stationarity length nor the correlation length. This
means that we cannot rule out some confusion between local intrinsic
variability and large-scale variations in the statistical properties of the
clouds. However, we believe that, for lack of a better solution, this does not
undermine the relevance of this study.

More precisely, we make multiple versions of the original (18 ″ resolution, 3
″pixel size) column density maps at different pixel sizes (6 ″, 9 ″, 12 ″, and
18 ″). The down-sampling is done using a bivariate spline approximation of
order 3 of the original maps. Then, we cut from these maps patches of size in
pixels $512\times 512$, with a step size of 256 pixels, such that two
neighboring patches will have 50% of their pixels in common. In the following,
these patches will be themselves tiled into 4 sub patches of size $256\times
256$ to perform a variability estimation. The choice of 6 ″ for the finest
sampling allows us to exploit the 18 ″ resolution of these maps with mitigated
sampling artifacts. The 9 ″, 12 ″, and 18 ″pixel sizes respectively correspond
to a relative shrinking of the 6 ″ pixel maps with ratios of 1.5, 2, and 3,
which allows to accommodate for the significant and quite uncertain range of
distances of the different MCs in the HGBS. This entire procedure, which is
illustrated in Fig. 5, and summarized in Table 1, leads to a total of $551$
patches with size $512\times 512$ pixels, each of which is then subdivided
into four $256\times 256$ sub-patches.

Figure 5: Illustration of the pre-processing of observational data. Details are given in Sec. 2.2. Region | Dist. | Coord. | # patches with resolution | References
---|---|---|---|---
|  $d$ (pc) |  $l$ (°) |  $b$ (°) | Tot. | 6” | 9” | 12” | 18” |
Aquila | 436 | 28 | 4 | 57 | 36 | 16 | 4 | 1 | Ortiz-León et al. (2018)
Chamaeleon I | 210 | 296 | -16 | 28 | 20 | 6 | 2 | 0 | Zucker et al. (2020)
Chamaeleon II III | 160 | 302 | -16 | 43 | 32 | 8 | 3 | 0 | Zucker et al. (2020)
Corona Australis | 160 | 0 | -19 | 44 | 34 | 7 | 3 | 0 | Zucker et al. (2020)
IC 5146 | 750 | 94 | -5 | 4 | 4 | 0 | 0 | 0 | Zucker et al. (2020)
Lupus I | 182 | 339 | 17 | 17 | 12 | 4 | 1 | 0 | Galli et al. (2013)
Musca | 200 | 301 | -9 | 6 | 6 | 0 | 0 | 0 | Knude & Hog (1998)
Ophiuchus | 139 | 353 | 17 | 61 | 42 | 12 | 6 | 1 | Mamajek (2008)
Orion B | 450 | 205 | -14 | 88 | 61 | 19 | 7 | 1 | Schlafly et al. (2014); Zucker et al. (2020)
Perseus | 285 | 159 | -20 | 44 | 32 | 10 | 2 | 0 | Zucker et al. (2020)
Pipe Nebula | 180 | 359 | 6 | 23 | 18 | 5 | 0 | 0 | Zucker et al. (2020)
Polaris Flare | 400 | 123 | 26 | 50 | 36 | 9 | 4 | 1 | Schlafly et al. (2014); Zucker et al. (2020)
Serpens | 436 | 31 | 4 | 46 | 33 | 9 | 4 | 0 | Ortiz-León et al. (2018)
Taurus | 145 | 170 | -16 | 40 | 30 | 6 | 4 | 0 | Yan et al. (2019)
Table 1: Summary of the properties of the different MCs studied in this paper,
as well as their division into patches at different sampling. References in
the rightmost column are for the distance estimates, which can be quite
uncertain (uncertainties are given in the references quoted).

##  3 Quantifying informative power of summary statistics on an unlabeled
dataset

###  3.1 General methodology

We want to quantify the amount of information we can exploit from a given set
of summary statistics on a given dataset. However, because we work with
unlabeled data, we cannot use supervised frameworks, such as Fisher analysis,
that have a label-based approach to quantify information. For instance, when
working with simulations, each sampled data ${x_{i}}$ can be labeled by its
corresponding physical parameters $\theta_{i}$. In our case, we have to deal
with an ensemble of unlabeled maps, for which we have a priori no notion of
distance between pairs.

In this paper, we choose to rely on a notion of compatibility between patches,
that can be estimated for any set of summary statistics. This approach allows
to be quantitative, even in this unsupervised setting. Nevertheless, without
supervision, it remains difficult to interpret a compatibility result between
two patches for a given set of summary statistics: do the two patches actually
have very similar properties, or do they have very distinct properties whose
differences are not effectively caught by the statistics?

This difficulty can be partially overcome by using complementary sets of
statistics. Indeed, if two patches are distinguished by a first set of
statistics, this is sufficient to assess that they are in practice different,
and therefore to highlight the degeneracies of another set of statistics. It
is this comparative approach, which is all the more relevant when the panel of
statistics compared is comprehensive, that we use in this paper.

The informative power of a set of statistics strongly depends on the family of
processes studied. For instance, we know that the empirical mean and power
spectrum are sufficient statistics for stationary Gaussian fields (Cover,
1999), which is not the case for physical processes in general. In this paper,
we will apply our approach to the dataset of observations defined above, as
well as to numerical simulations and synthetic logFBM models. In addition, we
note that this diagnostic can also be used to compare maps from two different
datasets.

In the rest of this section, we introduce the compatibility criterion, and
explain how we extend it from a pair level to a dataset level. The set of
statistics used in this paper will be presented in the following section.

Figure 6: Illustration of the proposed test to confront two sets of summary
statistics, $\phi_{A}$ vs $\phi_{B}$, on their degeneracy level for a given
dataset. Each star represents a pair of patches. Panel a: the presence of
stars in the bottom right region, i.e., where $d_{\phi_{A}}\gg 1$ and
$d_{\phi_{B}}\lessapprox 1$, reveals that some pairs of this dataset are
identified by $\phi_{A}$ as incompatible but not by $\phi_{B}$: such pairs
thus evidence degeneracies of $\phi_{B}$ lifted by $\phi_{A}$. Conversely, the
presence of stars in the top left region evidences degeneracies of $\phi_{A}$
lifted by $\phi_{B}$. Hence, if all the points land mainly in the sub-diagonal
part (panel b), this evidences that $\phi_{A}$ is better suited than
$\phi_{B}$ to compare the pairs of this dataset. If on the contrary the points
are spread both in the upper left and bottom right regions of the plot (panel
c), this shows that both $\phi_{A}$ and $\phi_{B}$ are individually not
sufficient to describe the processes of this dataset.

###  3.2 Statistical compatibility for a pair of patches

We want to measure a notion of compatibility, according to a given set of
summary statistics $\phi$, between the two processes that generated the
patches $(x_{i},x_{j})$. To estimate this $\phi$-compatibility, we need to
make a number of simplifications, given the low-data regime, which will bring
us back to a simplified case of statistical hypothesis testing. To do so, we
tile each patch $x_{i}$ into 4 sub patches $\\{x_{i}^{(l)}\\}_{1\leq l\leq
4}$, as illustrated in Fig. 5. We then compute the statistics at this sub
patch level: $\\{\phi(x_{i}^{(l)})\\}_{1\leq l\leq 4}$ and we assume that
these random variables can be considered as independent samples of the same
distribution, that we furthermore model as multivariate normal
distribution555This assumption is not far from being true when the law of
large numbers can be applied to the distribution of $\phi(x)$. This is for
instance the case when $\phi(x)$ is defined as an average over the image
pixels of a certain local distortion (such as filtering, possibly non linear)
$\phi_{loc}(x)$ of $x$:
$\phi(x)=\langle\phi_{loc}(x)[\@vec{u}]\rangle_{\@vec{u}}$ and when the
process $\phi_{loc}(x)$ is stationary with a correlation length that is small
compared to the image length. of mean $\mu_{i}$ and variance $\Sigma_{i}$:

$\phi(x_{i}^{(l)})\sim\mathcal{N}(\mu_{i},\Sigma_{i}).$ (1)

Under these assumptions, the problem boils down to testing the compatibility
between the two normal distributions, i.e., to test the hypothesis:
$\mu_{i}=\mu_{j}$ and $\Sigma_{i}=\Sigma_{j}$. In our case, however, we have
to estimate this compatibility from very few samples, most of the time fewer
than the dimension of the vector of statistics $\phi$. We thus choose to focus
only on testing if the means $\mu_{i}$ and $\mu_{j}$ are statistically
compatible, but not $\Sigma_{i}$ and $\Sigma_{j}$, a problem known as the
_multivariate two-sample mean test_. The most widespread test statistic for
this problem is the Hotelling’s two-sample $T^{2}$-statistic, a multivariate
extension of Student’s $t$-test (Hotelling, 1931). However, this test
statistic requires to invert an estimation $S$ of the full covariance matrix
$\Sigma_{i}+\Sigma_{j}$, which is usually intractable in our low data regime.

To overcome this, Srivastava & Du (2008) proposed a test statistic based only
on the diagonal $D_{S}$ of the covariance estimator $S$, and on the trace of
the square of its associated correlation matrix $R\equiv
D_{S}^{-1/2}SD_{S}^{-1/2}$:

$d_{\phi}^{2}(x_{i},x_{j})\equiv\alpha\Big{[}\big{(}\hat{\mu}_{i}-\hat{\mu}_{j}\big{)}^{T}D_{S}^{-1}\big{(}\hat{\mu}_{i}-\hat{\mu}_{j}\big{)}-\beta\Big{]},$
(2)

where $\hat{\mu}_{i}$ is the estimator of $\mu_{i}$ obtained through an
average over the $4$ sub-patches $x_{i}^{(l)}$, $\alpha$ is an overall factor
to normalize the variance of $d^{2}_{\phi}$, and $\beta$ is a debiasing term.
We emphasize that, through a dependence on $tr(R^{2})$, the $\alpha$ factor
accounts, at least partially, for the correlation structure of $\phi$. More
details about these terms are given in appendix C.

Under the assumption that $\Sigma_{i}=\Sigma_{j}$, the
$d_{\phi}^{2}(x_{i},x_{j})$ test statistic has a variance of order unity.
Thus, when it is much larger than 1, $\mu_{i}$ & $\mu_{j}$ cannot be
considered compatible:

$d_{\phi}^{2}(x_{i},x_{j})\gg 1\implies\mu_{i}\text{ \& }\mu_{j}\text{
incompatible},$ (3)

whereas when it is of the order of 1 or less, it is not possible to detect a
discrepancy between the two means $\mu_{i}$ and $\mu_{j}$ with the available
amount of data:

$\displaystyle d_{\phi}^{2}(x_{i},x_{j})\lessapprox 1\implies\mu_{i}\text{ \&
}\mu_{j}\text{ not incompatible based on}$ (4) the available amount of data.

###  3.3 Comparing summary statistics on a dataset

By extending the $\phi$-compatibility test introduced above, we set up a
comparison between two sets of summary statistics $\phi_{A}$ and $\phi_{B}$ on
a given dataset. This comparison consists in studying whether each set of
summary statistics has degeneracies that the other set can lift. These
degeneracies are evidenced by the presence in the dataset of pairs of maps
that are clearly incompatible for $\phi_{A}$, but for which $\phi_{B}$ detects
no incompatibility, or vice versa. To be relevant, this comparison requires
that the dataset contain maps that actually have different properties, which
will be the case for the data sets studied later.

In practice, we suggest the following algorithm:

  * •

Step 1: for every patch $x_{i}$ of the dataset, compute the statistics
$\phi_{A}$ and $\phi_{B}$ of its corresponding sub-patches
$\\{x_{i}^{(l)}\\}_{l}$.

  * •

Step 2: for every pair $\\{x_{i},x_{j}\\}$ with $i\neq j$, compute
$d_{\phi_{A}}^{2}(x_{i},x_{j})$ and $d_{\phi_{B}}^{2}(x_{i},x_{j})$ from the
quantities derived in step 1.

  * •

Step 3: place every pair as a point on a 2D scatter plot showing
$d_{\phi_{A}}^{2}(x_{i},x_{j})$ against $d_{\phi_{B}}^{2}(x_{i},x_{j})$, such
as Fig. 6.

The resulting scatter plot can then be used to detect whether there are some
pairs of processes in the dataset that are degenerate for $\phi_{A}$ but
lifted by $\phi_{B}$ (or vice versa), as explained in Fig. 6. In such plots,
the presence of points in the bottom right region, i.e., where
$d_{\phi_{A}}\gg 1$ and $d_{\phi_{B}}\lessapprox 1$, reveals that some pairs
of this dataset are identified by $\phi_{A}$ as incompatible but not by
$\phi_{B}$: such pairs thus evidence degeneracies of $\phi_{B}$ lifted by
$\phi_{A}$ (panel a). Conversely, the presence of points in the top left
region evidences degeneracies of $\phi_{A}$ lifted by $\phi_{B}$. Hence, if
all the points land mainly in the sub-diagonal part (panel b), this evidences
that $\phi_{A}$ is better suited than $\phi_{B}$ to compare the pairs of this
dataset. If on the contrary the points are spread both in the upper left and
bottom right regions of the plot (panel c), this shows that both $\phi_{A}$
and $\phi_{B}$ are individually not sufficient to describe the processes of
this dataset.

##  4 Summary statistics

We present below the different summary statistics that are used in this paper.
In the following, $x$ is an image, $x(\@vec{u})$ is the value of the image at
pixel $\@vec{u}$, $\tilde{x}(\@vec{k})$ is the discrete Fourier transform of
image $x$ evaluated at wave-vector $\@vec{k}$, $\star$ stands for the
convolution operator, and $\langle\rangle_{\@vec{u}}$ for the averaging over
pixels. When $\phi(x)$ is multivariate, $\phi(x)[i]$ is the value of its
$i$-th dimension. $\bar{x}$ designates the following normalization of $x$:
$\bar{x}\equiv x/\text{std}(x)$ where
$\text{std}(x)\equiv\sqrt{\big{\langle}\big{[}x-\langle
x\rangle_{\@vec{u}}\big{]}^{2}\big{\rangle}_{\@vec{u}}}$.

Note that for some usual statistics, we consider their logarithms instead of
their raw values. For instance we use $\log\,\langle x\rangle_{\@vec{u}}$
instead of $\langle x\rangle_{\@vec{u}}$. This is possible when we work with
positive-valued statistics. We made this choice to better fit the Gaussianity
assumption given Eq. (1). In addition, we found that logarithms expressed more
discriminatory power with our diagnostic. See appendix D for more details. If
no precision is made, $\log$ designates the logarithm in base 10.

###  4.1 One-point based statistics

We list below the one-point statistics used in this paper. Even though some of
these statistics can probe non Gaussian information such as sparsity, they are
all point wise statistics. This means in particular that they cannot capture
spatial arrangement in the maps.

  * •

The mean:

$\boxed{\phi_{\text{mean}}(x)\equiv\log\,\langle<EMAIL_ADDRESS>(5)

  * •

The variance:

$\boxed{\phi_{\text{var}}(x)\equiv\log\,\big{\langle}\big{[}x-\langle
x\rangle_{\<EMAIL_ADDRESS>(6)

  * •

The mean of the logarithm:

$\boxed{\phi_{\text{mean of log}}(x)\equiv\langle\log<EMAIL_ADDRESS>(7)

  * •

The variance of the logarithm:

$\boxed{\phi_{\text{var of log}}(x)\equiv\log\,\big{\langle}\big{[}\log
x-\langle\log x\rangle_{\<EMAIL_ADDRESS>(8)

  * •

The quantile functions666We consider quantile functions rather than
probability distribution functions, as the difficulty of defining a unique
binning range makes comparisons less efficient. Indeed, to define a range in
terms of quantiles allows for a natural adaptation to each MCs’ regions. (QF)
normalized by the median:

$\boxed{\phi_{\text{QF}}(x)[i]\equiv\log\,\big{[}q_{\alpha_{i}}(x)/q_{1/2}(x)\big{]},}$
(9)

where $q_{\alpha}(x)$ designates the $\alpha$-quantile of the distribution of
values $\\{x(\@vec{u})\\}_{\@vec{u}}$. Hence $q_{1/2}$ stands for the median
operator. We consider 10 quantiles $\\{\alpha_{i}\\}_{i}$ such that
$1-\alpha_{i}$ are logarithmically spaced between $10^{-4}$ and $0.4$. Low
quantiles are not considered as, for observational data, they are contaminated
(Ossenkopf-Okada et al., 2016). This logarithmic binning of the high column
density values is motivated by the log-normal to power-law behavior of the
tail of MCs’ PDFs. See for instance Pouteau et al. (2023) and references
therein.

###  4.2 Two-point based statistics

A very popular way to describe spatial properties of a process is through the
power spectrum (PS). It is defined as the Fourier transform of the auto-
correlation function, and describes the energy distribution scale by scale in
the process studied. In this paper, we consider the isotropic PS, defined as:

$PS[k]\equiv\langle\mathbb{E}[|\tilde{x}(\<EMAIL_ADDRESS>
A power-law behavior $PS[k]\propto k^{-\beta}$ is expected for fields arising
in turbulent MHD (Schekochihin, 2022). Thus we adopt a log-log representation,
to linearize it and for stability purposes (Bruna & Mallat, 2013):

$\boxed{\phi_{\text{PS}}(x)[i]\equiv\log\,\big{\langle}|x_{\text{apo}}\star\psi_{i}|^{2}\big{\rangle}_{\@vec{u}},}$
(10)

where $x_{\text{apo}}$ is apodized to mitigate non periodic boundary
conditions (PBC) as explained in appendix B. We use 6 band-pass filters
$\\{\psi_{i}\\}_{1\leq i\leq 6}$ defined as:

$\tilde{\psi}_{i}[\@vec{k}]=\mathrm{1}_{\|\@vec{k}\|\in[k_{i},k_{i+1}[},$

where $k_{i}$ is logarithmically spaced between
$k_{\min}=1/256\,\text{pix}^{-1}$ and $k_{\max}=1/4\,\text{pix}^{-1}$. For the
finest maps with pixel size $6\arcsec$, this $k_{\max}$ corresponds to a
smallest angular scale of $24\arcsec$, which remains above the $18\arcsec$
resolution of the observations.

Because some MCs tend to have log-normal behavior in their one-point
statistics, we also consider PS statistics of the logarithm of the maps:

$\boxed{\phi_{\text{PS of log}}(x)[i]\equiv\log\,\big{\langle}|(\log
<EMAIL_ADDRESS>(11)

###  4.3 Scattering statistics

The Wavelet Scattering Transform (WST) is a set of non-Gaussian descriptors
with a hierarchical and multi-scale structure (Bruna & Mallat, 2013). It has
been shown to be very efficient at describing astrophysical fields (see, for
instance, Allys et al., 2019; Saydjari et al., 2021). The usual WST consists
in two layers of statistics, the first of which depends on a single scale of
length $\simeq 2^{j}$ pixels with orientation $\theta$:

<EMAIL_ADDRESS>
The second layer probes a coupling between two oriented scales
$(j_{1},\theta_{1})$, and $(j_{2},\theta_{2})$, with $j_{1}<j_{2}$:

$S_{2}(x)[j_{1},j_{2},\theta_{1},\theta_{2}]\equiv\big{\langle}\big{|}|x\star\psi_{j_{1},\theta_{1}}|\star\psi_{j_{2},\theta_{2}}\big{|}\big{\rangle}_{\@vec{u}}\,/\,S_{1}(x)[j_{1},\theta_{1}].$

Then, observing that many processes of interest in astrophysics and cosmology
exhibit strong regularities in their angular dependencies, Allys et al. (2019)
proposed the Reduced WST (RWST): an angular compression of WST statistics for
2D data. We are going to use three of the main descriptors they introduced (we
refer to the above reference for more details):

  * •

$S_{1}^{Iso}$, a scale-by-scale isotropic descriptor, that we normalize by
estimating it on $\bar{x}$ instead of $x$:

<EMAIL_ADDRESS>(12)

Note that it computes a $L^{1}$ norm of the filtered field, in contrast with
the power spectrum that probes a $L^{2}$ norm.

  * •

$S_{2}^{Iso1}$ and $S_{2}^{Iso2}$, that measure an isotropic coupling between
scales and are defined by fitting the following model:

$\log_{2}\,S_{2}(x)[j_{1},j_{2},\theta_{1},\theta_{2}]\simeq
S_{2}^{Iso1}(x)[j_{1},j_{2}]\,+\,\\\
S_{2}^{Iso2}(x)[j_{1},j_{2}]\cos\big{[}2(\theta_{2}-\theta_{1})\big{]}.$

They characterize, respectively, a coupling between non-oriented scales and
the relative coupling between parallel and perpendicular scales.

In this paper, we will furthermore compress both $S_{2}^{Iso1}$ and
$S_{2}^{Iso2}$ coefficients that depend on two scales $(2^{j_{1}},2^{j_{2}})$
by only keeping their $j_{2}-j_{1}$ dependency, i.e., a dependency on their
ratio, by considering the following average:

$\boxed{\langle S_{2}^{Iso1}\rangle_{j_{2}-j_{1}}(x)[\delta]\equiv\langle
S_{2}^{Iso1}(x)[j_{1},j_{2}]\rangle_{j_{2}-j_{1}=\delta},}$ (13)

and

$\boxed{\langle S_{2}^{Iso2}\rangle_{j_{2}-j_{1}}(x)[\delta]\equiv\langle
S_{2}^{Iso2}(x)[j_{1},j_{2}]\rangle_{j_{2}-j_{1}=\delta}.}$ (14)

In practice, we consider 4 scales between $j_{\min}=2$ and $j_{\max}=5$, which
leads to the 3 possible values $1\leq\delta\leq 3$ for our reduced $S_{2}$
coefficients. We also divide the half-plane $[0,\pi[$ in 8 different $\theta$
angles. For the finest maps of pixel size $6\arcsec$, this $j_{\min}$
corresponds to a scale of $24\arcsec$ that remains above the $18\arcsec$
resolution of the observations.

For RWST computations, convolutions are done with local wavelets. This allows
us to mitigate non PBC after convolution by suitably cropping the corrupted
edges, without using apodization. These computations are done with the
pywst777https://github.com/bregaldo/pywst Python package (Regaldo-Saint
Blancard et al., 2020).

$\phi$ | Dim. | Non- | Scale | Scale
---|---|---|---|---
|  | Gaussian | description | interactions
mean | 1 | ✕ | ✕ | ✕
var | 1 | ✕ | ✕ | ✕
mean of log | 1 | log-Gauss. | ✕ | ✕
var of log | 1 | log-Gauss. | ✕ | ✕
QF | 10 |  | ✕ | ✕
PS | 6 | ✕ |  | ✕
PS of log | 6 | log-Gauss. |  | ✕
$S_{1}^{Iso}$ | 4 |  |  | ✕
$\langle S_{2}^{Iso1}\rangle_{j_{2}-j_{1}}$ | 3 |  | ✕ |
$\langle S_{2}^{Iso2}\rangle_{j_{2}-j_{1}}$ | 3 |  | ✕ |

Table 2: Overview of the sets of summary statistics, their dimensions and
properties (whether they are sensitive to non-Gaussianity, provide a
multiscale description, and probe couplings between scales).

###  4.4 Overview

The statistics introduced above are compiled in Table 2. In the following, in
addition to the sets of summary statistics defined above, we will also
consider some aggregated sets of summary statistics that are made from groups
of these building blocks. For instance, we will refer to
$\phi_{\text{Gaussian}}$ statistics to mean the joint set of
$\\{\phi_{\text{mean}},\phi_{\text{var}},\phi_{\text{PS}}\\}$. These groups
are defined in Table 3.

Set of stats. |  Composed of | Dim.
---|---|---
Gaussian |  mean + PS | 7
log-Gaussian |  mean of log + PS of log | 7
RWST |  $S_{1}^{Iso}$ \+ $\langle S_{2}^{Iso1}\rangle_{j_{2}-j_{1}}$ \+ $\langle S_{2}^{Iso2}\rangle_{j_{2}-j_{1}}$ | 10
final |  log-Gaussian, RWST | 17

Table 3: Aggregated sets of summary statistics used in this work.

##  5 Towards a low-degeneracy set of statistics

We want to construct a set of informative summary statistics for the
observational dataset. To do so, we will follow a bottom-up approach, starting
with usual low-order statistics, and increasingly trying to improve on them,
by exhibiting and lifting potential degeneracies. The underlying idea is that
low-order statistics are good candidates to concentrate most of the
informative power into a few coefficients. Such a concentration property is of
interest, especially as our low data regime prevents us from learning many
features. Another benefit of this approach is to first promote simple
statistics, and use more elaborate ones only if they add a significant
contribution.

###  5.1 Molecular clouds have Gaussian degeneracies

To begin with, we investigate whether Gaussian statistics (i.e. the set made
of mean and PS statistics) are degenerate for the observational data. To point
out such potential degeneracies, we confront these Gaussian statistics,
starting with a low-order non-Gaussian set of statistics, the quantile
function QF, using the methodology introduced in Sec. 3.3. As shown in Fig. 7,
this first confrontation on observational data already underlines a strong
degeneracy level for both Gaussian and QF features, showing that neither
statistics is by itself sufficient, according to the compatibility diagnostic
we introduced. In particular, the quantile statistics extract a significant
amount of information that cannot be efficiently captured by Gaussian
statistics.

Figure 7: Confronting Gaussian statistics with QF statistics on observational
data, based on the test presented in Fig. 6. Each set of statistics has strong
degeneracies lifted by the other set. To investigate the Gaussian confusions,
we pick six pairs of 512 $\times$ 512 patches, whose locations on the scatter
plot are given by the red stars. These pairs are shown in Fig. 8. Figure 8:
Examples of Gaussian confusions. Six pairs of 512 $\times$ 512 patches are
chosen, whose locations on the scatter plot of Fig. 7 are given there by the
red stars. The column density maps are shown in units of
$10^{20}\text{cm}^{-2}$. For each patch, we report: $s$ the sub-sampling
factor from the original 3″/pix map, $d$ and $b$ the approximated distance and
Galactic latitude of the cloud. The pixel size (in mpc) of a patch is thus
proportional to $s\times d$. If a pair has patches $(i,j)$ with incompatible
pixel sizes, that we define according to the following criterion
$\max\\{\frac{s_{i}d_{i}}{s_{j}d_{j}},\frac{s_{j}d_{j}}{s_{i}d_{i}}\\}\geq
3/2$, we color the labels in red. Figure 9: Statistics for the examples of
Gaussian confusions shown in Fig. 8. In each row, the orange filled line
(resp. band) corresponds to the mean (resp. std) of the statistics computed
over the four 256 $\times$ 256 sub-patches of the top patch of each pair of
Fig. 8, and the corresponding blue lines and areas refer to the bottom patch
of the pair. The top row corresponds to the statistics used in the $x$-axis of
the scatter plot of Fig. 7, i.e., QF statistics, plotted with 10 increasing
quantile values, while the bottom row corresponds to the $y$-axis, i.e.,
Gaussian statistics, plotted starting with mean and followed by the binned PS
with six decreasing scales. To better highlight the discrepancies between the
two patches of a given pair, we report in second and third rows the offsets of
the orange and blue filled lines with respect to their common mean.

In addition to this dataset-wide diagnostic, we display in Fig. 8 six randomly
picked pairs of 512 $\times$ 512 patches degenerate for Gaussian statistics
but with increasing QF incompatibility. The locations of these pairs on the
scatter plot of Fig. 7 are given there by the red stars. For these six pairs
of patches, we show in Fig. 9 the two sets of statistics used in the
diagnostic (QF statistics and Gaussian statistics), where the orange filled
lines (resp. bands) correspond to the mean (resp. standard deviation) of the
statistics computed over the four 256 $\times$ 256 sub-patches of the top
patch of each pair of Fig. 8, and the corresponding blue lines and areas refer
to the bottom patch of the pair. To better highlight the discrepancies between
the two patches of a given pair, the offsets of the orange and blue filled
lines with respect to their common mean are also shown in Fig. 9. These
results illustrate the ability of QF statistics to distinguish between images
with the same Gaussian statistics, which confirms that our compatibility
diagnostic works as expected. We emphasize here that the $\phi_{\text{QF}}$
statistics are complementary to the one-point properties probed by Gaussian
statistics.

This behavior is not particularly surprising as one-point properties of column
density maps of MCs are expected to be at least log-normal, if not with a
power-law tail. Hence, estimating such properties directly from the logarithm
of the maps, as does $\phi_{\text{QF}}$ but not $\phi_{\text{Gaussian}}$,
might enhance the discriminative power. To test this idea, we study in the
first column of Fig. 10 the degeneracy level of the set of mean and variance
statistics estimated respectively on the raw maps and their logarithm. More
precisely, top row plots confront QF statistics to
$\\{\phi_{\text{mean}},\,\phi_{\text{var}}\\}$, while bottom row plots
confront QF statistics to $\\{\phi_{\text{mean of log}},\,\phi_{\text{var of
log}}\\}$. We see in this figure that the pairs of observations are clearly
better discriminated when using the mean and variance computed on the
logarithm of the maps (d) rather than directly (a). This effect is very well
reproduced by logFBM data (middle column), but does not hold in general, as
shown with the Describable Textures Dataset (DTD) (right column).

These results suggest that the specific Gaussian degeneracies evidenced in
Fig. 7 are mainly explained by the inefficiency of Gaussian one-point
statistics $\\{\phi_{\text{mean}},\,\phi_{\text{var}}\\}$ to characterize one-
point properties of processes that have a log-normal (or heavier tail) nature.

Surprisingly, this analysis also shows that, for datasets such as observations
or logFBMs, probing the PDF properties through the prism of a mean and a
variance is more discriminative than probing its shape. Here, in both d) and
e), our suitably constructed set $\\{\phi_{\text{mean of
log}},\,\phi_{\text{var of log}}\\}$, of dimension 2, performs almost always a
better discrimination than the set $\phi_{\text{QF}}$, of dimension 10. This
emphasizes the importance of suitably constructed low dimensional descriptions
in such analysis. The performative power of such low-dimensional features is
usually due to the exploitation of underlying regularities of the studied
processes that allow these compression, without much loss of information. For
deeper insights on the connection between regularity, approximation and
sparsity, we recommend Pr. Mallat’s lecture888https://www.college-de-
france.fr/fr/agenda/cours/representations-parcimonieuses/le-triangle-
regularite-approximation-parcimonie. However, when dealing with datasets made
of a wide diversity of irregular processes, such as everyday life textures
like DTD, it is hard to identify a compression that does not lead, for some
pairs, to poorer performances (c, f).

Figure 10: Top row plots confront QF statistics to
$\\{\phi_{\text{mean}},\,\phi_{\text{var}}\\}$, while bottom row plots
confront QF statistics to $\\{\phi_{\text{mean of log}},\,\phi_{\text{var of
log}}\\}$. These results evidence that taking the logarithm of the map
enhances the discriminative power of the mean and variance statistics on both
observation and logFBM data (left and middle columns) but not on DTD (right
column).

###  5.2 Molecular clouds have log-Gaussian degeneracies

The ability of Gaussian one-point statistics to grasp efficiently one-point
properties of observations from the logarithm of the column density maps,
shown in Fig.10, suggests to shift towards log-Gaussian statistics, i.e., mean
and PS estimated on the logarithms of the maps. We thus investigate now
whether we may point out some degeneracies of these statistics on
observational data. However, if we confront this set with QF statistics, as we
did previously for Gaussian statistics, we do not expect to lift significant
log-Gaussian degeneracies. Instead, we suggest to search for such degeneracies
using a higher order set: the RWST. To make sure that, on other datasets, the
potential degeneracies lifted by these higher order statistics could not be
lifted by one-point statistics, we will confront the log-Gaussian set to the
RWST one.

We show in Fig. 11 that log-Gaussian statistics have degeneracies on
observational data lifted by RWST. We show six examples of degenerate pairs in
Fig. 12 and the statistics of these patches in Fig. 13. In most of these
pairs, the $S_{1}^{Iso}$ coefficients heavily contribute to lift the log-
Gaussian degeneracy. In the first and last pairs, $\langle
S_{2}^{Iso2}\rangle_{j_{2}-j_{1}}[\delta]$ and respectively $\langle
S_{2}^{Iso1}\rangle_{j_{2}-j_{1}}[\delta]$ are also very discriminative.

When applying the same analysis on logFBM data (Fig. 14.a), this finding
however does not hold. The RWST diagnostic brings no additional information
that was not already probed by log-Gaussian statistics. Indeed, as expected
for such a dataset, log-Gaussian statistics are sufficient (Cover, 1999).
Finally, simulations yield yet another result (Fig. 14.b). Although this
latter case leads to the same qualitative conclusions as observations, namely
underlining the insufficiency of log-Gaussian statistics, it differs from it
quantitatively. This supports the caveat of simulation-based inference
according to which observations and simulations of such processes lie in sub-
manifolds having different geometries, and motivates the observation-based
approach of this paper.

###  5.3 Final set of statistics

In this paper, we choose to limit our study to the log-Gaussian + RWST set,
which we call $\phi_{\text{final}}$, that seems well suited to describe this
observational dataset. To ensure this suitability any further, we could
continue to study the degeneracies of $\phi_{\text{final}}$ by confronting it
with other sets of statistics. In any case, it is difficult to guarantee with
certainty at any stage that we have not missed any remaining degeneracy of the
resulting set that would be obtained by such a construction, even if each
confrontation with a complementary set reinforces our confidence in its
completeness.

Once a set of statistics has been fixed, however, it is also possible to use a
complementary diagnostic to study it, which consists in checking the visual
proximity of maps identified as close. If they are clearly visually different,
this indicates complementary limitations of the set of statistics used. This
is done with $\phi_{\text{final}}$ in the following section, where we
introduce a distance to evaluate a notion of proximity rather than
compatibility between pairs of maps.

Finally, we want to emphasize that low-order statistics generally constitute a
precious and easily-accessed source of information. They should not be
underestimated, and attempting to tailor them to the type of data considered
can be fruitful, as shown here when shifting from a Gaussian to a log-Gaussian
description.

Figure 11: Same as Fig. 7, but using the log-Gaussian and RWST sets of
statistics on the observational dataset. To investigate the degeneracies of
the log-Gaussian statistics, six pairs of patches are selected, corresponding
to the six red stars on this plot. Figure 12: Examples of log-Gaussian
degeneracies. Six pairs of 512 $\times$ 512 patches are chosen, whose
locations on the scatter plot of Fig. 11 are given there by the red stars. The
column density maps are shown in units of $10^{20}\text{cm}^{-2}$. For each
patch, we report: $s$ the sub-sampling factor from the original 3″/pix map,
$d$ and $b$ the approximated distance and Galactic latitude of the cloud. If a
pair has patches $(i,j)$ with incompatible pixel sizes according to the
following criterion
$\max\\{\frac{s_{i}d_{i}}{s_{j}d_{j}},\frac{s_{j}d_{j}}{s_{i}d_{i}}\\}\geq
3/2$, we color the labels in red. Figure 13: Statistics for the examples of
log-Gaussian confusions shown in Fig. 12. In each row, the orange filled line
(resp. band) corresponds to the mean (resp. std) of the statistics computed
over the four 256 $\times$ 256 sub-patches of the top patch of each pair of
Fig. 12, and the corresponding blue lines and areas refer to the bottom patch
of the pair. The top row corresponds to the RWST, starting with
$S_{1}^{Iso}[j]$ coefficients with four increasing scales $j$, then $\langle
S_{2}^{Iso1}\rangle_{j_{2}-j_{1}}[\delta]$ coefficients with three increasing
scale ratios $\delta$ and finally $\langle
S_{2}^{Iso2}\rangle_{j_{2}-j_{1}}[\delta]$ coefficients with the same three
scale ratios. The bottom row represents in the following order: mean of log
followed by PS of log with 6 decreasing scales. The second and third rows show
the offsets of these statistics with respect to the mean of the two.

Figure 14: log-Gaussian statistics have degeneracies lifted by RWST for
simulations in (b), but not for logFBM processes in a). This is expected
because log-Gaussian statistics are sufficient for families made of logFBM
processes (Cover, 1999).

##  6 Comparing pairs and datasets

The previous section confronted multiple sets of statistics to assess their
information content through their level of degeneracy. In this section, we fix
the set of statistics: $\phi=\phi_{\text{final}}$, composed of 17 coefficients
(7 log-Gaussian descriptors and 10 RWST statistics), and use it to define a
distance between maps. We illustrate this distance by evidencing closest pairs
of images in a dataset, as well as between different datasets such as
observations and simulations.

###  6.1 Defining a morphological distance

Our objective is to define a distance between two patches based on
$\phi_{\text{final}}$. One of our requirements is to enable a comparison
between distance values for different pairs, so that these pairs can be
ordered according to the morphological proximity of their patches. This
requirement prevents us from using the statistical compatibility diagnostic
introduced earlier. Indeed, since it is weighted by the local variance of each
patch’s statistics, this can lead, for example, to some patches in a dataset
being compatible with almost all the others simply because their spatial
variance is very high. This could also encourage a simulation to approach an
observation, following this criterion, by arbitrarily increasing its variance
$\operatorname{Var}\phi(x_{SIM})$ without focusing on reducing the discrepancy
of its average statistical properties
$\hat{\phi}(x_{SIM})-\hat{\phi}(x_{OBS})$. This property was purposely used to
act as a penalization when confronting different sets of statistics in the
previous section, but is no longer desired for pairs’ ordering, once $\phi$ is
fixed.

To build a distance that avoids this drawback, we choose instead to normalize
it by the variability of the statistics evaluated over the entire
observational dataset. We thus use the following distance:

$d^{2}_{\mathcal{D}}(x_{i},x_{j})\equiv(\hat{\mu}_{i}-\hat{\mu}_{j})^{T}(\operatorname{diag}M_{\mathcal{D}})^{-1}(\hat{\mu}_{i}-\hat{\mu}_{j}),$
(15)

which is normalized by the spanning of the estimated values $\hat{\mu}_{i}$ of
$\phi_{\text{final}}$ over all maps of a given dataset $\mathcal{D}$:

$M_{\mathcal{D}}\equiv\big{\langle}(\hat{\mu}_{i}-\langle\hat{\mu}_{j}\rangle_{j})(\hat{\mu}_{i}-\langle\hat{\mu}_{j}\rangle_{j})^{T}\big{\rangle}_{i},$
(16)

where the brackets indicate an average over $\mathcal{D}$. Note, however, that
it is difficult to interpret the value of $d^{2}_{\mathcal{D}}$ in absolute
terms. Indeed, the $M_{\mathcal{D}}$ term does not describe a typical variance
for a given process, but describes the variety of morphologies encountered in
the entire dataset, that can be wide, as for the MC data investigated here.
Unlike statistical compatibility diagnostics, the $d^{2}_{\mathcal{D}}$
distance is therefore modified by the addition or removal of maps in the
$\mathcal{D}$ dataset, and can be affected by the presence of outliers.

In the following, we work with different datasets. For instance, we aim at
comparing the minimum distance between observations and simulations to the
typical distance of the closest pairs of observations. To do so, we use in
this paper a common metric $M_{\text{obs}}$ computed on the observation
dataset, as well as the associated $d^{2}_{\text{obs}}$ distance.

###  6.2 Closest pairs

We use the $d^{2}_{\text{obs}}$ distance to identify the closest pairs of
patches that can be found between two datasets. To do so, we report in Fig. 15
the cumulative distributions of distances associated to pairs of observation
patches (in gray), as well as to pairs made of one observation patch and one
logFBM, DTD, and simulation patch (in red, purple, and blue, respectively).
For each case, we pick up six of the typical closest pairs999We do not
investigate the $\sim 1\%$ closest pairs of observations, since the comparison
between observations is slightly biased with respect to the other ones. The
reason is that, in that case only, we have pairs made of non independent
patches, such as neighboring patches in the sky or two slightly different
scaled versions of the same region. This corresponds to approximately $1\%$ of
the pairs, that we therefore exclude., as well as six relatively distant
pairs, that we show in Fig. 16. In dashed lines with the same colors, we also
report the cumulative distributions of distances associated to pairs of
patches from a common dataset (logFBM, DTD101010Note that the $\sim 10$
closest pairs of DTD images, located at $d^{2}_{\text{obs}}\lessapprox 1$ are
artifacts in this dataset corresponding in practice to almost identical
images., and simulation patches, respectively).

Figure 15: Cumulative distributions of $d^{2}_{\text{obs}}$ distances between
pairs extracted from multiple datasets. The different curves correspond to
different choices of datasets from which the two patches of a given pair are
extracted. This distance is based on $\phi_{\text{final}}$, a set of 17
coefficients (7 log-Gaussian descriptors and 10 RWST statistics). The same
metric $M_{\text{obs}}$ is used for all distances and is defined based on the
observational dataset. The stars in the bottom sub-panel correspond to ”close”
and ”distant” pairs shown in Fig. 16. An horizontal line at the value $1\%$ is
drawn to represent the proportion of pairs of observations that are
neighboring patches in the sky. For these pairs, the distance is
underestimated and thus the gray curve below this line is not meaningful. Same
applies for pairs of simulations (blue dashed curve). Figure 16: Closest
(left) and more distant (right) pairs extracted from distributions of
distances reported in Fig. 15. We see that the closest (OBS, SIM) pairs are
much more distant ($d^{2}_{\text{obs}}\sim 7$) than the closests (OBS, OBS)
pairs ($d^{2}_{\text{obs}}\sim.7$). See Sec. 6.3 for more detailed
interpretation on these results. Some logFBM models end up quite close to the
most diffuse regions observed: Polaris Flare and Corona Australis. Such
regions are highly contaminated by CIB emission. Many (OBS, DTD) pairs are
found to be close whereas the textures look very different. This shows that
the set of summary statistics developed in this paper is tuned for ISM
observations but is far from being sufficient for any kind of data. This also
illustrates that MCs have much more regularity in terms of morphology than DTD
textures.

Closest (OBS, OBS) pairs and (OBS, logFBM) pairs are visually rather similar,
while distant pairs look very different. This relative agreement between
morphological proximity as probed by our statistical distance as well as by
human vision is encouraging, because it should be satisfied by an ideal
distance. However, it remains far from being an exhaustive diagnostic.

The closest (OBS, logFBM) pairs have a distance $d^{2}_{\text{obs}}\sim 2$,
which marks a high agreement. Indeed, only $\sim 0.5\%$ of (OBS, OBS) or
(logFBM, logFBM) pairs exhibit smaller distance values, and we have already
mentioned that the closest 1% of (OBS, OBS) pairs are very similar by
construction. Comparatively, simulations are more distant to observations: the
closest pairs between observations and simulations have a distance
$d^{2}_{\text{obs}}\sim 8$, four times larger than the closest (OBS, logFBM)
pairs, when already $\sim 10\%$ of pairs of observations exhibit smaller
distances (see Sec. 6.3 for more details on this interpretation). As expected,
these closest (OBS, SIM) pairs are visually less similar than the closest
(OBS, OBS) or (OBS, logFBM) pairs (left panels of Fig. 16), but are still not
so different, for instance with respect to the more distant pairs, that can be
seen in the right panels of Fig. 16.

The closest (OBS, logFBM) samples obtained here are very diffuse regions of
MCs such as Corona Australis, Polaris Flare, and Ophiuchus. This is not
surprising because for MCs, diffuse regions are closer to logFBM models than
dense regions, whose PDFs are known to deviate from log-normality. However,
these diffuse regions are still supposed to exhibit coherent structures that
should induce deviations from logFBM models. Here, such deviations are found
to be small. We note, however, that this proximity only means that certain
diffuse regions are close to logFBM processes relative to the total
variability of the observational dataset, and not necessarily that they are
well described by such models, or could not be distinguished by the previous
compatibility diagnostic. Moreover, this proximity to logFBM models is also
partly due to contamination by the Cosmic Infrared Background (CIB), that has
non negligible power at such low levels of column density, and is expected to
Gaussianize the data, including their RWST statistics (Auclair et al., 2024).
Incidentally, we remind that the observational dataset has some artifacts that
make it deviate from an ”ideal” MC column density dataset, even if we try to
limit as much as possible their effects (noise, finite resolution), as
discussed in Sec. 4.

On the contrary, the observations that are closest to the simulations of dense
MCs correspond, unsurprisingly, to dense regions of MCs such as Ophiuchus,
Aquila, and Serpens. Note that Ophiuchus exhibits patches that are close to
these dense simulations, but also at least one patch that is close to a logFBM
model, underlining the spatial variability of molecular cloud morphologies.
This illustrates the difficulty to treat a MC as a single entity, and
emphasizes the relevance of our local approach.

The comparison between observations and DTD shows the limitations of our
distance diagnostic. Indeed, the closest (OBS, DTD) pairs are found at a
distance $d^{2}_{\text{obs}}\sim 2$ that is the typical distance between close
(OBS, OBS) pairs or close (OBS, logFBM) pairs, although they are visually very
different. This illustrates that, for a diagnostic based on a low-dimensional
set of statistics, it is difficult to probe a distance over a family of
processes that has such a wide variety of textures as DTD. On the contrary,
because the simulations that are closest to observations are more distant
($d^{2}_{\text{obs}}\sim 8$), this suggests that the set of simulations does
not intersect the set of observations and that $\phi_{\text{final}}$ is able
to pinpoint this discrepancy, as discussed in the following subsection.

Finally, it is quite impressive to see that we can build a distance diagnostic
from a representation of dimension 17 only that still manages to identify
morphological similarity between maps quite satisfactorily. This illustrates
the possibility of constructing a highly informative but low-dimensional
description tailored to a family of processes from an ensemble of compressed
sets of usual statistics. It should be stressed, however, that this study
remains partial, notably because it is based mainly on the observation of a
few close pairs in our dataset. Yet, the confrontation diagnostics studied in
Sec. 5 showed that degenerate counterexamples remain largely in the minority.
In addition, we lack solid baselines since it is inherently difficult to
quantify visual impressions of morphological proximity, although some work has
been done in this direction (Peek & White, 2021).

###  6.3 Interpreting the minimal distance between observations and
simulations

A last question we tackle is whether the relatively high value of the minimal
distance $d^{2}_{\text{obs}}$ that we get between observations and simulations
indeed indicates distinct statistical properties of these two sets. More
precisely, we question whether the set of observations’ moments
$\\{\bar{\mu}_{i}\\}_{i\in\text{obs}}$ overlaps that of simulations
$\\{\bar{\mu}_{j}\\}_{j\in\text{sim}}$, where
$\bar{\mu}_{i}\equiv\mathbb{E}[\phi_{\text{final}}(x_{i})]$ stands for the
expected value over a given process $i$. In practice, we aim at retrieving the
minimal distance

$\bar{d}^{2}_{ij}\equiv(\bar{\mu}_{i}-\bar{\mu}_{j})^{T}(\operatorname{diag}M_{\text{obs}})^{-1}(\bar{\mu}_{i}-\bar{\mu}_{j})$
(17)

between $\bar{\mu}_{i}$ and $\bar{\mu}_{j}$ over (OBS, SIM) pairs $(i,j)$.
However, the distance $d^{2}_{\text{obs}}$ introduced previously is a
statistical estimator, so that in general

$d^{2}_{ij}\neq\bar{d}^{2}_{ij},$

and in particular: $\mathbb{E}[d^{2}_{ij}]>\bar{d}^{2}_{ij}$. Indeed, the
variance of the $\hat{\mu}_{i}$ estimator biases $d^{2}_{ij}$ with respect to
$\bar{d}^{2}_{ij}$:

$\mathbb{E}[d^{2}_{ij}]\equiv\bar{d}^{2}_{ij}+b_{ij}.$ (18)

This non-negative bias, which boils down to:

$b_{ij}=\mathrm{tr}\\{\operatorname{cov}\big{[}(\operatorname{diag}M_{\text{obs}})^{-1/2}(\hat{\mu}_{i}-\hat{\mu}_{j})\big{]}\\},$
(19)

prevents us from interpreting directly the value of $d^{2}_{ij}$ as
$\bar{d}^{2}_{ij}$. Furthermore, its value increases with the amplitude of the
fluctuations of $\hat{\mu}_{i}-\hat{\mu}_{j}$, and can thus change depending
on the pair $(i,j)$ considered. This dependency on the pair also prevents us,
without further check, from actually comparing $\bar{d}^{2}$ between different
pairs based on the estimations $d^{2}$.

To probe and interpret the minimal distance between observations and
simulations, we propose the following strategy:

  * •

find an observation $i$ and a simulation $j$ that minimize the biased
estimation $d^{2}_{ij}$. Such pairs are already reported in Fig.16 and are
good candidates to minimize $\bar{d}^{2}_{ij}$. The goal then becomes to
compare $\bar{d}^{2}_{ij}$ to the minimal value $\bar{d}^{2}_{ii^{\prime}}$,
for $i^{\prime}$ another observation, independent from $i$. But $\bar{d}^{2}$
is unknown. Instead:

  * •

find an observation $i^{\prime}$, that minimizes the biased distance
$d^{2}_{ii^{\prime}}$, while checking it does not overlap $i$ in the sky to
assume them independent.

  * •

Since
$\bar{d}^{2}_{ij}-\bar{d}^{2}_{ii^{\prime}}=\mathbb{E}[d^{2}_{ij}]-\mathbb{E}[d^{2}_{ii^{\prime}}]+b_{ii^{\prime}}-b_{ij}$,
the discrepancy between $\bar{d}^{2}_{ij}-\bar{d}^{2}_{ii^{\prime}}$ can be
estimated based on the measured discrepancy $d^{2}_{ij}-d^{2}_{ii^{\prime}}$
up to the unknown bias shift $b_{ii^{\prime}}-b_{ij}$.

  * •

If, in addition the observation patch $x_{i^{\prime}}$ is such that spatial
fluctuations of $\phi_{\text{final}}$ over its sub-patches are at least of the
order of those estimated from the sub-patches of the simulation patch $x_{j}$,
then based on Eq. 19, we have $b_{ii^{\prime}}\gtrapprox b_{ij}$. This allows
then to use $d^{2}_{ij}-d^{2}_{ii^{\prime}}$ as an estimated lower bound for
$\bar{d}^{2}_{ij}-\bar{d}^{2}_{ii^{\prime}}$.

We report in Fig. 17 such independent111111Indeed, only one of these pairs is
made of patches retrieved from the same MC (Ophiuchus), but these patches do
not overlap. pairs $(i,i^{\prime})$. In the first row, we fix $i$ associated
to the patch of Ophiuchus $x_{i}$ that is the closest observation to
simulations ($d^{2}_{ij}=7$, Fig.16). For that choice of patch, the reported
distances $d^{2}_{ii^{\prime}}$ of neighboring observations $i^{\prime}$ are
rather similar to $d^{2}_{ij}$. However, in the second row, where $i$ is
associated to the patch of Aquila that is the second closest observation to
simulations (as shown in Fig.16), we see other observations $i^{\prime}$, such
as Serpens or Orion B, that are closer than the closest simulation:
$d^{2}_{ii^{\prime}}=3<d^{2}_{ij}=8$. We checked that the fluctuations of
$i^{\prime}$ are comparable to the one of $j$ so that
$b_{ii^{\prime}}\gtrapprox b_{ij}$. This implies then
$\bar{d}^{2}_{ij}-\bar{d}^{2}_{ii^{\prime}}\gtrapprox 5$ while
$\bar{d}^{2}_{ii^{\prime}}\leq d^{2}_{ii^{\prime}}=3$. Hence,
$\bar{d}^{2}_{ij}/\bar{d}^{2}_{ii^{\prime}}\gtrapprox 8/3$. This second case
shows that one of the closest observation/simulation pair is definitely
further apart than this specific observation is with other observations, at
least according to $\phi_{\text{final}}$. This example also supports that the
potential bias on the $d^{2}_{\text{obs}}$ estimators remains moderate for the
study of (OBS, SIM) pairs ($\bar{d}^{2}_{ij}\gtrapprox 5$ and $d^{2}_{ij}=8$
implies $b_{ij}\lessapprox 3$).

We conclude that the minimum (OBS, SIM) distance value, obtained on all
possible pairs between those datasets, evidences a meaningful but moderated
distinction between these datasets. Note that the $\phi_{\text{final}}$
statistics set may still have some degeneracies, particularly for the (OBS,
SIM) comparison, and that accounting for these should likely deepen this gap.
We however believe that current analysis illustrates the usefulness of such a
distance, and leave a more detailed study for future work.

Figure 17: Comparing closest (OBS, SIM) pairs $(i,j)$ (gray boxes) to
neighboring (OBS, OBS) pairs $(i,i^{\prime})$. Top row focuses on the closest
observation to simulations, that is a patch of Ophiuchus, while bottom row
focuses on the second closest patch, that is in Aquila cloud. In the latter
case, neighboring (but independent) observations of Aquila are closer than the
closest simulations. The color-bars can change from a pair to another.

##  7 Conclusions

In this paper, we aim at studying the diversity of morphologies of observed
molecular clouds. To do so we construct a set of $\sim 500$ patches of size
$512\times 512$ pixels, extracted at different resolutions from column density
maps of 14 nearby clouds derived from the HGBS dust emission observations. We
compute several sets of statistics (mean, variance, quantile function, power
spectrum and scattering transform) from these maps, and we compare their
informative power. To do so, we introduce a new methodology (Fig. 6) that
allows to confront, without any supervision, two sets of summary statistics on
their respective abilities to detect statistical incompatibility between pairs
of processes in a given dataset.

Applying this methodology to this set of observations, we find that Gaussian
statistics have degeneracies for the observational dataset, some of which can
be lifted by one-point statistics (Fig. 7). We then show that even if log-
Gaussian statistics are much less degenerate, a compressed set of scattering
statistics still succeeds to demonstrate further degeneracies (Fig. 11). This
confirms that the diversity of morphologies arising in these observed clouds
cannot be sufficiently described by Gaussian nor log-Gaussian statistics. This
means that using such descriptions, typically to compare numerical simulations
with observational data, can lead to misleading projections: in addition to
missing potential absolute discrepancies, observations with different
properties can still be matched to the same simulation.

We apply the same diagnostic to simulations and to the set of logFBMs (Fig.
14), and find deviations from observations’ geometry. In particular, the
outcome of this diagnostic strikingly supports the sufficiency of log-Gaussian
statistics to discriminate between logFBM patches. Regarding observations and
simulations, there remain deviations between their respective geometries, even
though great care has been taken in this work to design robust and low
dimensional sets of summary statistics. This supports the difficulty of
transferring simulation-based priors121212Such as the outcome of a Fisher
analysis, or a trained neural network. to observations, especially high-order
information content, and also buttresses the supervision-detached approach
developed in this paper, along with the choice to work with compressed and
robust summary statistics.

From these results, we introduce a morphological distance $d^{2}_{\text{obs}}$
based on a set of summary statistics $\phi_{\text{final}}$ composed of 7 log-
Gaussian and 10 RWST coefficients. The similarity probed by this distance is
in agreement with visual impression when comparing observations with
themselves, with logFBMs and with simulations (Fig. 16). It remains however
insufficient to operate on datasets made of a wider diversity of textures such
as DTD.

This work opens multiple perspectives:

  * •

the methodology we developed to confront summary statistics requires very few
assumptions: it can operate in an unsupervised and very low data regime. Hence
such methodology can easily be assimilated by the ISM community and applied on
a wider set of statistics and physical tracers of data (velocities,
polarization, temperature).

  * •

The low dimensional and analytical morphological embedding
$\phi_{\text{final}}$ developed here allows for a directly interpretable
comparison, for instance between different observed clouds, between
observations and simulations or statistical models, but it also paves the way
for the use of more sophisticated unsupervised learning techniques.

  * •

Leveraging saliency maps
$\nabla_{\text{pixels}}d_{\text{obs}}^{2}(x_{i},x_{j})$, the distance
$d^{2}_{\text{obs}}$ can be used to highlight the main areas responsible for
morphological discrepancies between two patches $(x_{i},x_{j})$, broadening
the scope of the work initiated by Peek & Burkhart (2019) to the unsupervised
world of observations.

  * •

The confrontation methodology can be used to make a feature selection
algorithm, in the spirit of the FRAME model developed by Zhu et al. (1998),
but designed to optimize the comparison task of a non-parametric collection of
processes $\\{p_{i}\\}_{i}$, instead of modeling a single process.

The following improvements could also be of great benefit:

  * •

reduce the uncertainty in the compatibility diagnostic $d^{2}_{\phi}$ due to
the precision matrix estimation. A promising idea is to use a maximum entropy
model conditioned on the data on which to perform the precision estimation.

  * •

Reduce the overall dependency of the distance $d^{2}_{\mathcal{D}}$ on the
dataset $\mathcal{D}$, and in particular on its outliers. This can be done by
building a localized metric, based on local estimations of the geometry,
restricting $\mathcal{D}$ to $\mathcal{D}_{\text{loc}}[i]$, the neighbors of a
given process $i$.

###### Acknowledgements.

This research has made use of data from the Herschel Gould Belt survey (HGBS)
project131313http://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 work reused datasets available on the Galactica simulations
database141414http://www.galactica-simulations.eu.
This work reused the Describable Textures
Dataset151515https://www.robots.ox.ac.uk/~vgg/data/dtd/.

## References

  * Allys et al. (2019) Allys, E., Levrier, F., Zhang, S., et al. 2019, A&A, 629, A115, publisher: EDP Sciences
  * André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, Protostars and Planets VI, 27
  * André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102
  * Appel et al. (2022) Appel, S. M., Burkhart, B., Semenov, V. A., Federrath, C., & Rosen, A. L. 2022, The Astrophysical Journal, 927, 75
  * Arzoumanian et al. (2011) Arzoumanian, D., André, P., Didelon, P., et al. 2011, Astronomy & Astrophysics, 529, L6
  * Auclair et al. (2024) Auclair, C., Allys, E., Boulanger, F., et al. 2024, Astronomy & Astrophysics, 681, A1
  * Bruna & Mallat (2013) Bruna, J. & Mallat, S. 2013, IEEE transactions on pattern analysis and machine intelligence, 35, 1872
  * Brunt & Heyer (2002) Brunt, C. M. & Heyer, M. H. 2002, The Astrophysical Journal, 566, 276
  * Burkhart et al. (2009) Burkhart, B., Falceta-Gonçalves, D., Kowal, G., & Lazarian, A. 2009, The Astrophysical Journal, 693, 250
  * Burkhart & Lazarian (2016) Burkhart, B. & Lazarian, A. 2016, The Astrophysical Journal, 827, 26
  * Cimpoi et al. (2014) Cimpoi, M., Maji, S., Kokkinos, I., et al. 2014, in Proceedings of the IEEE Conf. on Computer Vision and Pattern Recognition (CVPR)
  * Cover (1999) Cover, T. M. 1999, Elements of information theory (John Wiley & Sons)
  * Elmegreen (2002) Elmegreen, B. G. 2002, The Astrophysical Journal, 564, 773
  * Falgarone et al. (2004) Falgarone, E., Hily-Blant, P., & Levrier, F. 2004, Ap&SS, 292, 89
  * Federrath & Klessen (2013) Federrath, C. & Klessen, R. S. 2013, The Astrophysical Journal, 763, 51
  * Fromang et al. (2006) Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371
  * Galli et al. (2013) Galli, P. A. B., Bertout, C., Teixeira, R., & Ducourant, C. 2013, A&A, 558, A77
  * Goodman et al. (2009) Goodman, A. A., Rosolowsky, E. W., Borkin, M. A., et al. 2009, Nature, 457, 63
  * Heithausen & Thaddeus (1990) Heithausen, A. & Thaddeus, P. 1990, ApJ, 353, L49
  * Hennebelle & Chabrier (2008) Hennebelle, P. & Chabrier, G. 2008, The Astrophysical Journal, 684, 395
  * Heyer & Brunt (2004) Heyer, M. H. & Brunt, C. M. 2004, The Astrophysical Journal, 615, L45
  * Hotelling (1931) Hotelling, H. 1931
  * Kainulainen et al. (2009) Kainulainen, J., Beuther, H., Henning, T., & Plume, R. 2009, Astronomy & Astrophysics, 508, L35
  * Knude & Hog (1998) Knude, J. & Hog, E. 1998, A&A, 338, 897
  * Könyves et al. (2015) Könyves, V., André, P., Men’shchikov, A., et al. 2015, Astronomy & Astrophysics, 584, A91, publisher: EDP Sciences
  * Levrier et al. (2006) Levrier, F., Falgarone, E., & Viallefond, F. 2006, Astronomy & Astrophysics, 456, 205
  * Levrier et al. (2018) Levrier, F., Neveu, J., Falgarone, E., et al. 2018, Astronomy & Astrophysics, 614, A124
  * Mamajek (2008) Mamajek, E. E. 2008, Astronomische Nachrichten, 329, 10
  * McKee & Ostriker (2007) McKee, C. F. & Ostriker, E. C. 2007, ARA&A, 45, 565
  * Miville-Deschênes et al. (2007) Miville-Deschênes, M.-A., Lagache, G., Boulanger, F., & Puget, J.-L. 2007, Astronomy & Astrophysics, 469, 595
  * Miville-Deschênes et al. (2010) Miville-Deschênes, M. A., Martin, P. G., Abergel, A., et al. 2010, A&A, 518, L104
  * Ntormousi & Hennebelle (2019) Ntormousi, E. & Hennebelle, P. 2019, Astronomy & Astrophysics, 625, A82
  * Ortiz-León et al. (2018) Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJ, 869, L33
  * Ossenkopf-Okada et al. (2016) Ossenkopf-Okada, V., Csengeri, T., Schneider, N., Federrath, C., & Klessen, R. S. 2016, Astronomy & Astrophysics, 590, A104
  * Palmeirim et al. (2013) Palmeirim, P. a., André, P., Kirk, J., et al. 2013, Astronomy & Astrophysics, 550, A38
  * Park et al. (2023) Park, C. F., Allys, E., Villaescusa-Navarro, F., & Finkbeiner, D. 2023, The Astrophysical Journal, 946, 107
  * Peek & Burkhart (2019) Peek, J. & Burkhart, B. 2019, The Astrophysical Journal Letters, 882, L12
  * Peek & White (2021) Peek, J. & White, R. 2021, Bulletin of the AAS, 53, https://baas.aas.org/pub/2021n6i301p06
  * Pilbratt et al. (2010) Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1
  * Planck Collaboration (2020) Planck Collaboration. 2020, A&A, 641, A12
  * Pouteau et al. (2023) Pouteau, Y., Motte, F., Nony, T., et al. 2023, A&A, 674, A76
  * Regaldo-Saint Blancard et al. (2020) Regaldo-Saint Blancard, B., Levrier, F., Allys, E., Bellomi, E., & Boulanger, F. 2020, Astronomy & Astrophysics, 642, A217
  * Remazeilles et al. (2011) Remazeilles, M., Delabrouille, J., & Cardoso, J.-F. 2011, MNRAS, 418, 467
  * Robitaille et al. (2019) Robitaille, J.-F., Motte, F., Schneider, N., Elia, D., & Bontemps, S. 2019, Astronomy & Astrophysics, 628, A33
  * Saydjari et al. (2021) Saydjari, A. K., Portillo, S. K., Slepian, Z., et al. 2021, The Astrophysical Journal, 910, 122
  * Schekochihin (2022) Schekochihin, A. A. 2022, Journal of Plasma Physics, 88, 155880501
  * Schlafly et al. (2014) Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29
  * Schneider et al. (2013) Schneider, N., André, P., Könyves, V., et al. 2013, ApJ, 766, L17
  * Schneider et al. (2022) Schneider, N., Ossenkopf-Okada, V., Clarke, S., et al. 2022, Astronomy & Astrophysics, 666, A165
  * Srivastava & Du (2008) Srivastava, M. S. & Du, M. 2008, Journal of Multivariate Analysis, 99, 386
  * Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
  * Vázquez-Semadeni et al. (1997) Vázquez-Semadeni, E., Ballesteros-Paredes, J., & Rodríguez, L. F. 1997, The Astrophysical Journal, 474, 292
  * Yan et al. (2019) Yan, Q.-Z., Zhang, B., Xu, Y., et al. 2019, A&A, 624, A6
  * Zavagno et al. (2023) Zavagno, A., Dupé, F.-X., Bensaid, S., et al. 2023, Astronomy & Astrophysics, 669, A120
  * Zhu et al. (1998) Zhu, S. C., Wu, Y., & Mumford, D. 1998, International Journal of Computer Vision, 27, 107
  * Zucker et al. (2020) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51

##  Appendix A Other datasets

###  A.1 Numerical simulations

Figure 18: Overview of the simulations. One snapshot is shown for each model
(”Hydro high res”, ”MHD”, and ”MHD high B”), and the maps show the central
$33\,\mathrm{pc}\times 33\,\mathrm{pc}$ field for gas column density
integrated along the $y$ axis.

The set of numerical simulations used in this paper as an example of state-of-
the-art attempts to reproduce the physics of the ISM in silico is taken from
the ORION project of the Galactica database161616http://www.galactica-
simulations.eu/db/STAR_FORM/ORION/. Using the adaptive mesh refinement (AMR)
code RAMSES (Teyssier 2002; Fromang et al. 2006), the data simulates the
collapse of a dense molecular cloud under self-gravity, including decaying MHD
turbulence but without stellar feedback. The focus of these simulations is to
study the early stages of the star formation process in a molecular cloud
($10^{5}\,\mathrm{M}_{\odot}$ in a 66 pc cubic box), and so fits in well with
the HGBS observational data used in this paper.

Three classes of simulations are considered, one without magnetic field
(”Hydro high res”), one with a typical magnetic field of order
$5-10\,\mu\mathrm{G}$ (”MHD”), and one with a higher value of the magnetic
field of order $10-25\,\mu\mathrm{G}$ (”MHD high B”). With the adaptive mesh
refinement scheme, the spatial resolution of the models can reach down to 1
mpc ($\sim 200\,$au), and even $100\,$au for the ”Hydro high res” model. For
more details about these simulations, we refer the reader to Ntormousi &
Hennebelle (2019).

For each of these models, we take two snapshots, integrate along the three
axes, and crop the resulting column density images to keep the central
$33\,\mathrm{pc}\times 33\,\mathrm{pc}$ field, on a regular $4096\times 4096$
grid. Examples of such images are shown in Fig. 18. As we did for
observations, we cut these panels into $512\times 512$ patches, keeping only
the ones where the effective resolution is fine enough, to avoid artifacts due
to the AMR scheme affecting our morphological analysis. The pixels in the
resulting patches have sizes $8\,\mathrm{mpc}$. This is in accordance with the
typical spatial sampling found in our observational dataset. Indeed, this
corresponds to the finest pixels (6 ″) of intermediate-distance clouds $d\sim
300\,$pc, but also to the closest clouds $d\sim 150\,$pc sampled at 12 ″.

###  A.2 logFBM models

To perform our analysis of observational data, we also consider purely
synthetic, parametric models of column density maps, derived from
exponentiations of fractional Brownian motions (fBm). Such fields have been
previously studied by the ISM community (Elmegreen 2002; Brunt & Heyer 2002;
Miville-Deschênes et al. 2007; Levrier et al. 2018). It is indeed a convenient
class of models that allows to simultaneously reproduce exactly the log-normal
one-point properties (encountered in quiescent regions) and power-law power
spectra to a good approximation. An example is given in Fig. 1.

In this paper, we consider the following model:

$X\equiv e^{\sigma F_{\beta}\star
Z+\mu},\,\,\,\textrm{with}\leavevmode\nobreak\ Z\sim\mathcal{N}(\@vec{0},I),$
(20)

parametrized by $\mu$, $\sigma$ and $\beta$, corresponding respectively to the
mean, standard deviation and spectral index of the fBm. The latter is
generated by sampling a $512\times 512$ Gaussian white noise map $Z$, which is
then filtered with $F_{\beta}$, that is defined in Fourier space as:

$\tilde{F}_{\beta}(\@vec{k})\equiv\begin{cases}\propto k^{-\beta/2}&\text{ if
}\@vec{k}\neq\@vec{0},\\\ 1&\text{ if }\@vec{k}=\@vec{0},\end{cases}$ (21)

such that $\sigma F_{\beta}\star Z+\mu$ is a real valued Gaussian process with
mean $\mu$, variance $\sigma^{2}$ and $\beta$-decaying power spectrum power-
law.

The parameters $(\mu,\sigma,\beta)$ are fitted from the observations. More
precisely, for each $512\times 512$ column density patch that has no negative
pixel (that is about $87\%$ of the 550 observational patches), we estimate the
mean $\mu$, variance $\sigma^{2}$ and spectral index $\beta$ on the logarithm
of the column density patch, expressed in $10^{22}\,\mathrm{cm^{-2}}$ units.
This procedure leads to a distribution of $\sim 480$ estimated parameters
triplets $(\mu,\sigma,\beta)$ that is reported in Fig. 19. For each triplet,
we then sample one $512\times 512$ logFBM patch. Contrary to the other patches
considered in this paper (observations, simulations, everyday textures), these
samples have periodic boundary conditions, by construction. However, their
$256\times 256$ sub-patches, that are the only ones on which the various
summary statistics are applied, do not, just as the rest of the data.

We dub these synthetic maps logFBM models, to recall that the statistics of
their logarithms are those of the Gaussian, fBm random fields. In the
literature, emission maps are often produced by integrating a 3D generated
field. In this paper, for simplicity and computational efficiency, we directly
generate 2D fields. Even though integrated 3D fields could be expected to
better model MCs’ column density maps, they do not differ much in morphology
from 2D models and both are far from reproducing MCs’ morphology. Thus, the 2D
logFBM models used here are suited to the scope of this paper.

Figure 19: log-Gaussian parameters fitted on the observations used to sample
logFBM data.

###  A.3 Describable Textures Dataset (DTD)

To better understand the peculiarities of MCs’ statistical properties, we
compare them to a very wide set of everyday textures: the Describable Textures
Dataset (DTD) (Cimpoi et al. 2014). This dataset is in open
access171717https://www.robots.ox.ac.uk/~vgg/data/dtd/. In this paper, we keep
only the images with smaller edge size larger than 512 pixels. We then crop
each one of these $\sim 1050$ resulting pictures around its center into a
$512\times 512$ patch. Then each patch is converted to gray-scale by summing
its three RGB channels. These maps have integer values that range from 0 to
$3\times 255=765$ that we finally divide by a factor $10^{3}$ to match the
typical values of the other datasets.

##  Appendix B Apodization

The $256\times 256$ sub-patches on which we perform 2D Fourier transform do
not have periodic boundary conditions (PBC). Thus, before applying such
transform during power spectrum estimation on a given sub-patch $x$, we first
apodize it as follows:

$x\mapsto rw\cdot(x-\langle x\rangle_{\@vec{u}}),$

where $w\cdot$ denotes the pixel-wise multiplication by the window $w$,
reported in Fig. 20 and $r$ is a scalar factor depending on the input map $x$
such that the output apodized map has same variance as the input. The
resulting map does not share the same mean with the input but this property is
not regarded by the Fourier-based statistics used here.

Figure 20: $256\times 256$ apodization window $w$ used. The orange line is a
1D slice of the middle of the window.

##  Appendix C Srivastava & Du test statistic

The test statistic $d^{2}_{\phi}(x_{i},x_{j})$ we introduce in Eq. 2 is
derived from Srivastava & Du (2008). We apply this test on the two collections
of summary statistics $\\{\phi(x_{i}^{(l)})\\}_{1\leq l\leq N_{i}}$ and
$\\{\phi(x_{j}^{(l)})\\}_{1\leq l\leq N_{j}}$ computed on the sub-patches
(indexed by $l$) of the patch $x_{i}$ and the patch $x_{j}$. In our case,
$N_{i}=N_{j}=4$. We recall Eq. 2 here:

$d_{\phi}^{2}(x_{i},x_{j})\equiv\alpha\Big{[}\big{(}\hat{\mu}_{i}-\hat{\mu}_{j}\big{)}^{T}D_{S}^{-1}\big{(}\hat{\mu}_{i}-\hat{\mu}_{j}\big{)}-\beta\Big{]}.$
(22)

This equation is based on the following quantities:

$\hat{\mu}_{i}\equiv\langle\phi(x_{i}^{(l)})\rangle_{1\leq l\leq N_{i}},$
$p=\text{dim }\phi,\text{ and }n=N_{i}+N_{j}-2,$
$\begin{cases}S_{i}\equiv\langle\big{[}\phi(x_{i}^{(l)})-\hat{\mu}_{i}\big{]}\big{[}\phi(x_{i}^{(l)})-\hat{\mu}_{i}\big{]}^{T}\rangle_{1\leq
l\leq N_{i}},\\\ S\equiv\frac{1}{n}(N_{i}S_{i}+N_{j}S_{j}),\\\
D_{S}\equiv\operatorname{diag}S,\end{cases}$ $\begin{cases}R\equiv
D_{S}^{-1/2}SD_{S}^{-1/2},\\\ c_{p,n}:=1+trR^{2}/p^{3/2},\\\
\alpha\equiv\frac{N_{i}N_{j}}{N_{i}+N_{j}}\frac{1}{\sqrt{2(trR^{2}-p^{2}/n)c_{p,n}}},\\\
\beta\equiv\frac{N_{i}+N_{j}}{N_{i}N_{j}}\frac{np}{n-2}.\end{cases}$

##  Appendix D Why taking the logarithm of some standard statistics?

In theory, because the mapping $\psi\mapsto\log\psi$ is invertible, the same
amount of information should be contained in $\phi(x)$ and in $\log\phi(x)$
(when the latter is well defined). However, for the compatibility diagnostic
we chose, these two options might perform differently. Such issues are not
specific to this test and also occur with standard diagnostics such as Fisher
analysis Park et al. (2023).

First, this non-linear mapping can change the fulfillment level of the test
assumptions. In the Srivastava & Du test used here, the two input
distributions that are to be compared are assumed to be normally distributed
and have equal variance. In many cases, typically when $\phi(x)>0$, taking the
logarithm of $\phi(x)$ decreases the deviation to normality and decreases the
relative deviations between the variances of $\phi(x_{i})$ and $\phi(x_{j})$.

Second, this could lead to dramatic improvement in incompatibility detection,
as illustrated in Fig. 21 for $\phi=\phi_{\text{var}}$ (provided that the
assumptions of the test are not too discarded in both cases for these results
to be appropriately interpreted).

Figure 21: Effect of the log rescaling after computing variance statistics on
observation data. Provided that the assumptions of the test are not too
discarded in both cases for these results to be appropriately interpreted on
both axes, this shows that $\big{\langle}\big{[}x-\langle
x\rangle_{\@vec{u}}\big{]}^{2}\big{\rangle}_{\@vec{u}}$ is strictly more
degenerate than $\log\,\big{\langle}\big{[}x-\langle
x\rangle_{\@vec{u}}\big{]}^{2}\big{\rangle}_{\@vec{u}}$.