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

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 parameterised 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$ characterising 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 HI emission observations of the Small Magellanic Cloud (SMC) from the GASKAP-HI survey. We find that the turbulence driving parameter varies between $b\sim 0.3$ and $b\sim 1.0$ within the main body of the SMC, but the median value converges to $b\sim0.51$, suggesting that the turbulence is overall driven more compressively ($b>0.38$). We observe no correlation between the $b$ parameter and HI or H$\alpha$ intensity, indicating that compressive driving of HI turbulence cannot be determined solely by observing HI or H$\alpha$ 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.


INTRODUCTION
The interstellar medium (ISM) is turbulent and mostly composed of multi-phase 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 (Mac Low & Klessen 2004;Elmegreen & Scalo 2004;Scalo & Elmegreen 2004;McKee & Ostriker 2007;Hennebelle & Falgarone 2012;Federrath & Klessen 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 timescales (Stone et al. 1998;Mac Low et al. 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 characterise 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. 2009Burkhart et al. , 2010;;Patra et al. 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 et al. 1997;Passot & Vázquez-Semadeni 1998;Federrath et al. 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;Padoan & Nordlund 2011;Hennebelle & Chabrier 2011;Federrath & Klessen 2012;Burkhart & Mocz 2019).Furthermore, studies by Price et al. (2011); Konstandin et al. (2012); Molina et al. (2012); Nolan et al. (2015); Federrath & Banerjee (2015); 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 σ ρ/ρ 0 is the three-dimensional (3D) standard deviation of density (ρ) scaled by the mean density (ρ0), 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(Federrath et al. , 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 eq. 9 in Federrath et al. 2010).The b parameter has been derived in theoretical models, and tested in simulations (Federrath et al. 2008(Federrath et al. , 2010;;Price et al. 2011;Molina et al. 2012;Nolan et al. 2015;Federrath & Banerjee 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 (σ ρ/ρ 0 and M), and the respective quantities that are accessible in observations, i.e., the column density and the intensity-weighted velocity centroid (Brunt et al. 2010;Kainulainen et al. 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 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 et al. (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), 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 Eq. (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 photoionising radiation.There has also been one extra-galactic 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 Hi 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 Hi 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, multi-phase atomic Hi 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 et al. 2020;Seta & Federrath 2022).Here we present a new method for mapping the turbulence driving parameter over large scales, and apply it to Hi 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 Sec. 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 Hi 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, 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 organised as follows.Sec. 2 introduced the data used in this work.Sec. 3 introduces our method for mapping the turbulence driving parameter, describing how we reconstruct the 3D density and velocity dispersions in large-scale observational datasets.Our results are presented in Sec. 4, and we discuss potential correlations of the b parameter with Hi and Hα emission in Sec. 5. Sec.6 discusses our SMC measurements of the b parameter in relation to other environments presented in the literature.We summarise our conclusions and pathways to future work in Sec. 8.

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 Hi 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 Hi data cube1 , resulting from 20.9 hours of total integration with ASKAP, used for our analysis is the most sensitive (rms brightness temperature σT = 1.1 K per 0.98 km s −1 spectral channel) and detailed view (30 ′′ restoring beam or ∼ 10 pc at the distance of the SMC) of the Hi 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 Sec. 3 of Pingel et al. (2022).

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 quantitiesb, as defined in Eq. ( 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 dv = 0.98 km s −1 , and integrate from v = 60.5 km s −1 to v = 235.0km 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 a signalto-noise ratio ≥ 10).We discuss the effect of the choice of signal-to-noise threshold on our results in Appendix A.

Influence of the multi-phase Hi
Eq. ( 1) was derived for isothermal, hydrodynamic gas.Hi emission does not originate from isothermal gas.The temperature structure of multi-phase neutral Hi can be roughly categorised into cold neutral medium (CNM; T ∼ 100 K), lukewarm neutral medium (LNM; 100 ≲ T /K ≲ 6000) and warm neutral medium (WNM; T ≳ 6000 K) (Wolfire et al.  2018), the overlays show the approximate location of the Bar (dashed box) and the Wing (dot-dashed box), and arrows indicate the directions towards the Magellanic Stream and the Magellanic Bridge linking to the LMC.The beam size is 30 ′′ , which is too small to show on this map.Each pixel in the map is 7 ′′ .The orange circle in the top right corner represents the full width at half maximum (FWHM) of the kernel used throughout the analysis in this work; it has a diameter of 10 instrument beams.Right-hand panel: Same as the left-hand panel, but for the moment-1 map, which shows the intensity-weighted velocity centroid, v. 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 et al. 2019).As such, the emission from Hi originates from a mixture of CNM, LNM and WNM.In order to apply Eq. ( 1) to the Hi 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 Hi 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 Hi 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 intra-clump 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 intra-clump velocity dispersion also traces the WNM velocity dispersion (Mohapatra et al. 2022;Kobayashi 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 LOS (i.e., the centroid velocity) for the velocity analysis (below in Sec.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% (Dempsey et al. 2022), meaning that the majority of the neutral Hi 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.

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 positionposition-velocity (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 Eq. ( 1) is constructed from the volume density dispersion σ ρ/ρ 0 and the turbulent sonic Mach number 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 (Sec.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 (Sec.3.4).To compute M, we then divide the 3D velocity dispersion by the sound speed of the gas (Sec.3.5).In order to isolate the true density and velocity fluctuations without contamination from nonturbulent 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 (Sec.3.2).This sequence of calculations is performed inside a Gaussian-weighted roving kernel (Sec.3.1), so as to build spatial maps of the 3D density variance, Mach number, and b parameter.

Roving kernel
The density contrast and velocity fluctuations are scaledependent 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 Sec.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.

Gradient subtraction
In calculating both the column density dispersion and the centroid velocity dispersion, we apply the gradientsubtraction 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 summarise 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. (2018Sharda et al. ( , 2019)); 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 normalised logarithmic column density log 10 (NHI/N0), where N0 is the mean column density inside the kernel, or the intensity-weighted velocity centroid v.

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 normalise the column density by the mean column density inside the kernel (left-hand panel of Fig. 2) and then fit a gradient to the normalised map, using Eq. ( 5) (middle panel of Fig. 2).We then subtract the gradient from the original normalised map via Eq.( 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, σ N 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 does not follow the anticipated Gaussian distribution, but after applying the gradient subtraction we see in the righthand panel that the resultant distribution well approximates  5); right: gradient-subtracted column density map via Eq.( 6) ).We can see a distinct gradient in the original map, while the gradient-subtracted map shows a clearer density contrast.The black circles in each of the top panels shows the size of the instrument beam, while the dashed circles represent the kernel's FWHM (10 beams).
The lower panels show the PDF in each case (shaded histogram), fitted with a Gaussian (solid line).The original data (left) has strong non-Gaussian components, as a result of the large-scale gradient.By contrast, the gradient-corrected data (right) is well-approximated by a Gaussian distribution in logarithmic column density, a universal feature of compressible turbulence.
a Gaussian, which is a hallmark of compressible isothermal turbulence2 .

Conversion from 2D to 3D density dispersion
To convert the 2D density dispersion σ N HI /N 0 to the 3D density dispersion, σ ρ/ρ 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 P2D(k) and the 3D density power spectrum P3D(k) is given by (Federrath & Klessen 2013), where k is the wave number.We exploit this relation to recover P3D(k).
We first compute P2D(k) of the gradient-subtracted column density, (NHI/N0 − 1 shown in Fig. 3), which immediately gives us P3D(k) of the quantity ρ/ρ0 as per Eq. ( 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, σ ρ/ρ 0 .Eq. ( 8) relies on the assumption that the input field is isotropic in k-space, and that therefore P2D(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 P2D(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 (c.f., Sec. 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 P2D(k), such that the distribution of the density variations in real 2D space may be similarly distributed in 3D space.

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 gradientsubtracted 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 factor3 of C any (c−grad) = 3.3 ± 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 Sec.7.2.

Mach number
The sonic Mach number (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 cs.The last step in the pipeline is to divide σv,3D by the sound speed to produce the quantity 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, kB is the Boltzmann constant, T is the gas temperature, µ is the mean particle weight, and mH is the mass of a hydrogen atom.Because we do not have a way to estimate the temperature of the combined Hi 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 ± 1.0) × 10 3 K, which gives a sound speed of cs = (8.3± 0.6) km s −1 .Estimates of the WNM temperature in the MW range from about 4000 K to 104 K, depending on the method used (Wolfire et al. 2003

RESULTS
In this section, we present and discuss the output of the analysis pipeline described in Sec. 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.

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 σ ρ/ρ 0 , following the methods described in Sec.3.3.The quantity σ ρ/ρ 0 measures the turbulent density variations in Eq. ( 1). .Same as Fig. 2, but for the intensity-weighted velocity v along the LOS (i.e., the moment-1).Similar to the column density, we find a significant gradient in the original data, leading to non-Gaussian components in the PDF of v.After gradient-correction, the PDF of v follows a Gaussian distribution, which is a hallmark of turbulent flows.Thus, the gradient-subtraction method (c.f.Sec.3.2) successfully filters out non-turbulent contributions and therefore isolates the turbulent velocity components in the data.
Table 1.Summary of key quantities.The measured quantities are averaged within the closed contour, and therefore do not represent an actual global value for the entire SMC.The error bars for the quantities found in this work represent the range between the 16 th and 84 th percentile, or the variation around the median within the main contour, not actual uncertainties in those values.The quantities noted as 'converged' are the best fit parameters found in Fig. 9, and their error bars are those produced by the fitting process.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 Sec. 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 Eq. ( 1).

Symbol
The right-hand panel of Fig. 5 shows the map of the turbulent sonic Mach number, M, following the methods in Sec.3.5.The variations within the analysis contour are ∼ 3, similar to the variations of σ ρ/ρ 0 .The Mach number distribution is also relatively uniform within the analysis contour, with notable regions of high M at the top of the Wing region, where it intersects with the Bar region, similar to σ ρ/ρ 0 .Comparing the 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 M ∼ 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 M ∼ 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 Hi column density, and found that 90% of the Hi 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  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% lower than inside the contour), which is why we have chosen to analyse only the data in the first closed contour.

Spatial distribution of the turbulence driving parameter
Combining the maps in Fig. 5 via Eq.( 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 Sec. 5, by studying the turbulence driving parameter as a function of the Hi and Hα emission in the SMC.

Influence of the kernel size
The 3D velocity dispersion and the volume density dispersion are scale-dependent quantities (Kim & Ryu 2005;Kritsuk et al. 2007;Kowal & Lazarian 2007;Federrath et al. 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, & 40.0.The panels of Fig. 7 show the driving parameter for these 5 values of the kernel size, and highlight the different levels of spatial granularity resulting from each kernel size.All 5 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 4 main quantities as a function of the kernel size (from left to right): the Brunt factor R 1/2 , the density dispersion σ ρ/ρ 0 , the Mach number M, and the driving parameter b.Fig. 8 demonstrates that the width and peak of the M and σ ρ/ρ 0 PDFs increase with increasing kernel size, and that the Brunt factor (R 1/2 ) decreases with increasing kernel size.Because the width of the M and σ ρ/ρ 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 panel of Fig. 9) is R 1/2 = 0.31 ± 0.01 + [N bpk /(0.16 ± 0.10)] −0.88±0.21, (12) such that in the limit of an infinitely large kernel (N bpk → ∞) the Brunt factor converges to R 1/2 = 0.31 ± 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 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 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 σ ρ/ρ 0 = (0.18 ± 0.01) (N bpk /10) 0.76±0.01 .( 14) and for the Mach number we find M = (0.37 ± 0.04) (N bpk /10) 0.70±0.05 . (15) The exponents of these two power laws agree within the uncertainties, which shows that σ ρ/ρ 0 and 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 M and σ ρ/ρ 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%.Fig. 10 shows the PDFs of R 1/2 , σ ρ/ρ 0 , M, and b inside the main analysis contour for the standard kernel size of N bpk = 10.We find R 1/2 = 0.34 ± 0.01, σ ρ/ρ 0 = 0.18 +0.06 −0.04 and 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 16 th and 84 th percentile range).Nevertheless, even the 1-sigma lower limit of b, i.e., the 16 th percentile value is with b16 = 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 Hi gas is subject to predominantly compressive turbulence driving mechanisms in the SMC.

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 Hi 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 Hi intensity, instead the turbulence driving seems to be in the mixed-to-compressive (b > 0.4) regime regardless of the Hi 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 Hi.
Similarly, we find that b is slightly compressive across the For the smallest kernel size (N bpk = 2.5), the peak of both σ ρ/ρ 0 and M are shifted towards smaller values, while the largest kernel size (N bpk = 40.0)shifts the peak of the distributions to higher values.On the other hand, the Brunt factor, R 1/2 , exhibits the opposite trend.While both σ ρ/ρ 0 and M depend on the kernel size, the ratio of the two, i.e., the turbulence driving parameter (Eq. 1) is independent of the kernel size for sufficiently large kernel size (N bpk ≳ 10).At the distance of the SMC, the physical sizes of the kernels are approximately 25 pc, 50 pc, 100 pc, 200 pc and 400 pc, respectively.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 Hii 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 dis-entangle 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 Hi gas in the SMC is not dominated by the star-formation activity.Same as left-hand panel, but for Hα data from MCELS (Smith et al. 1999;Winkler et al. 2015).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 Sec.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 M and σ ρ/ρ 0 change with the kernel or region size observed.We show the median values for the five kernel sizes investigated in Sec.4.2.1, cor-responding to a physical size of 25 pc, 50 pc, 100 pc, 200 pc 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; c.f., Fig. 9).

COMPARISON TO OTHER MEASUREMENTS OF THE TURBULENCE DRIVING PARAMETER IN DIFFERENT ENVIRONMENTS
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 spe- cific 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 σ ρ/ρ 0 = 0.6 ± 0.2 and M = 0.87 ± 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 discus-sion; 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 ∼ 100 pc.In Marchal & Miville-Deschênes ( 2021), the spatial scale on which they measure σ ρ/ρ 0 and M is ∼ 63 pc.Referring to Fig. 9, we can see that had we used a kernel N bpk = 5 (∼ 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 ± 1000 K, which is about 1000 K higher than the temperature used in Marchal & Miville-Deschênes (2021), but lower than the 10 4 K observed by Murray et al. (2018).There is a wide range of temperature estimates for the WNM in Galactic Hi, 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 M values to be about 10% 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 line of sight, it is possible that we have underestimated the volume density dispersion via the Brunt et al. (2010) method (Sec. 3.3.1).We discuss this issue further in Sec.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 line of sight, and therefore σ ρ/ρ 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.

Magnetic fields
The ISM is ubiquitously magnetised, 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 b where plasma beta, β = P th /Pmag = 2c 2 s /v 2 A , is the ratio of the thermal to magnetic pressure, and the Alfvén speed vA = B turb /(4πρ0) 1/2 for the turbulent component of the magnetic field, and ρ0 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 line of sight magnetic field strength (−0.3±0.1 µG) and the random component of the field (5 +3 −2 µG).For the WNM the number density is ∼ 0.2 − 0.5 cm −3 (Ferrière 2001).This gives us an Alfvén speed between vA ∼ 5 − 30 km s −1 .Combining this we estimate β ∼ 0.1 − 2. Using this in Eq. ( 16), the factor (1+β −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.

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 multi-component 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 any (c−grad) = 3.3 ± 0.5) to recover the LOS velocity fluctuations from the moment-1 map only, to within their quoted 10% accuracy, which then gives us the 3D turbulent velocity dispersion.

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 et al. 1954;Hindman et al. 1963;Stanimirović et al. 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 inter-actions between the SMC and LMC are required to consistently reproduce the Stream and the Bridge.(Besla et al. 2012;Diaz & Bekki 2012;Lucchini et al. 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 (Sec. 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 line of sight is comparable to the density dispersion in the plane of the sky.We further find that the dispersion in the planeof-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 line of sight, 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 Eq. 1) may be relatively robust against LOS averaging effects (similar to how it is independent of the kernel size on sufficiently large scales; c.f., 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.

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 (∼ 10%) and the error associated with our estimate for the WNM temperature and therefore the sound speed (∼ 10%).The error introduced by our conversion of the velocity variance via the Stewart & Federrath (2022) method is ∼ 10%, 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 ∼ 20% 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.

CONCLUSIONS
In this work we presented a new generalised 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 (Sec.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 (Sec.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 Eq.(1).In summary, the main steps of the pipeline performed in each instance of the roving kernel are: (i) To isolate only turbulent density variations, fit and subtract a linear gradient from the normalised column density map (see Fig. 2).
(ii) Via the "Brunt Method" (Sec.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, R 1/2 .
(iii) Using R 1/2 , reconstruct the volume density dispersion from the column density dispersion (see Fig. 5).
(iv) Fit and subtract a linear gradient from the centroid velocity map, thus isolating only turbulent velocity fluctuations (see Fig. 4).
(v) Find the dispersion of the gradient-corrected centroid velocity, and convert to 3D velocity dispersion via the "Stewart" method (see Sec. 3.4. (vi) Divide the 3D velocity dispersion by the sound speed to recover the turbulent sonic Mach number, M (see Fig. 5).
(vii) Finally, divide the volume density dispersion map by the Mach number map to create a map of the turbulence driving parameter, b, via Eq. 1 (see Fig. 6).
We have applied this method to high-resolution GASKAP-HI observations of 21 cm emission from the SMC (Sec.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 ≳ 100 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 Hi 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.opment of the ideas behind this work.I.A.G. would like to thank the Australian Government and the financial support provided by the Australian Postgraduate Award.C.F. acknowledges funding by the Australian Research Council (Future Fellowship FT180100495 and Discovery Projects DP230102280), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD).C.F. 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.S.S. and N.M.P. 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.  .Same as Fig. 10, but comparing the influence of correcting the Hi column density for optical-depth effects.We use correction relationships from Dickey et al. (2000) and Dempsey et al. (2022), and find that σ ρ/ρ 0 and b increase slightly for both correction cases by 12% and 17%, respectively, compared to no HISA correction.The Brunt factor is largely independent of the correction and the Mach number is completely independent of the correction (by definition, as it does not depend on the column density).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 Hi intensity.

Figure 1 .
Figure 1.Left-hand panel: Moment-0 map, which shows the integrated intensity.Because Hi is (mostly) optically thin (See Appendix B for further discussion on this topic), this quantity represents the column density N HI , which has been normalised by the mean of the emission inside the black contour, N 0,SMC = 4.5 × 10 21 cm −2 .The black contour (at 2.0 × 10 21 cm −2 ) denotes the first closed contour that encompasses the main body of the SMC.Following McClure-Griffiths et al. (2018), the overlays show the approximate location of the Bar (dashed box) and the Wing (dot-dashed box), and arrows indicate the directions towards the Magellanic Stream and the Magellanic Bridge linking to the LMC.The beam size is 30 ′′ , which is too small to show on this map.Each pixel in the map is 7 ′′ .The orange circle in the top right corner represents the full width at half maximum (FWHM) of the kernel used throughout the analysis in this work; it has a diameter of 10 instrument beams.Right-hand panel: Same as the left-hand panel, but for the moment-1 map, which shows the intensity-weighted velocity centroid, v.

Figure 2 .
Figure 2.An example of a single kernel, illustrating the process of fitting a gradient to the column density map.In this example the column density (N HI ) is normalised by the mean column density (N 0 ) in the kernel.The upper panels show the kernel-weighted maps (left: original column density map; middle: fitted gradient via Eq.(5); right: gradient-subtracted column density map via Eq.(6) ).We can see a distinct gradient in the original map, while the gradient-subtracted map shows a clearer density contrast.The black circles in each of the top panels shows the size of the instrument beam, while the dashed circles represent the kernel's FWHM (10 beams).The lower panels show the PDF in each case (shaded histogram), fitted with a Gaussian (solid line).The original data (left) has strong non-Gaussian components, as a result of the large-scale gradient.By contrast, the gradient-corrected data (right) is well-approximated by a Gaussian distribution in logarithmic column density, a universal feature of compressible turbulence.

Figure 3 .
Figure 3.An example of the Fourier image of the gradientsubtracted column density shown in Fig. 2. The wavenumber k = 1 corresponds to 3 × D kernel in real space, and the subscripts 'RA' and 'DEC' denote the orientation in real space.Given the high level of point symmetry in this image, especially for small wavenumbers, which correspond to large length scales inside the kernel, the conversion from 2D to 3D density dispersion via the Brunt et al. (2010) method is expected to introduce only a ∼ 10% uncertainty.
Figure 4. Same as Fig.2, but for the intensity-weighted velocity v along the LOS (i.e., the moment-1).Similar to the column density, we find a significant gradient in the original data, leading to non-Gaussian components in the PDF of v.After gradient-correction, the PDF of v follows a Gaussian distribution, which is a hallmark of turbulent flows.Thus, the gradient-subtraction method (c.f.Sec.3.2) successfully filters out non-turbulent contributions and therefore isolates the turbulent velocity components in the data.

Figure 5 .
Figure 5. Left-hand panel: 3D volume density dispersion σ ρ/ρ 0 .The black contour is the first closed contour at 2.0 × 10 21 cm −2 .The orange circle shows the FWHM of the kernel with N bpk = 10 telescope beams per kernel.The white regions around the edges are a result of our SNR cut (see Sec. 2).Right-hand panel: Same as left-hand panel, but for the turbulent Mach number.
of this work and Fig.8inBurkhart 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 inBurkhart et al. (2010)) increases the Mach number range we recover in this work (discussed further in Sec.4.2.1 below), and pushes the values up into the trans-sonic/supersonic regime.

Figure 6 .
Figure 6.Map of the turbulence driving parameter, b, calculated via Eq.(1), by combining the information in the maps shown in Fig. 5.The main contour (black line) and kernel (shown in the bottom left corner) are the same as in Fig. 5.The pink contours denote purely solenoidal driving (dark pink, b ∼ 0.3), naturally mixed driving (medium pink, b ∼ 0.38) and purely compressive driving (light pink, b ∼ 1.0).These lines are also shown on the colourbar.We see strong spatial variations in the turbulence driving parameter, ranging from purely solenoidal to purely compressive driving across the SMC.

Figure 7 .
Figure 7.Comparison the driving parameter map using different kernel sizes.From top to bottom, each panel uses a kernel with N bpk = 2.5, 5.0, 10.0, 20.0, 40.0.The orange circle in the bottom left-hand corner of each panel shows the kernel FWHM, and the black contour is the first closed contour as throughout this work.

NFigure 8 .
Figure8.PDFs of the four main analysis quantities (R 1/2 (blue panel), σ ρ/ρ 0 (pink panel), M (green panel) and b (purple panel)) for five different kernel sizes.For the smallest kernel size (N bpk = 2.5), the peak of both σ ρ/ρ 0 and M are shifted towards smaller values, while the largest kernel size (N bpk = 40.0)shifts the peak of the distributions to higher values.On the other hand, the Brunt factor, R 1/2 , exhibits the opposite trend.While both σ ρ/ρ 0 and M depend on the kernel size, the ratio of the two, i.e., the turbulence driving parameter (Eq. 1) is independent of the kernel size for sufficiently large kernel size (N bpk ≳ 10).At the distance of the SMC, the physical sizes of the kernels are approximately 25 pc, 50 pc, 100 pc, 200 pc and 400 pc, respectively.

Figure 9 .
Figure 9.The median values of the four main analysis quantities as a function of N bpk , with fitted functions.Top left: The Brunt factor R 1/2 , Top right: The density variance σ ρ/ρ 0 .Bottom left: The Mach number M. Bottom right: The driving parameter b.The error bars on each data point represent the 16 th and 84 th percentiles, and were used in the fits, with the errors on the fit parameters displayed in the legend of each panel.

Figure 10 .Figure 11 .
Figure10.Same as Fig.8but for the N bpk = 10 case.The solid lines in each panel represent the 50 th percentile, while the dashed lines show the 16 th and 84 th percentiles, respectively.These values are also reported in the legend.

Fig. 12
Fig.12concludes 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 Sec.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 M and σ ρ/ρ 0 change with the kernel or region size observed.We show the median values for the five kernel sizes investigated in Sec.4.2.1, cor-

Figure 12 .
Figure 12.A summary of the available observational estimates for the turbulence driving parameter b of several sources, contextualising our results for the SMC Hi gas.The y-axis shows the 3D (volume) density dispersion (σ ρ/ρ 0 ), and the x-axis shows the turbulent Mach number (M), including a factor involving plasma β (ratio of thermal to magnetic pressure), as some of the literature values shown (Taurus and G0.253+0.016)have been calculated using the magnetised version of the density dispersion -Mach number relation (see Sec. 7.1 for further discussion on this point).The three diagonal lines show the theoretical limits for compressive (b = 1.0, dotted), mixed (b = 0.38, solid) and solenoidal (b = 0.33, dashed) driving of the turbulence (Federrath et al. 2010).The SMC results of this work are shown in the lower left-hand corner, with the colourbar denoting the probability density of points in our SMC maps of σ ρ/ρ 0 and M for N bpk = 10.The hexagons show the median values for the SMC quantities in the five different kernel sizes we investigated in Sec.4.2.1, which correspond to roughly 25 pc, 50 pc, 100 pc, 200 pc and 400 pc.The error bars on these points show the 16 th to 84 th percentile on each axis.For context we include a variety of sources from the literature: Taurus (dark blue star) (Brunt 2010), which includes magnetic field estimates and revised Mach number estimations from Kainulainen & Tan (2013), using 13 CO line imaging observations; IC5146 (blue cross) (Padoan et al. 1997), using 12 CO and 13 CO observations; GRSMC 43.30-0.33(aqua plus) (Ginsburg et al. 2013), observed in H 2 CO absorption and 13 CO emission; 'The Brick' (G0.253+0.016,teal square) (Federrath et al. 2016), using HNCO observations; 'The Pillars of Creation' (NGC 3372 pillars, magenta diamonds) (Menon et al. 2021), from 12 CO, 13 CO and C 18 O; 'The Papillon Nebula' (LMC N159E, pink circle) (Sharda et al. 2022), again in 12 CO, 13 CO and C 18 O; and the WNM in the MW (lilac triangle) (Marchal & Miville-Deschênes 2021) (Hi observations).

Figure A1 .Figure A2 .
Figure A1.PDFs of the same quantities shown in Fig.10, but for three different signal-to-noise ratios.The data used to construct these PDFs are from inside the contour only.We can see that using a contour of 10σ provides a robust dataset, which is largely insensitive to the per-channel noise threshold applied to the raw PPV cube.The R 1/2 , σ ρ/ρ 0 , M, and b values are nearly invariant across the three noise levels we compare here.This gives us confidence that all analysis inside this contour is unaffected by any per-pixel noise.
Figure B1.Same as Fig.10, but comparing the influence of correcting the Hi column density for optical-depth effects.We use correction relationships fromDickey et al. (2000) andDempsey et al. (2022), and find that σ ρ/ρ 0 and b increase slightly for both correction cases by 12% and 17%, respectively, compared to no HISA correction.The Brunt factor is largely independent of the correction and the Mach number is completely independent of the correction (by definition, as it does not depend on the column density).

Fig
Fig. B1 shows a comparison of the Dickey et al. (2000), Dempsey et al. (2022), and no HISA correction, on the 4 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