Classification of Magnetohydrodynamic Simulations Using Wavelet Scattering Transforms

2.1.1. Three-point Correlation Function (3PCF)

2.1.2. 2D Wavelet Scattering Transform (2D-WST)

In 2D, a WST is specified by (J, L, M), where J is the maximum scale, L is the number of angular bins (i.e., the angular spacing is 180°/L), and M is the maximum order of the transform. The maximum scale and number of angular bins (J, L) define the set of Morlet wavelets ψ_j,l used in the transform. ¹⁴ The parent Morlet wavelet is the product of a plane wave of unit wavenumber with a Gaussian window of width σ. In one dimension, the wavelet is given by

where α and β are normalization factors (Ashmead 2010). We choose α so that the wavelets cover the Fourier domain as uniformly as possible. We choose β such that ψ has a null average.

In 2D and above, the inverse variance, 1/σ², of the Gaussian window becomes an inverse covariance matrix. We choose a diagonal covariance matrix, with σ_yy = 2σ_xx (where x is the direction of the plane wave), to enhance the directionality of the wavelets. We choose σ_xx = 0.8 for our 2D-WST, which limits the number of oscillations of the plane wave, thereby enhancing filament identification. Each wavelet used in computing the WST coefficients is said to represent an oriented scale, indexed by (j, l) ∈ ([0, …, J − 1], [0, …, L − 1]), and is obtained from the parent wavelet by rotation and scaling of the argument (by 2^j ) and frequency (by 2^−j ). ¹⁵ While we use l as an index, we denote the corresponding angles as θ ∈ [π/2 − π/L,..., − (L − 1)π/2L, − π/2]. The wavelet, ψ_j,l , has a bandwidth of ∼2^−j , and is centered at a wavevector of magnitude ∼2^−j and angle θ, corresponding to l in the Fourier domain (Figure 1).

Zoom In Zoom Out Reset image size

The WST can be represented as a convolutional network, with layers labeled m from 0 to M. We will outline the computation of coefficients for each layer for a network where M = 2. Coefficients for the mth layer are derived from m convolutions of the image field I( x ), with Morlet wavelets under the usual modulus for complex numbers. Normalization of the coefficients is relative to the response to a Dirac delta distribution, δ_2D( x ):

Here, is the indicator function, which is 1 everywhere in the image. For μ₀ and μ₁, the convolution simplifies to the area of the image, and the sum of over the image, respectively. ¹⁶ Note that first-order coefficients depend on one oriented scale (j₁, l₁), while second-order coefficients depend on two oriented scales, (j₁, l₁) and (j₂, l₂). Since the first convolution obscures information at scales below pixels, only coefficients with j₂ > j₁ are computed, such that one is then only looking at scales . We have one coefficient for m = 0, JL coefficients for m = 1, and J(J − 1)L²/2 coefficients for m = 2, for a total of J(J − 1)L²/2 + JL + 1. The WST coefficients at the mth layer depend on correlation functions up to order 2^m (proof in Mallat 2012). Thus, M > 1 is required to capture non-Gaussianity, because for m = 1, the WST coefficients depend predominantly on the 2PCF. For many applications, the spectral energy falls off rapidly with m; as such, M = 2 provides sufficient sensitivity to non-Gaussianity with a tractable number of coefficients.

We normalize the 2D-WST coefficients following Allys et al. (2019), and Bruna et al. (2015), as

This normalization makes the first- and second-order coefficients invariant under multiplication of the field by a constant factor, and makes the second-order coefficients invariant to modification of the spectrum of the field by the action of a class of linear filters (Bruna & Mallat 2018).

2.1.3. 3D Wavelet Scattering Transform (3D-WST)

In 3D, the WST coefficients are defined similarly, except that solid harmonic wavelets (denoted by ) are used in place of Morlet wavelets, with m ∈ [ − l,…,l − 1, l].

Here, Y_l ^m are the Laplacian spherical harmonics, with (j, l) ∈ ([0,…,J − 1, J], [0, …, L − 1, L]). This corresponds to the choice of σ = 1 for the Gaussian window; the j-dependent wavelets are obtained by scaling the argument, as for 2D. We follow the implementation in kymatio ¹⁷ (Eickenberg et al. 2018), averaging over m in computing both the first- and second-order coefficients. For first-order coefficients, this is explicitly

where q = 1. For the 3D-WST, S₂(j₁, j₂, l₁) is a function of one angular scale only, we take l₂ = l₁. This means that the second-order coefficients are obtained between scales at the same angles only (rather than the Cartesian product of all angles). We compute the convolution coefficients for q = 1/2 and q = 2, in addition to q = 1, which is the exclusive power used in 2D-WST, although we comment only briefly on the q ≠ 1 cases. ¹⁸ There are (L + 1)(J + 1) coefficients for m = 1, and (L + 1)J(J + 1)/2 coefficients for m = 2, for a total of (L + 1)(J² + 3J + 2)/2 (for each power q computed).

2.1.4. Reduced Wavelet Scattering Transform (RWST)

The RWST reduction (which is only defined in 2D) takes advantage of the periodicity observed in the WST coefficients to remove the angular dependence (see Appendix C). Explicitly, the RWST coefficients are obtained from the WST coefficients by a least-squares fit of the first- and second-order coefficients. For the first order,

where , , and are fit coefficients. For the second order,

where , , , , and are fit coefficients. There are 3J coefficients for m = 1, and 5J(J − 1)/2 coefficients for m = 2, for a total of J(5J + 1)/2. ¹⁹

In order to interpret these high-dimensional (>100 coefficients), non-Gaussian descriptors, we employ dimensional reduction techniques, including principal component analysis (PCA), linear discriminant analysis (LDA), and quadratic discriminant analysis (QDA). ²⁰ PCA is an unsupervised technique, in which the ith principal component is orthogonal to the first i − 1 principal components, and is in the direction of maximal remaining variance (Jolliffe & Cadima 2016). In contrast, LDA and QDA are supervised techniques, meaning that they take class labels (here, hydrodynamic parameters) in addition to the high-dimensional descriptor of a data point (here, the WST, RWST, or 3PCF coefficients). LDA assumes that each component of the high-dimensional descriptor is normally distributed, and that all classes share the same covariance matrix; this latter assumption is known as homoscedasticity (Bandos et al. 2009). LDA then finds the hyperplanes that best separate the K classes in the K − 1 dimensional space containing the class means. The hyperplanes are found using eigenvalue decomposition, involving the inter- and intra-class scatter matrices. This is equivalent to evaluating the Gaussian likelihood that each input is in each class, and assigning the input to the class with the greatest likelihood. Note that the LDA components (which maximize the separation between classes) need not be orthogonal. The LDA components are in units of σ_ii , the on-diagonal components of the shared covariance matrix. QDA is similar to LDA, except that it drops the assumption that the covariance matrix must be the same for all classes, and thus finds conic sections that best separate the classes, rather than hyperplanes. The boundary between two neighboring classes occurs where the Gaussian likelihoods of each class are equal, which is in general a quadratic function (as in QDA), but simplifies to a linear function when the class covariances are assumed equal (as in LDA).

Interstellar dust is an indirect tracer of interstellar turbulence, and provides an observational point of comparison to the MHD simulations discussed below. Dust grains are heated to ∼20K by ambient starlight, and produce thermal emissions in the far-infrared (Low et al. 1984). We use the SFD dust map (Schlegel et al. 1998), which is based on this thermal emission. The map derives dust morphology from 100 μm data from the Infrared Astronomical Satellite (IRAS) Sky Survey Atlas (Neugebauer et al. 1984; Wheelock et al. 1994), corrected for temperature variation using data from the Diffuse InfraRed Background Experiment (DIRBE) on the Cosmic Background Explorer (COBE) (Mather 1982; Boggess et al. 1992).

To obtain a representative sample of 2D dust images, we reproject the SFD dust map with bilinear interpolation to a TAN (gnomonic) projection, centered on random locations on the celestial sphere. Each 256² pixel image has a random position angle, a pixel scale of 30, and a point-spread function (PSF) full width at half maximum (FWHM) of 61. The SFD bitmask flags regions with missing IRAS data (bit 128), pixels in the Magellanic clouds or M31 (bit 64), and regions with no point-source subtraction (mostly at low Galactic latitude, ∣b∣ < 5°; bit 32). An image is discarded if any pixel has any of these mask bits set, leaving a sample of 10,222 images (Figure 2). ²¹ By rejecting images near the Galactic plane, we focus on the morphology of relatively nearby dust (on the order a few 100 pc to a few kpc), where each image is dominated by dust at one, or a few, distances. The statistics of the dust distribution near b = 0 might be substantially different, owing to the large number of clouds along each line of sight, and the steep dependence on latitude. The images are normalized to have [min,max] = [0,1] for input to the WST, and zero mean with unit variance for input to the 3PCF, consistent with other images used in the paper.

Zoom In Zoom Out Reset image size

We use MHD simulations of turbulent gas flows published in previous works (Cho & Lazarian 2003; Burkhart et al. 2009; Portillo et al. 2018; Bialy & Burkhart 2020; Burkhart et al. 2020), and review only key properties. The simulations solve the ideal MHD equations constraining the divergence to zero (∇ · B = 0), assuming an isothermal equation of state, and periodic boundary conditions. No gravity is included, and the simulations are “scale-free.” The simulation box measures 256³ voxels, and is driven isotropically by random, divergence-free forcing, with wavenumber k ≈ 2.5 (i.e., 1/2.5 the box size). No dissipation is assumed, other than numerical dissipation from the finite resolution of the simulation. The magnetic field is initialized uniformly in the x-direction.

Each simulation is characterized by its Alfvénic Mach number (M_A), and sonic Mach number (M_S). The Alfvénic Mach number is defined as M_A ≡ ∣v∣/〈v_A〉, where v is the flow velocity, and 〈v_A〉 the Alfvén (magnetic wave) speed averaged over the simulation box. The sonic Mach number is defined as M_S ≡ ∣v∣/c_s , where c_s is the isothermal speed of sound. We use M_A = 0.7, 2.0, and M_S = 0.7, 1.0, 2.0, 4.0 to include sub-Alfvénic (M_A = 0.7) and super-Alfvénic (M_A = 2.0) cases, as well as subsonic (M_S = 0.7), transonic (M_S = 1.0), and supersonic (M_S = 2.0, 4.0) cases (Figure 3). The images are normalized to have [min,max] = [0,1] for compatibility with the log normalization of the WST, but we have confirmed that our WST classification result is reproduced if one uses zero mean and unit variance normalization. The images are normalized to have zero mean and unit variance when computing the 3PCF, in order to exclude any contributions from the 2PCF (see Appendix B).

Zoom In Zoom Out Reset image size

We assess the performance of the WST, its reduction the RWST, and the 3PCF as high-dimensional descriptors for the characterization of non-Gaussianity in 2D images. Each of these descriptors is then used in conjunction with two dimension-reduction techniques, PCA, and LDA. PCA and LDA are used to probe the information contained in the high-dimensional parameter spaces, corresponding to the WST, RWST, and 3PCF coefficients. We begin by confirming that these descriptors are sensitive to non-Gaussianity by applying them to 2D dust map images. We compute the WST coefficients (J = 8, L = 8, M = 2) on 10,222 SFD images (“SFD”), ²² 2000 images where each pixel is a uniformly drawn pseudo-random number [0,1) (“Rand”), and those same random images, convolved with a σ = 4 pixel uniform Gaussian (“Smooth”). This latter set imposes nearest-neighbor correlations, which are present in empirical images due to the continuity of real density fields and instrumental point-spread functions. To further prove that the WST captures non-Gaussianity, we compute the WST coefficients for images where the Fourier power spectra matches the SFD images, but with small higher-order correlations (“NHC”). Specifically, for each image in the SFD set, we apodize the image (i.e., multiply by a generalized Gaussian, ²³ p = 6, σ = 100), compute the Fourier transform (FT), randomize the phase (with a unique pseudo-random seed), and compute the inverse FT. ²⁴ The RWST and 3PCF (N_bins, L) = (8, 9) coefficients are also computed on the same four data sets for comparison.

To investigate the variance of these coefficients with respect to prototypical 2D dust, we perform PCA ²⁵ on the coefficients from the SFD set and then project the results from the other three control image sets (“Rand,” “Smooth,” “NHC”) onto those principal components. Each point in Figure 4 corresponds to an image from one of the data sets above.

Zoom In Zoom Out Reset image size

In WST-PCA space, the SFD data set presents as an ellipsoid, with a population of outliers extending along the first PC axis in the positive direction. The “NHC,” “Smooth,” and “Rand” data sets present as progressively narrower ellipsoids overlapping the main mass of “SFD.” While slight, the different orientation of these ellipsoids may encode more detailed information differentiating the data sets. We observe that, to first order, the diameter of the ellipsoids decreases with the progressively decreasing pixel – pixel correlations in the “SFD” > “NHC” > “Smooth” > “Rand” series. This suggests that the variance in the WST coefficients on “SFD” images primarily arises from higher-order correlations, as desired. As shown in the cumulative variance plot, most of the variance of “SFD” is described by the first few PCA components, suggesting that WST-PCA admits a low-dimensional representation of SFD dust. By cross-referencing the Galactic coordinates of the images associated with the outliers in the “SFD” data set, using glue, ²⁶ we associated most of the outliers with images of the Orion nebula. The formation of massive stars in the Orion nebula modifies the dust distribution, and is likely to give rise to outliers in the WST-PCA space. This was also true for the outliers in the RWST-PCA and 3PCF-PCA bases. The sensitivity of WST-PCA to outliers may suggest its potential application in the area of outlier detection and flagging in large survey pipelines.

When the input coefficients to PCA are not scaled, the mean and variance of the WST coefficients are sufficient to partially separate all four data sets in WST-PCA space (see the interactive version of Figure 4). This is in contrast to what is expected for a simple N-point correlation function, which should have an expectation value of zero for N > 2 on “NHC,” “Smooth,” and “Rand.” While exclusion of the S₀ coefficient does not visibly inhibit class separation, we cannot rule out that information regarding the mean of the density distribution does enter the higher-order coefficients, owing to the [min,max] = [0,1] normalization we chose. We therefore believe that a more careful study of the distribution and expectation values for the WST coefficients is called for in future work.

The distributions of the four data sets under RWST-PCA are similar to the distribution under WST-PCA, except that the ellipsoidal distributions have lower eccentricity. The unscaled RWST-PCA basis also partially separated the four classes, and further improved the separation of “Rand” to have no overlap with the other three classes. The cumulative variance rises much more slowly as a function of the number of principal components, suggesting that the RWST does not admit as low-dimensional a representation, and that the variance in “SFD” under the RWST is similar for many directions.

We select the real part of the 3PCF coefficients for analysis by PCA. The distribution of the four data sets under 3PCF-PCA is similar to that under WST-PCA, except that “Smooth” and “Rand” have higher and lower eccentricities, respectively. The unscaled 3PCF-PCA basis separates “SFD” from the other three data sets, with little overlap. However, the expectation value of the 3PCF for “NHC,” “Smooth,” and “Rand” must be zero. As such, the centroids for all three data sets are coincident. Despite overlapping centroids, the major axes of “NHC” and “Smooth” are distinct, allowing partial separation. The distribution of “Rand” is so tight and circular that it overlaps entirely with “NHC” and “Smooth.” The cumulative variance of the 3PCF is between that of the WST and RWST, suggesting the 3PCF admits an intermediate-dimensional representation of SFD dust. We repeated this analysis, taking the modulus of the 3PCF coefficients instead. While the exact shape of the distribution for each data set changed, and the centroid of “NHC,” “Smooth,” and “Rand” shifted, the separation was not strongly affected. ²⁷

The lack of complete separation in the PC basis does not mean that these descriptors do not contain the information necessary to distinguish the four classes. PCA is an unsupervised technique, and as such cannot optimize for the separation of classes, as it has no knowledge of them. To determine whether the computed coefficients can distinguish the four classes, we turn to a supervised learning technique, LDA. We split the data sets into test/train sets (80%/20%) and compute the confusion matrix for the test set between the true labels and the assignment predicted by the LDA model, which was trained on the first set. ²⁸ In the confusion matrix, we show the actual number of images having a true class label that were predicted to have a certain class label. The grayscale coloring of the matrix indicates the true positive percentage relative to the number of test images, with the corresponding true label (the row sum). In this way, the row sum is always the test fraction (∼0.2) times the total number of images in each class. To quantify the predictive power of LDA on this multi-class problem, we report the precision score, which is the total fraction of correctly classified images (true positive percentage).

Figure 5 illustrates that the WST contains enough information to completely distinguish the four classes, and can thus at least partially characterize the non-Gaussianity of SFD dust. The RWST performs similarly well, apart from the presence of an outlier in the RWST-LDA space that is not associated with any of the clusters. These outliers tend to result from undertraining, i.e., failing to fully constrain non-discriminatory dimensions. While not shown, we also analyzed an RWST model, where the frequency in Equations (6) and (7) were fit in addition to the usual RWST parameters. This unfortunately created a large degeneracy between the reference angle and frequency, which caused fitting instability and massive variability that were common to all data sets (not discriminatory). Outliers may also arise when the cosinusoidal dependence is so small in amplitude that the angle is undefined, or when the cosinusoidal frequency differs from the model. These points demonstrate that more work remains to be done to regularize the RWST technique.

Zoom In Zoom Out Reset image size

For improved classification, we take the modulus of the 3PCF coefficients as inputs to LDA. Under these conditions, 3PCF-LDA almost completely distinguishes the four classes, with only a few images of “SFD” confused for “NHC”, and one extreme outlier. Viewed in the interactive version of Figure 5, these confused images between “NHC” and “SFD” are associated with the “SFD” cluster, suggesting that the misclassification derives from the strictly linear boundaries to which LDA is limited. In addition, the assumption of equal class covariances (homoscedasticity) is strongly violated for the 3PCF, because the variance is strongly class-dependent (Figure 5, bottom right). To attempt to further improve classification with 3PCF, we use QDA, which relaxes the homoscedasticity condition, and removes the confusion between “NHC” and “SFD.” In all other sections of the manuscript, LDA and QDA produce nearly identical results, and so only the simpler LDA results are presented. When attempting classification with only the real part of the 3PCF coefficients, there was a large confusion between “NHC,” “Smooth,” and “Rand,” ²⁹ which are forced to have exactly overlapping centroids. This condition is lifted upon taking the modulus, since the expectation value of the modulus of the 3PCF coefficients need not be zero on images drawn from Gaussian distributions.

In summary, all three non-Gaussian descriptors contain enough information to distinguish “SFD,” “NHC,” “Smooth,” and “Rand,” but with the WST, the classes can be distinguished using a simple linear classifier (in contrast to classification with the 3PCF), and without outliers (in contrast to classification with the RWST).

Since WST, RWST, and 3PCF distinguish Gaussian from non-Gaussian fields, we seek to use them to classify different physical processes that can contribute to non-Gaussianity. One application with implications for the ISM is inferring characteristic MHD parameters, such as the sonic and Alfvénic Mach number. We demonstrate the ability of WST-LDA to do this on idealized MHD simulations (Figure 3), as a proof of concept. For each of the eight conditions (sonic/Alfvénic Mach number pairs) we have the full density cube (256³ pixels) at nine time steps. Because the separation between these time steps is longer than the eddy turnover time, each can be treated as an independent sample for the given condition. For each cube, we take 2D slices (256² pixels) along each axis (x, y, z), spaced every eight pixels so as to reduce the sampling redundancy. ³⁰ Therefore, we computed the WST coefficients for 6,912 (8 × 9 × 32 × 3) images, or 864 images per class. To ensure that the test set is not correlated with the training set, all images in the test set are selected from two independent time steps, withheld from the training set. Since each time step is an independent sample, the test and train sets are thus independent.

We achieve a precision of 97% for WST-LDA (J = 8, L = 8, M = 2; N_coeff = 1, 857), 98% for RWST-LDA (N_coeff = 164), and 87% for 3PCF-LDA (N_bins = 8, L = 9; N_coeff = 360) (Figure 6). The precision for WST-LDA is comparable to that obtained previously for differentiating sub/super-Alfvénic simulations after training a neural net (Peek & Burkhart 2019), but notably required no training except for a deterministic LDA step. As shown in the interactive version of Figure 6, the first three LDA components separate each of the classes. There is a small amount of crosstalk between adjacent classes, due to scatter of the WST coefficients of images in the same class.

Zoom In Zoom Out Reset image size

Since these MHD parameters are actually continuous variables, the variance of each class can be viewed as an indicator of the uncertainty of an inference of these MHD parameters from a single image. To obtain a rough estimate of this uncertainty, we assume M_S varies linearly as a function of separation between clusters of the same M_A value. Then, we find that M_S only varies ∼10% within 1σ of the centroid of each cluster. Thus, for these simulations, we can approximately infer M_S ± 10%. Note that the sonic Mach numbers studied here are lower than the average estimated for molecular clouds in the ISM from CO line widths (∼10–20, Kainulainen & Tan 2013; Kainulainen & Federrath 2017), but comparable to measurements for 21-cm cold neutral medium (CNM) (Heiles & Troland 2003; Burkhart et al. 2010). However, whether this approximately linear mapping of M_S to LDA space holds at larger M_S remains an open question. We hold off on a comparison between observations and simulations until we can characterize a wider range of Mach numbers and other physical parameters.

In this WST-LDA space, the classes of images with different Alfvénic Mach number form a helical structure (Figure 3, interactive version). Along the vertical axis of the helix, the sonic Mach number varies monotonically while the separation in the plane perpendicular to this axis determines the Alfvénic Mach number. We can also see this behavior in the 2D projections, where LD₀ is the helical axis (Figure 6). There is little structure as a function of Mach numbers in the 4th and higher LDA components, suggesting the classification occurs primarily in a 3D subspace of the WST-LDA basis.

The increased precision of RWST-LDA as compared to WST-LDA appears to derive from reduced scatter around the centroids of each cluster, although the centroids are separated by the same distance. Since LDA selects the optimal linear combination of the WST coefficients, the success of RWST is not simply in selecting an optimal subset of the coefficients containing the most information; the RWST is not simply performing dimensional reduction. The nonlinear combination of WST coefficients used by the RWST preferentially selects information relevant to differentiating the MHD classes.

For improved classification, we take the modulus of the 3PCF coefficients as inputs to LDA, and find partial separation of the classes. We find that the variance is strongly class-dependent, with higher sonic Mach numbers having larger variance, and lower sonic Mach numbers having smaller variance. This suggests that the variance of the 3PCF is sensitive to sharp density fluctuations. Furthermore, the spacing of the centroids is highly nonlinear. The small spacing between these centroids for M_S = 0.7 is the source of most of the confused images. ³¹ This small spacing suggests the contribution of an important influence from M_A in the fourth-order moments, which is present in the WST, but not in the 3PCF. As aforementioned, for cold gas in the ISM, M_S is likely to be >0.7; as such, this lack of discriminatory power for low M_S may not be a problem for some applications. The fact that gas near young hot stars has lower M_S suggests that it may be more difficult to apply WST-LDA to active star-forming regions (Gallegos-Garcia et al. 2020).

Since dust maps such as SFD project a 3D volume along the line of sight down to 2D, we investigate how WST-LDA, RWST-LDA, and 3PCF-LDA perform with respect to projected images by repeating our sampling procedure on the cumulative sum of density along a given coordinate axis. This is equivalent to sampling the column density along the line of sight in the optically thin limit. We computed the WST coefficients for the same number of images as before, i.e., 6912. Note that the line of sight sums introduce further correlations, thereby reducing the effective size of our data set, and can be seen as clusters or lines within a given class. We find a precision of 57% for WST-LDA, 83% for RWST-LDA, and 44% for 3PCF-LDA (see Appendix D, and the figure therein). The improvement provided by the nonlinear transforms in the RWST becomes even clearer in the context of cumulative sums, and we anticipate that this will be instrumental in applications featuring multiple clouds along the line of sight.

Overall, the WST and RWST seem most suited to the continuous inference problem of estimating Mach numbers from measured density fields, since M_S maps approximately linearly into the LDA classification space, in contrast to the highly nonlinear mapping for 3PCF-LDA. While the nonlinear transformations of WST coefficients made by the RWST improve the classification, as well as rendering the classification more robust to cumulative sums, future work should seek an optimal nonlinear transformation of WST coefficients that does not also create outliers. Since we do not yet have such an optimal transformation, we focus here on characterizing the performance of WST-LDA on several test cases.

It is important to establish how many coefficients the WST needs to characterize a non-Gaussian process. However, this question is intertwined with that of the scale at which the process of interest is occurring, and how many different scales the WST needs to characterize a non-Gaussian process. We study this trade-off in our MHD test case by computing those scattering coefficients that include only a subset of the scales in the problem. Firstly, we compute the WST coefficients for L = 8, M = 2, and J = 3, 4, 5, 6, 7, and 8 to investigate the effect of the maximum scale (). ³² We can follow a similar process to find by subsampling the 256² image down to a image (i.e., binning down by a factor of ), and computing the WST with L = 8, M = 2, and J = 8, 7, 6, 5, 4, and 3. Both results are shown in Figure 7. Note that the (var) and (var) curves cross at j = 4. The j = 0 (smallest spatial) scale is likely to contain little information if the image is well sampled. Above the j = 0 scale, the linearity of each of the data sets suggests that there is approximately the same information content in each scale.

Zoom In Zoom Out Reset image size

The number of coefficients (the dimension of the WST vector) decreases as we narrow the range of scales in the above analysis. To test whether the loss of information results explicitly from reduction in the number of coefficients versus loss of information in the image due to binning effectively eliminating certain scales, we generate a series of images of fixed size but with increasing redundancy. We leverage the fact that these simulations have periodic boundary conditions, and can therefore repeat a given image in a 2ⁿ × 2ⁿ matrix, downsample the image to 256 × 256 resolution, and always compute the WST for J = 8, L = 8, M = 2. These results ( “constant”) are always within 6% of the “variable”, so we conclude that the trend observed in the “variable” is predominantly the result of the scales excluded, and is not due the number of coefficients.

While the highest-resolution maps of the ISM are 2D, the resolution of 3D maps is improving (Green et al. 2019; Leike et al. 2020). Owing to the fact that the physical processes driving the structure of the ISM are happening in 3D, one may expect a measurement of non-Gaussianity that is computed in 3D to capture more information. We show that 3D-WST-LDA (J = 6, L = 8, M = 2) is comparable to 2D-WST-LDA in terms of separating the eight classes in our set of MHD simulations when the same scales are included (Figure 7). ³³ 3D-WST-LDA specifically outperforms 2D-WST-LDA by ∼5% for j = 5, 6. Note that 3D-WST is a strictly lower dimensional representation for all j studied, as we do not compute second-order coefficients between different angular bins in the 3D-WST procedure. ³⁴ 2D-WST and 3D-WST also use different wavelets (see Equations (1) and (4)) and so we leave a more detailed comparison of the different bases to future work. 3D-WST-LDA can clearly provide similar classifications to 2D-WST-LDA, though the increased computational expense means it is only likely to be used when 3D correlations are necessary (see Appendix A for computational costs).

4.4.1. Using Velocity Information

While few 3D maps of the ISM exist in position–position–position (PPP) space, and those that do exist are of relatively low resolution (1 pc; Leike et al. 2020), high-resolution maps of gas tracers (H i, ¹²CO, ¹³CO) in position–position–radial velocity (PPV) space are available (Dame et al. 2001; Ridge et al. 2006; Pineda et al. 2008; Peek et al. 2018; Clark & Hensley 2019). These PPV maps are also used to constrain the sonic Mach numbers of molecular and atomic clouds (Heiles & Troland 2003; Burkhart et al. 2010; Kainulainen & Tan 2013). We convert our MHD simulations to PPV space by interpolation, using the density field and x-velocity field. ³⁵ Representative 2D slices in PP and PV space are shown in Figure 8. Classification of 2D slices was difficult (average precision 63%), but by utilizing the 3D-WST (J = 6, L = 8, M = 2), we achieved a precision of 71%. ³⁶ We suspect that the reason for the difficulty in terms of 2D classification is that there is little correlation between the density and velocity fields observed by Burkhart et al. (2009). However, those authors noted a correlation between Mach numbers and velocity dispersion, which may explain why the 3D-WST in PPV space, which captures a large range of velocities, performs well. The spherical harmonic-based wavelets used here treat all three axes identically. While the velocity axis is not interchangeable with the two position axes, we use these simple wavelets to demonstrate the benefit of the correlation between velocity and position channels. ³⁷ This ability of WST-LDA to operate in PPV space has potential application both to a wider variety of data, and to higher-resolution 3D data. This higher-resolution PPV data may prove invaluable for mapping regions with overlapping clouds observed in dust and gas tracers. Given that a large range of velocities was required to determine the Mach number, resolving overlapping clouds may be difficult, unless the clouds are well-separated in velocity space.

Zoom In Zoom Out Reset image size

4.4.2. Artifacts and Missing Data

Pattern noise is a common feature in imaging data sets. Long-term sensitivity drift in scanning detectors produced stripe artifacts in IRAS (Wheelock et al. 1994) and AKARI (Doi et al. 2015). A ground loop in the readout electronics of a charge-coupled device (CCD) can imprint nearby radio transmissions on the data. ³⁸ Even highly sophisticated instruments on the Hubble Space Telescope can experience electronics problems that leave subtle “herringbone” patterns in the data (Jansen et al. 2010). A desideratum of a classifier is that it be robust to such pattern noise, or at least be able to give warning of its presence.

For testing purposes, we mimic pattern noise by adding a sinusoidal pattern with a random wavevector to the MHD slices. Specifically, we generate for each MHD slice (6,912 draws) a random amplitude (10%–30% of signal max), phase, period width (2–25 pixels), and angle relative to the x and y axis, and add it to the MHD slice prior to computing the WST coefficients (Figure 9). In some cases, large contiguous regions of an image may be lost due to satellite transits, diffraction spikes, or ghosts of bright stars. We model this by introducing a cutoff for eliminating pixels, based on the sinusoidal field used above (in this case with periods in the range of 43-256 pixels, and amplitudes of 10%–20%). The sinusoidal field was not added to the data, but the values in the image were set to zero for all indices where the sinusoidal field exceeded 10% (Figure 9). In both cases, we observe at most a slight reduction in precision (96% and 94%, respectively, compared to 97%) and increased variance for each class. However, the quantitative classification is evidently robust to these artifacts, and the qualitative structure of the clusters in the space of the first three LDA components was also unaffected.

Zoom In Zoom Out Reset image size

4.4.3. Dependence on Point-spread Function

In observational images, the true signal is also convolved with a point-spread function (PSF). We model this as a Gaussian PSF, with a FWHM of 2.355σ (in pixels). Moreover, we do not assume that the PSF is rotationally symmetric, and allow for nonzero aspect ratios AR = σ_x /σ_y . We investigated the following PSFs in detail: (σ, AR) = (2,1), (2,2), (4,1), (4,2), (4,3). The MHD classification by WST-LDA for samples from a single PSF was qualitatively unchanged, with a precision of ≥94%. Further, when all the WST coefficients from the PSF-blurred images and original images were combined, LDA was able to classify the images on Mach numbers with a precision of 92%.

However, we must emphasize that the optimal LDA space for each individual PSF is different. For instance, the LDA space for (2,1), and the LDA space for the original images are not the same. In the interactive version of Figure 10, we demonstrate how the WST coefficients for (2,1) look in the LDA space for the original images, and vice versa. While these spaces are not identical, the helix that we observed in the native LDA space is still present, but shifted and compressed, or expanded. We observe that for small σ and AR, the old LDA space can be approximately continuously deformed into the new LDA space (not shown).

Zoom In Zoom Out Reset image size

Ideally, we would find a classification space which is as insensitive as possible to the PSF. To approximate this, we take the original images, in addition to five blurred copies ((2,1), (2,2), (4,1), (4,2), (4,3)), and define the WST-LDA space for the Mach number classification task on that entire set. To better understand how the position of an image’s WST coefficients in the LDA space depend on the σ and AR of the PSFs, we take four random slices from each time step of each MHD class, and compute the WST coefficients for a Cartesian product of σ and AR. We then plot the WST coefficients corresponding to each each image with a given PSF in the space defined above in Figure 10. ³⁹ We observe that increasing the PSF radius moves images to lower LD₂, as compared to the original helical cluster. On this tail induced by the large PSF radii, the images are separated along LD₁ by AR. As such, both σ and AR are encoded in the WST-LDA space. Similarly, the WST-LDA spaces for = 0 and = 1 are related, but not identical, suggesting that the scale of a feature relative to the minimum pixel scale will also deform the WST-LDA space. We foresee the creation of a classifier which is invariant to, or deforms in a known way in response to, changes in the PSF and pixel scale. This is likely to be the next step before WST-LDA can be applied to classify the relevant hydrodynamic parameters (or any other parameter of interest) for observational data. It is primarily for this reason that we hold off regarding any analysis of the SFD data set, in terms of the classification space derived from the MHD simulations.

In the MHD turbulence simulations, the original magnetic field direction in the simulation box is topologically constrained by the periodic boundary conditions of the simulation. This introduces a preferred direction in the simulation, because the characteristics and scaling of turbulence, both parallel and perpendicular to the magnetic field, depend on the Alfvénic state of the gas (Goldreich & Sridhar 1995). The mean magnetic field remains relatively unperturbed in the sub-Alfvénic regime, in contrast to the super-Alfvénic regime, where the direction of field lines can wander (Lazarian & Vishniac 1999; Cho et al. 2002). In the simulations studied here, the mean magnetic field is initialized along the x-axis. We therefore pooled for all sub-Alfvénic MHD classes (M_A = 0.7) density slices taken perpendicular to the x, y, and z axes, and show that WST-LDA picks out the x-axis, as distinct from the y- and z-axes (Figure 11). While the precision for separating x, y, and z is shown in the Figure, the fact that we find a precision of 99% for differentiating y and z from x is a more meaningful conclusion, in light of the degeneracy of y and z. For the super-Alfvénic MHD classes (M_A = 2.0), WST-LDA returns a precision of only 53%. This reduced precision, as compared to the sub-Alfvénic classes, indicates that WST-LDA is sensitive to the effect of magnetic field wandering; the signatures of the magnetic field occurring in the density field for the sub-Alfvénic classes were no longer present. Thus, WST-LDA may provide an insight into the direction of homogeneous magnetic fields in the ISM.

Zoom In Zoom Out Reset image size