CHARACTERIZING MAGNETOHYDRODYNAMIC TURBULENCE IN THE SMALL MAGELLANIC CLOUD

Higher order statistical moments of density fluctuations have been studied extensively. For example, the variance of density fluctuations has been shown to increase with the sonic Mach number, ${\cal M}_s$ (Nordlund & Padoan 1999; Ostriker et al. 2001; Cho & Lazarian 2003). Therefore, if turbulence is the dominant structuring mechanism, ${\cal M}_s$ can be estimated from the variance of density fluctuations. However, the observable that is the most easily available from observations for various ISM tracers is the column density. While this is a less direct measure of turbulence compared to velocity, a comparison between observations and simulations is the most straightforward in the case of column densities.

Only very recently, Kowal et al. (2007), henceforth referred to as KLB, have investigated how variance, skewness, and kurtosis, the second-, third-, and fourth-order moments respectively, depend on ${\cal M}_s$. They found strong correlations: as the sonic Mach number increases, so does the Gaussian asymmetry of the column density (and density) PDFs due to gas compression via shocks, resulting in the increase of variance, skewness, and kurtosis. KLB used limited resolution models of 128³ in their study, while Burkhart et al. (2009), henceforth referred to as BFKL, saw the same trends using high-resolution isothermal simulations. In both studies, the moments had little dependence on the Alfvén Mach number $({\cal M}_A)$, or the line-of-sight (LOS) orientation used for integrating three-dimensional simulated data cubes. These studies are motivational and imply that the sonic Mach number of turbulence in interstellar clouds could be characterized by variance, skewness, and kurtosis of observed column density distributions.

The two-dimensional spatial power spectrum characterizes the energy distribution over spatial scales. SX99 found that the spatial power spectrum of individual H i velocity slices is well fit by a power law, with an average slope of ≈−3. Stanimirović & Lazarian (2001) showed that the power-law slope steepens when several channels are integrated together, and used this to estimate the density power-law slope of −3.3 and the velocity slope, which was estimated to be −3.4 by applying the Lazarian & Pogosyan (2000) VCA technique. No obvious breaks or preferred scales were found over the range of 30 pc to 4 kpc. The density and velocity spectral slopes are similar, and in the case of incompressible MHD turbulence, density behaves as a passive scalar and thus scales in the same way as velocity (Monin & Yaglom 1967; Lithwick & Goldreich 2001; Cho & Lazarian 2003). However, the estimated slopes are slightly more shallow than the predictions for the Kolmogorov spectrum (k^−11/3), which is expected for incompressible fluids.

However, different types of turbulence are expected to have different spectral slopes. For example, in incompressible fluids with a strong magnetic field, the spectrum is expected to be even steeper and scale as k^−13/3 (Biskamp 2003).

When the medium is supersonic (as we will see later applies to parts of the SMC), these relations are no longer valid due to shocks forming highly asymmetric density structures. KLB demonstrated that the spectral slope becomes more shallow with increasing ${\cal M}_{s}$. This can be understood as shocks in supersonic turbulence create more small-scale structure in density (Beresnyak et al. 2005). This behavior was found to be weakly dependent of the Alfvén Mach number.

While the spatial power spectrum has long been applied on both observations and simulations, it essentially uses only the amplitude of the Fourier transform of the initial signal, while the phase information is totally ignored. The bispectrum, however, is a three point statistical measure which incorporates both the amplitude and phase of the correlation of signals in Fourier space. Because of this, it can be used to search for nonlinear wave–wave interactions and characterize how shocks affect turbulent properties of the ISM. The bispectrum has been used in cosmology and gravitational wave studies to detect departures from Gaussianity (Fry 1998; Scoccimarro 2000; Liguori et al. 2006), and for the characterization of wave–wave interactions in laboratory plasmas (Intrator et al. 1989; Tynan et al. 2001). The first application of the bispectrum on synthetic astrophysical MHD turbulence was in BFKL.

BFKL found a general correlation between the bispectrum of two-dimensional column density and three-dimensional density maps for simulated data cubes of 512³ resolution. Also, supersonic models showed a much greater degree of correlation between structures of different scales than subsonic models. While comparing models with the same sonic Mach number, models with a stronger magnetic field (i.e., sub-Alfvénic) showed enhanced correlations. These results are encouraging and suggest that the bispectrum can be also used to provide insight into the nature of the turbulence cascade.

The small-scale H i structure of the SMC was observed with the ATCA, a radio interferometer, in a mosaicking mode (Staveley-Smith et al. 1997). Observations of the same area were obtained also with the 64 m Parkes telescope, in order to map the distribution of large-scale features. Both sets of observations were then combined (see SX99), resulting in the final H i data cube with angular resolution of 98'', velocity resolution of 1.65 km s⁻¹, and 1σ brightness temperature sensitivity of 1.3 K, to the continuous range of spatial scales between 30 pc and 4.4 kpc. The velocity range covered with these observations is 90–215 km s⁻¹. For details about the ATCA and Parkes observations, data processing, and data combination (short spacings correction), see Staveley-Smith et al. (1997) and SX99.

The H i column density image is shown in Figure 1. We corrected the original image by multiplying the H i column density of each pixel (N_HI in atom cm⁻²) with the correction factor f_c derived by SX99:

$$\begin{eqnarray*} f_c = \left\lbrace \begin{array}{@{}lr} 1+0.667(\log N_{{\rm HI}}-21.4) & : \log N_{{\rm HI}} > 21.4 \\ \ms 1 & : \log N_{{\rm HI}} \le 21.4. \end{array} \right. \end{eqnarray*}$$

Zoom In Zoom Out Reset image size

As the original position–position–velocity (PPV) SMC cube has 578 × 610 × 78 pixels, the resulting column density image is a two-dimensional array with 578 × 610 pixels.

We generate a database of 13 three-dimensional numerical simulations of isothermal compressible (MHD) turbulence by using the code of KLB and varying the input values for the sonic and Alfvénic Mach number. The sonic Mach number is defined as ${\cal M}_s \equiv \langle |{\bf v}|/C_s \rangle$, where v is the local velocity, C_s is the sound speed, and the averaging is done over the whole simulation. Similarly, the Alfvénic Mach number is ${\cal M}_A\equiv \langle |{\bf v}|/v_A \rangle$, where $v_A = |{\bf B}|/\sqrt{\rho }$ is the Alfvénic velocity, B is magnetic field, and ρ is density. KLB used resolution of 128³, while we use 512³. Details about the code used in KLB were published (see Cho et al. 2002) and the code was used in several studies. We briefly outline the major points of their numerical setup.

The code is a second-order-accurate hybrid essentially nonoscillatory (ENO) scheme (see Cho & Lazarian 2002) which solves the ideal MHD equations in a periodic box:

$$\begin{eqnarray} &&\frac{\partial \rho }{\partial t} + \nabla \cdot (\rho {\bf v}) = 0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&\frac{\partial \rho {\bf v}}{\partial t} + \nabla \cdot \left[ \rho {\bf v} {\bf v} + \left(p + \frac{B^2}{8 \pi } \right) {\bf I} - \frac{1}{4 \pi }{\bf B}{\bf B} \right] = {\bf f}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\frac{\partial {\bf B}}{\partial t} - \nabla \times ({\bf v} \times {\bf B}) = 0, \end{eqnarray} \tag{ 3 }$$

with zero-divergence condition ∇ · B = 0, and an isothermal equation of state p = C²_sρ, where p is the gas pressure. On the right-hand side, the source term $\bf {f}$ is a random large-scale driving force.³ We drive turbulence solenoidally, at wave scale k equal to about 2.5 (2.5 times smaller than the size of the box). This scale defines the injection scale in our models in Fourier space to minimize the influence of the forcing on the generation of density structures. Density fluctuations are generated later on by the interaction of MHD waves. The time t is in units of the large eddy turnover time (∼L/δV) and the length in units of L, the scale of energy injection. The scale of energy dissipation is defined by the numerical diffusivity of the scheme.⁴ The magnetic field consists of the uniform background field and a fluctuating field: B = B_ext + b. Initially b = 0.

We divided our models into two groups corresponding to sub-Alfvénic (B_ext = 1.0) and super-Alfvénic (B_ext = 0.1) turbulence. For each group we computed several models with different values of gas pressure (see Table 1). We ran compressible MHD turbulent models, with 512³ resolution, for t ∼ 5 crossing times, to guarantee full development of energy cascade. Since the saturation level is similar for all models and we solve the isothermal MHD equations, the sonic Mach number is fully determined by the value of the isothermal sound speed, which is our control parameter. The models are listed and described in Table 1.

Table 1. Description of the Simulations: MHD, 512³

Model	pgas	Bext	${\cal M}_s$	${\cal M}_A$	Description
1	2.00	1.00	0.1	0.7	Subsonic & sub-Alfvénic
2	1.00	1.00	0.7	0.7	Subsonic & sub-Alfvénic
3	0.10	1.00	2.0	0.7	Supersonic & sub-Alfvénic
4	0.025	1.00	4.4	0.7	Supersonic & sub-Alfvénic
5	0.01	1.00	7.0	0.7	Supersonic & sub-Alfvénic
6	0.0077	1.00	8.4	0.7	Supersonic & sub-Alfvénic
7	0.0049	1.00	10	0.7	Supersonic & sub-Alfvénic
8	1.00	0.10	0.7	2.0	Subsonic & super-Alfvénic
9	0.10	0.10	2.0	2.0	Supersonic & super-Alfvénic
10	0.025	0.10	4.4	2.0	Supersonic & super-Alfvénic
11	0.01	0.10	7.0	2.0	Supersonic & super-Alfvénic
12	0.0077	0.10	8.4	2.0	Supersonic & super-Alfvénic
13	0.0049	0.10	10	2.0	Supersonic & super-Alfvénic

Download table as:  ASCIITypeset image

For each model, the results of MHD simulations are the isothermal three-dimensional density field, three components of velocity (V_x, V_y, V_z), and three components of magnetic field. As an example, Figure 2 (left) shows a density field for a sub-Alfvénic, subsonic simulation. To calculate the column density distribution we integrate the three-dimensional density fields along a given LOS. We introduce the following nomenclature throughout the paper: "x column density" or "column density in the x-direction" refers to the density cube being integrated along the x-direction (parallel to the mean magnetic field). This description is similar for the y and z-directions (perpendicular to B_ext). In the case of Figure 2 (left), the LOS is along the z-axis, and the magnetic field is oriented along the x-axis.

Zoom In Zoom Out Reset image size

3.2.1. Scaling of the Column Density

In order to compare simulations with the SMC data, we apply a simple scaling to the simulated column density distributions:

$$\begin{equation} N_{\rm norm}=\frac{N-\langle N \rangle }{\sigma (N)}, \end{equation} \tag{ 4 }$$

where σ(N) denotes the standard deviation and 〈N〉 denotes the arithmetic mean of the column density image. This method, often referred to as the standard score, standardizes all data sets used in this study. After the scaling, the new data set has values which represent the difference between the original data values and the sample mean, in units of the standard deviation. The same scaling was applied on the SMC H i column density image.

3.2.2. Simulating Cloud Boundaries

While it is useful to know the Mach number of the ISM in a global sense, even more interesting is to see how it varies spatially and whether it correlates with local star forming regions where high turbulence is expected. To characterize small-scale departures from Gaussian distributions within the column density distributions, we create moment maps of the simulated column densities using a circular moving kernel. We do this by moving a circle of a given radius r across the image and calculating the skewness and kurtosis at all points.

First, we must decide on a kernel size that will enclose enough pixels to provide reliable statistics. According to Tabachnick & Fidell (1996), the standard error in skewness (SES) can be estimated by $\sqrt{6/N}$ and the standard error in kurtosis (SEK) can be estimated by $\sqrt{24/N}$, where N is the number of points used to calculate the skewness/kurtosis. Generally, if a value of skewness/kurtosis falls within ±2×SES/SEK, then the distribution is considered to be normal, while values of skewness/kurtosis larger than twice the absolute value of SES/SEK imply significant non-Gaussianity. Figure 4 shows 2× SES and SEK as a function of the kernel radius r. Clearly, SES/SEK is proportional to 1/r; the smaller the radius, the higher the absolute value of SES/SEK. Thus, we select r = 35 pixels, which corresponds to the point where 2× SES/SEK starts to deviate from zero and ensures that Gaussianity is represented by 2× SES/SEK ∼0.

Zoom In Zoom Out Reset image size

Using a circular moving kernel we generate third and fourth moment maps of simulated column density distributions. To ensure there are no resolution artifacts, we briefly explore correlations between skewness and kurtosis of the derived moment maps on a pixel-to-pixel basis. Based on the work by KLB and BFKL, and also our Figure 3, we expect for supersonic models a correlation of both skewness and kurtosis with the sonic Mach number, and that skewness and kurtosis agree in sign and relative value. However, as subsonic models show unexpected behavior at low resolutions (from the KLB study), the dependence of skewness and kurtosis on the sonic Mach number is not strong.

Figure 5 shows pixel-to-pixel comparison between skewness and kurtosis for two example models: a supersonic with ${\cal M}_{s}=7.0$ (left panel) and a subsonic with ${\cal M}_{s}=0.7$ (right panel). It is evident that the supersonic model shows a good correlation between kurtosis and skewness, while the subsonic model does not. This is exactly as expected. This result proves that our approach of deriving moment maps is not resolution limited although we use a smaller number statistics within each circular kernel (πr² = 1380 pixels). This also highlights another interesting property: regions of moment maps with a good correlation between skewness and kurtosis could be interpreted by supersonic MHD turbulence, while a poor correlation may be caused by subsonic turbulence. We present in the Appendix another technique that could be used to further identify subsonic regions.

Zoom In Zoom Out Reset image size

As shown above, with our limited angular resolution of 30' we see strong hints that different regions in the SMC have different turbulent properties, suggesting spatial variations of ${\cal M}_{s}$. We would now like to combine all this information into an image of ${\cal M}_{s}$ for the SMC. As already discussed, based on the KLB study, kurtosis is the least prone to effects caused by cloud boundaries and as Figure 3 suggests has an almost linear shape for ${\cal M}_{s}\sim 1\hbox{--}8$. Therefore, we fit a linear function for all values of kurtosis > ∼ −1. This function can be inverted to estimate ${\cal M}_{s}$: ${\cal M}_{s} \sim ({\rm kurtosis} + 1.44)/1.05$. Pixels in the ${\cal M}_{s}$ image with corresponding kurtosis <−1 are set to zero and most likely mark subsonic regions, or could be caused by additional processes. The resultant map is shown in Figure 8 (left).

Zoom In Zoom Out Reset image size

As expected already, most of the area of the SMC bar and the eastern wing could be interpreted as having subsonic or transonic Mach numbers, ${\cal M}_{s}\sim 0\hbox{--}2$. Several concentrated regions across the SMC indicate very quiescent environments; most of them unfortunately have a size close to our angular resolution. Regions with the highest sonic Mach number, ${\cal M}_{s}\sim 4$, are found around the bar and correspond to compressed H i contours.

We can quantify the fraction of H i with different ${\cal M}_{s}$. The most likely quiescent regions (set to zero in our map) comprise 8% of the mapped area. Regions with $0<{\cal M}_{s}\le 1$ comprise about 40%, while about the same fraction is contained in transonic regions with ${\cal M}_{s}=1\hbox{--}2$. Regions with higher turbulence, ${\cal M}_{s}>2$ constitute about 10% of the area.

Figure 8 (right) shows contours of the Hα emission (obtained by Kennicutt et al. 1995) smoothed to 30' and overlaid on the ${\cal M}_{s}$ map. Sites of the most recent star formation in the bar and the eastern wing have ${\cal M}_{s}\sim 1$, suggesting that the most turbulent regions are not associated with star formation. The most turbulent regions appear to trace shearing and tidal motions. We note though that our resolution of 30' is too low to trace individual star-forming regions.

As an estimate of the velocity dispersion along the LOS is often used as a possible measure of H i turbulence, we investigated how the H i velocity dispersion compares to our ${\cal M}_{s}$ map. Stanimirović et al. (2004) found that regions with the highest dispersion (∼30 km s⁻¹) appear to be associated with the positions of the three largest super-giant shells. We find no obvious correlation between the H i velocity dispersion and our ${\cal M}_{s}$ map.

This lack of correlation could be due to projection effects, and/or large LOS depth of the SMC. In addition, velocity dispersion (similar to velocity centroids) is subject to the interplay of several statistical effects (Lazarian & Esquivel 2003; Esquivel et al. 2003, 2007; Esquivel & Lazarian 2005). The most obvious is the contribution of both velocities and densities to the resulting measure. More subtle, but still important, are effects of velocity–density correlations and non-Gaussianity. These effects make centroids unreliable when dealing with supersonic turbulence (Esquivel et al. 2007). Effects of phase transitions resulting in the existence of cold and warm H i and the pliable equation of state corresponding to interstellar H i are likely to act in a similar way. Thus, the velocity dispersions obtained from spectral lines are substantially affected by densities, which are the quantities that we deal with in this study. They should be distinguished from the true velocity dispersions, which, unfortunately, are not available through spectral line observations.

The bispectrum is closely related to the power spectrum. In a discrete system, the power spectrum is defined by Equation (11). In a similar way, the bispectrum is defined as

$$\begin{equation} B(k_{1},k_{2})=\sum _{\vec{k_{1}}=k_{1}}\sum _{\vec{k_{2}}=k_{2}} A(\vec{k_{1}}) \cdot A(\vec{k_{2}}) \cdot {A}^{*}(\vec{k_{1}}+\vec{k_{2}}), \end{equation} \tag{ 12 }$$

where $\vec{k_{1}}$ and $\vec{k_{2}}$ are the wavenumbers of two interacting waves, and $A(\vec{k})$ is the original discrete time series data with finite number of elements with $A^{*}(\vec{k})$ representing the complex conjugate of $A(\vec{k})$.

The bispectrum is a complex function which measures both phase and magnitude information between different wave modes. As this is the first application of the bispectrum on the H i data we show in Figure 11(a) visual example of the difference between the power spectrum and the bispectrum. This figure shows the original H i column density image of the SMC (top left), and the same image but with manipulated phases (top right). To obtain the top right image we Fourier transformed the SMC image and randomized the phases with a Gaussian random distribution but left the amplitudes the same. Even though the phases of the top images are very different, the power spectrum uses only amplitude information and is identical for both images, as shown in the bottom panels of Figure 11. However, the bispectrum looks very different for the two images, as expected, and offers an insight into the phase information.

Zoom In Zoom Out Reset image size

In practice, our calculation of the bispectrum involves performing a Fast Fourier Transform (FFT) of the column density functions and the application of Equation (12). We randomly choose wavevectors and their directions, k₁ and k₂ and iterate over them, calculating k₃ = k₁ + k₂. We limit the maximum length of the wavevectors to half of the box size. We normalize direction vectors to unity, calculate positions in Fourier space, and finally, compute the bispectrum, which yields a complex two-dimensional array. When plotting the bispectrum (Figures 12–15) we plot bispectral amplitudes, which give information about the degree of mode correlations in the original column density distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to interpret the bispectrum of the SMC H i column density we first compute the bispectrum of the MHD simulations (scaled by the standard score method and without applying cloud boundaries) and look for dependence of wave–wave correlations on the sonic and Alfvén Mach numbers. Figure 12 shows the bispectrum of simulated column density distributions for the case of fixed ${\cal M}_{A}=0.7$ and two extreme cases of ${\cal M}_{s}$: 0.7 and 7.0. In Figure 13, we fix ${\cal M}_{s}=2.0$ and look at two cases of ${\cal M}_{A}$: 0.7 and 2.0.

The first thing to notice in the bispectral contour maps is that the k₁ = k₂ case always shows high correlation since this is a trivial case of two wavenumbers being the same. For k₁ ≠ k₂ the bispectral amplitude and isocontour shape vary with turbulent properties but generally amplitude decreases gradually radially from k₁ = k₂ = 0. Models with circularly/broad-shaped isocontours have high amplitudes and wave–wave correlations, while models with more narrow isocontours have lower amplitudes and therefore weaker correlations. For a fixed ${\cal M}_{A}$, supersonic models show a higher degree of wave–wave correlations over subsonic models. For a fixed ${\cal M}_{s}$, models with a higher magnetic field (sub-Alfvénic, e.g., ${\cal M}_{A}$ = 0.7) show somewhat stronger correlations than the models with a weaker magnetic field (super-Alfvénic) although this difference is not striking.

The bispectrum of the SMC H i column density is shown in Figure 14. The top row (first column) shows the bispectrum of the column density derived by integrating all 78 velocity channels. The middle and bottom rows show the bispectrum of the first (for channels 1–42) and the last half (for channels 43–78) of the data cube. To facilitate interpretation the wavenumbers k₁ and k₂ are shown in terms of linear size⁵. While there are small differences in the bispectral amplitudes, there is essentially no significant difference in the contour shape for the three bispectra.

To properly compare the SMC bispectrum with simulations we apply a windowing function on the SMC column density image to simulate the effect of periodic boundaries. We use the Hanning window w(n) function, defined as

$$\begin{equation} w(n)=0.5 \left(1- \cos \frac{2\pi n}{N-1} \right), \end{equation} \tag{ 13 }$$

where N is the number of pixels along the x- or y-axis, and n ranges from 0 to N − 1. This windowing function makes the map periodic in the Fourier space. The second column of Figure 14 shows the resultant bispectra for all three cases of H i column density integrations. The windowing seems to slightly change bispectral amplitudes on the large scales, but leaves the general structure of the bispectral contours intact, although some smoothing effects are observed. The increase in large-scale correlations is at a level of ∼10%.

Before interpreting the SMC bispectrum we briefly investigate how the presence of noise in the observational data affects the bispectrum. We made the noise image of the H i column density by taking a single velocity channel from the SMC H i data cube without H i emission and integrating this image over the whole velocity range. The bispectrum of the noise image (normalized using the standard score method) is shown in Figure 15 (right). In the case of purely Gaussian noise, we expect a random distribution over the whole space of wavelengths. As a result, the bispectrum will have a gradual increase in amplitude, with the highest amplitude being at large k (or small scale). Figure 15 (left) illustrates this effect for a pure Gaussian distribution. The bispectrum of the SMC noise image shows essentially the same trend, suggesting that we are dealing mainly with Gaussian (random) noise. As is obvious from Figures 12– 14, the bispectral amplitude increases in the opposite direction (from large to small scales). The effect of noise is therefore the most important at the smallest spatial scales, and we keep this in mind when interpreting the bispectrum.

Figure 14 shows the degree of wave–wave correlations in the H i distribution. Modes are obviously strongly correlated at the largest scales and show weaker correlation at small scale. This trend is similar to what was found for various MHD simulations and obviously very different from the case of Gaussian noise. In addition, the k₁ = k₂ line shows the strongest correlations, as is also seen in the simulations. However, contrary to simulations where the level of correlation gradually decreases along this line, we see significant small-scale variations. Several local peaks are noticeable and there is a strong break at mid-scales. The break occurs around k₁ = k₂ ≈ 160 pc and is seen in all three integrations of the H i data cube. This could be due to a lack of interactions between turbulent eddies at the scale of ∼160 pc, caused possibly by numerous expanding shells in the SMC. Hatzidimitriou et al. (2005) and Stanimirović (2007) showed that shells in the SMC range in size from 30 to 800 kpc, however, the radius distribution peaks at ∼60 pc. The break in the bispectrum may be signifying the lack of correlations on scales close to the typical shell diameter, or some additional physical processes.

We can also compare the bispectrum of the SMC with the bispectra of various simulations. By visual inspection of Figures 12 and 13, the closest simulation to the SMC in terms of contour shapes is the one for ${\cal M}_s= 2.0$ and ${\cal M}_A= 2.0$. Further work in this area is required to derive a real measure to quantitatively compare observed and simulated bispectra.

Using third and fourth statistical moments of the H i column density image and bootstrapping turbulent information from a database of isothermal MHD simulations, we have mapped spatial variations of the sonic Mach number across the SMC. While most of the H i seen in emission in the SMC appears to be subsonic or transonic, several supersonic regions have emerged from our study. It is interesting that these regions do not correlate well with the most recent sites of star formation and seem to point to large scale shearing or tidal flows. Commonly, it is believed that supernovae and superbubbles are the main drivers of galactic turbulence (McCray & Snow 1979), with a typical size of ∼100 pc. While we do not have high enough resolution to see changes on such small scales (∼10') in our derived map of ${\cal M}_s$, most of the star-forming bar of the SMC appears to have subsonic or transonic properties when viewed at resolution of 30'.

The most turbulent regions in the SMC may be tracing some kind of shearing flows between the SMC bar and the surrounding H i. This suggests that SMC's chaotic history with the LMC and our own Milky Way has probably left strong turbulent imprints on the H i gas. The lack of a turnover in the H i spatial power spectrum on the largest observed scales is also indicative of the fact that turbulent energy injection happens largely on scales larger than the size of the SMC (Stanimirović & Lazarian 2001). Similarly, Goldman (2000) suggested that the H i turbulence in the SMC was induced by large-scale flows from tidal interactions with the Milky Way and the LMC about 2 × 10⁸ yr ago. Such large-scale bulk flows could have generated turbulence through shear instabilities, and this turbulence has not have had enough time to decay.

Most of the H i in the SMC has a sonic Mach number of 1–2. This is on average at least 2 times smaller than what we inferred from H i absorption observations for the CNM in the SMC, ${\cal M}_s\sim$ 3.5–4. Similarly, for the CNM in the Milky Way Heiles & Troland (2003) found ${\cal M}_s \sim 3$ with a large dispersion. A sonic Mach number of about 4–5 is commonly assumed for cold gas in molecular clouds (Federrath et al. 2009). For example, Heyer et al. (2006) measured from CO observations ${\cal M}_s=4.2\pm 0.17$ for the Rosette molecular cloud, and ${\cal M}_s=4.7\pm 0.12$ for G216-2.5. On the other hand, Hill et al. (2008) found that the distribution of the warm ionized medium (WIM) in the Milky Way can be best fit by models for mildly supersonic turbulence with ${\cal M}_s\sim 1.4\hbox{--}2.4$. Turbulent properties of the H i in emission in the SMC are therefore closer to properties of the WIM in the Milky Way than properties of the CNM. This may suggest a large fraction of warm relative to cold H i being traced in H i emission.

In the Milky Way, the H i is known to consist of at least two components with different temperatures: the WNM with T_warm = 6000 K and the CNM with T_cold = 70 K. In addition, there could be a substantial amount of gas at intermediate temperatures (Heiles & Troland 2003). Due to its lower metallicity, the H i in the SMC has different properties. Dickey et al. (2000) found T_cold = 40 K, in agreement with theoretical expectations by Wolfire et al. (1995) whereby the existence of the two-phase medium is possible only at higher pressures compared with the range that applies for solar neighborhood conditions. They also estimated the fraction of cold H i in the SMC to be ≲15%. This is lower than ∼25% found for the Milky Way.

As our simulations are isothermal it is obvious to wonder how the multi-phase H i affects our statistics and conclusions. We investigate this in Figure 16 by producing a simulated data cube from a weighted combination of two cubes, one subsonic $({\cal M}_{s}= 0.7)$ and one supersonic. The subsonic cube represents contribution from warmer gas, while the supersonic cube represents the cold gas. We combined the two cubes with different emphasis on warm versus cold gas, obtained the column density image of the resultant cube, and calculated its moments.

Zoom In Zoom Out Reset image size

Figure 16 shows that for the case when the supersonic cube has ${\cal M}_{s}$ = 2.0 skewness of the final cube with up to 50% of subsonic gas will be biased toward supersonic gas and will appear dominated by cold gas. If we increase the sonic Mach number of the supersonic cube to 4.0, the dominance of the higher turbulence is even more pronounced. A cube with up to 60%–70% of subsonic gas and 25% of supersonic gas, will still have high skewness biased by the supersonic contribution. Considering that the H i column density image results in the mean ${\cal M}_{s}=1\hbox{--}2$, significantly lower than what is expected for the cold H i, the fraction of the CNM along any LOS is most likely ≲25%. This supports the Dickey et al. (2000) estimate of the fraction of cold H i in the SMC being about 15%.

Two different statistical approaches in our study suggest that the H i gas in the SMC seen in emission is super-Alfvénic. As we have shown in Figure 10, in addition to the sonic Mach number the spectral slope of the spatial power spectrum is sensitive to the Alfv´enic Mach number for ${\cal M}_{s}<2$. The sub-Alfvénic models generally show steeper slopes due to large-scale influence of magnetic fields. Thus, if one independently knows the sonic Mach number, it is possible to estimate the Alfv´enic one using just the column density data. While the dependence of the spectral slope on the Alfv´enic Mach number has not received much attention in the past, it is somewhat expected. Essentially, magnetization decreases compression in the shocks. Strong magnetic forces mix up density clumps preventing formation of isolated peaks, which results in a steeper spectrum. In addition, in the sub-Alfvénic case we expect oblique shocks to be disrupted by Alfvén shearing, which in turn, produces more small-scale shocks Beresnyak et al. (2005).

Another indication that the H i in the SMC is super-Alfvénic comes from the bispectrum. The very sharp decrease in the bispectral amplitudes from large to small scales observed for the SMC is the closest to the trend found for simulated data for the case of ${\cal M}_{s}\sim 2$ and ${\cal M}_{A}\sim 2$. Detailed comparison between simulated and observed bispectra awaits future work; however, this qualitative comparison is certainly encouraging.

Assuming on average ${\cal M}_{s}\sim 1\hbox{--}2$, the power spectrum slope suggests a super-Alfvénic H i in the SMC with ${\cal M}_{A}\sim 1\hbox{--}3$. This is generally in agreement with the observationally inferred strength of the magnetic field by Mao et al. (2008). Using their estimate for B_ext = 2 μG, a radius of the SMC of 2 kpc, the total hydrogen mass of 4.2 × 10⁹ M_☉, and a typical velocity dispersion of 20 km s⁻¹ (Stanimirović et al. 2004), we estimate ${\cal M}_{A}\sim 3$. As the Alfv´enic Mach number shows the nature of the interplay between gas pressure and magnetic fields, it appears that the gas pressure in the SMC dominates over the magnetic pressure.

Our bispectrum analysis of the SMC H i data was the first attempt to apply bispectrum on observed astrophysical data. While more detailed comparison between observations and simulations awaits future work, we clearly see trends in the bispectral amplitudes similar to what was found for simulations of supersonic MHD turbulence. The most interesting finding is, however, the effect of small-scale variations in the k₁ = k₂ correlations and a strong break in correlations at a scale of 160 pc. Such small-scale variations, or jumps, have not been seen in the bispectrum of simulated data. We can speculate about several possible scenarios that could explain their existence. The jumps could be caused by the energy injection due to processes other than turbulence affecting specific spatial scales. Alternatively, the jumps may be marking the presence of colder or multi-phase gas. Similarly, the observed break in the bispectrum at about 160 pc is intriguing. As we already pointed out, it is interesting that most expanding shells in the SMC (more than 500 were cataloged so far) have a diameter of ∼120 pc. The break could be due to the lack of correlations on scales similar to the distance between two shell centers. Obviously this will require further studies.

A natural question to ask is how results presented in this paper depend on the resolution of numerical simulations. For example, Kritsuk et al. (2007) investigated how resolution of numerical simulations affects the power spectrum of density. These authors found that the spectral index estimates based on low resolution simulations bear large uncertainties due to the bottle neck contamination, and that the power spectra of 512³ simulations are substantially shallower than models with resolution of 1024³. However, while Kritsuk et al. (2007) only examined hydrodynamic turbulence, Beresnyak & Lazarian (2009) showed that the slopes were very different between the MHD and pure hydrodynamic cases. For instance, the slopes for hydrodynamic simulations showed a pronounced and well-defined bottleneck effect, while the MHD slopes were much less affected. This is indicative of MHD turbulence being less local than the Kolmogorov turbulence, and suggests that our simulations will be less affected by resolution. In addition, Kritsuk et al. (2007) found a difference in the slope between hydrodynamic 512³ and 1024³ simulations to be 0.17. This would result in a change of $d{\cal M}_{s}\sim 0.5$ only and will not change our interpretation. We also add that in the case of higher statistical moments BFKL confirmed trends noticed by KLB at lower resolution of 128³.

Another issue that should be further addressed and that could affect our results is the type of numerical forcing of turbulence. Federrath et al. (2009) recently investigated the effects of solenoidal versus compressive (divergence-free versus curl-free) forcing on a variety of statistics including PDFs and higher order moments. They found that both types of driving mechanisms are compatible with observations of molecular clouds; however, depending on the data studied, one type could be superior to the other in terms of the statistics and reproduced observables. This implies that different regions in the SMC may exhibit statistical signatures of either compressive or solenoidally driven turbulence. In addition, Federrath et al. (2009) compared power spectra at 256³, 512³, and 1024³ for solenoidal and compressive forcing, and demonstrated that the change in the slope is only 3% between 512³ and 1024³. Thus, it is safe to say that our resolution of 512³ is sufficient for this type of study.

We have investigated the possibility of constraining turbulent properties of the ISM, specifically the sonic Mach number, by using the H i column density image and a database of numerical simulations with a range of sonic and Alfv´enic Mach numbers. By applying the third and fourth statistical moments on both observed and simulated data we have derived the spatial distribution of the sonic Mach number across the SMC with angular resolution of 30'. To provide an estimate of the Alfv´enic Mach number we used two approaches: the spatial power spectrum and the bispectrum. Using the database of numerical simulations we have confirmed that the spatial power spectrum varies with both the sonic and Alfv´enic Mach numbers. If the sonic number is known the Alfv´enic number can be constrained from this dependence. The bispectrum shows the level of correlation between turbulent eddies of different size and depends greatly on the sonic Mach number, and somewhat on the Alfv´enic Mach number. By comparing the bispectra of observations and simulations we have gauged the importance of magnetic fields relative to the gas pressure in the SMC. The following results were discussed in the paper.

• 1.  
  Skewness and kurtosis of the H i column density generally correlate well and are within the range expected from MHD simulations. This suggests that departures from Gaussianity could be interpreted as being governed by MHD turbulence.
• 2.  
  Most of the H i in the SMC bar and the eastern wing is subsonic or transonic with ${\cal M}_{s}\sim 0\hbox{--}2$. Sites of most recent star formation have ${\cal M}_{s}\sim 1$. Regions with the highest skewness and kurtosis, which could be interpreted as having ${\cal M}_{s}\sim 4$, correspond to the edges of the bar. The most turbulent regions are most likely tracing tidal or shearing flows. The fraction of the SMC with different turbulent properties is 10% with ${\cal M}_{s}>2$, 80% with $0<{\cal M}_{s}<2$, and about 10% with very low values of ${\cal M}_{s}$.
• 3.  
  Using H i absorption profiles from Dickey et al. (2000) we have estimated that the CNM in the SMC is highly supersonic with ${\cal M}_{s}=3.5\hbox{--}4$. This is at least a factor of 2 higher than what we measured from the higher statistical moments for the H i gas seen in emission. One possible reason for this discrepancy could be that H i emission is dominated by warm gas and the fraction of the CNM in the SMC is ≲20%.
• 4.  
  The slope of the spatial power spectrum and the bispectrum is suggestive that the H i in the SMC is mildly super-Alfv´enic with ${\cal M}_{A}\sim 1\hbox{--}3$. This implies that the gas pressure dominates over the magnetic pressure.
• 5.  
  The bispectrum of the H i column density shows large-scale wave correlations suggesting a large-scale energy injection mechanism. Contrary to simulations which show a smooth decrease of wave–wave correlations from large to small scales, the SMC bispectrum shows localized enhancements of correlations and at least one prominent break at ∼160 pc. We speculate that the multi-phase medium, and/or energy injection by processes other than turbulence, could be responsible for the correlation jumps. The break on the other hand appears at a scale similar to the diameter of the majority of expanding shells in the SMC.