A new method for spatially resolving the turbulence-driving mixture in the ISM with application to the Small Magellanic Cloud

ABSTRACT

Turbulence plays a crucial role in shaping the structure of the interstellar medium. The ratio of the three-dimensional density contrast (⁠|$\sigma _{\rho /\rho _0}$|⁠) to the turbulent sonic Mach number (⁠|$\mathcal {M}$|⁠) of an isothermal, compressible gas describes the ratio of solenoidal to compressive modes in the turbulent acceleration field of the gas, and is parameterized by the turbulence driving parameter: |$b=\sigma _{\rho /\rho _0}/\mathcal {M}$|⁠. The turbulence driving parameter ranges from b = 1/3 (purely solenoidal) to b = 1 (purely compressive), with b = 0.38 characterizing the natural mixture (1/3 compressive, 2/3 solenoidal) of the two driving modes. Here, we present a new method for recovering |$\sigma _{\rho /\rho _0}$|⁠, |$\mathcal {M}$|⁠, and b, from observations on galactic scales, using a roving kernel to produce maps of these quantities from column density and centroid velocity maps. We apply our method to high-resolution |${\rm H}\,\rm{\small I}$| emission observations of the Small Magellanic Cloud (SMC) from the GASKAP-HI survey. We find that the turbulence driving parameter varies between b ∼ 0.3 and 1.0 within the main body of the SMC, but the median value converges to b ∼ 0.51, suggesting that the turbulence is overall driven more compressively (b > 0.38). We observe no correlation between the b parameter and |${\rm H}\,\rm{\small I}$| or H α intensity, indicating that compressive driving of |${\rm H}\,\rm{\small I}$| turbulence cannot be determined solely by observing |${\rm H}\,\rm{\small I}$| or H α emission density, and that velocity information must also be considered. Further investigation is required to link our findings to potential driving mechanisms such as star-formation feedback, gravitational collapse, or cloud–cloud collisions.

1 INTRODUCTION

The interstellar medium (ISM) is turbulent and mostly composed of multiphase hydrogen gas. The density distribution in cold, dense, self-gravitating molecular clouds embedded in the diffuse ISM is strongly influenced by the turbulence and magnetic field of the warm, diffuse medium from which they condense and are embedded in. The star formation rate and efficiency are functions of the density distribution of molecular clouds (Elmegreen & Scalo 2004; Mac Low & Klessen 2004; Scalo & Elmegreen 2004; McKee & Ostriker 2007; Federrath & Klessen 2012; Hennebelle & Falgarone 2012; Padoan et al. 2014). Since the ISM is observed to be ubiquitously turbulent at all densities and temperatures, and since turbulence decays on short time-scales (Mac Low et al. 1998; Stone, Ostriker & Gammie 1998), ISM turbulence must be sustained by continuous energy injection into all phases of the ISM over a large range of scales.

Understanding the processes that drive the structure and evolution of the diffuse ISM, and the condensation of the warm neutral medium (WNM) to colder, denser gas phases from which stars are formed, is imperative to our continual quest for a coherent theory of ISM structure and star formation. Supernovae, stellar winds, protostellar outflows, radiative feedback, large-scale galactic dynamics, and galaxy interactions will all impact the evolution of density and velocity structures of the diffuse ISM in different ways, but they do so primarily through their ability to precipitate and drive turbulence (Elmegreen 2009; Federrath et al. 2017).

In order to characterize the turbulence in the ISM through comparison of observations, 3D simulations, and analytic models, various statistical methods have been developed that are applicable to the information available in observational data (see review by Burkhart 2021). Measuring the spatial power spectrum is one commonly used method, which reveals the energy injection and/or dissipation scale and the energy cascade in the turbulent ISM (e.g. Stanimirovic et al. 1999; Stanimirović & Lazarian 2001; Kowal & Lazarian 2007; Heyer et al. 2009; Chepurnov et al. 2015; Nestingen-Palm et al. 2017; Pingel et al. 2018; Szotkowski et al. 2019). In addition, investigating the higher order statistical moments, such as the skewness and kurtosis, or the probability distribution function (PDF) of column density is useful given the highly non-Gaussian behaviour of the density field of cold gas (Kowal & Lazarian 2007; Burkhart et al. 2009, 2010; Patra, Chengalur & Begum 2013; Bertram et al. 2015; Maier et al. 2017).

In (magneto)hydrodynamic simulations of isothermal, supersonic gas it has been shown that the PDF of the gas density can be described by a log-normal function, meaning that the logarithm of the density follows a normal Gaussian distribution (Vazquez-Semadeni 1994; Padoan, Jones & Nordlund 1997; Passot & Vázquez-Semadeni 1998; Federrath, Klessen & Schmidt 2008; Hopkins 2013; Squire & Hopkins 2017; Beattie et al. 2021). It is the width (standard deviation) of this density PDF that describes the density fluctuations of the gas, and is a key parameter in theoretical models of star formation (Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012; Burkhart & Mocz 2019). Furthermore, studies by Price, Federrath & Brunt (2011), Konstandin et al. (2012), Molina et al. (2012), Federrath & Banerjee (2015), Nolan, Federrath & Sutherland (2015), and Kainulainen & Federrath (2017) have shown that the width of the density PDF is proportional to the turbulent Mach number. Formally, this is described by

where |$\sigma _{\rho /\rho _0}$| is the three-dimensional (3D) standard deviation of density (ρ) scaled by the mean density (ρ₀), |$\mathcal {M}$| is the standard deviation of the (3D) turbulent velocity dispersion divided by the sound speed, and b is a constant of proportionality known as the turbulence driving parameter (Federrath et al. 2008, 2010). This constant of proportionality can be used to quantify how turbulence is being driven in the gas: the acceleration field that drives the turbulence can have purely solenoidal modes (divergence-free; b = 1/3) or purely compressive modes (curl-free; b = 1), or any combination of the two extremes. A value of b = 0.38 is the natural mixture, which corresponds to 1/3 of the power in the driving field in compressive modes and 2/3 of the power in solenoidal modes (see equation 9 in Federrath et al. 2010). The b parameter has been derived in theoretical models, and tested in simulations (Federrath et al. 2008, 2010; Price et al. 2011; Molina et al. 2012; Federrath & Banerjee 2015; Nolan et al. 2015; Beattie et al. 2021). In addition, crucial methods have been developed to allow for a mapping between the intrinsically 3D properties of the gas (⁠|$\sigma _{\rho /\rho _0}$| and |$\mathcal {M}$|⁠), and the respective quantities that are accessible in observations, i.e. the column density and the intensity-weighted velocity centroid (Brunt, Federrath & Price 2010; Kainulainen, Federrath & Henning 2014; Federrath et al. 2016; Stewart & Federrath 2022). Thus, we can recover the turbulence driving parameter b from observations, which allows us to learn more about what kind of physical processes are driving turbulence in those regions.

There have been several studies (primarily of molecular star-forming regions) in the Milky Way (MW) that have measured the b parameter from observations of the variance in the density and velocity. In Taurus, Brunt (2010) used |$\mathrm{^{13}CO\, J=1-0}$| observations to derive |$b=0.48^{+0.15}_{-0.11}$|⁠, which lies in the mixed-to-compressive regime, possibly a result of active star formation in that region. Ginsburg, Federrath & Darling (2013) study a non-star-forming giant molecular cloud (GMC) towards W49A and find a lower limit for b ≳ 0.4, concluding that it is being driven compressively relative to the natural mixture case. The authors posit that because this GMC is representative of all GMCs in the solar neighbourhood, compressive driving may be a common feature of GMCs in general. Federrath et al. (2016) find that the turbulence in a molecular cloud in the central molecular zone (CMZ), G0.250+0.016 (aka ‘The Brick’), is being driven solenoidally (b ∼ 0.22 ± 0.12) due to the strong shearing motions in the CMZ, which may account for inefficient star formation in that region of the MW. In agreement with Ginsburg et al. (2013) and Kainulainen & Federrath (2017) find that the lower limit of b for 15 solar neighbourhood clouds is ≳ 0.4 and rule out that any of these clouds can be dominated by solenoidal driving, although the authors point out that their density and velocity variance measurements do not always follow the standard relation described by equation (1). Menon et al. (2021) use ALMA observations of several CO isotopologues to measure the turbulence driving parameter in the ‘Pillars of Creation’ in the Carina Nebula, and find that all the pillars are being compressively driven, with b ∼ 0.7−1. The compressive driving in these star-forming regions could be a result of compression induced by photoionizing radiation. There has also been one extragalactic measurement of the b parameter in the Papillon Nebula in the Large Magellanic Cloud (LMC) by Sharda et al. (2022), which is also a molecular star-forming region. They find that this region is being driven almost purely compressively, with b ∼ 0.9, and speculate that this filamentary molecular cloud could have been formed by large-scale |${\rm H}\,\rm{\small I}$| flows due to tidal interactions between the LMC and the Small Magellanic Cloud (SMC). Marchal & Miville-Deschênes (2021) derive a driving parameter for the warm atomic gas at high galactic latitudes in the MW using |${\rm H}\,\rm{\small I}$| emission data, and found that the turbulence was driven mostly compressively in that region, with b ∼ 0.68. The scale on which all of these previous measurements of b have been made is typically on the molecular cloud scale or smaller, and they all derive a single estimate for the turbulence driving parameter of those clouds/regions. This is largely because the resolution required to robustly estimate the density and velocity variance is available only for nearby regions of the ISM, and the star-forming molecular clouds are the most well studied.

However, multiphase atomic |${\rm H}\,\rm{\small I}$| comprises the majority of the ISM, and understanding how turbulence driving affects its structure and evolution as it condenses into molecular hydrogen is crucial in understanding how molecular clouds form (Gazol et al. 2001; Koyama & Inutsuka 2002; Audit & Hennebelle 2005; Vázquez-Semadeni et al. 2006; Mandal, Federrath & Körtgen 2020; Seta & Federrath 2022). Here, we present a new method for mapping the turbulence driving parameter over large scales, and apply it to |${\rm H}\,\rm{\small I}$| emission from the SMC, observed with the Australian Square Kilometre Array Pathfinder (ASKAP) (Johnston et al. 2008; Hotan et al. 2021). These data are unprecedented in both spectral and angular resolution (see Section 2), and provide the ideal test bed for our new method of mapping the b parameter, particularly because 21-cm observations of the ISM provide a spatially coherent field in which we can probe the variations in density and velocity. Hitherto, no attempt has been made to spatially map the turbulence driving parameter over a large region, such as an entire galaxy that may contain an array of different physical turbulence driving mechanisms. Mapping the turbulence driving parameter requires high spatial and spectral resolution in order to accurately recover the density and velocity statistics, which makes |${\rm H}\,\rm{\small I}$| particularly useful for this large-scale mapping technique.

The SMC is of particular interest because it has a rich dynamic history of interactions with the LMC and the MW. It is a disrupted, star-forming, low-metallicity, and irregular dwarf galaxy (Kallivayalil et al. 2013). There are, therefore, many physical mechanisms present by which to drive turbulence in the ISM across the SMC, providing a rich environment to test our method for mapping the b parameter, and correlating our findings with known turbulence drivers.

The rest of this work is organized as follows. Section 2 introduced the data used in this work. Section 3 introduces our method for mapping the turbulence driving parameter, describing how we reconstruct the 3D density and velocity dispersions in large-scale observational data sets. Our results are presented in Section 4, and we discuss potential correlations of the b parameter with |${\rm H}\,\rm{\small I}$| and Hα emission in Section 5. Section 6 discusses our SMC measurements of the b parameter in relation to other environments presented in the literature. We summarize our conclusions and pathways to future work in Section 8.

2 OBSERVATIONS

The observations of the SMC were obtained as part of the pilot phase of the atomic neutral hydrogen component of the Galactic ASKAP survey (GASKAP-HI), which aims to reveal the structure, kinematics, and thermodynamics of |${\rm H}\,\rm{\small I}$| in the MW and Magellanic System at angular resolutions a factor of 3 to 30 times finer than existing 21-cm surveys of the Southern sky. The |${\rm H}\,\rm{\small I}$| data cube,¹ resulting from 20.9 h of total integration with ASKAP, used for our analysis is the most sensitive (rms brightness temperature |$\sigma _{T} = 1.1\, \mathrm{K}$| per |$0.98\, \mathrm{km\, s}^{-1}$| spectral channel) and detailed view (30 arcsec restoring beam or |$\sim 10\, \mathrm{pc}$| at the distance of the SMC) of the |${\rm H}\,\rm{\small I}$| associated with the SMC made to date. The details on the GASKAP-HI imaging pipeline, including the data validation, calibration, imaging, and combination with observations from the 64-m Parkes single dish telescope (Murriyang) to fill in the low spatial frequencies filtered out by ASKAP, are discussed in section 3 of Pingel et al. (2022).

2.1 Moment maps

Our analysis pipeline uses the moment-0 and moment-1 maps of some given position–position–velocity (PPV) cube as input, and processes them in such a way as to recover the 3D (volume) density dispersion and the Mach number, and thus the ratio of those two quantities – b, as defined in equation (1). The moment-0 (M0) is the integrated intensity, given by

where T_b(v) is the brightness temperature (in Kelvin) in the channel, and dv is the channel spacing. The moment-1 (M1) is the intensity-weighted centroid velocity, and is given by

where v is the velocity of each channel. We have |${\rm d}v=0.98\, \mathrm{km\, s}^{-1}$|⁠, and integrate from |$v=60.5$| to |$235.0\, \mathrm{km\, s}^{-1}$|⁠.

By excluding any pixels below some multiple of the rms brightness temperature per channel prior to integration to create M0 and M1, we can be confident that the data we include in our analysis are robust and we do not introduce any spurious effects to our calculation of the turbulence driving parameter due to noisy pixels. For our data, we choose a 10σ_T cut per channel, i.e. we only consider data with an signal-to-noise ratio (SNR) ≥10. We discuss the effect of the choice of signal-to-noise threshold on our results in Appendix  A.

2.1.1 Influence of the multiphase |${\rm H}\,\rm{\small I}$|

Equation(1) was derived for isothermal, hydrodynamic gas. |${\rm H}\,\rm{\small I}$| emission does not originate from isothermal gas. The temperature structure of multiphase neutral |${\rm H}\,\rm{\small I}$| can be roughly categorized into cold neutral medium (CNM; |$T \sim 100\, \mathrm{K}$|⁠), lukewarm neutral medium (LNM; 100 ≲ T/K ≲ 6000), and WNM (⁠|$T \gtrsim 6000\, \mathrm{K}$|⁠) (Wolfire et al. 1995, 2003; Kalberla & Haud 2018). The WNM is the most diffuse phase and has the largest volume-filling factor, while the CNM is expected to exist as smaller structures (such as clumps, sheets or filaments) dispersed throughout it (Clark, Peek & Miville-Deschênes 2019). As such, the emission from |${\rm H}\,\rm{\small I}$| originates from a mixture of CNM, LNM, and WNM. In order to apply equation (1) to the |${\rm H}\,\rm{\small I}$| emission data, we must make some careful assumptions about which phase dominates the moment-0 and moment-1 maps derived from the emission, which are the inputs to our analysis pipeline.

The relative mass fractions of CNM, LNM, and WNM will impact column density and the centroid velocity. For instance, narrow peaks in the spectrum due to the presence of CNM will shift the velocity centroid towards the centroid of the narrow peak, rather than measuring the centroid velocity of only the WNM. The strength of this effect, however, depends on the relative intensities of these CNM peaks with respect to the WNM component(s). The column density will also be affected by the presence of CNM, as the emission may be underestimated due to |${\rm H}\,\rm{\small I}$| self-absorption (HISA) by the cold gas (discussed further in Appendix  B). Furthermore, the statistics associated with these quantities are also affected by the ratio of WNM/LNM to CNM. If we imagine that the |${\rm H}\,\rm{\small I}$| structure is such that cold clumps/clouds of CNM are embedded in the diffuse WNM, then the dispersion of the centroid velocity will recover the intraclump velocity dispersion, rather than the internal dispersion of the cold clumps themselves. Because the clumps are moving within the WNM, we can assume that this intraclump velocity dispersion also traces the WNM velocity dispersion (Kobayashi et al. 2022; Mohapatra et al. 2022), which is what we primarily aim to measure. We also note that while the CNM may introduce a distortion in the individual line profiles, we only require a reasonably good estimate of the intensity-weighted average velocity along the line of sight (LOS, i.e. the centroid velocity) for the velocity analysis (below in Section 3.4, which is not expected to be strongly affected by individual features in the line profiles. The effect that the CNM has on the column density and centroid velocity and associated statistics is a function of the relative CNM mass fraction. In the SMC, the CNM mass fraction is about 11 per cent (Dempsey et al. 2022), meaning that the majority of the neutral |${\rm H}\,\rm{\small I}$| is WNM/LNM. Furthermore, because the SMC is relatively metal-poor (Russell & Dopita 1992), we expect that the portion of the gas that is not CNM is mostly comprised of WNM, as metal-line cooling will be less efficient in low-metallicity environments, as will photoelectric heating (Bialy & Sternberg 2019), and therefore the amount of LNM is likely much smaller than the WNM. At such low CNM percentages, the effects described above are relatively small, so we can assume that the emission data we use throughout this work is representative of the WNM, while acknowledging that the column density and centroid velocity maps are a product of a mixture of WNM, LNM, and CNM.

3 METHODS

Here, we present the analysis pipeline developed for producing a map of the turbulence driving parameter on galactic scales using the moment-0 and moment-1 maps of PPV data as input. In this section, we will describe each step in detail, but we begin with a brief summary of the steps in the pipeline:

The turbulence driving parameter b in equation (1) is constructed from the volume density dispersion |$\sigma _{\rho /\rho _0}$| and the turbulent sonic Mach number |$\mathcal {M}$|⁠. To recover the volume density dispersion, we first compute the column (2D) density dispersion from the moment-0 map and use it to reconstruct the volume (3D) density dispersion via the Brunt method (Section 3.3). The Mach number is a function of the 3D velocity dispersion and the sound speed, and we again start by finding the velocity centroid dispersion using the moment-1 map, and then extrapolating it to 3D (Section 3.4). To compute |$\mathcal {M}$|⁠, we then divide the 3D velocity dispersion by the sound speed of the gas (Section 3.5). In order to isolate the true density and velocity fluctuations without contamination from non-turbulent motions due to bulk rotation and/or large-scale hierarchical density structure, we employ a gradient-correction method in our calculations of the density and velocity dispersion (Section 3.2). This sequence of calculations is performed inside a Gaussian-weighted roving kernel (Section 3.1), so as to build spatial maps of the 3D density variance, Mach number, and b parameter.

3.1 Roving kernel

The density contrast and velocity fluctuations are scale-dependent quantities. Measuring the standard deviation of these quantities across the plane-of-the-sky requires some minimum number of spatial resolution elements. Our approach is to move a circular kernel across the input maps, pixel-by-pixel, calculate the dispersion within the kernel and assign that value to the central pixel, thus building up a map of the dispersion of the quantity from the input map (e.g. moment-0 or moment-1). The pixels of the resultant map are therefore not independent measurements and are correlated with one another on the scale of the diameter of the roving kernel.

In this study, we choose a Gaussian kernel, defined by

where i and j denote the central pixel of the kernel, x and y are the position of each pixel inside the kernel, and σ is the full width at half-maximum (FWHM) of the kernel divided by 2.355. We choose the number instrument beams per kernel FWHM to be 10, a quantity that is referred to throughout this text as N_bpk = D_kernel/D_beams, where D_kernel is the FWHM of the kernel, and D_beams is the diameter of the instrument beam. This gives a sufficiently large number of independent data points of the centroid velocity and column density to reliably calculate the dispersion in those quantities (see Sharda et al. 2018, for a discussion on the choice of the number of resolution elements per kernel; and a detailed study on the effects of the kernel size is presented below in Section 4.2.1). We apply the Gaussian kernel to a cut-out of a square region three times the FWHM of the kernel, which truncates the Gaussian. The size of the kernel is a flexible parameter in our pipeline and should be chosen as a function of the beam and pixel size of the input data. The kernel is also weighted by the primary beam response of the mosaic – i.e. the sensitivity across the field of view – so that pixels with lower sensitivity contribute less to the calculations of the dispersion. This weighting is an optional input to the pipeline, and must be the same shape as the moment-0 and moment-1 images, and should be a grid of values between zero and one, which is then multiplied with the Gaussian to give a grid of weights that can be used to weight the dispersion calculations.

3.2 Gradient subtraction

In calculating both the column density dispersion and the centroid velocity dispersion, we apply the gradient-subtraction method described in Federrath et al. (2016) and Stewart & Federrath (2022). We refer the reader to those works for a complete discussion on the details of the method, which we summarize in relation to this work here. Our aim is to compute the turbulent dispersion of the centroid velocity and column density maps, and we therefore must eliminate any systematic gradients in motion or intensity which are caused by other factors besides turbulence, primarily the hierarchical column density structure (e.g. galaxies tend to be denser in their centres), and any bulk rotation of the gas. This method was previously applied in Federrath et al. (2016), Sharda et al. (2018, 2019), and Menon et al. (2021), but only to the velocity field. Here, we apply the method to both density and velocity as large-scale gradients are present in both quantities (seen in Fig. 1).

To do this, we fit a linear gradient to the field, on the scale of the kernel in which we perform our calculations, then subtract it from the original field. Formally, for some quantity q(i, j) with pixel coordinates i and j, we define a gradient

where a, b, c are fit parameters and (i, j) is a point on the plane of the sky. After fitting for a, b, and c, we can compute the gradient-subtracted field,

For our purposes q is either the normalized logarithmic column density log₁₀(N_HI/N₀), where N₀ is the mean column density inside the kernel, or the intensity-weighted velocity centroid v.

3.3 Density dispersion

To calculate the 2D column density dispersion, we take the moment-0 map as input (in this case, the left-hand panel of Fig. 1), and move the Gaussian kernel across it, computing the following steps for each pixel. We first normalize the column density by the mean column density inside the kernel (left-hand panel of Fig. 2) and then fit a gradient to the normalized map, using equation (5) (middle panel of Fig. 2). We then subtract the gradient from the original normalized map via equation (6), so that we are left with the turbulent density variations (right-hand panel of Fig. 2). We then calculate the standard deviation of the resultant map, which gives the 2D column density dispersion, |$\sigma _{N_{\mathrm{HI}}/N_0}$|⁠.

In Fig. 2, we have fitted Gaussians to the PDFs in the lower panels. The raw data, in the left-hand panel, clearly do not follow the anticipated Gaussian distribution, but after applying the gradient subtraction we see in the right-hand panel that the resultant distribution well approximates a Gaussian, which is a hallmark of compressible isothermal turbulence.²

3.3.1 Conversion from 2D to 3D density dispersion

To convert the 2D density dispersion |$\sigma _{N_{\mathrm{HI}}/N_0}$| to the 3D density dispersion, |$\sigma _{\rho /\rho _0}$|⁠, we follow the method outlined in Brunt et al. (2010). This method is based on reconstructing the power spectrum of the volume density contrast from the power spectrum of the column density contrast, by measuring the column density power spectrum and extending its 2D power into 3D space. The relation between the 2D density power spectrum P_2D(k) and the 3D density power spectrum P_3D(k) is given by (Federrath & Klessen 2013),

where k is the wave number. We exploit this relation to recover P_3D(k).

We first compute P_2D(k) of the gradient-subtracted column density, (N_HI/N₀ − 1 shown in Fig. 3), which immediately gives us P_3D(k) of the quantity ρ/ρ₀ as per equation (7). The ratio of the sums over k-space of these two quantities gives the density dispersion ratio (Brunt et al. 2010),

known as the ‘Brunt Factor’, which then allows us to recover the volume density dispersion, |$\sigma _{\rho /\rho _0}$|⁠.

Equation (8) relies on the assumption that the input field is isotropic in k-space, and that therefore P_2D(k) is also isotropic. This does not assume that the input field is isotropic in real space, only that its power spectrum is symmetric and isotropic, so we must check the Fourier image to make sure that this assumption holds, and we are not introducing further uncertainty in our measurements. The P_2D(k) image is shown in Fig. 3, where we can see that the power spectrum has good point symmetry. There are some angle-dependent features in the Fourier image. However, overall the distribution of power can be approximated as point symmetric around the k = 0 mode (centre). The sampling of the power spectrum follows the spatial sampling of the data. Noise is accounted for through the SNR threshold applied (cf .Section 2). The effect of the beam size is accounted for in that it is primarily the low-k modes (i.e. scales larger than the beam size) that contribute to the total power of the spectrum. We have also checked other regions (kernels) of the SMC and find qualitatively similar results. This implies that spherical symmetry is a good approximation for P_2D(k), such that the distribution of the density variations in real 2D space may be similarly distributed in 3D space.

3.4 Velocity dispersion

The 3D turbulent velocity dispersion is defined as

which cannot be measured in PPV space as we do not have access to the two velocity components in the plane of the sky. To recover σ_{v, 3D} in a given kernel, we follow the methods developed by Stewart & Federrath (2022), who show that σ_{v, 3D} can be recovered from PPV space using the standard deviation of the gradient-corrected moment-1 map together with a correction (1D/2D to 3D conversion) factor. As with our method for computing the density contrast, we apply a gradient correction to the moment-1 map that captures large-scale systematic motions (e.g. rotation) before computing the centroid velocity variance. First, we apply the kernel to the moment-1 map (right-hand panel of Fig. 1). Next, we fit and subtract the gradient, and then take the standard deviation within the kernel to get the gradient-subtracted centroid dispersion σ_{v, 1D}. This process is illustrated in Fig. 4.

Subtracting the large-scale gradients from the velocity map isolates the turbulent motions, however, the variance of the gradient-corrected map is not a true representation of the 3D turbulent velocity dispersion, because it does not take into account the line-of-sight dispersion. Therefore, we multiply by a correction factor³ of |$C_\mathrm{(c-grad)}^\mathrm{\, any}=3.3\pm 0.5$|⁠, which was determined based on synthetic observations of 3D hydrodynamical simulations of rotating, turbulent clouds (Stewart & Federrath 2022) to convert the centroid velocity dispersion into a 3D turbulent velocity dispersion σ_{v, 3D}. The choice of using only the moment-1 map as opposed to using the moment-2 (or the ‘parent’ velocity dispersion, which is a combination of the moment-1 and moment-2; see Stewart & Federrath 2022) is discussed further in Section 7.2.

3.5 Mach number

The sonic Mach number (⁠|$\mathcal {M}$|⁠) of the turbulent component of the velocity field is given by

To construct a map of this quantity, we need the velocity dispersion σ_{v, 3D} (described above) and an estimate of the sound speed c_s. The last step in the pipeline is to divide σ_{v, 3D} by the sound speed to produce the quantity |$\mathcal {M}$|⁠. Depending on available data, the sound speed could be input to our analysis pipeline as either a constant or as a spatially varying parameter (if spatially resolved information about the temperature or sound speed is available). As we do not have access to spatially varying data for the sound speed, we choose a constant speed defined as

where γ = 5/3 is the adiabatic index, k_B is the Boltzmann constant, T is the gas temperature, μ is the mean particle weight, and m_H is the mass of a hydrogen atom. Because we do not have a way to estimate the temperature of the combined |${\rm H}\,\rm{\small I}$| phases present in our data, we will assume that the phase which dominates the temperature of the neutral gas is the WNM, and so we use its molecular weight of μ = 1.4 (Kauffmann et al. 2008). In this study, we adopt |$T = (7.0 \pm 1.0) \times 10^3\, \mathrm{K}$|⁠, which gives a sound speed of |$c_\mathrm{s}= (8.3 \pm 0.6) \, \mathrm{km}\, \mathrm{s}^{-1}$|⁠. Estimates of the WNM temperature in the MW range from about 4000 K to |$10^4\, \mathrm{K}$|⁠, depending on the method used (Wolfire et al. 2003; Murray et al. 2014, 2018; Marchal & Miville-Deschênes 2021). Bialy & Sternberg (2019) use simulations to investigate temperature and pressure structures in atomic neutral gas for varying metallicities. They find that the temperature structure changes for metallicity values smaller than |$0.1\, Z_\odot$|⁠. The SMC is more metal poor than the MW (Z ∼ 0.2Z_⊙, Russell & Dopita 1992), and therefore choosing a temperature in line with solar metallicity values in accordance with Bialy & Sternberg (2019) seems a reasonable assumption. However, ultimately the temperature only enters with a square-root dependence in the sound speed, and therefore in the Mach number; hence, these assumptions do not introduce any major source of uncertainty in our calculations (discussed further in Sections 6 and 7).

4 RESULTS

In this section, we present and discuss the output of the analysis pipeline described in Section 3 as applied to the GASKAP-HI emission cube of the SMC. A summary of all the relevant physical parameters and measurements for the SMC is provided in Table 1.

4.1 Spatial distribution of volume density dispersion and turbulent Mach number

In the left-hand panel of Fig. 5, we present the volume density dispersion map, i.e. the plane-of-the-sky variations of |$\sigma _{\rho /\rho _0}$|⁠, following the methods described in Section 3.3. The quantity |$\sigma _{\rho /\rho _0}$| measures the turbulent density variations in equation (1). We see that the density dispersion varies by about an order of magnitude across the map. Within the analysis contour, the variations are about a factor of ∼3. We see a slight tendency of higher dispersion values towards the edges of the main contour (where the emission density begins to drop off), but overall the density dispersion does not show distinct regions of low or high values, with exception of the region at the top of the Wing and Bar regions. The relatively small variation and the lack of any large-scale global gradients in density dispersion within the main body of the SMC is a result of the gradient-subtraction method, which, although calculated on the kernel scale, successfully accounts for large-scale gradients (also) on the global scale. This is clear when visually comparing with the left-hand panel of Fig. 1, where column density gradients towards the centre and bar of the SMC are visible. Given that these gradients would be expressed in the column density dispersion, and the volume density dispersion is directly related to that quantity (see Section 3.3.1), we would expect to see large-scale gradients in the volume density dispersion if the gradient subtraction method had not accounted for them. In summary, using the gradient-subtraction method, we have isolated the overall turbulent density fluctuations, which enter equation (1).

The right-hand panel of Fig. 5 shows the map of the turbulent sonic Mach number, |$\mathcal {M}$|⁠, following the methods in Section 3.5. The variations within the analysis contour are ∼3, similar to the variations of |$\sigma _{\rho /\rho _0}$|⁠. The Mach number distribution is also relatively uniform within the analysis contour, with notable regions of high |$\mathcal {M}$| at the top of the Wing region, where it intersects with the Bar region, similar to |$\sigma _{\rho /\rho _0}$|⁠. Comparing the |$\mathcal {M}$| map to the right-hand panel of Fig. 1, we can see that the gradient-subtraction method has successfully removed the large-scale rotation of the SMC. We find that the Mach numbers are distinctly subsonic within the main analysis contour, with typical values of |$\mathcal {M}\sim 0.1$|–0.7. Turbulent Mach numbers of this magnitude are expected for the WNM (Burkhart et al. 2010), where the gas remains largely subsonic, while the CNM usually exhibits trans-sonic to mildly supersonic turbulence with |$\mathcal {M}\sim 1$|–2 (Vázquez-Semadeni et al. 2006; Hennebelle & Audit 2007; Heitsch et al. 2008). Burkhart et al. (2010) derived a spatial map of the sonic Mach number based on lower resolution |${\rm H}\,\rm{\small I}$| column density, and found that 90 per cent of the |${\rm H}\,\rm{\small I}$| in the SMC was sub- or trans-sonic. While the primary beam of the data used in that work was much larger than the primary beam in the GASKAP-HI data, or the kernel FWHM we use in this work, the spatial distribution of the variations in the Mach number is in agreement. Higher Mach numbers towards the edges of the Wing, and surrounding the lower portion of the Bar are features in Fig. 5 of this work and fig. 8 in Burkhart et al. (2010). The difference in absolute values of the Mach number is likely due to a difference in the methods used, as well as the difference in resolution of the two data sets. Increasing the kernel size (and thus mimicking the lower resolution in Burkhart et al. 2010) increases the Mach number range we recover in this work (discussed further in Section 4.2.1), and pushes the values up into the trans-sonic/supersonic regime.

Both maps in Fig. 5 show some gradients at the edges of the fields, which are a result of lower emission density in those regions, and therefore lower SNR per-channel (below SNR = 10), as well as lower instrument sensitivity (40–80 per cent lower than inside the contour), which is why we have chosen to analyse only the data in the first closed contour.

4.2 Spatial distribution of the turbulence driving parameter

Combining the maps in Fig. 5 via equation (1), we compute the turbulence driving parameter, b, for the whole SMC, which is shown in Fig. 6. It is immediately clear that there are spatial variations in b. Within the main contour, b varies between ∼0.3 and 1, i.e. between purely solenoidal and purely compressive driving of the turbulence. Regions towards the Stream and the upper parts of the Bar seem to exhibit more compressive driving, while the central regions of the SMC tend towards more solenoidal and mixed (1/3 < b < 0.38) driving. The key result from this map is that we clearly find spatial variations in the turbulence driving mode. The exact cause of these variations is difficult to determine and requires detailed matching of the new b-parameter map obtained here, with information about potential physical drivers of turbulence (Elmegreen 2009; Federrath et al. 2017), such as (1) feedback, i.e. from star formation/evolution activity, including supernovae, winds, jets/outflows, and radiation, or (2) dynamics of the Magellanic system, which includes accretion onto the SMC, large-scale flows inside the SMC, such as induced by rotation or shear, or tidal interactions with the hot Galactic halo causing ram pressure stripping. Disentangling and cross-matching all of these effects is challenging and will be subject of future work. However, we start investigating some of these potential correlations in Section 5, by studying the turbulence driving parameter as a function of the |${\rm H}\,\rm{\small I}$| and Hα emission in the SMC.

4.2.1 Influence of the kernel size

The 3D velocity dispersion and the volume density dispersion are scale-dependent quantities (Kim & Ryu 2005; Kowal & Lazarian 2007; Kritsuk et al. 2007; Federrath, Klessen & Schmidt 2009; Beattie et al. 2019; Federrath et al. 2021). In order to investigate the influence of the size of the analysis kernel, we compute the four main analysis quantities for 5 different kernel FWHMs, such that N_bpk = 2.5, 5.0, 10.0, 20.0, and 40.0. The panels of Fig. 7 show the driving parameter for these five values of the kernel size, and highlight the different levels of spatial granularity resulting from each kernel size. All five kernel sizes, we test are small compared to the size of the SMC. This should always be the case when applying this method to observations of entire galaxies, as the method is designed to probe the (small-scale) 3D turbulence properties.

In Fig. 8, we show PDFs of the four main quantities as a function of the kernel size (from the left- to right-hand panel): the Brunt factor |$\mathcal {R}^{1/2}$|⁠, the density dispersion |$\sigma _{\rho /\rho _0}$|⁠, the Mach number |$\mathcal {M}$|⁠, and the driving parameter b. Fig. 8 demonstrates that the width and peak of the |$\mathcal {M}$| and |$\sigma _{\rho /\rho _0}$| PDFs increase with increasing kernel size, and that the Brunt factor (⁠|$\mathcal {R}^{1/2}$|⁠) decreases with increasing kernel size. Because the width of the |$\mathcal {M}$| and |$\sigma _{\rho /\rho _0}$| PDFs scale with kernel size in roughly the same way, the turbulence driving parameter b is not a strong function of the kernel size (see right-hand panel of Fig. 8). However, it is clear that the kernel size plays a role in the distribution and peak of all the analysis quantities, so in Fig. 9, we plot the median value of the PDFs shown in Fig. 8 for each analysis quantity as a function of N_bpk, to study the convergence behaviour of each of the relevant quantities.

Fig. 9 shows how the median value of each analysis quantity behaves as a function of N_bpk. The top two panels show the Brunt factor and the driving parameter, which both converge to a constant value as the kernel size increases. We find that the best fit for the Brunt factor (top left-hand panel of Fig. 9) is

such that in the limit of an infinitely large kernel (N_bpk → ∞) the Brunt factor converges to |$\mathcal {R}^{1/2}= 0.31\pm 0.01$|⁠. This value of roughly 0.3 is consistent with previous findings (Brunt 2010; Ginsburg et al. 2013; Federrath et al. 2016; Menon et al. 2021; Sharda et al. 2022).

The driving parameter (bottom right-hand panel of Fig. 9) converges to a value of b = 0.51 ± 0.01, in the limit of an infinitely large kernel, which tends towards the compressive driving end of the scale (b > 0.38). The functional form of the fit is given by

We take these two converged values to be the overall median Brunt factor and driving parameter in the SMC. The Mach number and volume density dispersion are not expected to converge, as they are scale-dependent quantities. We therefore fit power-law functions to each of these quantities. For the volume density dispersion, we derive

and for the Mach number we find

The exponents of these two power laws agree within the uncertainties, which shows that |$\sigma _{\rho /\rho _0}$| and |$\mathcal {M}$| do indeed change in the same way with increasing kernel size, which can also be seen in Fig. 8. Thus, b converges to a constant value with increasing kernel size.

Because |$\mathcal {M}$| and |$\sigma _{\rho /\rho _0}$| do not converge, we must necessarily choose one value of N_bpk when presenting the maps of these quantities and b. As previously outlined, we chose N_bpk = 10, because it provided a sufficient level of granularity and had enough resolution elements per kernel to accurately recover the main analysis quantities (in particular the driving parameter b). We can see in Fig. 9 that this choice of N_bpk recovers the converged value of b within about 2.5  per cent.

Fig. 10 shows the PDFs of |$\mathcal {R}^{1/2}$|⁠, |$\sigma _{\rho /\rho _0}$|⁠, |$\mathcal {M}$|⁠, and b inside the main analysis contour for the standard kernel size of N_bpk = 10. We find |$\mathcal {R}^{1/2}= 0.34\pm 0.01$|⁠, |$\sigma _{\rho /\rho _0}= 0.18^{+0.06}_{-0.04}$|⁠, and |$\mathcal {M}= 0.35^{+0.16}_{-0.09}$|⁠. As expected for the WNM, the Mach numbers lie in the subsonic regime. For the turbulence driving parameter across the SMC, we find |$b=0.49^{+0.22}_{-0.13}$|⁠. However, there are substantial variations of b from region to region (as quantified by the 16th and 84th percentile range). Nevertheless, even the 1σμ lower limit of b, i.e. the 16th percentile value is with b₁₆ = 0.36 still in the regime of a natural mixture of compressive and solenoidal driving. This is an interesting result as it may indicate that the |${\rm H}\,\rm{\small I}$| gas is subject to predominantly compressive turbulence driving mechanisms in the SMC.

5 CORRELATIONS BETWEEN THE TURBULENCE DRIVING PARAMETER AND THE GAS DENSITY AND/OR STAR FORMATION ACTIVITY

In the left-hand panel of Fig. 11, we investigate whether there is any correlation between |${\rm H}\,\rm{\small I}$| intensity and the turbulence driving parameter b within the main contour region of the SMC. We do not find evidence of a correlation between b and |${\rm H}\,\rm{\small I}$| intensity, instead the turbulence driving seems to be in the mixed-to-compressive (b > 0.4) regime regardless of the |${\rm H}\,\rm{\small I}$| emission density. This is expected since we found b to be compressive overall in Fig. 10. It also shows that b is not simply a function of the column density, so one cannot simply assume that the turbulence is more compressively driven in regions of higher column density, in the case of |${\rm H}\,\rm{\small I}$|⁠. Similarly, we find that b is slightly compressive across the entire range of H α (MCELS Smith et al. 1999; Winkler et al. 2015) intensity, as shown in the right-hand panel of Fig. 11. To first order, the presence of Hα signifies active, massive star formation, which may provide a feedback mechanism by which to drive turbulence compressively, through winds, expanding |$\rm{H{\small II}}$| regions (Menon et al. 2021), outflows and ultimately supernovae (Elmegreen 2009; Federrath et al. 2017). Compressive driving also positively influences star formation rates (Federrath & Klessen 2012), and so it is difficult to disentangle whether compressive driving where H α is present is causing star formation, or whether star formation is causing compressive driving. Because we do not find any correlation with increased H α intensity and more compressive driving, it is possible that the overall compressive nature of the driving in the |${\rm H}\,\rm{\small I}$| gas in the SMC is not dominated by the star-formation activity.

6 COMPARISON TO OTHER MEASUREMENTS OF THE TURBULENCE DRIVING PARAMETER IN DIFFERENT ENVIRONMENTS

Fig. 12 concludes the discussion of our results by presenting a selection of other observational measurements of the density dispersion – Mach number relation in regions of the MW, as well as one molecular star-forming region in the LMC. Here, we can see that our data for the SMC lie between the lines denoting mixed and fully compressive driving, and that the points lie in the subsonic and low-density variance regime. With reference to our discussion in Section 4.2.1, however, it should be noted that the points we present on this figure (both our work and previous studies) are all a function of the scale (and selected region) on which they have been measured, which means that |$\mathcal {M}$| and |$\sigma _{\rho /\rho _0}$| change with the kernel or region size observed. We show the median values for the five kernel sizes investigated in Section 4.2.1, corresponding to a physical size of 25, 50, 100, 200, and 400 pc. We can see that the smallest kernel (25 pc) is squarely in the solenoidal regime (however, this value is not converged; see Fig. 9), and as the kernel size increases, compressive driving becomes more dominant (and leads to a converged value of b; cf. Fig. 9).

Included in Fig. 12 are star-forming molecular clouds such as the Pillars of Creation (Menon et al. 2021), Taurus (Brunt 2010; Kainulainen & Tan 2013), IC5416 (Padoan et al. 1997), and the Papillon Nebula (Sharda et al. 2022), which all exhibit supersonic, compressively driven turbulence, likely driven by star formation feedback. There are also two molecular clouds that are not star-forming; one which exhibits supersonic mixed/compressively driven turbulence – GRSMC 43.30-0.33 (Ginsburg et al. 2013), and another, which is supersonic and solenoidally driven – ‘The Brick’ (Federrath et al. 2016). The kind of turbulence driving in each of these molecular clouds depend on the specific physical mechanisms at play in each region, and while it is interesting to note that a variety of b values are recovered and can be correlated to the environmental conditions (for instance, strong shearing motions giving rise to predominantly solenoidal driving in the Brick Federrath et al. 2016), it is difficult to do so for the present SMC measurements, as these contain contributions from the entire galaxy in the data shown in this plot. Correlating b with environmental conditions will ultimately involve studying the spatial variation of b as shown in the map of Fig. 6, and linking it to other measurements that provide information about the dynamics and potential physical drivers of turbulence, as discussed in the preceding sections.

The most direct comparison we can make with our data is the MW WNM datapoint (Marchal & Miville-Deschênes 2021) (lilac triangle in Fig. 12). Marchal & Miville-Deschênes (2021) found |$\sigma _{\rho /\rho _0}= 0.6 \pm 0.2$| and |$\mathcal {M}= 0.87 \pm 0.15$|⁠, both somewhat higher than our measurements for the SMC. Because the density variation and the Mach number are scale-dependent quantities (see previous discussion; and Figs 7–9), we must first consider whether the spatial scale on which we measure our quantities is comparable to the scales on which the quantities in the MW were measured. We chose to use a kernel N_bpk = 10, which at the distance to the SMC is |$\sim 100\,$| pc. In Marchal & Miville-Deschênes (2021), the spatial scale on which they measure |$\sigma _{\rho /\rho _0}$| and |$\mathcal {M}$| is |$\sim 63\,$| pc. Referring to Fig. 9, we can see that had we used a kernel N_bpk = 5 (⁠|$\sim 50\,$| pc), we would have derived even smaller median values for these two quantities. It is therefore likely that the difference between the measured SMC and MW values is not a result of the scale on which they were measured.

We have assumed that the WNM temperature in the SMC is |$T = 7000\pm 1000\, \mathrm{K}$|⁠, which is about 1000 K higher than the temperature used in Marchal & Miville-Deschênes (2021), but lower than the 10⁴ K observed by Murray et al. (2018). There is a wide range of temperature estimates for the WNM in Galactic |${\rm H}\,\rm{\small I}$|⁠, and we have chosen a temperature range in keeping with observations from Murray et al. (2014), following results from Bialy & Sternberg (2019), who show that at the metallicity of the SMC, the temperature structure should be similar to that of the MW. This influences our |$\mathcal {M}$| values to be about 10 per cent lower than if we had used the MW temperature estimate from Marchal & Miville-Deschênes (2021). However, taking our lowest estimate for the temperature in the SMC and the highest estimate for the temperature in the MW, our median Mach number value still does not result in any overlap with the Mach number estimate for the WNM in the MW. We conclude that this cannot account completely for the difference between our Mach number values and those estimated in Marchal & Miville-Deschênes (2021). It is likely that because we do not perform a phase separation as in Marchal & Miville-Deschênes (2021), and assume that the emission is dominated by WNM, this is a contributing factor in our differing result, but still may not account for it entirely.

The difference between our median value for the volume density dispersion as compared to the reported MW value could be explained by significant variation in the depth of the SMC. If the SMC is highly extended along the LOS, it is possible that we have underestimated the volume density dispersion via the Brunt et al. (2010) method (Section 3.3.1). We discuss this issue further in Section 7.3, but considering that we do not have a robust estimate for the extent of the SMC in the third dimension, we can only assume that our measured column density dispersion is a reasonable representation of the dispersion along the LOS, and therefore |$\sigma _{\rho /\rho _0}$| is truly smaller than in the MW WNM.

In summary, our analysis of the SMC WNM in comparison to the MW WNM region studied by Marchal & Miville-Deschênes (2021) may indicate that the values of density dispersion and Mach number are indeed physically smaller by a factor of ∼2–3 in the SMC compared to the MW, but the average turbulence driving parameter of b ∼ 0.7 for the MW and b ∼ 0.5 for the SMC, indicates a similar dominance of compressive turbulence driving in the WNM of both the MW and the SMC.

7 CAVEATS

7.1 Magnetic fields

The ISM is ubiquitously magnetized, and the influence of magnetic pressure on the density dispersion – Mach number relation has been derived theoretically and investigated in simulations (Padoan & Nordlund 2011; Molina et al. 2012; Körtgen 2020) and observations (Federrath et al. 2016; Kainulainen & Federrath 2017). We have assumed the purely hydrodynamical relation in this work, as we are unable to map the magnetic field strength of the SMC in a way that would allow us to meaningfully incorporate it into our calculation of b. From Molina et al. (2012), we know that for cases in which the magnetic field strength is proportional to the square root of the density, b is given by

where plasma beta, |$\beta = P_{\mathrm{th}}/P_{\mathrm{mag}} = 2c_\mathrm{s}^2/v_{\mathrm{A}}^2$|⁠, is the ratio of the thermal to magnetic pressure, and the Alfvén speed v_A = B_turb/(4|$\pi$|ρ₀)^1/2 for the turbulent component of the magnetic field, and ρ₀ is the mean volume density (e.g. Federrath & Klessen 2012; Federrath 2016). The best we can currently do is estimate a single value for β in the SMC, and therefore quantify its bulk effect on our spread of b values.

Our best estimate for the magnetic field strength across the SMC comes from Livingston et al. (2022), who study Faraday rotation towards 80 sources across the SMC to estimate the LOS magnetic field strength (⁠|$-0.3 \pm 0.1\, \mu$|G) and the random component of the field (⁠|$5^{+3}_{-2}\, \mu$|G). For the WNM, the number density is |$\sim 0.2 - 0.5\, \mathrm{cm^{-3}}$| (Ferrière 2001). This gives us an Alfvén speed between |$v_{\mathrm{A}} \sim 5 - 30 \, \mathrm{km\, s^{-1}}$|⁠. Combining this, we estimate β ∼ 0.1−2. Using this in equation (16), the factor (1 + β⁻¹)^1/2 ∼ 1−3, which means that the turbulence driving parameter could increase by up to a factor of three. Hassani et al. (2022) estimate that β < 1 in the WIM and HIM, and although they make no estimate for β in the WNM, it is in-keeping with our estimate. Using the magnetohydrodynamical version of the density dispersion – Mach number relation does not provide a particularly useful constraint in this instance, but would generally increase our values of b, pushing our results towards the compressive end of the driving spectrum. However, given that we do not have a map of the random magnetic field strength across the SMC at this time, the hydrodynamic approach used in our analysis is robust enough to provide a lower limit of the turbulence driving parameter in the SMC.

7.2 Calculating 3D velocity dispersion from the centroid velocity map

Stewart & Federrath (2022) discuss three methods to determine the 3D turbulent velocity dispersion of a cloud. They find that the gradient-corrected parent velocity dispersion, the sum in quadrature of the gradient-corrected moment-1 and moment-2 maps, is the most robust way of recovering the 3D turbulent velocity dispersion of a cloud. We initially attempted to use this method, but were unable to reliably disentangle the various components along the LOS in a given pixel, causing an overestimation of the contribution of moment-2 to the sum. In future work, we plan to explore new methods for decomposing multicomponent spectra, at which time we will revisit this aspect of the pipeline. However, the unprecedented resolution of the PPV cube in both velocity and position allows us to use the correction factor found in Stewart & Federrath (2022) (⁠|$C_\mathrm{(c-grad)}^\mathrm{\, any}=3.3\pm 0.5$|⁠) to recover the LOS velocity fluctuations from the moment-1 map only, to within their quoted 10  per cent accuracy, which then gives us the 3D turbulent velocity dispersion.

7.3 Depth of the SMC

The Magellanic system is highly dynamic and complex, and the 3D structure of the gas in the SMC is, for all intents and purposes, an unknown quantity. The robust gradient in the integrated velocity (see right-hand panel of Fig. 1) has in the past been interpreted as evidence that the SMC has a disc-like structure (Kerr, Hindman & Robinson 1954; Hindman, Kerr & McGee 1963; Stanimirović, Staveley-Smith & Jones 2004; Di Teodoro et al. 2019). However, the kinematics of gas-tracing stars, mapped using proper motions from Gaia, are inconsistent with a rotational disc model (Murray et al. 2019). 3D hydrodynamical simulations that attempt to reconstruct the dynamic history of the Magellanic system show that either one or two previous interactions between the SMC and LMC are required to consistently reproduce the Stream and the Bridge (Besla et al. 2012; Diaz & Bekki 2012; Lucchini, D’Onghia & Fox 2021).

It is likely that the SMC is not a rotating disc, but a torn and extended structure, with an estimated depth of 10–30 kpc (Tatton et al. 2021). This is not, however, a particular problem for our method. Because we use the Brunt et al. (2010) method (Section 3.3.1) to recover only the volume density dispersion, rather than any estimate for the volume density itself, we do not need to directly account for the depth over which we integrate. The assumption behind this method is that the gas density dispersion along the LOS is comparable to the density dispersion in the plane of the sky. We further find that the dispersion in the plane-of-the-sky is relatively isotropic (see Fig. 3), supporting the assumption in the Brunt et al. (2010) method. We may be estimating a lower limit on the volume density fluctuations present in each kernel if the SMC is highly extended along the LOS, given that integration over a large distance tends to wash out variance in the column density. However, the moment-1 map, used to quantify the turbulent velocity fluctuations in the plane-of-the-sky, is subject to similar smoothing along the LOS. Thus, both the moment-0 and the moment-1 maps are expected to be affected by LOS averaging in a similar way, and therefore, the ratio of the two (i.e. the turbulence driving parameter b; see equation 1) may be relatively robust against LOS averaging effects (similar to how it is independent of the kernel size on sufficiently large scales; cf. Figs 8 and 9). Without further exploration which is beyond the scope of this work and the available data, we cannot investigate this effect further at this time and therefore leave such an analysis to future work.

7.4 Relative uncertainties

As we have outlined in the above sections, there are several major assumptions made in this work, which likely dominate any errors associated with our results. However, under the assumptions we have made, we can quantify the relative error in our calculation of the turbulence driving parameter. We must account for the error introduced through the Brunt et al. (2010) method (⁠|$\sim 10~{{\ \rm per\ cent}}$|⁠) and the error associated with our estimate for the WNM temperature and therefore the sound speed (⁠|$\sim 10~{{\ \rm per\ cent}}$|⁠). The error introduced by our conversion of the velocity variance via the Stewart & Federrath (2022) method is |$\sim 10~{{\ \rm per\ cent}}$|⁠, which is similar to the range of applicable values, so the choice of conversion factor does not have a significant effect on our results. Combined, this gives a relative uncertainty of |$\sim 20~{{\ \rm per\ cent}}$| associated with our b values. Referring to Fig. 10, we can see that this error is smaller than the variance of b itself, and as such, we are primarily measuring the spatial variations of the turbulence driving parameter across the SMC, rather than variations introduced due to noise or uncertainties in the method.

8 CONCLUSIONS

In this work, we presented a new generalized method for creating a map of the turbulence driving parameter in the ISM of galaxies using column density and centroid velocity information. Our method can be applied on scales from molecular clouds to entire galaxies, provided there are enough spatially coherent resolution elements available (Section 4.2.1), and the data are sufficiently sensitive to provide reliable moment-0 and moment-1 maps. We use a Gaussian-weighted roving kernel (Section 3.1) to recover the volume density dispersion and turbulent Mach number, and construct maps of these quantities (Fig. 5) which we then use to create a map of the turbulence driving parameter (Fig. 6) via equation (1). In summary, the main steps of the pipeline performed in each instance of the roving kernel are as follows:

• to isolate only turbulent density variations, fit and subtract a linear gradient from the normalized column density map (see Fig. 2);.

• via the ‘Brunt Method’ (Section 3.3.1), use the 2D power spectrum of the gradient-corrected column density to construct the 3D power power spectrum, the ratio of which is the Brunt factor, |$\mathcal {R}^{1/2}$|⁠;

• using |$\mathcal {R}^{1/2}$|⁠, reconstruct the volume density dispersion from the column density dispersion (see Fig. 5);

• fit and subtract a linear gradient from the centroid velocity map, thus isolating only turbulent velocity fluctuations (see Fig. 4);

• find the dispersion of the gradient-corrected centroid velocity, and convert to 3D velocity dispersion via the ‘Stewart’ method (see Section 3.4);

• divide the 3D velocity dispersion by the sound speed to recover the turbulent sonic Mach number, |$\mathcal {M}$| (see Fig. 5); and

• finally, divide the volume density dispersion map by the Mach number map to create a map of the turbulence driving parameter, b, via equation (1) (see Fig. 6).

We have applied this method to high-resolution GASKAP-HI observations of 21 cm emission from the SMC (Section 2). We find that spatial variations of the turbulence driving parameter span the entire driving range from purely solenoidal to purely compressive driving, with some regions exhibiting very compressive driving (b ∼ 1, e.g. towards the Bridge and the Stream), while other regions are driven more solenoidally (b ∼ 1/3, e.g. in the centre). We find that the driving parameter is a weak function of the scale on which it is measured (see Fig. 9), but that it converges on kernel scales |$\gtrsim 100\, \mathrm{pc}$|⁠) to a constant value of b ∼ 0.5, which is towards the compressive end of the spectrum (i.e. b > 0.38, which defines the natural driving mixture), and with 16th to 84th percentile variations in b between ∼0.3 and ∼0.6 across the SMC. In the context of other measurements of the turbulence driving parameter in the WNM, specifically Marchal & Miville-Deschênes (2021), we find that while both the volume density dispersion and the Mach number are significantly lower than MW values, b is similarly compressive overall (∼0.7 in the MW). We do not find evidence for a correlation of b with either |${\rm H}\,\rm{\small I}$| or H α emission intensity. In future work, we will delve deeper into specific regions of the SMC and correlate variations in the b parameter with known physical turbulence driving mechanisms.

Acknowledgement

The authors acknowledge Interstellar Institute’s program ‘With Two Eyes’ and the Paris-Saclay University’s Institut Pascal for hosting discussions that nourished the development of the ideas behind this work. IAG would like to thank the Australian Government and the financial support provided by the Australian Postgraduate Award. CF acknowledges funding by the Australian Research Council (Future Fellowship FT180100495 and Discovery Projects DP230102280), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). CF further acknowledges high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi, and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grant ek9) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme, through which the data analyses presented in this paper were performed. SS and NMP acknowledge support provided by the University of Wisconsin – Madison Office of the Vice Chancellor for Research and Graduate Education with funding from the Wisconsin Alumni Research Foundation, and the NSF Award AST-2108370.

DATA AVAILABILITY

The code used in the analysis is freely available on GitHub. The GASKAP-HI emission cube used in this study is freely available to download at https://doi.org/10.25919/www0-4p48.

Footnotes

References

APPENDIX A: INFLUENCE OF THE SNR

As we previously outlined, we choose a high SNR cut per velocity channel of 10 to ensure the robustness of our data analysis. As defined in Section 2, the rms noise level is 1.1 K per 0.98 km s⁻¹ spectral channel. We apply multiples of this as a per-channel SNR threshold. We investigated three different SNR cuts: 3, 5, and 10. Fig. A1 demonstrates that the choice of SNR cut does not impact the analysis quantities particularly, because the data within our analysis contour has a very high SNR.

We can see this via visual inspection in Fig. A2, where the nonviable kernels that contain low SNR pixels are shown in white. Clearly, all the kernels inside the analysis contour contain pixels that have velocity channels above the SNR =10 threshold, which accounts for how similar the PDFs of the analysis quantities are. The slight variations are caused by some specific velocity channels in a given pixel being excluded by the per-channel SNR cut, resulting in less channels contributing to the integration that creates the moment-0 and moment-1 maps. This is the case in a minority of pixels towards the contour edge, which, in turn, result in the slight variations in the PDFs of the analysis quantities. Pixels which contain no velocity channels above the SNR threshold are not considered, and kernels containing such pixels are ignored entirely. This results in the lower SNR thresholds having more viable kernels than the higher SNR thresholds, as shown in Fig. A2.

APPENDIX B: CORRECTION FOR OPTICALLY THICK |${\rm H}\, \rm {\small I}$|

In this work, we chose to account for HISA so as to avoid underestimating the true column density. Using absorption data, the true optical depth of the gas can be estimated and a correction factor derived, such that

under the assumption that the gas is isothermal. In Fig. B1, we present PDFs of our primary analysis quantities with two different correction factors for the column density applied. The first is from Dickey et al. (2000), and takes the form,

such that the correction factor is applied to column densities above |$10^{21.4}\, \mathrm{cm}^{-2}$|⁠. An updated version of this correction factor was derived by Dempsey et al. (2022) using new ASKAP absorption data, and takes the form,

Despite the fact that the Dempsey et al. (2022) correction is only fitted to absorption data in the Bar of the SMC, we chose to apply it to our column density map none the less. While the correction factor in the Bar is not directly applicable to data from the Wing (and these are the two main regions enclosed by our analysis contour), a Wing-only correction is not available. By inspecting the data presented in fig. 10 of Dempsey et al. (2022), it seems likely that the fit parameters in such a correction would not be too dissimilar to those in the Bar correction, as there is a substantial amount of overlap between the data from the two regions.

Fig. B1 shows a comparison of the Dickey et al. (2000), Dempsey et al. (2022), and no HISA correction, on the four main quantities studied in this work. Without any HISA correction, we find |$b = 0.42^{+0.11}_{-0.19}$|⁠, but with the Dickey et al. (2000) correction factor we find |$b = 0.47^{+0.41}_{-0.16}$| and with the Dempsey et al. (2022) correction we find |$b = 0.49^{+0.22}_{-0.13}$|⁠. The Mach number is not affected by the correction factor because it is not a function of the column density, but interestingly, the Brunt factor also remains largely invariant, suggesting that the ratio of the column density dispersion to the volume density dispersion is not particularly sensitive to corrections of the |${\rm H}\,\rm{\small I}$| intensity.

Fig. B2 shows the respective maps of b without HISA correction (left-hand panel), and using the Dickey et al. (2000) (equation B2, middle panel) and Dempsey et al. (2022) (equation B3, right-hand panel) corrections, respectively.