© The Authors 2025

1 Introduction

The interstellar medium (ISM) in galaxies comprises a variety of components with distinct densities, kinematics, chemical compositions, and thermal phases – from cold and warm neutral medium (CNM and WNM) to hot ionized gas, and from dense molecular clouds near the Galactic plane to diffuse structures at high latitudes. Interwoven in this multiphase ISM is a pervasive population of interstellar dust grains, which plays a crucial role in regulating star formation and chemical enrichment (e.g., Dwek 1998; Draine & Li 2007; Klessen & Glover 2016). Dust is involved in multiple radiative transfer processes; therefore, dust emission has been widely used as a diagnostic tool for probing the physical conditions within the ISM (e.g., Draine & Li 2001; Draine 2003a; Zubko et al. 2004; Compiègne et al. 2011).

Diffuse radiation from dust in the Milky Way, in the optical historically referred to as (the reflected portion of) the diffuse Galactic light (DGL), is a prominent manifestation of the diffuse ISM. All-sky infrared (IR) mapping (Low et al. 1984) revealed characteristically wispy, filamentary structures in dust emission prominently at high Galactic latitudes, which became referred to as Galactic cirrus. Similar structures had already been seen in optical imaging as high latitude reflection nebula (Sandage 1976). In this work we use the common term “cirrus” to denote this spatially resolved form of the DGL, recognizing that the radiative mechanisms are distinct at different wavelengths.

Infrared cirrus has been extensively investigated using observations from infrared missions such as the IR Astronomical Satellite (IRAS), the Diffuse Infrared Background Experiment (DIRBE) aboard the Cosmic Background Explorer (COBE), the Wide-field Infrared Survey Explorer (WISE), the Herschel Space Observatory (e.g., Laureijs et al. 1987; Boulanger & Perault 1988; Arendt et al. 1998; Zagury et al. 1999; Miville-Deschênes et al. 2007; Martin et al. 2010; Bracco et al. 2011; Sano et al. 2015; Schisano et al. 2020). Far-infrared (FIR) observations probe thermal emission from characteristically larger dust grains, whereas mid-infrared (MIR) observations detect primarily nonequilibrium emission from small grains. These observations have provided characterizations of the spatial distribution and spectral energy distribution (SED) of the IR DGL, thereby informing our understanding of the variations in dust properties, including dust temperature, composition, and emissivity (Planck Collaboration XXIV 2011; Planck Collaboration Int. XVII 2014; Planck Collaboration XI 2014).

The optical cirrus, or DGL, first revealed through early observations before the 90s (e.g., Elvey & Roach 1937; Lillie & Witt 1976; Sandage 1976; Mattila 1979; Guhathakurta & Tyson 1989), and more recently through deep charge-coupled device (CCD) imaging surveys (e.g., Witt et al. 2008; Ienaka et al. 2013; Miville-Deschênes et al. 2016; Román et al. 2020; Zhang et al. 2023; Zhao et al. 2024), originates predominantly from the scattering of the interstellar radiation field (ISRF) by large (a > 0.1 μm) dust grains (e.g., Weingartner & Draine 2001; Draine 2003b; Compiègne et al. 2011). Optical observations can thus offer unique insights into the physical properties of dust grains, scattering anisotropy, and the characteristics of the incident ISRF (Gordon 2004; Brandt & Draine 2012; Ienaka et al. 2013; Zhang et al. 2023; Zhao et al. 2024). In addition, the higher angular resolution achievable makes observations of the optical cirrus particularly valuable for investigating the turbulent cascade processes in the ISM (Miville-Deschênes et al. 2016).

Observations of optical cirrus have not yet been utilized extensively as a tool to study the diffuse ISM because, in the optically thin conditions most easily modeled, the cirrus is very faint, typically only a few percent as bright as the night sky. Similar to many other low surface brightness phenomena, the photometry of optical cirrus is susceptible to various systematic effects, such as stray light from off-axis diffraction and the extended point-spread function (PSF) of sources in the image, improper sky subtraction, and flat-fielding errors. These systematics substantially suppress, if not completely remove, the cirrus signals. Only recently have advancements in instrumental design and in data reduction tailored for low surface brightness science (e.g., Abraham & van Dokkum 2014; Trujillo & Fliri 2016; Liu et al. 2023; Watkins et al. 2024; Cuillandre et al. 2025) substantially overcome these challenges. These improvements have facilitated the imaging of optical Galactic cirrus with optimized observing and data reduction strategies (e.g., Miville-Deschênes et al. 2016; Román et al. 2020; Zhang et al. 2023; Liu et al. 2025), opening new opportunities to study the diffuse ISM.

Existing characterizations of optical Galactic cirrus have focused primarily on their photometric characteristics: surface brightness, colors, and correlations with their IR counterparts. “Pixel-by-pixel” correlations, however, necessitate downgrading high-resolution data to the lower resolution of the pair and may be nonlinear. Although such correlation encodes important insights on dust physical properties, the spatial coherence information about the hierarchical structures of cirrus is lost, because pixels of the same intensity are binned regardless of whether they are from structures with very different shapes and scales.

Statistical approaches to characterize the cirrus morphology remain rather limited but offer an intriguing path forward. Miville-Deschênes et al. (2016) found that the angular power spectrum of optical cirrus in deep imaging data from the Canada France Hawaii Telescope (CFHT) follows a power law with an index of γ = −2.9, consistent with Planck and WISE results in the same field. They found no break in the power spectrum, as an indication of energy dissipation, down to a physical scale of 0.01 pc. Marchuk et al. (2021) studied the fractal properties of optical and IR cirrus in the Sloan Digital Sky Survey (SDSS) Stripe82 field and found a mean 2D fractal dimension of 1.69 and 1.38, respectively. They concluded that this difference cannot be attributed solely to differing angular resolutions of optical and IR data and might reflect intrinsic physical differences such as imposed by the scattering phase function.

In this work, we provide characterizations of the morphology of optical Galactic cirrus using a suite of statistical approaches. We place our investigation of optical cirrus in a multiwavelength context by comparing the morphology of optical cirrus with those of other dust tracers. Although cirrus maps from various tracers may appear visually similar, to what extent are they quantitatively the same (or different) remains an under-explored question. We address this by seeking “distribution-to-distribution” approaches that preserve and extract structural information across a range of scales. Furthermore, observations in different wavelengths have varying beam sizes, sensitivity, reduction processes, and systematics. We corrected these effects where feasible, and where corrections are not feasible, we discuss how these effects could result in different statistics at certain scales. The main objectives of this paper are therefore to:

• Quantify the spatial coherence of optical Galactic cirrus across a range of angular scales.

• Investigate the morphological similarities and differences among cirrus maps at different wavelengths, corresponding to different dust tracers.

• Develop and evaluate statistical tools for exploratory and diagnostic data analysis in forthcoming deep imaging surveys.

The goals focus on developing a better understanding of cirrus for its own sake, but we have an additional objective. Galactic cirrus is often a source of foreground contamination on deep images, masking out other low surface brightness phenomena – including ultra-diffuse galaxies (Zaritsky et al. 2021), intracluster light (ICL) (Mihos et al. 2017; Kluge et al. 2025), and tidal features (Bílek et al. 2020; Martínez-Delgado et al. 2023) – which can profoundly impact their detection and measurement. Existing techniques to distinguish are based on colors (Román et al. 2020; Mattila et al. 2023; Smirnov et al. 2023) or supplementary IR data (Davies et al. 2010; Besla et al. 2016; Mihos et al. 2017) More recently Liu et al. (2025) proposed a method combining morphological filtering with color constraints to distinguish faint diffuse galaxies from cirrus. However, this approach might not be optimal for extragalactic sources with angular scales comparable to cirrus (e.g., ICL) or with visually similar morphology (e.g., tidal tails), and many imaging surveys may not have high resolution IR coverage or multiband photometry available at the same depth and/or resolution. Therefore, a single-band approach would be highly valuable to facilitate low surface brightness analyses. Another important objective of this paper is therefore to:

• Identify morphological clues of optical cirrus that can differentiate cirrus from other faint diffuse extragalactic emission.

The paper is organized as follows. Section 2 describes the datasets employed, supplemented by Appendix A on the Dragonfly Telephoto Array. Section 3 outlines statistical methods that quantify the cirrus spatial coherence using local probability density functions (PDFs), angular power spectra (supplemented by Appendix B) and its variants, and wavelet scattering transforms (WSTs, supplemented by Appendix C). Sections 4–6 present the principal results from application of these methods. Section 7 explores distinguishing between extragalactic light (from tidal tails) and cirrus based on a single band and discusses the principles of the methodology, with supplementary material in Appendix D. Finally, Section 8 presents a summary and prospects for application in deep wide-field imaging surveys.

2 Datasets

The subsections below describe the datasets used in this work. The main results are focused on optical imaging data obtained from the Dragonfly Telephoto Array (hereafter Dragonfly for short). To characterize the similarity and difference of dust morphology with different tracers, we used FIR data from Herschel, MIR data from WISE, and the radiance image from thermal dust modeling of Planck observations. We also used H I observations from DHIGLS as supplementary data. The beam widths and pixel sizes of the datasets described below are summarized in Table 1.

2.1 The Spider field

The “Spider” field was chosen as the example for this work. This diffuse region of the ISM at intermediate Galactic latitude, centered at (l, b) ~ (135°, 40°) and located at a distance of ~320 pc and a height of ~205 pc above the Galactic plane (Zucker et al. 2020; Marchal & Martin 2023), is part of the North Celestial Pole Loop (NCPL; Meyerdierks et al. 1991; Martin et al. 2015; Marchal & Martin 2023). For the most part, the field is faint and optically thin (Zhang et al. 2023). It has been observed by several facilities at a range of wavelengths, making it an interesting target for ISM investigations and dust modeling.

2.2 Dragonfly optical imaging

At visible wavelengths, Galactic cirrus radiation mainly originates from the scattering of starlight by larger-size (a > 0.1 μm) dust grains. Our primary dataset is deep optical imaging of Galactic cirrus obtained from Dragonfly. Dragonfly is a mosaic aperture telescope optimized for imaging diffuse extended emission (Abraham & van Dokkum 2014). Details about the telescope and data reduction are summarized in Appendix A.1. In particular, we adopt the background modeling recipes in Liu et al. (2023) to preserve the faint diffuse cirrus signal, which prevent it from being suppressed or smeared out by conventional sky subtraction methods.

This work uses the same images of the Spider field presented in Liu et al. (2023). The final co-add was a mosaic of six Dragonfly fields, with an average of 53 frames in g and 64 frames in r for each field passing quality checks. In the g (r) band, the 3σ surface brightness limit of the co-add is 28.5 (28.3) mag arcsec⁻² on 10″ × 10″ scales, and 29.6 (29.3) mag arcsec⁻² on 60″ × 60″ scales¹. The total field-of-view of the co-add is $4_{\cdot }^{\circ }1\times 3_{\cdot }^{\circ }8$. An RGB composite image constructed from the g and r-band co-adds (R:r; G:(g+r)/2; B:g) is displayed in the upper left panel of Figure 1.

We performed source removal following the procedures in Liu et al. (2025), followed by a component separation process to yield the diffuse light from cirrus. This process, which aims to eliminate the contamination from faint extragalactic sources and residuals in the source removal, is based on cirrus morphology and SED constraints calibrated using the Planck thermal dust model (Planck Collaboration XI 2014). The compact source removal and cirrus component separation are summarized in Appendix A.2. During this process, the images were converted to physical units in kJy sr⁻¹ with zero-points from Planck, assuming that the optical cirrus is the counterpart of Planck thermal emission. We then combined the decomposed g and r cirrus images into a brightness image approximating the V band² and resampled to a pixel size of 6″ to enhance the signal-to-noise ratio (S/N). The image size is [2460x2272] pix². This optical image of the Spider field is displayed in the top middle panel of Fig. 1.

2.3 Herschel 250 μm

FIR (sub-mm) observations probe thermal emission from interstellar dust. We used data from the Spectral and Photometric Imaging Receiver (SPIRE) instrument on the Herschel Space Observatory (Griffin et al. 2010). SPIRE produces images in three bands: 250 μm, 350 μm, and 500 μm. The 250 μm image was chosen because it has the highest S/N and spatial resolution. We retrieved the image from the Herschel Science Archive (HSA)³. We used the Level 2.5 products created using HIPE12 pipeline version 14.0 (Ott 2010). The SPIRE images have been zero-level corrected and calibrated in units of MJy sr⁻¹, with a beam FWHM of 18″ on 6″ pixels. The image size is the same as the optical image. The Herschel 250 μm image is shown in the top right panel of Fig 1.

To reduce the impact of cosmic infrared background (CIB) anisotropies on the morphology of dust emission, we first masked all bright point sources using the Herschel/SPIRE Point Source Catalog (HSPSC; Schulz et al. 2017)⁴. Next, we removed relatively faint sources. This was done by fitting a correlation between the intensities of Dragonfly g + r data and Herschel 250 μm data. The g + r image was scaled to remove the diffuse background to create a residual image for the detection of faint compact sources. We then masked all sources above S/N of 3, and replaced masked pixels with median filtering. To mitigate the impact of fainter sources below the detection limit, we further performed a low-pass filtering through multiplying the Fourier transformed image by a Gaussian filter whose FWHM corresponds to twice that of the Herschel SPIRE 250 μm beam. Effectively, this low-pass filtering smooths small-scale structures below this scale.

It should be noted that these conventional approaches are not optimal because the CIB anisotropies (CIBA) still has a residual contribution with power at low spatial frequencies (Viero et al. 2013; Matsuura et al. 2017; Singh & Martin 2022). More advanced decomposition techniques, such as by Auclair et al. (2024)⁵, ought to be preferred but are deferred to future work.

Fig. 1 Dragonfly RGB mosaic image of the Spider field at the top left and multiwavelength images from different dust tracers used in this work. Top middle: Dragonfly combination of g and r after source removal, approximating diffuse radiation in the V band. Top right: Herschel 250 μm. Bottom left: WISE 12 μm. Bottom middle: Planck radiance. Bottom right: DHIGLS H I LVC column density. See text for details. The images are displayed on a linear scale with the same contrast between 0 and the 99.99% quantile.

2.4 WISE 12 μm

The emission at MIR bands is dominated by the nonequilibrium emission from the smallest (a < 0.01 μm) dust grains such as polycyclic aromatic hydrocarbons (PAHs; Tielens 2008). We used the WISE 12 μm (W3 band) map from the all-sky WISE 12 μm atlas reprocessed by Meisner & Finkbeiner (2014)⁶, which optimized the removal of compact sources and artifacts. The WISE intensities were converted from counts to MJy sr⁻¹ following the prescription of Cutri et al. (2012), Sect. 4.4h.

We further performed a source residual cleaning. First, we did a first visual inspection and masked the negative holes and bright sources in the map. We then followed the same approach as for the FIR data, i.e., we subtracted a diffuse component scaled based on the g + r image and masked sources detected above S/N > 3 in the residual image. The masked pixels were replaced with median filtering.

The native FWHM of the WISE W3 band is $6_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}5$ but the reprocessed all-sky WISE 12 μm atlas has been smoothed to 15″. The cleaned WISE map was resampled to the finer pixel resolution of the Dragonfly image for comparison of structures on the same grid, which of course does not gain back the original WISE W3 resolution. The image size is again the same as the optical image. The WISE 12 μm map of the Spider field is shown in the lower left panel of Fig. 1.

2.5 Planck radiance

The all-sky thermal dust model derived from Planck observations (Planck Collaboration XI 2014) was used as complementary data (see also Planck Collaboration Int. XLVIII 2016). In particular, we used the product dust radiance map, defined as the integral of the model thermal emission: ℛ = ∫ I_ν dν. The ℛ map retrieved from the Planck Legacy Archive⁷ has a beam width of 5′ in a HEALPix representation with pixel size approximately $1_{\cdot }^{\mathrm{\prime }}72$ and for the Spider field was resampled to a pixel size of 1′. The Planck image size is [246 × 227] pix², amounting to the same angular extent as the optical image. The map was visually inspected and two point sources were masked and refilled by smoothing. The radiance map is shown in the lower middle panel of Fig. 1 in units of 10⁻⁷ W m⁻² sr⁻¹.

One major advantage of using the Planck thermal dust model is the reduction of the impact of the CIBA. In Liu et al. (2025), we also demonstrated the advantage of using ℛ over other quantities like optical depth as the dust surrogate (see also discussion in Liu et al. 2023). In summary, ℛ is a better representation of the total thermal emission observed by Planck and is less affected by optical depth effects in the illuminating radiation.

2.6 Different couplings to the interstellar radiation field

The morphology in each of the products discussed above obviously depends on the amount of dust from which the radiation arises along the line of sight. But for each product there is also a different coupling to the ISRF. For high latitude fields like the Spider, the incident ISRF outside of the cirrus structures can be treated as uniform across the field. However, this ISRF could be affected internally by the optical depth of the cirrus itself, which is strongly wavelength dependent, thus potentially impacting the morphology.

The cirrus in optical images arises from the scattering of light by the classical “large” grain component of the dust size distribution at particular optical wavelengths where the optical depth is small (though probably not completely negligible). The FIR radiation considered is thermal radiation from the same large dust grains responsible for the scattered light. The SED can be modeled as a modified blackbody. In equilibrium, the radiance, the integral over this SED, is equal to the integral of the radiation absorbed from the local ISRF. Most of the energy in the ISRF is in the optical and near infrared, where optical depth effects are not large. For more in-depth discussion, see Zhang et al. (2023). To the extent that the ISRF decreases, so does the radiance and the SED also shifts to slightly lower temperature. The emission in the 250 μm band, chosen to sample the SED at wavelengths longward of the peak of the SED, would decrease but not nonlinearly as for bands at the peak of the SED or shortward.

Thus it can be anticipated a priori that the optical-FIR coupling of our morphology-tracing products will be fairly robust. But any changes in properties of the large dust grains with environment across the field could induce some decoupling.

The MIR is nonequilibrium emission from small grains responding directly to the intensity of incident UV photons in the ISRF. Therefore, the morphology traced could be different not only because of changes in the grain size distribution (large vs. small) across the field but also because of changes in the shape of the ISRF arising from higher attenuation in the UV. For these fundamental reasons, some optical-MIR decoupling can be expected.

2.7 DHIGLS HI

Gas and dust are known to be spatially correlated, especially at high Galactic latitudes (Planck Collaboration XXIV 2011). It is therefore interesting to investigate the morphology of neutral hydrogen (H I) column density as an indirect tracer of dust. We used data for the DF field of the DHIGLS H I survey (Blagrave et al. 2017)⁸, which targeted intermediate-to-high Galactic latitude fields using the Synthesis Telescope at the Dominion Radio Astrophysical Observatory (DRAO). The H I emission can be distinguished by a variety of gas velocity components (VCs). When the emission is optically thin, the brightness temperature cube can be integrated over these component velocity ranges to obtain maps of the column density of H I (N_H I) for low-, intermediate-, and high-VCs (LVC, IVC, and HVC, respectively). Because IVC of the Spider field is faint and the HVC is negligible, we used the diffuse emission map of the LVC in units of 10¹⁹ cm⁻². The velocity range defining the H I LVC emission is 30 > v > −15 km s⁻¹ (row 2 of Table 2 in Blagrave et al. 2017). The H I image has a beam width of 55″ and pixel size 18″. The image size is [820x758] pix², again amounting to the same angular extent as the optical image. The H I map is shown in the lower right panel of Fig. 1.

3 Methods

We employ complementary families of statistical methods for structural analysis on Galactic cirrus. The local intensity statistics approach (Sect. 3.1) is focused on local non-Gaussianity, while Fourier-domain approaches (Sect. 3.2) quantify global coherence across scales. The scattering transform approach (Sect. 3.3) incorporates both local and global structural information; although physical interpretation of its wealth of outcomes can be subtle, benefits can arise from using summary statistics derived from this approach to build on insights obtained from other methods. As a concrete example of how all these approaches can be applied together to interpret a panchromatic dataset, they are applied to investigate the Spider field in Sects. 4, 5, and 6.

3.1 Local intensity statistics

A general approach to characterizing the ISM is to analyze the intensity statistics of the map or cube (e.g., 2D column densities/velocities, 3D densities). In particular, the probability density function (PDF) and its derived statistics encapsulate critical information regarding the physical processes that shape the ISM (Nordlund & Padoan 1999; Kowal et al. 2007; Burkhart et al. 2009). Previous studies have identified correlations between the PDF and turbulence (Burkhart et al. 2009), self-gravity (Burkhart et al. 2015; Corbelli et al. 2018), stellar feedback (Boyden et al. 2016), and astrochemistry (Boyden et al. 2018), the first seemingly most relevant in the Spider field. In this analysis, we compute local PDFs of intensity maps across the field of view for various dust tracers and examine their properties. We derive statistics from the PDFs to characterize cirrus morphology and quantify the similarities and differences among maps.

Fig. 2 Illustration of the extraction of local PDFs and PDF statistics. (a) Local PDF extracted from the intensity map within a circular region moving across the field. This example shows an extraction with radius of r = 3′ at the same position on the optical (left), FIR (middle), and MIR (right) map. The lower right panel shows the KDE-smoothed PDFs of the logarithm of the normalized intensity in the circular regions ($z=\log {\tilde{I}}_{\nu }$). The PDFs are used for computing statistics and distance metrics. DH indicates the distance between the PDF of the FIR/MIR data and that of the optical data within the particular subregion. (b) Skewness, kurtosis, and Gini coefficient of the local PDFs at the same position measured on the optical, FIR, and MIR cirrus data with varying kernel radius r. The rightmost panel shows the PDF distance of the FIR and MIR data relative to the optical as a function of kernel radius.

3.1.1 Probability density function

Local PDFs encode structural information on specific physical scales, assuming uniform distance to the density field (Burkhart et al. 2010; Corbelli et al. 2018; Boyden et al. 2018). Because cirrus structures span a wide range of scales, we employed a series of radii to extract coherent morphological information at different angular scales and examine its variation as a function of scale. Given the broad dynamical range of intensities in the maps, we computed the PDFs and their associated statistics using the logarithm of the intensity. This approach is also motivated by the fact that diffuse ISM morphology is governed largely by turbulence, where density fields driven by compressible isothermal turbulence typically follow a log-normal distribution (Passot & Vázquez-Semadeni 1998).

At a given position (x, y) on the I_ν map, we extracted a representation of the underlying PDF using the pixel intensity values I_ν,i within a circular region centered at (x, y) with radius of r. The PDFs are computed over on a grid of positions {x_i, y_i} by moving the circular kernel across the map, where the grid spacing is equal to one-third of the kernel radius. To compute local PDFs and compare among different tracers, we normalized the intensities in the region, {I_ν,i(x, y|r)} by the local mean value so that the mean of the normalized intensities ${\tilde{I}}_{\nu }$ is equal to one. We then took the logarithm of ${\tilde{I}}_{\nu }$ as the input map $z=\log {\tilde{I}}_{\nu }$ (these PDFs are referred to as p(z)dz below⁹.

To mitigate the influence of extreme outliers (e.g., due to inadequate removal of foreground/background sources), we applied a kernel density estimation (KDE) to the discrete values to produce a smooth empirical PDF $p(z=\log {\tilde{I}}_v\mid x,y,r)$. The bandwidth of the estimator h is determined using the empirical Scott’s rule: h = n^−1/(d+4), with n the number of pixels used for evaluation and d the number of dimensions.

Figure 2a illustrates the extraction of local PDFs using a circular kernel with a radius of 3′. The right-hand circular panels display the zoomed-in regions at a fixed position in the cirrus maps I_ν for optical, FIR, and MIR. The lower right panel shows the KDE-smoothed PDFs corresponding to the intensity distributions within the circular region. The PDFs show clear deviation from a single Gaussian and display multiple components that represent substructures that likely correspond to those identified by the dendrogram (Goodman et al. 2009). Both the 2D and 1D distributions indicate that, although the optical, FIR, and MIR maps appear visually similar on large scales (see Fig. 1), they exhibit local differences. In particular, optical cirrus more closely resembles the FIR emission than the MIR emission. This can be explained by the fact that dust scattering in the optical is mainly contributed by the same dust population that emits thermally in the FIR, while MIR emission originates from stochastic heating of ultrasmall dust grains (Draine & Li 2007). We discuss this further in Section 4. The greater similarity between optical and FIR in this local region is further corroborated by the much lower Hellinger distance D_H between the PDFs (see Section 3.1.2).

3.1.2 PDF statistics

While the local PDFs derived from maps encode comprehensive information, it is often more efficient to characterize the shape of the PDF and quantify differences using summary statistics and distance metrics. These measures are scale-invariant and their spatial variations across the field can be visualized. We computed the following statistics and metrics based on the PDF p(z):

• Skewness and Kurtosis: skewness and kurtosis are the third- and fourth-order moments of the PDF. They serve to characterize the shape of the PDFs. Skewness measures the symmetry of the distribution around the center value, with positive skewness indicating an excess of high values and negative skewness indicating an excess of low values. Kurtosis is a measure of the “peakiness” of the distribution; a distribution that is flatter than a Gaussian yields positive kurtosis, whereas a more concentrated distribution yields negative kurtosis. Skewness and kurtosis have been found to be correlated with turbulent properties, particularly the sonic Mach number (Burkhart et al. 2010). The PDF-weighted skewness and kurtosis are computed using the following equations: $$skew=\frac{1}{{\sigma }_z^3}\int (z-\overline{z}{)}^3\cdot p(z)dz,$$(1) $$kurt=\frac{1}{{\sigma }_z^4}\int (z-\overline{z}{)}^4\cdot p(z)dz-3,$$(2)

  where ${\sigma }_z^2$ and $\overline{z}$ are the PDF-weighted variance and mean given by: ${\sigma }_z^2=\int (z-\overline{z}{)}^2\cdot p(z)dz$ and $\overline{z}=\int z\cdot p(z)dz$.

• Gini coefficient: the Gini coefficient (“Gini”) is a measure of inequality in a given set of data values¹⁰. Gini ranges from 0 to 1, with a higher value indicating greater imbalance. For a continuous PDF, Gini can be computed using the following equation (Gastwirth 1972)¹¹: $$Gini=\frac{1}{2\overline{w}}\iint |w-w^{\mathrm{\prime }}|\cdot p(w)p\left(w^{\mathrm{\prime }}\right)dwd{w}^{\mathrm{\prime }},$$(3)

  where p(w) dw is first derived from p(z) dz by a min-max normalization along z after clipping extreme outliers so that w ∈ [0, 1]. The calculation of Gini does not require a center, which makes it well-suited for characterizing complex morphologies. We computed Gini for each map as a function of scale using the extracted $p(\log {\tilde{I}}_{\nu })$ at each grid position.

• Hellinger Distance: when comparing PDFs from different maps with the same kernel at a given position, it is useful to have a measure that quantifies their similarities. One distance metric commonly used for PDFs is the Hellinger distance (D_H). Compared to other distance measures such as the Kullback-Leibler Divergence, D_H is less sensitive to outliers. The Hellinger distance between two PDFs p(z) and q(z) is defined as (Hellinger 1909): $$D_{\mathrm{H}}=\frac{1}{\sqrt{2}}{\{\int [\sqrt{p(z)}-\sqrt{q(z)}{]}^{2}dz\}}^{\frac{1}{2}}.$$(4)

  By definition, D_H ranges from 0 to 1, with lower values indicating greater similarity between p(z) and q(z). However, it is noteworthy that perfect similarity in the PDFs does not necessarily imply a one-to-one correspondence of the maps at the pixel level nor in frequency space. In our analysis, we choose the optical map as the baseline, i.e., we compute the distance of the local intensity PDFs extracted from each map relative to that of the optical map at corresponding positions.

  Figure 2b illustrates the skewness, kurtosis, and Gini calculated from the local PDFs of the optical/FIR/MIR maps at the same position in Fig. 2a as a function of kernel radius r from 1′ to 6′. The right panel shows the Hellinger distance of the FIR/MIR data relative to the optical data. In this example region, the skewness and kurtosis between optical and FIR are marginally consistent above r = 2′; the small-scale differences are probably due to the presence of CIBA in FIR data, beam effects, and residuals from source removal in both bands. On nearly all scales r > 2′, Gini in the optical and FIR are remarkably consistent. However, statistics of MIR data, in particular, skewness and Gini, show very different trends from those of optical and FIR. The difference is also revealed by its larger PDF distance D_H to the optical data than FIR at all scales.

3.2 Fourier statistics

While local intensity statistics serve as useful diagnostics, they are often influenced by local anomalies and instrumental effects. Thus, a widely adopted approach in ISM studies is to characterize the structures in Fourier space by analyzing how structures correlate across different scales. It is also inherently connected to turbulent processes in the ISM, as theoretical models of the turbulent cascade and energy transfer (e.g., the Kolmogorov theory) are formulated in Fourier space.

We employed three widely used statistical tools: (angular) power spectrum, Δ-variance, and cross-power spectrum. It is worth noting that several other techniques – such as bispectrum/bicoherence (Burkhart et al. 2009) and multi-fractal analysis (Elia et al. 2018) – also characterize structures in Fourier space and provide additional phase information. However, these methods could be computationally demanding for high-resolution, wide-field data. Therefore, in this work, we adopt the stated three methods owing to their clarity and computational simplicity.

3.2.1 Power spectrum analysis

The power spectrum is a standard tool for studying the statistical properties of the diffuse ISM (e.g., Miville-Deschênes et al. 2002). Both observations and simulations indicate that the power spectrum of the ISM column density is closely related to the underlying turbulent flow and dissipation processes in the ISM (Miville-Deschênes et al. 2007). For an optically thin ISM, the 2D column density power spectrum is equivalent to that of the 3D density (Miville-Deschênes et al. 2002). The slope of the power spectrum, as well as the presence of any characteristic scale, could therefore provide valuable insights into how turbulence regulates the density structures of the ISM.

We computed the 1D power spectrum as follows. The 2D power spectrum is calculated as the square of the modulus of the 2D Fast Fourier Transform (FFT). To mitigate edge effects (i.e., Gibbs phenomena) resulting from the FFT applied to nonperiodic distributions, we performed an apodization using a split cosine-bell function. The shape parameters of the apodization kernel were tuned until the edge effects effectively disappear. This choice primarily affects the power at only the largest couple of scales. The 2D power spectrum was then azimuthally averaged in annuli at k to yield the 1D power spectrum, P(k), where k denotes the spatial frequency (wavenumber).

The 1D power spectrum, P(k) is modeled using the following formula: $$P(k)=B(k)\times [A{k}^{\gamma }+C{k}^{\beta }]+D,$$(5)

which consists of three components:

• (1)

  The power term of the dust radiation, Ak^γ;

• (2)

  The contribution from other astrophysical sources and noise in the foreground or background sky signals, Ck^β;

• (3)

  Instrumental noise and systematics at the CCD level, D.

The first two terms are convolved with a Gaussian beam-related factor, B(k) that can be finely adjusted in the fit. This model is adapted from the modeling in Miville-Deschênes et al. (2016), with the addition of the D term to account for instrumental effects that are not scale-dependent.

Previous studies have reported values of γ ranging between −2 and −5 for ISM emission in various environments (Kiss et al. 2003; Hennebelle & Falgarone 2012). Accordingly, we adopt a broad fitting range for γ between −2 and −5. The second term is typically assumed to represent Poisson noise (i.e., β = 0) but β could be negative due to residuals in source removal and/or non-white noise characteristics (e.g., 1/f noise signal from clustering of sources, for which β ≃ −1). Therefore we set the fitting range of β between 0 and −2. The constant term D is constrained with an upper limit set by the P(k) value at the highest frequency in the fitting range.

For the Planck radiance, the beam size is fixed at 5′. The single channel maps used for the modified black body fit are each contaminated by CIBA, which in principle would be characterized by different values of β because of a partially uncorrelated set of galaxies dominating in each map. It is not clear how this contamination would be propagated through the nonlinear operator, integration over the fitted model modified black body SED. We simply adopted β = 0 given that it is a high-order effect. The free parameter D is still present, because although the integral of the SED model might seem noise-free, another set of observations with an independent set of noises would produce a slightly different radiance (D is a measure of the reproducibility of the radiance over repeated observations).

For the H I DHIGLS LVC data, there is no cosmic background and so C is zero. Following Martin et al. (2015) and Blagrave et al. (2017), the noise term D is replaced by a noise template constructed from the power spectra of a set of emission-free channels and multiplied by a free parameter η that should be near unity. The beam is asymmetical but is taken to be cirular with FWHM ${56}_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}8$, which is an adequate approximation (see Blagrave et al. 2017 for further details).

The fitting is performed after excluding the P(k) values at the largest scale (i.e., smallest k) and at scales much smaller than the beam size (i.e., the FWHM of the PSF). The smaller the upper value of k is, the less sensitive the value of γ is to other details of the model. The fitting range is thus restricted as for each map.

3.2.2 Δ-Variance analysis

The Δ-variance method is a useful tool for characterizing the structure of the ISM. Originally introduced by Stutzki et al. (1998) and refined by Ossenkopf et al. (2008), this technique has been used to study the turbulence of molecular clouds regulated by various physical processes – such as magnetic fields, shock waves, and gravity – in both simulations (Bertram et al. 2015) and observational data (Boyden et al. 2016; Dib 2023), as well as H I gas in nearby galaxies (Dib et al. 2021). In Liu et al. (2025), this approach was used to quantify the amount of structure before and after the separation of optical cirrus.

The Δ-variance method can be regarded as a variant of the power spectrum method; it measures the power of structures across a range of spatial scales by convolving the image with a series of kernels of increasing width. Compared to power spectrum analysis, Δ-variance is more robust when handling observational data with irregular boundaries or missing pixels.

We adopt the implementation of Δ-variance calculation provided by the publicly available TurbuStat package (Koch et al. 2019), which utilizes the formulation and kernel separation outlined in Ossenkopf et al. (2008). Briefly, the implementation generates a series of Mexican hat (Ricker) kernels, separated into their core and outer components, convolves the intensity map with these kernels, and computes the weighted variance of the resulting map. As the kernel width varies, the Δ-variance is a function of the spatial scale L: $${\sigma }_{\mathrm{\Delta }}^2(L)={⟨{(I(x,y)\ast {\odot }_L)}^2⟩}_{x,y},$$(6)

where the average is computed over the entire map, * denotes convolution, and ⊙_L represents the kernel at scale L. The inverse variance map of the intensity map is used as a weight map in the spatial integration. By convention, L is referred to as the “lag.”

If the underlying density field exhibits a power-law behavior in its power spectrum, i.e., P(k) ∝ k^γ, then ${\sigma }_{\mathrm{\Delta }}^2$(L) is expected to also show a power-law scaling: $${\sigma }_{\mathrm{\Delta }}^2(L)\sim L^{\alpha },$$(7)

where the power index α is related to that of the power spectrum via α = −γ − 2 (Stutzki et al. 1998). We fitted a power-law function to the main portions of the Δ-variance spectra and examine this relationship.

3.2.3 Cross-power spectrum analysis

While the power spectrum measures the statistical properties of intensity fluctuations in the field, comparisons between datasets need to account for systematics, contamination by other sky signals, and instrumental noise effects first. The cross-power spectrum provides a complementary means to characterize the correlation of structures across spatial scales in two fields (Tristram et al. 2005). By construction, systematics, contamination, and noise exclusive to a single dataset are suppressed in the cross-spectrum, thereby enhancing the S/N of genuine correlation. This method has been adopted widely in analysis of various astrophysical maps in cosmology, such as the cosmic microwave backgrounds, 21 cm intensity maps, and weak gravitational lensing fields (e.g., Harnois-Déraps et al. 2016; Planck Collaboration XI 2020), yielding more robust constraints to theoretical models than auto-power spectra alone.

For two given maps, I_a and I_b, the 2D cross-power spectrum is defined as: $$C{P}_{a\times b}(\mathbf{k})=⟨\mathrm{Re}[{\mathcal{I}}_a(\mathbf{k})\times {\mathcal{I}}_b^{\ast }(\mathbf{k})]⟩,$$(8)

where ℐ_a and ℐ_b denote the discrete Fourier transforms of I_a and I_b, respectively (* denoting the complex conjugate), k is the wavevector, and ⟨·⟩ represents averaging over radial binning of k. Because the input maps are in real values and we aim to measure the correlations between common structures with the same spatial phase, we retain only the real part in the product of their transforms. In practice, the 1D cross-power spectrum was then calculated by azimuthally averaging the 2D spectrum over radial bins (annuli) in Fourier space. Following the prescription in Sect. 3.2.1, we applied an apodization using a split cosine-bell function to both maps¹².

We further calculate the cross-correlation ratio, which is the 1D cross-power spectrum normalized by the auto-power spectra of each map: $${\xi }_{a\times b}(k)=\frac{C{P}_{a\times b}(k)}{\sqrt{P_a(k)P_b(k)}},$$(9)

where P_a(k) and P_b(k) are the auto-power spectra of I_a and I_b, respectively. At a given spatial frequency k, one would expect ξ_a×b=1 for perfect correlation, 0 for no correlation, and negative values in the case of anticorrelation.

3.3 Wavelet scattering transform (WST)

Although power spectrum-based approaches provide useful information about the energy partition across scales, they are insufficient to extract the non-Gaussian information. This is particularly important because many physical processes in the ISM are intrinsically non-Gaussian. Consequently, two fields may exhibit identical first- and second-order statistics in Fourier space while displaying markedly different morphologies. One novel statistical tool for quantifying the non-Gaussianity and morphological characteristics of a physical field is the wavelet scattering transform. The WST has been successfully applied to cosmology (e.g., Cheng et al. 2020; Greig et al. 2022; Jiang et al. 2025) and ISM astrophysics (e.g., Allys et al. 2019; Saydjari et al. 2021; Lei & Clark 2023). Compared to conventional high-order statistics, the WST estimator offers advantages in robustness, rapid convergence, and stability against additive noise and deformations in observational data.

3.3.1 Formalism to generate WST fields

Here, we briefly summarize the formalism of the WST approach, which is illustrated in Figure 3. For a 2D input field I₀ = I(x, y), the scattering transform computes a set of first-order fields I₁ = I₁(x, y) by convolving I₀ with a family of Morlet wavelets {ψ^j,θ(x, y)} and applying a modulus operation: $$I_0\to I_1^{j_1,{\theta }_1}\equiv |I_0\ast {\psi }^{j_1,{\theta }_1}|,$$(10)

where $I_1^{j_1,{\theta }_1}$ represents a set of fields {I₁} characterized by the scale index j₁ and the orientation index θ₁ of the wavelet. The second-order fields I₂ = I₂(x, y) can be generated by applying the same process to I₁: $$\left\{I_1\right\}\to I_2^{j_1,j_2,{\theta }_1,{\theta }_2}\equiv |I_1\ast {\psi }^{j_2,{\theta }_2}|=||I_0\ast {\psi }^{j_1,{\theta }_1}|\ast {\psi }^{j_2,{\theta }_2}|,$$(11)

where $I_2^{j_1,j_2,{\theta }_1,{\theta }_2}$ represents a set of fields {I₂} characterized by the wavelet indices (j₁, θ₁) and (j₂, θ₂). By recursion, the n-th order fields are generated by: $$\left\{I_{n-1}\right\}\to I_n^{j_1,...,j_n,{\theta }_1,...,{\theta }_n}\equiv |I_{n-1}\ast {\psi }^{j_n,{\theta }_n}|.$$(12)

By convention, the number of scales J is defined in a dyadic sequence of 2^j for 0 ≤ j < J, with the largest scale 2^J being smaller than the size of the field¹³. Because the wavelets scatter the energy into larger scales, only coefficients with j₂ > j₁ are physically meaningful.

The number of orientations, Θ, defines wavelets with position angles ϕ = πθ/Θ in the range 0 ≤ ϕ < π generated by integers θ over the range 0 ≤ θ < Θ.

Fig. 3 Illustration of wavelet scattering transform applied to cirrus images, adapted from Figure 4 of Cheng & Ménard (2021). The input field is “scattered” through convolution with a bank of Morlet wavelets at different scales j and orientations θ (not displayed for clarity), followed by a nonlinear operation (modulus). The output fields are shown up to j2 = 3. The scattering coefficient Sn is computed by taking the spatial average of the output field. Only j2 > j1 coefficients are used.

3.3.2 Formalism to generate WST coefficients

The set of nth-order WST coefficients, S_n, is computed by taking the mean amplitude of the nth-order fields: $$S_n={⟨I_n⟩}_{x,y},$$(13)

where ⟨·⟩_x,y denotes spatial averaging over the field. These coefficients form a compact set of translation-invariant descriptors that capture the coherence of the structures among different scales and orientations. In practice, most of the information in a physical field is contained in the leading orders of the scattering coefficients, as the norm of the coefficients decays exponentially with increasing order. Therefore, as widely adopted in the literature (Lei & Clark 2023; Auclair et al. 2024), we focus on WST coefficients up to the second order.

Because the fields I_n averaged in the coefficients S_n are computed from the previous-order fields I_n–1, S_n is correlated with S_n–1 under the same index family. Therefore, it is usually convenient to normalize the coefficients according to s_n = S_n/S_n–1, which renders the coefficients dimensionless and scale-invariant.

3.3.3 Intepretability

Here we briefly describe the intepretability of the WST coefficients. Readers with deeper interests can refer to Allys et al. (2019) and Cheng & Ménard (2021) for further exposition.

The first-order coefficients, S₁, provide information similar to that obtained from the power spectrum (or the two-point correlation function), as both quantify the power of density structures as a function of scale. The primary differences are that the power spectrum adopts an ℒ² norm and a Fourier kernel, while the WST uses an ℒ¹ norm and localized wavelet kernels.

It is often more interesting to investigate the second-order coefficients S₂. When averaged over orientations (yielding isotropic statistics), S₂ (j₁, j₂) characterizes the clustering of patterns at a scale j₁ as a function of scale j₂, thereby encoding information about the coherence and interaction between scales. A larger value of S₂ (j₁, j₂) indicates that structures at scale j₁ occur more frequently on a scale j₂.

In the case where anisotropy is also of interest, it is often useful to reduce the second-order coefficients by averaging only over certain orientations. In particular, we follow the reduction formalism in Lei & Clark (2023), which retains the angular difference dependence Δθ = θ₂ − θ₁ and averages over θ₁: $$S_2(j_1,j_2,\mathrm{\Delta }\theta )={⟨S_2(j_1,j_2,{\theta }_2-{\theta }_1,{\theta }_1)⟩}_{{\theta }_1},$$(14)

where ⟨·⟩_θ denotes averaging over θ. Large components with small Δθ ~ 0 delineate coherent structures with certain scales distributed along the same direction on larger scales (e.g., cirrus filaments), whereas large components with Δθ ~ 90° correspond to structures predominantly extending in the perpendicular direction (e.g., sheets and patches).

3.3.4 Adopted summary statistics

We examine only the two summary statistics proposed by Cheng & Ménard (2021), which further reduce the second-order coefficients:

• Sparsity s₂₁: the sparsity is defined as the ratio between S₂ and S₁ coefficients averaged over orientations: $$s_{21}\equiv {⟨S_2/S_1⟩}_{{\theta }_1,{\theta }_2}.$$(15)

  This diagnostic measures the degree of deviation from Gaussian random fields. A higher s₂₁ indicates the structures being more localized, i.e., concentrated at certain positions in the field, whereas a lower s₂₁ suggests a more widespread distribution.

• Linearity s₂₂: the linearity is defined as the ratio between the parallel (Δθ = 0°) and perpendicular (Δθ = 90°) components of the coefficients averaged over orientations: $$s_{22}\equiv {⟨S_2^{\parallel }/S_2^{\perp }⟩}_{\theta }={⟨S_2^{\mathrm{\Delta }\theta =0}/S_2^{\mathrm{\Delta }\theta ={90}^{\circ }}⟩}_{\theta },$$(16)

  which is a scale- and rotation-invariant measure of the shape of structures. Fields with s₂₂ < 1 tend to exhibit curvy, bubbly, or swirly patterns, while fields with s₂₂ > 1 display more linear, wispy, filamentary features.

Both summary statistics depend solely on the two scales j₁ and j₂. Because S₂ and S₁ at different scales are correlated, the sparsity and linearity may also exhibit correlations across various scale combinations. Consequently, in our analysis, we further compute the sparsity and linearity averaged over all scale combinations with 0 ≤ j₁ < j₂ < J, denoted as ${\tilde{s}}_{21}$ and ${\tilde{s}}_{22}$, which are treated as scale-averaged statistics.

3.3.5 Computation

The WST coefficients are computed using the publicly available code scattering (Cheng et al. 2020), which is developed based on the Kymatio package (Andreux et al. 2018). This implementation performs the wavelet convolution through FFT and supports GPU acceleration.

The choice of the number of scales (J) and the number of orientations (Θ) is a necessary input for WST analysis. In this work, we choose J = ⌊log₂N − 1⌋, where N is the minimum dimension (in pixels) of the field and ⌊·⌋ denotes rounding to the nearest integer. Thus, the statistical description is restricted to scales that are smaller than or comparable to half the size of the field to ensure reliable estimation. Following Regaldo-Saint Blancard et al. (2020) and Saydjari et al. (2021), we set Θ = 8 to decompose the half-circle into discrete directions, which provides as a reasonable trade-off between the compactness of the descriptor set¹⁴ and smooth sampling of the directions.

We did not apply apodization in the computation because in comparison with Fourier power spectra the use of localized wavelets in WST analysis makes it less sensitive to nonperiodic boundary conditions. However, we note that the cutoff of cirrus structures at the field edges might still bias the measured shape parameters (such as s₂₂) for the largest couple of scales. This effect, however, does not affect the relative comparisons between datasets in our analysis below.

Note our analysis applies to the intensity image as I₀ rather than the logarithm of intensity. Therefore, when applying such analysis for studying the cirrus morphology, bright interlopers – such as nearby galaxies or bright stars – should be carefully removed to minimize their impact on the non-Gaussian statistics of the field. This task, however, could be challenging in a crowded field (such as galaxy cluster or near the Galactic center) if the extended PSF wings dominate the sky background.

3.4 Summary

The statistical methods adopted in this work are summarized in Table 2. These methods cover a range of algorithmic complexity, describing the morphology of cirrus structures across angular scales in different angles, with their own strengths and limitations. In the following sections, we apply these methods to statistical characterizations of the morphology of optical cirrus. For a multi-wavelength context, we compare these results with results obtained using alternative dust tracers, highlighting the quantitative similarities and differences relative to optical cirrus. Where relevant, we discuss the connection between the observed structures and the interplay of physical mechanisms and instrumental effects that shape them.

4 Local intensity statistics

Here, we show analysis based on the local intensity statistics described in Section 3.1 using Dragonfly optical, Herschel FIR, and WISE MIR data. The motivation here is to look into the “neighborhoods” around each pixel as an extension of the pixel-by-pixel comparison. Planck and H I maps are not used for comparing statistics given their significantly larger beam widths.

However, the H I map is used here to exclude (mask) regions with very low ISM column densities, because, unlike simulations, PDFs in low column density regions in observational data have higher uncertainties and may suffer from incomplete sampling at the faint end (Lombardi et al. 2015). We produce a mean H I map by computing the local mean value using a moving circular region (with the same kernel scale). We then exclude low column density regions with N_H I below the first quartile of the N_H I distribution in the mean H I map. For reference, this corresponds to 2.6 × 10⁻¹⁸ cm⁻² for d_k = 5′.

4.1 PDF moment maps: a visualization of where morphology changes

Figure 4 shows the spatial distribution of skewness, kurtosis, and Gini extracted with a characteristic scale d_k (d_k = 2 × r) of 10′, which correspond to ~0.9 pc assuming a distance of 320 pc (Zucker et al. 2020). The maps were produced via bilinear interpolation using the grid values, as described in Sect. 3.1.

In general, regions with positive skewness occur at the boundaries of sharp filaments and extend alongside bright structures, indicating strong discontinuities. Given the correlation between higher-order moments and turbulence (Kowal et al. 2007; Burkhart et al. 2010), these discontinuities likely correspond to the fronts of ISM turbulence – possibly undergoing shearing or tidal flows. In the absence of energetic processes such as supernova explosions or stellar winds, the morphology is primarily shaped by intermittent turbulence and magnetic fields.

Positive kurtosis occurs more sporadically but is often cospatial with positive skewness. Conversely, negative kurtosis regions are mainly found along the ridgelines of filaments.

At small angular scales the skewness and kurtosis maps of FIR and MIR images exhibit several localized peaks. This is due to high-order PDF moments being vulnerable to extreme values or outliers, for example, residuals from imperfect source removal. The FIR image also tends to have higher values on average, which could be attributed to the presence of residuals of CIBA.

The spatial distribution of Gini is similar to that of skewness: the boundaries of cirrus filaments and patches exhibit higher Gini, indicating greater density inequality in these regions, while patches and ridgelines have lower Gini, reflecting a flatter intensity distribution. Indeed, Gini is correlated with skewness. However, it should be noted that Gini and skewness still trace different structural aspects because they respond differently to the same PDF. For instance, a local region containing two symmetric peaks relative to the mean may exhibit near-zero skewness, but Gini can be high given its significant deviation from a uniform distribution. The Gini maps of the FIR and MIR data present fewer local peaks at small angular scales, suggesting a much lower susceptibility to source residuals.

In summary, the local PDF departs strongly from a Gaussian or a uniform distribution at the edges of sharp structures such as filaments, revealing local intermittency. The PDF moments thus serve as good indicators of where cirrus morphology changes.

4.2 Which is closer to optical cirrus locally: FIR or MIR?

The left column of Figure 5 shows the histograms of the distributions of skewness, kurtosis, and Gini across the field for scales d_k ranging from 2′ to 20′. The lower end of the range of kernel scale is chosen to be several times larger than the beam widths to ensure reliable statistics. The upper end still retains good spatial resolution, while providing about sufficient points to sample the PDFs after masking the low ISM column density regions. Fainter histograms represent smaller d_k.

For skewness, most values lie between −1 and 1. The skewness distribution of optical data is centered around zero, whereas the FIR data tend to be more positive and the MIR tend to be more negative. Kurtosis values range between −1.5 to 1 and peaks around −0.5. The optical data display a broader distribution than the FIR and MIR data, likely owing to the smaller beam width, which allows it to capture finer structures. The distributions of Gini for the three dust tracers are similar to those of skewness in a relative sense.

The middle column of Fig. 5 shows the distributions of the differences between the statistics of the optical and FIR/MIR maps. For all statistics, the differences between optical and FIR (orange histograms) peak around zero. However, the differences in skewness and Gini between the MIR and optical maps (green histograms) are skewed toward positive values, indicating systematic offsets. The kurtosis difference relative to optical is broader for MIR than for FIR.

These trends are further demonstrated in the right column of Fig. 5, which shows the mean absolute difference as a function of scale. The mean of absolute difference quantifies the average discrepancy between the two maps. For d_k < 4′, the differences between optical and FIR statistics are comparable to those between optical and MIR, while both follow a decreasing trend with increasing scale. This is likely due to residual CIB. For d_k > 4′, the differences continue to drop at increasing scales for the FIR data, exhibiting more similar statistics to those of the optical data, whereas the MIR data present larger systematic offsets. These trends reinforce that the structures probed by optical and FIR are more similar than what MIR traces, as anticipated in Sect. 2.6.

The greater similarity between FIR and optical can be further corroborated by the Hellinger distance D_H. The left panel of Figure 6 shows the histograms of D_H across the map extracted using various kernel scales, where lighter shades represent smaller scales. The distributions for FIR relative to optical are significantly narrower and peak at smaller values than those for MIR. The right panel of Fig. 6 plots the mean D_H as a function of kernel scale, with the FIR distance being consistently smaller than the MIR distance at all scales. Because the PDF distance directly measures the proximity of local intensity distributions between two maps, these results indicate that the morphology of dust in optical imaging is closer to that traced by FIR than by MIR. This finding is consistent with the distinct origins of the diffuse light in the three bands, with the same characteristic dust size contributing to the optical and FIR but not the MIR (Sect. 2.6).

Furthermore, we note that the PDF distance approaches constant (~0.15 for MIR and ~0.09 for FIR) on larger scales. Because the beam widths in FIR and MIR are several times smaller than the kernel scales, we attribute the increase in PDF distances at small scales to the presence of CIB residuals. We infer that if CIBA were optimally removed from the images, the PDF distance for the FIR data would shrink to a constant down to the beam width. At scales comparable to the beam width, PSF blurring smooths the density structures, and consequently smooths the intensity PDFs.

Fig. 4 Spatial distributions of local skewness, kurtosis, and Gini extracted from the Dragonfly optical (left column), Herschel FIR (middle column), and WISE MIR (right column) data, using a moving circular kernel with a kernel width of dk = 10′. These maps trace where cirrus morphology changes.

Fig. 5 Left column: histograms of skewness, kurtosis, and Gini measured from the local PDFs in optical, FIR, and MIR data. Fainter histograms represent smaller kernel scales used for PDF extraction. Middle column: histograms showing the distributions of differences of the statistics between FIR or MIR and optical data. The difference is calculated by [optical-X], where X stands for the band. Right column: mean absolute difference of statistics between FIR or MIR and optical as a function of kernel scale.

4.3 ℳ_s in diffuse ISM inferred from dust morphology

The sonic Mach number, ℳ_s, is a measure of the compressibility of the medium. It is defined as the ratio of the local flow velocity to the sound speed of the medium. This parameter is crucial for hydrodynamics of astrophysical fluids, as the physics of compressible turbulence differs markedly from that of incompressible turbulence. Thus, investigating the variation of ℳ_s in the ISM density field is of great interest.

ℳ_s is known to correlate with statistical moments (variance, skewness, and kurtosis) of both observed column density distributions and simulated density fields (Kowal et al. 2007; Burkhart et al. 2009). The degree of non-Gaussian asymmetry in the PDF increases with ℳ_s. Although this correlation is strong for supersonic models, it is relatively weak for subsonic models where additional physical mechanisms may be at play; hence, caution is warranted when interpreting subsonic regions.

While higher-order moments generally correlate with ℳ_s, we follow Burkhart et al. (2010) in using the fourth moment (kurtosis) to estimate ℳ_s via the empirical relation: ℳ_s ~ (kurt +1.44)/1.05, which exhibits an almost linear correlation for ℳ_s in the range 1–8 and is less affected by boundary effects. This model is derived from simulations with Alfvén Mach number ℳ_A = 0.7 and the moments show little dependence on ℳ_A (Kowal et al. 2007). We similarly assigned regions with very low kurtosis (< −1) with ℳ_s = 0 to mark subsonic conditions whose behavior is not well characterized. We masked a high-kurtosis region in the lower right of the field with a circular aperture, as it likely results from a residual of source removal.

Figure 7 shows the ℳ_s map derived from the local intensity PDF of optical cirrus maps using a kernel scales of 10′. Most regions are subsonic (ℳ_s < 1 or kurt < −0.39) or transonic (ℳ_s ~ 1) – see the kurtosis histograms in the middle of the left column in Fig. 5. For a scale of 5/10′, ~55/45% of the map comprises regions with ℳ_s < 1, and only ~9/7% of the map has ℳ_s > 1.5 (kurt > 0.14). These results are qualitatively consistent with those reported for the H I in the Small Magellanic Cloud (Burkhart et al. 2010), where approximately 90% of the areas in the Small Magellanic Cloud exhibit ℳ_s < 2 based on PDF statistics. On average, ℳ_s in the diffuse ISM is considerably lower than that found in star-forming regions, where shocks from stellar feedback can compress the medium and efficiently produce supersonic turbulence. In the right panel we overlay contours of regions with ℳ_s > 1.5 and ℳ_s < 0.5 (in cyan and red, respectively) on the mean H I column density map. These regions tend to cluster along or near the edges of the cirrus filaments and patches, which suggest transition of cirrus morphology from more compact structures to diffuse ones.

Fig. 6 Left: histograms of the Hellinger distance (DH) between FIR or MIR and optical maps, which indicates the proximity of two cirrus maps. Smaller DH represents a closer local PDF. Histograms with lighter shades represent smaller kernel scales used for the PDF extraction. The PDF distance of FIR relative to optical is systematically lower than that of MIR. The dispersion is also much smaller. Right: mean DH between FIR (orange) or MIR (green) and optical data as a function of kernel scale. The mean DH of FIR relative to optical is lower than MIR in all scales. The dashed lines show the asymptotic values at large scales.

Fig. 7 Left: maps of sonic Mach number inferred from optical cirrus map illustrated with a kernel scale of dk = 10′. The circular kernel is illustrated with the purple disk. Right: mean H I column density fields extracted with the same kernel in the upper panels. Regions with ℳs < 0.5 and > 1.5 in the upper panels are indicated by red and cyan contours, respectively, which aggregate at edges of cirrus filaments or patches.

5 Fourier statistics

Although local intensity statistics probe where the field becomes non-Gaussian and provide a localized picture of cirrus structures, a global picture of how the structures are distributed across scales is missing. Structural analysis in Fourier space is better suited to extracting such information. Interpretation of the power-law slope can be directly linked to turbulence theories with supports from other investigations in the literature.

In this section, we show the results of Fourier-domain statistics, as described in Sect. 3.2. We present the power spectrum and Δ-variance of the optical cirrus in the example dataset. We also compute the results from other dust tracers for comparison with optical cirrus, including cross-correlations.

5.1 Power spectrum analysis

5.1.1 A unified power spectrum across dust tracers

The power spectra were computed following the procedures in Sect. 3.2.1. The fitting was performed using Capfit¹⁵, for robust nonlinear least square fitting. The power indices of the signal term and noise term from the fitting and their estimated uncertainties are listed in Table 3. The ISM terms for the optical, FIR, and MIR maps have power index values ranging from γ = −2.82 to −3.00, which is comparable to values found using IRAS 100 μm data (Gautier et al. 1992; Miville-Deschênes et al. 2002; Miville-Deschênes et al. 2007) and using optical data from CFHT (γ ~ −2.9; Miville-Deschênes et al. 2016). The power index for the Planck radiance lies near the lower (steeper) end of this range. For the H I LVC map, γ = −3.00, which is also in accordance with the values reported Blagrave et al. (2017).

The index of the sky contamination (noise) term, however, shows differences among the three maps in which it was included in Eq. (5). The optical map has a nearly flat noise term (β ~ −0.1), whereas the FIR has a gray noise closer to β ~ −1. This difference is likely due to the presence of CIBA residuals and/or imperfect source removal in the FIR maps. We find a flatter dependence for the MIR (β ~ −0.3). Miville-Deschênes et al. (2016) reported a flat noise term for WISE and a non-flat noise term for CFHT MegaCam. The difference in the optical might be attributed to differences in the sky area and source removal procedures. In any case, the variations in this noise term do not significantly affect the power index determined for the ISM component, which in our fitting is mainly dependent on the power at large scales.

Figure 8 shows the combined 1D power spectra of maps of various dust tracers. We only show the ISM component (the k^γ term in Eq. (5)), with other components subtracted from the total power spectra. The power spectra are vertically shifted to an arbitrary normalization. The overlapping ranges plotted are conservative relative to the available data (c.f., Fig. 10), ensuring that the ISM model is not sensitive to details of the multi-component modeling.

In general, the power of structures decreases at smaller scales, following power-laws that are remarkably similar, despite being derived from very different data products or wavelengths. The MIR map exhibits a similar power spectrum to optical and FIR even though it traces a different dust population. Also, the dust emission maps show power spectra consistent with that of H I gas. This correspondence suggests that the power of density fluctuations on different structural scales, when averaged into 1D in Fourier space, is highly coherent in diffuse areas such as Spider¹⁶. This is interesting considering that the tracer intensities are not always linearly correlated with each other, e.g., the optical and FIR correlation becomes nonlinear at high intensities.

The fact that the power spectrum obtained from the Dragonfly optical imaging follows the same well-defined power law as other tracers indicates success in preserving cirrus light on large scales. By contrast, in Appendix B.2, we show the power spectrum analysis on supplementary optical data from the DESI Legacy imaging survey (Dey et al. 2019), where the background systematics suppress or remove cirrus light on large scales and thus bias the power spectrum measurements.

Fig. 8 Combined 1D power spectra of different dust/gas tracer maps. We only show the ISM components of the power spectra, i.e., the contaminating sky noise and instrumental noise components are subtracted. The power spectra from different maps are similar, with a power index of γ ranging from −2.8 to −3.0 (Table 3). Each power spectrum is displayed for a conservative range at higher k, excluding where the contaminating noise is most prominent and/or the cutoff by the beam is significant. The normalization is in arbitrary units, with the power-laws being aligned in the range where they overlap at lower k. The power spectra are overplotted in order of increasing beam size (except FIR and MIR, guided by Fig. 10), so that it can be appreciated that the spectrum for the Dragonfly optical (g+r) map extends to the highest k. The dashed blue line shows the result of Miville-Deschênes et al. (2016) with P(k) ∝ k−2.9.

5.1.2 Insights on turbulence in the diffuse ISM

The shape of the power spectrum provides valuable information about turbulent flows. Previous studies revealed that the slope of the power spectrum correlates with ℳ_s (Kim & Ryu 2005; Burkhart et al. 2009). In Sect. 4.3, our results using local intensity statistics indicate that most regions in this field are subsonic or transonic (ℳ_s < 1.5). A subsonic isothermal compressible turbulent flow is expected to exhibit a power spectrum of γ = −11/3 ≈ −3.7 in the column density distribution (Saury et al. 2014), which is steeper than the γ ~ 2.8–3 slope we observed and those reported in other observations in the diffuse ISM (e.g., Miville-Deschênes et al. 2016). In contrast, supersonic compressible turbulence can produce much shallower power spectra than the subsonic regime (Kim & Ryu 2005). This discrepancy can be reconciled if the atomic gas in the local diffuse ISM is thermally bistable (i.e., non-isothermal), without requiring supersonic turbulence (Saury et al. 2014). Miville-Deschênes et al. (2016) found that the 21 cm line emission in their data indicated thermally bistable H I, and concluded that the γ ~ −2.9 spectra found in their results can be caused by thermal instability in H. The power spectra in our results, combined with the results in Sect. 4.3, support the scenario in which thermal instability could play a non-negligible role in shaping the cirrus structures at intermediate to high Galactic latitudes¹⁷.

In addition to the slope, the presence of any characteristic scale would provide useful information. One would expect a loss of power below certain scales if the turbulent energy is dissipated at small spatial scales (Ntormousi et al. 2016). Such an imprint of dissipation was not observed down to 0.01 pc in the analysis of Miville-Deschênes et al. (2016) using deep CFHT imaging of a high latitude region. Owing to the large pixel size of the Dragonfly, we are not capable of detecting such a break in the power spectrum at these small scales. However, the recently launched Euclid mission, with its unprecedented resolution of $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}1$ and superb sensitivity at low surface brightness, will shed new light on the turbulent nature of the ISM.

Fig. 9 Δ-variance spectra of different dust tracer maps. The “lag” (L) corresponds to the structure scale and ${\sigma }_{\mathrm{\Delta }}^2$ quantifies the variance of density fluctuations. WISE, Planck, and H I results are multiplied by a factor of 104, 5, and 5 × 10−3, respectively, for this visualization. A similar coherence, characterized by a power-law scaling ${\sigma }_{\mathrm{\Delta }}^2$ ∝ Lα, is revealed by the different dust tracers over a wide range of scales. A linear model is fit (dashed black lines) to find the slope, in the range α = 0.7–1. At very small and large scales, the coherence is altered by map-specific systematics.

5.2 Δ-variance analysis

5.2.1 An alternative view of coherence in cirrus structures across scales

Δ-variance is a variant of the power spectrum that is more robust at field edges and bad/missing pixels. We computed Δ-variance as a function of angular scale as described in Sect. 3.2.2 using Dragonfly optical, WISE MIR, Herschel FIR, Planck radiance, and DHIGLS N_H I data.

Figure 9 shows the Δ-variance spectra of the various cirrus maps. In general, Δ-variance increases with the structure scale over a wide range, largely following power laws. This is a hall-mark of the self-similarity of the density structures, indicative of the strong coherence between structures at different scales and the fractal nature of cirrus (Marchuk et al. 2021). Similar to power spectra, such power-law behavior is expected from interstellar turbulence in the ISM. The linear portions of the spectra were fit with a power law (Eq. (7)), as indicated by the dashed black lines. The fitted power-law slopes (α) are listed in Table 3, which ranges from 0.7 to 0.94. Recall that theoretically the power index of the power spectrum is related to the slope of the Δ-variance spectrum via γ = −α −2; hence, the structural coherence revealed by Δ-variance has indeed a quantitative likeness to that derived from the power spectrum analysis.

Fig. 10 Combined 1D cross-power spectra and power spectra of optical, MIR, and FIR cirrus maps. For visualization, the spectra are vertically offset, with the y-axis in arbitrary units. The colored ticks above the x-axis indicate twice the beam width of each map with the corresponding color. The median error bars are indicated in the upper right for each spectrum. At larger spatial scales (lower k), the cross-power spectra between optical and FIR or MIR follow a power law similar to the auto-power spectra, indicating strong structural correlation. At small scales comparable to or below the beam width, the power actually declines because the non-cirrus systematics are uncorrelated between the maps. Purple triangles show the cross-power spectrum between the CFIS-r image and the Herschel 250 μm image of a cirrus-rich field (“FLS”) measured by Lim et al. (2023) (see text).

5.2.2 What breaks the structural coherence at small and very large scales?

The Δ-variance spectra of various dust tracers follow a simple power law over a wide scale range, indicative of structural coherence across scales. However, they are affected by several factors at small and very large scales.

On small scales (≲1′), the optical, FIR, and MIR spectra show clear flattening. The flattening is more pronounced and occurs at slightly larger scales in FIR and MIR than in optical, likely due to the combined beam effects (including filtering and PSF blurring) and the presence of CIBA residuals in these datasets. In the simple case of white noise in a map, we would expect a spectrum with α = −2, and as the beam cuts off the real sky signals in the Δ-variance spectrum, a transition is made to this steep noise-related spectrum. The H I map also shows an up-bending at small scales comparable to its beam, because of the transition to the Δ-variance spectrum of the noise map from the cirrus-free end channels. In contrast, the Planck radiance map shows steepening on scales smaller than or comparable to the Planck beam; it is derived from thermal dust modeling and the reproducibility noise is very small.

On very large scales (≳$0_{\cdot }^{\circ }6$), Δ-variance spectra of molecular clouds often display a plateau or turnover, presumably caused by wind feedback (Ossenkopf et al. 2008). The spectra of the diffuse Spider region do not reveal any clear characteristic scale. However, the MIR and optical spectra exhibit a small degree of flattening compared to other tracers. The flattening in the MIR map might result from imperfect gradient correction during the stacking process used to subtract zodiacal light and faint off-axis moon-glow (Meisner & Finkbeiner 2014). Alternatively, it could be a result of WISE 12 μm tracing emission from ultrasmall, stochastically heated dust grains that are distributed slightly differently from large grains.

The flattening observed in the optical map cannot be explained by differences in dust compositions. But the flattening could stem from imperfect flat-fielding or subtle background mismatch in the mosaicking process. Boundary effects due to the finite size of the field could also influence the largest couple of scales (Bensch et al. 2001). There are other considerations with astrophysical origins. Optical cirrus can suffer from attenuation effects in the regime where it starts to transition from optically thin to optically thick (Mattila et al. 2023). Several studies have revealed nonlinearity in the optical-FIR correspondence in high column density regions (Román et al. 2020; Mattila et al. 2023; Zhang et al. 2023). However, it remains to be demonstrated that this would account for the flattening on large scales. This discrepancy between optical and FIR large-scale structures might relate to differences in the fractal dimension of optical cirrus in Stripe82 compared to its FIR counterpart reported by Marchuk et al. (2021). Finally, we cannot dismiss the possibility that astrophysical sources other than dust – in particular the cosmic optical background (COB) (e.g., Bernstein 2007; Zemcov et al. 2017; Postman et al. 2024) – also contribute faint emission on large scales. Further analysis with a larger dataset and refined data processing will be helpful in understanding the patterns observed in the optical Δ-variance spectra at these scales.

5.3 Cross-power spectra analysis

5.3.1 Structural correlation decoupled from systematics

A cross-power spectrum analysis enables the decoupling of map-specific systematics. We show the angular cross spectrum results computed following Sect. 3.2.3 using Dragonfly optical, WISE MIR, and Herschel FIR images. This cross-correlation analysis complements the comparisons of optical/FIR/MIR data based on local intensity statistics in Sect. 4.

We first plot the auto-power spectra for the optical, FIR, and MIR images as a consistency check; these spectra are in accordance with those shown in Fig. 8. The main difference here is that we show the gross measurements, i.e., we do not separate the ISM components (k^γ terms) by subtracting other components. At small scales (high k), below about three to four times the beam widths, the spectra deviate from a single power law due to CIBA residuals, beam effects, etc. The colored tick marks at the bottom of the plot denote twice the beam width for each map.

The 1D cross-power spectra between the optical image and the FIR or MIR images are shown as open markers in Fig. 10. These measurements can be well described by a power law with an index of γ = −2.87 and −2.95, respectively, close to that of the auto-power spectra, demonstrating a strong correlation between structures in the optical and FIR or MIR data. These correlations persist down to scales comparable to the FIR or MIR beam widths; below these scales (higher k), density fluctuations become de-correlated due to differing beam effects, pixel window functions, noise properties, and other map-specific systematics. As a result, the power index measured using this approach is less sensitive to systematics and decoupled from the degeneracies with other components specific to a single map (c.f., Eq. (5)).

Lim et al. (2023) conducted cross-correlation analyses between the Canada-France Imaging Survey (CFIS) optical images and the Herschel SPIRE images. In their work, the goal was to study the correlation between the COB and the CIB, where cirrus is a major contaminant. We plot their crosss-pectrum between the CFIS r-band image and the Herschel 250 μm map in the “FLS” field (their Fig. 13), which contains a substantial amount of cirrus (see their Fig. 5). Their cross-power spectrum is marginally shallower than ours, plausibly due to an inherent COB-CIB correlation. Note that the CIBA has a power-law dependence P(k) ∝ k^{−1−−1.5}, much shallower than the Galactic component (Viero et al. 2013), supported by their measurements in fields with minimal cirrus contamination. Therefore, it could be expected that our cross-power spectra focusing on the Galactic cirrus component, combined with efforts to mitigate extragalactic contributions by removal of point sources, yield a steeper power-law dependence.

Fig. 11 Correlation ratio between the optical and FIR (red)/MIR (blue) cirrus maps as a function of spatial frequencies. For both paired maps, the correlation approaches unity at large scales, indicating strong correlation of dust structures across wavelengths. At small scales, the correlation falls to zero because of uncorrelated contamination, noise, and map-specific systematics. Notably, the FIR-optical correlation exceeds the MIR-optical correlation over a broad range of angular scales, reflecting the greater similarity between dust morphology traced by thermal emission and that traced by scattered light.

5.3.2 Which is closer to optical cirrus globally: FIR or MIR?

In Figure 11, we show the correlation ratio ξ_a×b between the optical and FIR or MIR data computed using Eq. (9). For both paired maps, the correlation approaches unity at large spatial scales (low k), indicating the concordance of dust structures traced by different radiative mechanisms. At small scales, the correlation declines toward zero due to uncorrelated noise and systematics inherent to each map. It is worth noting that the FIR–optical correlation remains higher than the MIR–optical correlation across a wide range of scales. This result, from a global perspective, complements our findings in Sect. 4 at local scales (note that the smallest kernel scale used there is d_k = 2′) indicating that optical cirrus morphology more closely resembles that of FIR cirrus than MIR cirrus. This quantifies what was anticipated in Sect. 2.6.

6 Quantitative morphology of cirrus with wavelet scattering transforms

Local intensity statistics probe local non-Gaussianity but ignore multi-scale coupling, while the angular power-spectrum (and its variants) describe global coherence across scales but discard phase and non-Gaussian information. It is particularly useful to find a statistical approach that covers both strengths. With this motivation, we explore the potential of the wavelet scattering transform. One key hurdle is that the physical interpretation of the full set of WST coefficients is not straightforward. However, we take advantage of the results in Sects. 4 and 5 as a basis for extending our understanding, and use summary statistics that encapsulate the rich morphological information encoded in the coefficients.

Below we show as morphological measures the summary statistics computed from the second-order WST coefficients using optical and IR maps. The results for the Planck radiance map and the H I map are presented in Appendix C.1. We then explore the patterns of WST coefficients and summary statistics by dividing the field into subregions to study ensemble properties. Finally, as an application we show that mock cirrus can be generated using WST outputs from observations.

6.1 WST statistics as morphological measures

We start with the two summary statistics, sparsity s₂₁ and linearity s₂₂, defined by Eqs. (15) and (16) and derived from the WST coefficients computed from the Dragonfly optical, Herschel FIR, and WISE MIR maps using J = 10 and Θ = 8. These statistics are dimensionless and invariant under translation and rotation. Relative to the power spectrum and Δ-variance results in Sect. 5, the WST statistics incorporate additional information about the non-Gaussianity of the structures. Relative to the results based on local intensities in Sect. 4, they provide a more compact description of structure coherence across different scales.

Figure 12 illustrates the hierarchical patterns of s₂₂ (upper panel) and s₂₁ (lower panel) obtained from the three maps for different (j₁, j₂) combinations. Here j₁ represents the structure scale (indicated by different colors) and j₂ (plotted on the x-axis) represents the coherence scale that characterizes the assembly or clustering of j₁-scale structures separated by j₂-scale. We display only terms with j₂ > j₁ because the terms with j₂ ≤ j₁ are mainly determined by the properties of wavelets and thus do not carry significant physical information about the density field.

Across all scale combinations, the structures consistently exhibit s₂₂ > 1, indicating that the cirrus morphology is predominantly thin and filamentary. At a given coherence scale j₂, larger structures (i.e., those with higher j₁) generally exhibit higher linearity compared to smaller structures (with the exception of the largest scale, j₁ = 8). This suggests that cirrus morphology appears to be more elongated on larger scales, which can be interpreted as a consequence of the structure growth driven by the cascade of interstellar turbulence. Sparsity values are low, ranging from 0.01 to 0.12, indicating that cirrus is mostly diffuse and widespread across the field. At a given coherence scale j₂, larger structures exhibit higher sparsity (except for the largest scale). This suggests that large structures are more “localized,” i.e., occupying specific portions of the field.

Overall, both the FIR map and MIR map display trends highly consistent with the optical map, reinforcing the overall similarity in dust morphology. However, deviations from the optical data occur at specific scale combinations. For example, for both FIR and MIR maps, small-scale structures (j₁ ≤ 2) are slightly less filamentary at comparable coherence scales (j₂ = j₁ + 1), and are more evenly distributed (lower sparsity) at large coherence scales (j₂ ≥ 6). Additionally, for intermediate-large structures (j₁ = 5 and 6) the FIR data appear to have slightly higher sparsity on larger coherence scales, while the MIR data show more deviation to higher linearity than the optical at the largest coherence scales. Note that FIR and MIR data have larger beam widths (18″ and 15″, respectively, in Table 1). This would smooth structures with j₁ = 0 that correspond to 2^j₁+1 = 2 pixels, which would lower both linearity and sparsity. However, this factor alone could not explain the observed difference at larger scales, where the linearity is higher. We attribute these differences to residuals of extragalactic emission. For example, the CIBA can leave large-scale imprints across a wide range of structure scales (Auclair et al. 2024). The presence of COB and/or its distinction from CIB may also contribute to the differences in these hierarchical patterns between optical and IR. We note that in a different field, other diffuse extragalactic emission such as ICL and intragroup light may also contribute, if not well deblended from the cirrus foreground (see discussion in Sect. 7.2).

Comparisons with Planck and H I maps are provided in Appendix C.1. In general, the WST statistics computed from the Planck radiance map show trends similar to those of the optical map, while the H I map exhibits deviations. We speculate that the deviations may result from a combination of the differences in data processing and the possible decoupling between gas and large dust grains; further analysis with larger datasets is needed to determine whether these deviations persist and to assess their physical implications.

One application of this method is the evaluation of statistical similarity between two fields. To be considered realistic, mock cirrus images, for example produced by generative models, should exhibit WST statistics with patterns comparable to what is observed. Further applications are discussed in the following sections.

Fig. 12 Linearity s22 and sparsity s21 of optical, FIR, and MIR cirrus maps with different (j1, j2) combinations, which are summary statistics quantifying the non-Gaussianity of density fluctuations in the field. The structure scale is color-coded by j1. The coherence scale j2 characterizes the clustering on a scale j2 of j1-scale patterns. Only j2 > j1 terms are shown because j2 ≤ j1 terms are not physically significant. Values of s22 > 1 correspond to thin, filamentary structures, whereas values of s22 < 1 correspond to swirly, curvy structures. Lower s21 indicates that the structures are more widely spread, whereas s21 = 0 is the expectation for a Gaussian random field. In general, s22 and s21 of optical, FIR, and MIR cirrus maps have similar patterns except for small-scale structures.

6.2 A continuous representation of changing cirrus morphology

In this section, we explore the ensemble properties and variations of the cirrus morphology across the field. To achieve this, we divided the field into small subregions, each with a size of 45′ × 40′, as shown in Figure C.2 in Appendix C.2. This division respects a compromise between covering a sufficiently large area for structural analysis and obtaining a reasonable number of samples. We computed the WST coefficients for all subregions {I_k} (with k = 1 to 30) using J = 7 and Θ = 8.

Subsequently we calculated the sparsity s₂₁ and linearity s₂₂ for each subregion at different (j₁, j₂). The patterns are consistent with the patterns shown in Fig. 12. The statistics s₂₁ and s₂₂ calculated for each subregion also show correlations at several scale combinations (see below).

We further calculated the normalized statistics, denoted as ${\tilde{s}}_{21}$ and ${\tilde{s}}_{22}$, by averaging over all scale combinations with j₂ > j₁, with the average weighted by the inverse variance of the corresponding distribution among subregions. The left panel of Fig. 13 shows the distribution of ${\tilde{s}}_{21}$ and ${\tilde{s}}_{22}$ for all of the subregions, revealing a clear correlation between these statistics. The Pearson correlation coefficient R for ${\tilde{s}}_{21}$ versus ${\tilde{s}}_{22}$ is R = 0.81, indicating a strong correlation among subregions. We note that this scale-averaged correlation is mainly driven by specific scale combinations; in particular, the correlations are strong for small coherence scales (j₂ ≤ 2) or structures with comparable scales (j₂ − j₁ ≤ 3), but are weak or absent at larger scales (j₁ ≥ 4). For illustrative purposes, in the right panel of Fig. 13, we show the scatter plot of s₂₁ and s₂₂ for j₁ = 1, j₂ = 3, which has a Pearson correlation coefficient R of 0.89, indicating a strong correlation in play at this scale combination, stronger than for the averages plotted in the left panel. It can therefore be inferred that not all scales are relevant to describe the cirrus morphology with some scale combinations showing stronger correlations than others.

In Fig. 14, we display the subregions sorted by their positions on the left scatter plot of Fig. 13. Although the scatter plot layout is not strictly preserved, this representation effectively traces the correlation and visualizes the variation in cirrus morphology. In general, there is a discernible trend from the lower left to the upper right of the plot. This correlation between sparsity and linearity quantifies the visual impression that when cirrus structures are more concentrated, they tend to appear more filamentary, and conversely, when structures are more evenly distributed, they appear more diffuse. This morphological trend could be a result of the interplay between local ISM conditions and physical mechanisms. For example, Taank et al. (2022) found that CNM exhibits more small-scale structures than WNM using data from DHIGLS. Likewise, Lei & Clark (2023) found that the fraction of CNM in H I is correlated with small-scale linearity, indicating that CNM is more populated by H I filaments. Future studies with larger datasets and multiple tracers may further elucidate the connection between cirrus morphology and physical conditions in the ISM.

From another perspective, this representation from WST analysis provides a 2D characterization of the non-Gaussianity of each image, organizing various cirrus morphologies in a continuous and interpretable manner without requiring training or classification. Such representation is similar to the output of deep learning methods such as convolutional neural networks (CNNs). In fact, one interesting aspect of this approach is that it offers insights into the “blackbox” embedded in CNNs. As discussed in Cheng & Ménard (2021), many elements of the WST are analogous to, or simplification of, features characteristic of a CNN. In effect, the WST bridges the gap between traditional power spectrum methods and CNNs. Using predefined convolutional kernels, the WST can be viewed as a non-trainable shallow CNN (Cheng & Ménard 2021), which circumvents potential overfitting issues while simplifying the training process, architecture, and interpretable bottleneck associated with CNNs.

Fig. 13 Scatter plot of sparsity vs. linearity for subregions of cirrus in the example Spider field. Each marker represents a subregion within which the statistics have been evaluated. The left panel shows the distribution of reduced statistics, ${\tilde{s}}_{21}$ and ${\tilde{s}}_{22}$, after a weighted average over by all scale combinations with j2 > j1. The right panel shows s21 vs. s22 for one of the combinations with strong correlations (j1 = 1, j2 = 3). In both cases, the sparsity and linearity show strong correlations quantified by high Pearson correlation coefficients. Furthermore, some scale combinations are more relevant than others in describing the continuous changes in cirrus morphology.

Fig. 14 Illustration of sparsity vs. linearity for subregions in the example field showing the cutout images, averaged over by all scale combinations with j2 > j1. Note that the image cutouts are sorted by their positions on the left panel of Fig. 13, so the original layout does not strictly hold. In general, the patterns of the cirrus present a discernible trend varying with the two WST summary statistics, showing a regulation in cirrus morphology.

6.3 Application: statistical synthesis of mock cirrus

One interesting application is generation of mock cirrus fields. Such synthesis can be useful for testing detection and identification algorithms for extragalactic low surface brightness phenomena as well as providing estimates of bias on their measurements, where cirrus becomes one of the major sources of confusion.

We demonstrate using the same subregions of optical cirrus observations as shown in Fig. C.2, Appendix C.2. For computational efficiency, we binned the subregion [4x4] and cropped the binned image into the same dimensions in height and width. This leads to a [100 × 100] pix² image cutout corresponding to 40′ × 40′. To avoid artifacts caused by the FFT run on nonperiodic density fields, we folded the pixel values at the boundaries to make the image twice as large as the original cutout and become periodic at the boundaries. As described in detail in Appendix C.3, the image synthesis is performed to match the scattering covariance computed from WST run on these cutouts¹⁸. A synthesis can be performed with WST coefficients only, but we use scattering covariance so that cross-scale structure couplings are reproduced as well. In brief, the scattering covariance is obtained by doing a cross-correlation of WST map outputs (wavelet-convolved footprints), which further captures scale interactions across different wavelets. We adopt the calculation implemented in the scattering package (Cheng et al. 2020, 2024).

Here we demonstrate one-to-one synthesis by generating a mock cirrus field from a single cutout. An illustration of the image synthesis is presented in Figure 15, where the mock cirrus is displayed with the data from which it is generated side-by-side. We show three examples of realizations from top to bottom with different cirrus morphologies. The visual similarity between the mock and the data is striking. In addition, the right panels show the histograms of the intensity distributions of the data and the mock cirrus, indicating good correspondence in intensity PDFs. By construction, the non-Gaussian statistics of the field characterized by the WST match up to the second order, which can be extended to higher orders at the cost of computational time.

Although generative cirrus fields can be useful for image simulation tests, it should be noted that these synthesized images are purely phenomenological, i.e., no physical inference is involved and no new information is produced. The quality of the synthesis depends on several factors including the initial condition, the optimization, the parameterization, the data quality (e.g., surface brightness limits), and data processing effects (e.g., sky modeling, removal of foreground/background sources). Nevertheless, the fact that the mock cirrus fields look visually similar to the data strongly demonstrates the power of the WST approach in characterizing the cirrus morphology.

Fig. 15 Mock cirrus using scattering covariance computed from 40′ × 40′ subregions of the example dataset (from top to bottom: #13, #20, and #21 in Fig. C.2). The images are moderately filtered at high frequency to remove minor artifacts involved in the FFT process. The right panels show the histograms of the intensity distributions of target and mock cirrus.

7 Can extragalactic light be distinguished from cirrus via morphology?

We have demonstrated how a complementary suite of morphological diagnostics can robustly characterize the filamentary, non-Gaussian signatures of optical Galactic cirrus. An interesting application is whether these signatures can be used to distinguish extragalactic light from cirrus. Of the various phenomena, we choose tidal features for demonstration because cirrus is widely recognized as significant confusion to tidal features given their visually similar wispy morphology. If one of our statistical tools can reliably differentiate in this challenging regime, it will hold potential for other emissions such as galactic halos and ICL.

7.1 Identification of tidal features injected into cirrus-rich sky areas

Román et al. (2020) showed that it is possible to use color differences to distinguish extragalactic stellar populations from cirrus. It is interesting to explore the feasibility of distinguishing between tidal tails and cirrus based solely on their morphological characteristics. The underlying motivation lies in the different physics, tidal structures being molded by gravitational interactions, while cirrus structures are primarily shaped by turbulence, thermal instability, and (modulated by) magnetic fields. As tidal tails orbit the host galaxy, they tend to present different curvatures (e.g., arcs, loops, and pretzel-like patterns) compared to those manifested by cirrus, which can be distinguished, though subjectively, by trained human eyes (e.g., Martínez-Delgado et al. 2025). The following proof-of-concept experiment aims to identify the presence of tidal structures in cirrus-rich regions by leveraging differences in the statistical properties of the two types of features. Here we experiment with a single-band WST analysis using the g + r image, which does not require supplementary datasets.

We performed an injection test using realistic simulated tidal tail structures, which is illustrated in Fig. 16. The simulated tidal structures were generated by Miró-Carretero et al. (2025). Details about the generation of mock galaxies and the injection process are described in Appendix D.1. WST coefficients were computed for all subregions with tidal feature injections using the same scale indices (J = 7, Θ = 8). Instead of reducing the coefficients into summary statistics as in Sect. 3.3, we examined the full set of the second-order coefficients – after averaging over θ₁ (Eq. (14)) – yielding a total of J (J − 1) Θ/2 = 168 coefficients. We then employed dimensionality reduction techniques to extract the features that may distinguish subregions with tidal tails. The following investigations are motivated by Saydjari et al. 2021, who incorporated WST with dimensionality reduction to classify ISM structures in magneto-hydrodynamical simulations. Specifically, we applied principal component analysis (PCA) and linear discriminant analysis (LDA) on the WST coefficients.

We begin with PCA, a classic technique that finds the orthogonal components in a low-dimensional space that explain maximal variance. The input consists of the N = 30 subregions from the original cirrus field, each described by 168 variables¹⁹, which are scaled to have zero mean and unit variance. This WST-PCA analysis allows the WST coefficients to be described by the first few PCA components, with one, two, and three components accounting for 68%, 84%, and 90% of the variance, respectively, which indicates that the statistical properties of cirrus morphology can be effectively represented in a low-dimension space.

We then projected the WST coefficients of subregions with and without tidal tail injections onto those principal components. Figure 17 shows the projection onto the first two principal components, PC1 and PC2. Each data point represents a subregion – with or without tidal injection – labeled according to Fig. C.2. The two outliers (#25 and #30) correspond to regions with very little cirrus, which we designate as “cirrus-free.” In this WST-PCA space, subregions with tidal tails are shifted uniformly to the left in PC1 but up and down in PC2. This is not unexpected, because PCA does not aim to maximize separation between groups. The top histograms with respect to PC 1, which explains 68% of the variance, shows a distinct separation of peaks in an unsupervised manner. There is some overlap of the histograms (including #25 and #30 with tidal tails added). Increasing the surface brightness of the tidal features would shift the points further to the left, making them more distinguishable. One potential application, to search for tidal feature candidates entangled with cirrus in large surveys, is to examine the tail/outliers of this distribution in the WST-PCA space. In general, this supports the feasibility of using dimension reduction on WST coefficients to distinguish tidal features from cirrus, or at least to establish the presence of tidal features despite the contamination by cirrus.

We now turn to supervised learning to optimize the separation between patches with and without tidal injections. LDA is a classic machine learning technique that identifies hyperplanes in the N − d space that best separate the predefined classes. Unlike PCA, LDA components need not be orthogonal. Although more advanced techniques exist, we used LDA here for exploratory purposes, noting that future analyses may improve with alternative algorithms. The input data now consist of N = 60 subregions – with and without tidal tail injections²⁰ – each characterized by 168 WST coefficients. We split the input into training and test sets using five-fold cross-validation²¹. The input data were standard-scaled similarly as with WST-PCA. We combined multiple LDA classifiers across the folds through weighted voting based on accuracy to ensure robustness in the classification.

Figure 18 illustrates the WST-LDA results, which clearly distinguish the different groups. The accuracy and precision of the classification from cross-validation are 0.98 and 0.97, respectively. Note that this performance is based on a single injection case. We tested with different realizations of simulated tidal tails and obtained comparable results. Whether the performance is effective for a large number of tidal features remains to be tested, which requires more systematic investigation and is deferred to future work. We also tested the WST-LDA analysis using isotropic coefficients only (with J(J-1)=21 variables), which still yielded reasonably good performance with clear boundaries between categories. However, the ratio of intergroup and intragroup variance (i.e., the Calinski–Harabasz index) decreased remarkably. Therefore, retaining anisotropic information in the WST is preferable for the classification.

These results suggest that although cirrus and tidal features may appear visually quite similar, the statistical properties of their morphologies are subtly different. This distinction likely arises because the two types of structures are shaped by fundamentally different physical mechanisms. An intuitive explanation is that tidal features are more localized (as they are associated with the host galaxy) and become sparser on large scales, whereas cirrus appears more filamentary on large scales, because ISM physics operates over a broader range of scales than tidal stripping. We examined the patterns of WST statistics of pure tidal features and confirmed these trends (motivating this work; see Appendix D.2). Combining morphological measures with color information holds promise for identifying the presence of tidal features entangled with cirrus in deep wide-field surveys.

Fig. 16 Illustration of tidal tail injection into cirrus patches (see Appendix D.1). Top: (from left to right) (a) simulated galaxy from TNG 50 with tidal features (for DECam). (b) Same galaxy with satellites and density peaks removed, after convolution to Dragonfly PSF. (c) galaxy with tidal tails injected into the cirrus subregion. Bottom: (from left to right) (d) example subregion of cirrus (#22 in Fig. C.2). (e) tidal tails, excluding the central galaxy. (g) tidal tails injected into the cirrus subregion.

Fig. 17 Projection of the WST coefficients onto principal components PC1 and PC2 for subregions with (orange) and without (cyan) tidal tail injection. The top panel shows the histogram with respect to PC1. The PCA is trained with cirrus-only subregions. Even with only one component, the histograms for subregions with and without tidal tails show clustering with different peaks, although there is some overlap. This supports the feasibility of using dimension reduction on WST coefficients to distinguish tidal features from cirrus.

7.2 Separation of extragalactic light and DGL via morphology

Sect. 7.1 focuses on the detection challenge of extragalactic light in the presence of cirrus contamination. Component separation of extended emission for low surface brightness studies is a much more challenging task for advanced deep imaging surveys, but is essential for unbiased characterization of the target signal. Addressing this issue is beyond the scope of this paper. However, we outline possible morphology-based approaches for isolating the Galactic component from the integrated diffuse light.

Through training, a neural network can process subtle morphological cues that differentiate cirrus from extragalactic light, enabling an automated workflow scalable to large sky areas. However, one major obstacle is the lack of labeled training data. Smirnov et al. (2023) developed a deep learning framework for the segmentation of cirrus filaments using annotated SDSS Stripe82 data. Their method can greatly facilitate the cataloging of cirrus, which is very time-consuming using manual annotation. The separation task, however, requires pixel-level “ground truth” maps indicating the fractional contribution from each component, which is rarely the case for faint diffuse light such as cirrus or tidal tails where they are intertwined. In fact, optical imaging data of cirrus itself has been very limited. Wide-field imaging surveys optimized for low surface brightness covering the high Galactic latitude sky, such as the Dragonfly Ultra Wide Survey (DFUWS) and the Euclid Wide Survey, will provide valuable datasets for training. These can be combined with the synthesis method (see Sect. 6.3) for the purpose of data augmentation to mitigate overfitting.

An alternative approach is statistical modeling. Statistical modeling enables component separation in a robust, unsupervised manner (e.g., Siahkoohi et al. 2023; Lei et al. 2025). In this paper we present statistical characterizations of optical Galactic cirrus, and we experiment with those signatures to identify the presence of tidal features among cirrus in Sect. 7.1. Going one step forward, one could use them for component separation by assuming that the statistics of optical cirrus that is entangled with extragalactic light are identical to those without extragalactic light. A successful methodology can seen in the work by Auclair et al. (2024), who used WPH statistics to separate dust emission and CIB in FIR data. Similar algorithms using WPH or WST covariance (Appendix C.3) can be employed on optical images. Note that other diffuse light components, such as scattered light from the wide-angle PSF, should be carefully subtracted in advance. It is possible to also include COB in the modeling in this framework by using “flat-sky” images. The main advantage of such approach is that there is no need for labeled training data. Aside from greater interpretability, the methodology is also more robust than deep learning in terms of domain transfer (from survey to survey) or the risk of overfitting. The drawbacks, however, include the longer convergence time and the requirements in systematics control (i.e., pre-processing). Therefore, such methods are best applied on a case-by-case basis.

The fundamental basis for any separation arises from the difference in the morphologies of Galactic cirrus and extragalactic light, which are driven by distinct physical mechanisms. Figure 19 presents a schematic comparison of optical Galactic cirrus and two typical low surface brightness emissions of extragalactic origin, tidal features and ICL, summarizing their characteristic morphologies and the underlying mechanisms that shape them. We also indicate differences in the radiative sources. While the investigation in the present work has focused on morphology alone, different radiative processes would lead to distinct SEDs characterized by multiband photometry. Incorporating this complementary information with morphological signatures will further benefit the separation of cirrus and extragalactic light.

Fig. 18 LDA projection of the WST coefficients of subregions with and without tidal tail injections. Trained with labeled subregions, WST-LDA is able to find a good division in low-dimensional space between the categories. The good separation indicates that it is possible to identify the presence of tidal structures in cirrus-rich areas from their single-band morphology.

Fig. 19 Schematic summary of the typical morphology, underlying morphological drivers, and radiative mechanisms for optical Galactic cirrus, tidal features, and intracluster light. The differences in their morphologies and SED shapes together will provide critical information to break the confusion and separate cirrus and diffuse extragalactic components in the combined diffuse light.

8 Summary and prospects for application in deep wide-field imaging surveys

In this paper we applied multiple statistical approaches to characterize the morphology of optical Galactic cirrus. Although this morphological information is useful in itself, it is often most interesting to address the following:

• What astrophysical insights can cirrus morphology provide?

• Why does it matter to other low surface brightness science?

Here we summarize lessons learned from the characterization of cirrus morphology, and discuss prospects for its practical applications with modern deep wide-field imaging surveys in the context of both ISM studies and extragalactic investigations.

• (1)

  Local intensity statistics: PDF-based approaches have been widely applied to simulations and H I observations to study turbulence in the ISM. Moments and metrics such as skewness, kurtosis, and Gini can be used to probe where the morphology changes. However, few investigations have been extended to the optical regime. With next-generation deep imaging surveys, statistics of dust-scattered light derived from observations may be compared with those in state-of-the-art hydrodynamical ISM simulations at exceedingly high resolution to constrain dust prescriptions. Vice versa, if correlations between observables and ISM physical parameters (such as ℳ_S) exist and can be calibrated, they will provide a new avenue for inferring the latter.

  The greater similarity between optical and FIR cirrus shown in Sect. 4 is expected from their related origins involving larger dust particles, while MIR cirrus exhibits differences due to its distinct origin as nonequilibrium emission for smaller dust particles such as PAHs. Therefore, caution should be taken when using MIR as a dust surrogate to remove cirrus from optical images; FIR supplementary data would be preferred where available. Conversely, the PDF distance can be used as a metric to identify regions of strong PAH excitation versus thermal emission from cold dust.

  We also note that if extragalactic light is of interest at scales much larger than the Herschel beam width, such as in the case of ICL, it is possible to do a component separation by minimizing the local PDF distance of the optical cirrus to that of FIR over a range of spatial scales and across the sky. It is important that forefround/background systematics in both maps are removed, which could be demanding.

• (2)

  Fourier statistics: although power spectrum analysis, along with its variants, is a conventional tool, it has not been fully exploited over a wide sky area at visible wavelengths. The angular power spectrum of the optical cirrus in our example dataset follows a γ ≃ −2.9 power law, similar to other dust tracers, revealing coherent structures. Investigating the spatial variations of the cirrus power spectrum (and its variants) to study their dependence on the environment and ISM properties (e.g., dust temperature, CNM fraction), will be interesting.

  Second, the fact that the power spectra across dust tracers have consistent power indices indicates that the normalizations of power spectra can be used to calibrate the zero-point of the diffuse light background contributed by cirrus in optical images at the CCD chip level, with no need to dramatically downsample the optical data to, for example, the Planck beam width. Another useful application is that one can determine whether the cirrus signal is suppressed by sky background subtraction, and if so, at what scales. This is clear in the illustration shown in Appendix B.2, which presents the result for the same Spider sky area in the Legacy survey imaging.

  Furthermore, as discussed in Sect. 5.1, whether higher-resolution data would reveal a characteristic scale will be interesting. The superb resolution and sensitivity of Euclid and Roman will render them the ideal facilities for this kind of study in the future.

• (3)

  WST analysis: this approach enables the extraction of higher-order, non-Gaussian information that is missing in power spectrum analyses, thereby enhancing our understanding of the complex, multi-scale cirrus structures. As stated in Sect. 6, it is not yet well understood why some cirrus is highly filamentary and how its morphology correlates with ISM properties. Therefore, it will be interesting to explore the correlation between ISM characteristics and morphological signatures, for example, by extending the analysis in Lei & Clark (2023) to the optical regime.

  The WST statistics could also be useful when comparing the morphology in different observations and that in simulations, where co-spatial correspondence does not exist. Furthermore, mock cirrus can be generated using wavelet-based covariance, as shown in Sect. 6.3. This can be useful for estimating cirrus contamination of extragalactic signals, as well as algorithmic improvement of detection and measurement pipelines directed at targets of interest.

  In Section 7.1, we demonstrated the potential of using WST coefficients to distinguish extragalactic light from cirrus. This is particularly useful when color information is not available. For example, in the Euclid Wide Survey, cirrus in near-infrared (NIR) may not be comparably significant to cirrus in the optical due to a combination of the lower sensitivity in Euclid NIR bands (Cuillandre et al. 2025) and the fact that dust-scattered light is intrinsically fainter in the NIR in the absence of large grains (Zubko et al. 2004). Where colors are available, they will provide additional information for distinguishing extragalactic light from cirrus.

In summary, the morphology of Galactic cirrus in optical bands contains a wealth of information beyond measurements of intensity. Just as the morphology of a galaxy tells a tale of its physical properties and evolutionary history, the morphology of cirrus encodes insights about ISM properties and physical mechanisms that shape its structures, thereby rendering its morphology subtly different from extragalactic light. This information has been largely neglected. However, deep wide-field imaging surveys with low surface brightness optimizations should be capable of utilizing this information, advancing both ISM studies and extragalactic science.

Acknowledgements

Q.L. would like to thank Sihao Cheng, Wei Zhang, Yunning Zhao, Javier Roman, Jean-Charles Cuillandre, Nina Hatch, and Marc-Antoine Miville-Deschênes for helpful discussions. Q.L. is supported by the Oort Postdoctoral Fellowship. The research of R.G.A. and P.G.M. is supported by grants from the Natural Sciences and Engineering Research Council of Canada. The Dunlap Institute is funded through an endowment established by the David A. Dunlap family and the University of Toronto. This research has made use of data from the Planck Legacy Archive that provides online access to all official data products generated by the Planck mission. This research has made use of the APASS database, located at the AAVSO website. Funding for APASS has been provided by the Robert Martin Ayers Sciences Fund. The authors thank the excellent technical staff at the New Mexico Skies Observatory where Dragonfly is sited. This research has made use of data from the Herschel Space Observatory through the ESA Herschel Science Archive. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This research has made use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This research has made use of data from the Legacy Surveys. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 12120101003, # 11433005). The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO.

References

Appendix A Acquisition, reduction, and post-processing of optical imaging data

Appendix A.1 The Dragonfly telephoto array

Dragonfly is a distributed aperture telescope composed of 48 Canon 400 mm f /2.8 IS II USM-L telephoto lenses, together equivalent to a 1-m f/0.39 refractor. Each telephoto lens is equipped with a Santa Barbara Imaging Group (SBIG) CCD camera, with a field of view of $2_{\cdot }^{\circ }6\times 1_{\cdot }^{\circ }9$ and a pixel scale of $2_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}84$ pix⁻¹. The cameras take exposures with Sloan g- and r-filters. By design, Dragonfly is optimized for low surface brightness science. Readers are referred to Abraham & van Dokkum (2014) and Danieli et al. (2020) for detailed descriptions of the configuration of Dragonfly.

The general strategy of data acquisition and the data reduction workflow of Dragonfly is described in Danieli et al. (2020), summarized as follows. Briefly, Dragonfly takes 10-minute science exposures with the 48 lenses with large dithers and performs initial quality checks. The mean FWHM at the New Mexico Skies Observatory under good conditions is ≈ 5″ and the native pixel size is $2_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}85$. The exposures are reduced using the upgraded Dragonfly data reduction pipeline DFReduce (Bowman et al. in prep). The pipeline performs dark-subtraction, flat-correction (using twilight flats), astrometric solutions, photometric calibration, and quality checks on individual frames. Pixel area maps are used to correct the distortion by airmass difference given the large field-of-view. For cirrus-rich fields as used in this work, we adopt sky subtraction and combination techniques specifically developed for the preservation of the faint diffuse light from cirrus. This is critical, because improper background modeling would oversubtract cirrus, which leaves significant bias in the image co-add and dramatically alters the radiometric precision of the signal. Details about the principles and procedures are described in Liu et al. (2023). In brief, we adopt Planck thermal dust models (Planck Collaboration XI 2014) as priors for the background modeling of individual frames, and combine them with control in their background consistency. Demonstrations of improvement in the sky modeling can be found in Figure 2 and Appendix C of Liu et al. (2023).

Appendix A.2 Component separation of Galactic cirrus

We performed source modeling followed by a post-processing of the optical images to decompose the faint diffuse light from cirrus. The source models were built using the mrf package (van Dokkum et al. 2020) developed for modeling compact sources in low surface brightness imaging. In brief, mrf makes use of high-resolution imaging data and finds a matching kernel between the low-resolution image (here the Dragonfly data) and the high-resolution image. The flux models are built by convolving the high-resolution image with the kernel, which is then subtracted from the low-resolution image to leave out extended light in the image. We used imaging data from the DESI Legacy imaging survey DR9 (Dey et al. 2019) as the high-resolution source images. The scattered light in the extended PSF wings was incorporated into the modeling following the method described in Liu et al. (2022). Although an iterative process was implemented by removing the local cirrus signal, we caution that the presence of cirrus could still affect the determination of the PSF wing on large scales. The core parts of compact sources are masked and refilled using the algorithm in van Dokkum & Pasha (2024).

After the compact source removal, we applied a cirrus decomposition technique to the residual g and r images following the procedures described in Liu et al. (2025). The technique is based on two assumptions: (1) cirrus exhibits mostly filamentary or patchy morphology extending on scales larger than extragalactic sources, and (2) cirrus, which results from dust scattering in the Milky Way, is well constrained in its SED compared to that of extragalactic sources depending on the evolutionary tracks of stellar populations. This process involves a spatial filtering using a Rolling Hough Transform in combination with color constraints of cirrus calibrated using the Planck thermal dust model. This processing aims to remove any faint blob that does not agree with the cirrus SED (i.e., extragalactic sources), as well as residuals from compact source removal. Interested readers can refer to Liu et al. (2025) for further descriptions of the principles, procedures, and demonstrations.

It is noteworthy that the cirrus decomposition in this work adopted a simplified optical depth-dependent color model by using a piecewise linear model as described in Liu et al. (2025). In general, optically thick regions of cirrus appear redder than optically thin regions and this effect is nonlinear (Mattila et al. 2023). Although the piecewise linear model accounts for the overall reddening at high optical depths, it is a simplified model and does not account for the variation of optical depth at small scales. Therefore, the presumption of color constraint in the decomposition does not hold strictly in regions with potentially higher optical depth (see discussion in Liu et al. 2025), specifically, the central region of the Spider field (subregion #18, Fig. C.2). This could lead to a potential difference in the contrast of the cirrus structures at the center of the field relative to the rest of the regions, compared with other dust tracers. This is partially mitigated by using the g + r image to average out the effect, compared to using individual bands. Because of the limited area of potentially optically thick cirrus regions, we defer the investigation tackling the optical depth effect to future analysis with a larger dataset.

Appendix B Supplementary results of power spectrum analysis

Appendix B.1 Results of the HOTT radiance map

Thermal dust emission can also be modeled with modified blackbody parameterization using the high-quality Herschel PACS/SPIRE 160, 250, 350, and 500 μm images. This has been derived by the HOTT (Herschel Optimized Tau and Temperature) analysis (Singh & Martin 2022)²², which provides much higher-resolution (compared to Planck) dust optical depth, temperature, and radiance maps for a suite of fields including the Spider field. We made use of the HOTT maps for supplementary analysis. The HOTT radiance map has a beam width of 36″, which produces a larger dynamical range in spatial frequencies across this field compared to Planck.

The HOTT radiance map also include sources and noisy pixels in faint diffuse regions. Therefore, we performed a cleaning on the map to mitigate the impact from sources and the CIBA by masking sources detected above S/N threshold of 2 with mean fitted temperature below/above the 5%/95% quantile and pixels with χ² of fitting above 10. The criteria are empirical based on that the dust temperature is relatively uniform in these fields (Singh & Martin 2022) and emission other than thermal dust emission typically yields bad fitting. The masked map was convolved using a Gaussian kernel with a size of twice the HOTT beam width to fill in the masked pixels. In addition, we cropped out the map to avoid the diffuse noisy regions. This process mitigates the contribution of bright sources at large spatial frequencies, however, a more delicate component separation will be needed to fully remove those contaminations. The cropped and cleaned HOTT radiance map is displayed in Figure B.1, with the Planck radiance map shown in the background.

Fig. B.1 HOTT radiance map derived from Herschel data. The map is cropped to avoid the noisy regions and cleaned to remove bright sources. The Planck radiance map from thermal dust modeling is shown in the background.

We performed the power spectrum analysis on the HOTT radiance map following Sect. 5.1. We adopted Eq. 5 for the power spectrum modeling due to the presence of (faint) sources and CIBA, yielding γ = −2.9 and β = −1.1. The power spectrum of the HOTT radiance map is shown in Fig. B.2, with the contamination noise and instrumental components subtracted. The power index of the cirrus component is similar to that of Planck radiance and other dust tracers, while the contamination component has a similar power index to that of Herschel 250 μm FIR image, indicating that the same components exist consistently in Herschel data (Table 3). The uncertainties of the measured power spectrum are larger likely due to residuals in the source removal and CIBA. However, the HOTT map extends the power spectrum of the radiance to larger spatial frequencies than could be reached using Planck radiance given the large Planck beam.

Appendix B.2 Power spectrum measured in the DESI Legacy Imaging Survey

Here we show how sky background removal can oversubtract the cirrus signal in optical imaging and bias the structural statistics. We used the imaging data from the Data Release 10 of the DESI Legacy Imaging survey (Dey et al. 2019) obtained with the Dark Energy Camera (DECam). As pointed out by Zhao et al. (2024), it is necessary to carefully mask the foreground and background sources to separate their contribution from cirrus signals. Therefore, we first ran a local background subtraction to remove the diffuse light for compact source detection. All detected sources above signal-to-noise ratio of 2 were masked. In addition, saturated stars were masked using an empirical law according to their magnitudes, whose masks were manually enlarged after visual inspection to cover the extended PSF wings. The mask map was then applied to the original image and refilled following the recipe in Liu et al. (2025). The g and r images were then combined into a g+r image in a similar way as Dragonfly g+r and resampled to 6″ resolution. This reprocessed image is displayed in Figure B.3.

Fig. B.2 Combined 1D power spectra of different dust tracer maps as Fig 8 added with the result of HOTT radiance map (pink markers). Only the ISM components are shown. The normalization is in arbitrary unit. The power spectrum of the HOTT map follows the same power law as other dust tracers. The larger uncertainties at high frequencies are likely due to the CIBA.

It is clear by inspection that cirrus light in the Legacy survey image is suppressed. As discussed in Liu et al. (2023) (see their Figure 1 for a zoomed-in inspection), this is because the sky subtraction in the Legacy pipeline high-pass filters the image on spatial scales comparable to about half of the CCD chip size in DECam. This effectively eliminates the investigation of cirrus light on larger scales, in particular, its zero-point, and thus may bias the color measurements of cirrus. Such bias could partially explain the bluer mean color in the Legacy survey compared to other observations (Zhao et al. 2024). It is noteworthy that such bias could not be resolved by binning adjacent pixels.

In the context of this work, the structures of cirrus on large scales are significantly smeared out or removed. This can be confirmed by the measured power spectrum shown in Figure B.4, where its shape is flattened at low spatial frequencies by the sky background removal. We performed power spectrum modeling following Eq. 5 for k higher than where it flattens, yielding a power index of −2.6 for the cirrus signal. At even larger scales, the power spectrum rises again for reasons unknown. Our speculation is that some structural information remains or it could be some weird field edge effect.

For future analysis on cirrus using deep optical imaging surveys over a wide sky area, caution should be taken on the sky background modeling aiming for the preservation of the cirrus light (see further discussion in Liu et al. 2023). On the other hand, power spectrum analysis as a tool can be used to evaluate whether sky background removal has affected the cirrus structures. In this example, it indicates that cirrus light on scales above about 3′ in the Legacy imaging data is affected.

Fig. B.3 DESI Legacy g+r image of the field at 6″ resolution with sources removed. The intensity scale is linear and stretched to the 95% quantile to enhance the cirrus signal. Signals on large scales are largely removed by sky background subtraction.

Fig. B.4 1D power spectrum of the Legacy g+r image of the Spider field. The dashed blue line, dotted green line, and black line show the fitted cirrus component, the contamination noise component, and the combined model, respectively. The cirrus signal at large scales is significantly suppressed by sky background removal affecting the power spectrum at the low spatial frequencies in the shaded area.

Appendix C Supplementary results of WST analysis

Appendix C.1 Comparison of WST statistics with Planck and H I maps

In this section, we compare WST results of the optical map with Planck and H I maps. The Planck radiance map has a beam width of ~5′ with a pixel size of 1′, and the DHIGLS N_H I map has a beam width of ~1′ with a pixel size of 18″. To match the beam and pixel grid for comparison, we convolved the optical cirrus map to the corresponding beam widths of the two maps, and resampled them to the same grid. We then computed the summary statistics of the beam-matched, downsampled optical maps matching the Planck or H I map.

The left two panels of Figure C.1 show s₂₂ and s₂₁ of the Planck radiance map compared to the matched optical map. Because the first structure j₁ = 0 is much smaller than the beam, we only show terms with j₁ ≥ 1. However, we note j₁ = 1 terms also mildly suffer from smoothing by the large beam. In general, s₂₂ and s₂₁ of the Planck map follow trends similar to those of the optical map, suggesting that the dust morphology traced by radiance derived from FIR to sub-millimeter thermal emission is close to that traced by dust scattering. Optical cirrus appears to be more filamentary on larger coherence scales, and more evenly distributed (diffuse) for small-scale structures. The difference might be caused by the variation of the scattering phase function or dust temperature, however, further investigation on a much wider sky area will be needed to characterize the pattern differences and the possible causes.

The right two panels of Figure C.1 show the two summary statistics for the DHIGLS H I map and the matched optical map. Similarly, we do not plotj₁ = 0 terms due to the large beam width, and note that j₁ = 1 terms are comparable to the beam width. Values for the H I map shows clear deviation from the optical map at several scale combinations. At small scales (j₁ ≤ 2), the difference (especially in sparsity) can be caused by the complex noise properties of H I data (Blagrave et al. 2017), which is partially manifested in the complex noise template in the power spectrum that increases at high k. The reasons for the differences at large scales, however, remain in question. One possible explanation is that large dust grains traced by optical data become decoupled from gas flows and tend to be more concentrated than gas. This is qualitatively consistent with results from hydrodynamical simulations in denser regions (Hopkins & Lee 2016; Mattsson et al. 2019). Auclair et al. (2024) used a denoising algorithm on the H I velocity data cube prior to constructing the H I image to create a noise-free mock dust map. It would be useful to apply similar approaches in further analysis to investigate whether instrumental effects or the physical decoupling of dust and gas is the major cause of the observed difference in their morphologies as captured in the WST statistics.

Appendix C.2 Division of the field into subregions

Figure C.2 shows the division of the Spider field into subregions to study the ensemble properties of the WST coefficients. The two subregions, #25 and #30, correspond to the two rightmost points in Fig. 17 and the two “cirrus-free” subregions in Fig. 18.

Appendix C.3 Scattering covariance and ensemble synthesis of cirrus

The WST coefficients probe only the scale interactions within a given wavelet in Fourier space. A natural extension of the WST approach is to compute the cross-correlations of the different wavelet-convolved footprints to capture the coupling between wavelets, which is significant when the sequence is dyadic (Cheng & Ménard 2021). However, it is necessary to introduce nonlinearities in the operation for the characterization of interactions between different scales due to phase misalignment of the fluctuations (Allys et al. 2020). To compute this correlation, one can take the modulus of the footprint, and use a nonlinear Phase Harmonic operator to tune the phase information between footprints. This approach is called the Wavelet Phase Harmonic (WPH) approach. Details of the mathematical framework can be found in Allys et al. (2020). The WPH approach has been applied to ISM and cosmology (e.g., Allys et al. 2020; Regaldo-Saint Blancard et al. 2021; Mousset et al. 2024; Auclair et al. 2024). The cost of a richer description is a larger set of coefficients by one or two orders of magnitude than the WST coefficients.

Fig. C.1 WST statistics of Planck radiance map (left column) and DHIGLS NH I map (right column) compared to those computed using the optical map with matched beam widths. The x-axis describes the coherence scale j2. The structure scale is color-coded by j1. Only j2 > j1 ≥ 1 terms are shown because j1 = 0 structures are significantly smaller than the beams and j2 ≤ j1 terms are dependent on the wavelets. In general, s22 and s21 of Planck follow similar trends with those of optical, except for s22 at large coherence scales. However, H I shows deviation from the optical, which could be caused by different processing of H I data and possibly indicate the decoupling between gas and dust.

Fig. C.2 Subregion division to study the ensemble property of cirrus morphology. Each cutout has a size of 45′ × 40′. Subregion #18 might be affected by optical depth effect (see Appendix A.2). The two subregions at the upper right corner, #25 and #30, are marked as “cirrus-free” in the excise described in Section 7.1.

With similar objectives as the WPH approach, the WST covariance of an input field I down to the second order can be computed from the convolved footprints²³: $$\begin{array}{rl}& C_{01}^{j_1,j_2,{\theta }_1,{\theta }_2}=Cov(I\ast {\psi }^{j_2,{\theta }_2},|I\ast {\psi }^{j_1,{\theta }_1}|\ast {\psi }^{j_2,{\theta }_2}),\\ & C_{11}^{j_1,j_2,j_3,{\theta }_1,{\theta }_2,{\theta }_3}=Cov(|I\ast {\psi }^{j_1,{\theta }_1}|\ast {\psi }^{j_3,{\theta }_3},|I\ast {\psi }^{j_2,{\theta }_2}|\ast {\psi }^{j_3,{\theta }_3}).\end{array}$$(C.1)

The covariance depends on the amplitude of the power spectrum and so is normalized by the autocorrelation of the terms.

Figure C.3 illustrates how variations in WST covariance alter the morphology of structures through image synthesis based on a square cutout of the optical cirrus (shown in the middle) from the example dataset. The cutout is 25′ × 25′, binned [4x4], and is run with J = 5 and Θ = 4. The amplitudes of both the real and imaginary parts of C₀₁ and C₁₁ were multiplied by a scaling factor K increasing from left to right. The original image (not synthesized) is displayed in the middle panel. As the covariances increase, the cirrus morphology appears more filamentary and less diffuse. This corresponds to moving from the lower left to the upper right of Fig. 14.

The synthesis process aims to generate fields with similar texture to the input field(s) by reproducing the summary statistics (here the scattering covariances). For the synthesis, we employ the utilities of the scattering package (Cheng et al. 2020). This is done by using a gradient descent in pixel space to match the statistical estimators of the input field. The initial field is a random Gaussian field with the same power spectrum as the input image. The example of one-to-one synthesis has been displayed in Fig. 15, which was run with J = 6 and Θ = 4. on binned cutouts of 40′ × 40′ for computational efficiency. Such synthesis can be generalized to a multiple-to-multiple synthesis that aims to reproduce the average statistics of an ensemble of input images. Figure C.4 illustrates the ensemble synthesis for mock cirrus based on a set of subregions of the Spider data, showing the diversity of cirrus morphology that can be reproduced through the characterization by the WST approach.

Fig. C.3 Illustration of the morphological variation of cirrus structures by artificially increasing or decreasing scattering coefficients. The third panel is a 25′ × 25′ cutout of the optical cirrus map in the example dataset, which is used as the base image. The amplitudes of scattering covariance C01 and C11 are scaled by a factor K that is reduced in the left two panels and increased in the right two panels. From left to right, the sparsity and linearity of the field increase.

Fig. C.4 Illustration of ensemble synthesis using a set of cirrus subregions. The initial fields are random Gaussian fields following the same power spectra as the data in the left panel, which are the subregions in Fig. C.2 excluding faint regions #1, #11, #25, #30, and the optically thick region #18. The mock cirrus cutouts are generated by matching the averaged statistics of the input subregions.

Appendix D Injection of simulated tidal features into cirrus regions

Appendix D.1 Procedures of injecting tidal tail features

We inject realistic tidal features into the example cirrus field. The test galaxy with tidal tails was generated from the TNG 50 simulation of the IllustrisTNG project (Pillepich et al. 2019). Details about generating mock images from the simulation are referred to Miró-Carretero et al. (2025). In Miró-Carretero et al. (2025), the mock images are generated as appropriate for the DECam instrument at 70 Mpc. For illustrative purposes, below we describe the injection with one example mock (we have also performed experiments using a suite of mock galaxies and cirrus patches).

To match the optical data obtained by Dragonfly while covering a range of spatial scales of the typical cirrus structure characterized in this work, we convolved the mock DECam image to the Dragonfly PSF and resampled the image which equivalently placed the mock at a closer distance of 15 Mpc²⁴.

Because the central galaxy is much brighter than the tidal features and cirrus, it needs to be masked to avoid the power of the galaxy dominating over tidal structures. We masked the central galaxy using a combination of an aperture mask and surface brightness mask where pixels from the central galaxies with μ_r < 25.3 mag/arcsec² were masked. Through visual inspection, we further masked background sources and nearby faint satellite galaxies that have not been stripped or merged into the central galaxy. We then replaced the masked pixels using an iterative filling process. In brief, in each iteration the nearest N pixels in the unmasked portion of the image are refilled by convolving the previous image with a Gaussian kernel with increasing sizes until all the masked area is filled. This ensures the continuity of the infilled pixel values and smoothness in the center. The remaining faint tidal features have a mean surface brightness of ⟨μ_r⟩ = 26.2 mag/arcsec²²⁵. Their surface brightness is converted to the same unit of the cirrus map and then injected into subregions (patches) of the g + r intensity image following the division shown in Fig. C.2 and placed in the center of the patch.

Steps in the injection process are illustrated in Figure 16. The bottom right panel, which shows the injection after removing the central galaxy, is used as the input for the WST analysis.

In Sect. 7.1 the investigation was performed with this single mock. We have tested with several different mocks and obtained similar results. We noticed that some geometries and projections of features make them easier to distinguish than others. Since this is a proof-of-concept experiment, we defer a more comprehensive analysis to future work.

Appendix D.2 WST statistics of tidal tails

Here we show the comparison of the WST summary statistics derived from simulated tidal features with the results from optical Galactic cirrus. Fig. D.1 shows the linearity and sparsity of the tidal features in Fig. 16 (bottom middle panel) computed with J = 7 and Θ = 8. The host galaxy has been removed from that image. Both metrics exhibit trends that diverge markedly from those for optical cirrus. Tidal features appear more filamentary at small coherence scales, but their linearity consistently decreases as the scale increases. Conversely, sparsity rises substantially at larger scales. The difference trends can be ascribed to fundamentally different mechanisms governing the morphologies: tidal stripping is a localized gravitational process centering on the host galaxy, whereas turbulence and other physical processes in the ISM act over a wide range of angular scales.

Because these differences persist in the summary statistics, it can be inferred that they are also manifest in the original set of WST coefficients, which leads to the differentiation of the two phenomena in the WST-PCA and WST-LDA analyses in Sect. 7.1. Interestingly, such differences in morphology can be discerned by experienced observers. The quantitative descriptions introduced here now provide a rigorous basis for that perceptual knowledge.

Fig. D.1 WST summary statistics of the simulated tidal features (dashed lines) compared to results of optical cirrus shown in Fig. 12. From small to large coherence scales, s22 of tidal features consistently decreases, while s21 shows opposite trends. Note that the range of the y-axis of s21 here is larger than that in Fig. 12. The results here quantify the differences in the morphology of optical Galactic cirrus and tidal features, which lead to their differentiation shown in Figs. 17 and 18.

All Tables

All Figures

Fig. 1 Dragonfly RGB mosaic image of the Spider field at the top left and multiwavelength images from different dust tracers used in this work. Top middle: Dragonfly combination of g and r after source removal, approximating diffuse radiation in the V band. Top right: Herschel 250 μm. Bottom left: WISE 12 μm. Bottom middle: Planck radiance. Bottom right: DHIGLS H I LVC column density. See text for details. The images are displayed on a linear scale with the same contrast between 0 and the 99.99% quantile.

In the text

Fig. 2 Illustration of the extraction of local PDFs and PDF statistics. (a) Local PDF extracted from the intensity map within a circular region moving across the field. This example shows an extraction with radius of r = 3′ at the same position on the optical (left), FIR (middle), and MIR (right) map. The lower right panel shows the KDE-smoothed PDFs of the logarithm of the normalized intensity in the circular regions ($z=\log {\tilde{I}}_{\nu }$). The PDFs are used for computing statistics and distance metrics. DH indicates the distance between the PDF of the FIR/MIR data and that of the optical data within the particular subregion. (b) Skewness, kurtosis, and Gini coefficient of the local PDFs at the same position measured on the optical, FIR, and MIR cirrus data with varying kernel radius r. The rightmost panel shows the PDF distance of the FIR and MIR data relative to the optical as a function of kernel radius.

In the text

Fig. 3 Illustration of wavelet scattering transform applied to cirrus images, adapted from Figure 4 of Cheng & Ménard (2021). The input field is “scattered” through convolution with a bank of Morlet wavelets at different scales j and orientations θ (not displayed for clarity), followed by a nonlinear operation (modulus). The output fields are shown up to j2 = 3. The scattering coefficient Sn is computed by taking the spatial average of the output field. Only j2 > j1 coefficients are used.

In the text

	Fig. 4 Spatial distributions of local skewness, kurtosis, and Gini extracted from the Dragonfly optical (left column), Herschel FIR (middle column), and WISE MIR (right column) data, using a moving circular kernel with a kernel width of dk = 10′. These maps trace where cirrus morphology changes.
In the text

Fig. 5 Left column: histograms of skewness, kurtosis, and Gini measured from the local PDFs in optical, FIR, and MIR data. Fainter histograms represent smaller kernel scales used for PDF extraction. Middle column: histograms showing the distributions of differences of the statistics between FIR or MIR and optical data. The difference is calculated by [optical-X], where X stands for the band. Right column: mean absolute difference of statistics between FIR or MIR and optical as a function of kernel scale.

In the text

Fig. 6 Left: histograms of the Hellinger distance (DH) between FIR or MIR and optical maps, which indicates the proximity of two cirrus maps. Smaller DH represents a closer local PDF. Histograms with lighter shades represent smaller kernel scales used for the PDF extraction. The PDF distance of FIR relative to optical is systematically lower than that of MIR. The dispersion is also much smaller. Right: mean DH between FIR (orange) or MIR (green) and optical data as a function of kernel scale. The mean DH of FIR relative to optical is lower than MIR in all scales. The dashed lines show the asymptotic values at large scales.

In the text

Fig. 7 Left: maps of sonic Mach number inferred from optical cirrus map illustrated with a kernel scale of dk = 10′. The circular kernel is illustrated with the purple disk. Right: mean H I column density fields extracted with the same kernel in the upper panels. Regions with ℳs < 0.5 and > 1.5 in the upper panels are indicated by red and cyan contours, respectively, which aggregate at edges of cirrus filaments or patches.

In the text

Fig. 8 Combined 1D power spectra of different dust/gas tracer maps. We only show the ISM components of the power spectra, i.e., the contaminating sky noise and instrumental noise components are subtracted. The power spectra from different maps are similar, with a power index of γ ranging from −2.8 to −3.0 (Table 3). Each power spectrum is displayed for a conservative range at higher k, excluding where the contaminating noise is most prominent and/or the cutoff by the beam is significant. The normalization is in arbitrary units, with the power-laws being aligned in the range where they overlap at lower k. The power spectra are overplotted in order of increasing beam size (except FIR and MIR, guided by Fig. 10), so that it can be appreciated that the spectrum for the Dragonfly optical (g+r) map extends to the highest k. The dashed blue line shows the result of Miville-Deschênes et al. (2016) with P(k) ∝ k−2.9.

In the text

Fig. 9 Δ-variance spectra of different dust tracer maps. The “lag” (L) corresponds to the structure scale and ${\sigma }_{\mathrm{\Delta }}^2$ quantifies the variance of density fluctuations. WISE, Planck, and H I results are multiplied by a factor of 104, 5, and 5 × 10−3, respectively, for this visualization. A similar coherence, characterized by a power-law scaling ${\sigma }_{\mathrm{\Delta }}^2$ ∝ Lα, is revealed by the different dust tracers over a wide range of scales. A linear model is fit (dashed black lines) to find the slope, in the range α = 0.7–1. At very small and large scales, the coherence is altered by map-specific systematics.

In the text

Fig. 10 Combined 1D cross-power spectra and power spectra of optical, MIR, and FIR cirrus maps. For visualization, the spectra are vertically offset, with the y-axis in arbitrary units. The colored ticks above the x-axis indicate twice the beam width of each map with the corresponding color. The median error bars are indicated in the upper right for each spectrum. At larger spatial scales (lower k), the cross-power spectra between optical and FIR or MIR follow a power law similar to the auto-power spectra, indicating strong structural correlation. At small scales comparable to or below the beam width, the power actually declines because the non-cirrus systematics are uncorrelated between the maps. Purple triangles show the cross-power spectrum between the CFIS-r image and the Herschel 250 μm image of a cirrus-rich field (“FLS”) measured by Lim et al. (2023) (see text).

In the text

Fig. 11 Correlation ratio between the optical and FIR (red)/MIR (blue) cirrus maps as a function of spatial frequencies. For both paired maps, the correlation approaches unity at large scales, indicating strong correlation of dust structures across wavelengths. At small scales, the correlation falls to zero because of uncorrelated contamination, noise, and map-specific systematics. Notably, the FIR-optical correlation exceeds the MIR-optical correlation over a broad range of angular scales, reflecting the greater similarity between dust morphology traced by thermal emission and that traced by scattered light.

In the text

Fig. 12 Linearity s22 and sparsity s21 of optical, FIR, and MIR cirrus maps with different (j1, j2) combinations, which are summary statistics quantifying the non-Gaussianity of density fluctuations in the field. The structure scale is color-coded by j1. The coherence scale j2 characterizes the clustering on a scale j2 of j1-scale patterns. Only j2 > j1 terms are shown because j2 ≤ j1 terms are not physically significant. Values of s22 > 1 correspond to thin, filamentary structures, whereas values of s22 < 1 correspond to swirly, curvy structures. Lower s21 indicates that the structures are more widely spread, whereas s21 = 0 is the expectation for a Gaussian random field. In general, s22 and s21 of optical, FIR, and MIR cirrus maps have similar patterns except for small-scale structures.

In the text

Fig. 13 Scatter plot of sparsity vs. linearity for subregions of cirrus in the example Spider field. Each marker represents a subregion within which the statistics have been evaluated. The left panel shows the distribution of reduced statistics, ${\tilde{s}}_{21}$ and ${\tilde{s}}_{22}$, after a weighted average over by all scale combinations with j2 > j1. The right panel shows s21 vs. s22 for one of the combinations with strong correlations (j1 = 1, j2 = 3). In both cases, the sparsity and linearity show strong correlations quantified by high Pearson correlation coefficients. Furthermore, some scale combinations are more relevant than others in describing the continuous changes in cirrus morphology.

In the text

Fig. 14 Illustration of sparsity vs. linearity for subregions in the example field showing the cutout images, averaged over by all scale combinations with j2 > j1. Note that the image cutouts are sorted by their positions on the left panel of Fig. 13, so the original layout does not strictly hold. In general, the patterns of the cirrus present a discernible trend varying with the two WST summary statistics, showing a regulation in cirrus morphology.

In the text

	Fig. 15 Mock cirrus using scattering covariance computed from 40′ × 40′ subregions of the example dataset (from top to bottom: #13, #20, and #21 in Fig. C.2). The images are moderately filtered at high frequency to remove minor artifacts involved in the FFT process. The right panels show the histograms of the intensity distributions of target and mock cirrus.
In the text

Fig. 16 Illustration of tidal tail injection into cirrus patches (see Appendix D.1). Top: (from left to right) (a) simulated galaxy from TNG 50 with tidal features (for DECam). (b) Same galaxy with satellites and density peaks removed, after convolution to Dragonfly PSF. (c) galaxy with tidal tails injected into the cirrus subregion. Bottom: (from left to right) (d) example subregion of cirrus (#22 in Fig. C.2). (e) tidal tails, excluding the central galaxy. (g) tidal tails injected into the cirrus subregion.

In the text

Fig. 17 Projection of the WST coefficients onto principal components PC1 and PC2 for subregions with (orange) and without (cyan) tidal tail injection. The top panel shows the histogram with respect to PC1. The PCA is trained with cirrus-only subregions. Even with only one component, the histograms for subregions with and without tidal tails show clustering with different peaks, although there is some overlap. This supports the feasibility of using dimension reduction on WST coefficients to distinguish tidal features from cirrus.

In the text

	Fig. 18 LDA projection of the WST coefficients of subregions with and without tidal tail injections. Trained with labeled subregions, WST-LDA is able to find a good division in low-dimensional space between the categories. The good separation indicates that it is possible to identify the presence of tidal structures in cirrus-rich areas from their single-band morphology.
In the text

	Fig. 19 Schematic summary of the typical morphology, underlying morphological drivers, and radiative mechanisms for optical Galactic cirrus, tidal features, and intracluster light. The differences in their morphologies and SED shapes together will provide critical information to break the confusion and separate cirrus and diffuse extragalactic components in the combined diffuse light.
In the text

	Fig. B.1 HOTT radiance map derived from Herschel data. The map is cropped to avoid the noisy regions and cleaned to remove bright sources. The Planck radiance map from thermal dust modeling is shown in the background.
In the text

	Fig. B.2 Combined 1D power spectra of different dust tracer maps as Fig 8 added with the result of HOTT radiance map (pink markers). Only the ISM components are shown. The normalization is in arbitrary unit. The power spectrum of the HOTT map follows the same power law as other dust tracers. The larger uncertainties at high frequencies are likely due to the CIBA.
In the text

	Fig. B.3 DESI Legacy g+r image of the field at 6″ resolution with sources removed. The intensity scale is linear and stretched to the 95% quantile to enhance the cirrus signal. Signals on large scales are largely removed by sky background subtraction.
In the text

Fig. B.4 1D power spectrum of the Legacy g+r image of the Spider field. The dashed blue line, dotted green line, and black line show the fitted cirrus component, the contamination noise component, and the combined model, respectively. The cirrus signal at large scales is significantly suppressed by sky background removal affecting the power spectrum at the low spatial frequencies in the shaded area.

In the text

Fig. C.1 WST statistics of Planck radiance map (left column) and DHIGLS NH I map (right column) compared to those computed using the optical map with matched beam widths. The x-axis describes the coherence scale j2. The structure scale is color-coded by j1. Only j2 > j1 ≥ 1 terms are shown because j1 = 0 structures are significantly smaller than the beams and j2 ≤ j1 terms are dependent on the wavelets. In general, s22 and s21 of Planck follow similar trends with those of optical, except for s22 at large coherence scales. However, H I shows deviation from the optical, which could be caused by different processing of H I data and possibly indicate the decoupling between gas and dust.

In the text

	Fig. C.2 Subregion division to study the ensemble property of cirrus morphology. Each cutout has a size of 45′ × 40′. Subregion #18 might be affected by optical depth effect (see Appendix A.2). The two subregions at the upper right corner, #25 and #30, are marked as “cirrus-free” in the excise described in Section 7.1.
In the text

Fig. C.3 Illustration of the morphological variation of cirrus structures by artificially increasing or decreasing scattering coefficients. The third panel is a 25′ × 25′ cutout of the optical cirrus map in the example dataset, which is used as the base image. The amplitudes of scattering covariance C01 and C11 are scaled by a factor K that is reduced in the left two panels and increased in the right two panels. From left to right, the sparsity and linearity of the field increase.

In the text

Fig. C.4 Illustration of ensemble synthesis using a set of cirrus subregions. The initial fields are random Gaussian fields following the same power spectra as the data in the left panel, which are the subregions in Fig. C.2 excluding faint regions #1, #11, #25, #30, and the optically thick region #18. The mock cirrus cutouts are generated by matching the averaged statistics of the input subregions.

In the text

Fig. D.1 WST summary statistics of the simulated tidal features (dashed lines) compared to results of optical cirrus shown in Fig. 12. From small to large coherence scales, s22 of tidal features consistently decreases, while s21 shows opposite trends. Note that the range of the y-axis of s21 here is larger than that in Fig. 12. The results here quantify the differences in the morphology of optical Galactic cirrus and tidal features, which lead to their differentiation shown in Figs. 17 and 18.

In the text