# Multi-phase HI clouds in the Small Magellanic Cloud halo

F. Buckland-Willis<fnref n="1" />,<fnref n="2" />,*, M. A. Miville-Deschênes<fnref n="1" />,<fnref n="2" />, A. Marchal<fnref n="3" />, J. R. Dawson<fnref n="4" />,<fnref n="5" />, H. Dénes<fnref n="6" />, E. M. Di Teodoro<fnref n="7" />, J. M. Dickey<fnref n="8" />, S. J. Gibson<fnref n="9" />, I. P. Kemp<fnref n="10" />,<fnref n="11" />, C. Lynn<fnref n="3" />, Y. K. Ma<fnref n="3" />, N. M. McClure-Griffiths<fnref n="3" />, C. E. Murray<fnref n="12" />,<fnref n="13" />, N. M. Pingel<fnref n="14" />, S. Stanimirović<fnref n="14" />, and J. Th. Van Loon<fnref n="15" />

Received 17 October 2024 / Accepted 28 November 2024

## ABSTRACT

**Context.** The Galactic ASKAP collaboration (GASKAP) is undertaking an HI emission survey of the 21cm line to map the Magellanic system and the Galactic plane with the Australian Square Kilometre Array Pathfinder (ASKAP). One of the first areas observed in the Pilot Phase I of the survey was the Small Magellanic Cloud (SMC). Previous surveys of the SMC have uncovered new structures in the periphery of the SMC, along relatively low column density lines of sight.

**Aims.** In this work we aimed to uncover the phase distribution of three distinct structures in the periphery of the SMC. This work will add to the constraints we have on the existence and survival of the cold neutral medium (CNM) in the SMC.

**Methods.** We used ROHSA, a Gaussian decomposition algorithm, to model the emission across each cloud and classify the HI emission into their respective phases based on the linewidths of the fitted Gaussians. We created maps of velocity and column density of each phase of the HI across these three clouds. We measured the HI mass and CNM number density for each cloud. We also compared the HI results across the different phases with other gas tracers.

**Results.** We find that in two clouds, the ends of each cloud are almost completely CNM dominated. Analysis of these two clouds indicates they are experiencing a compressive force from the direction of the SMC main body. In the third cloud we find a uniform CNM distribution along one wall of what is likely a supershell structure. Comparison with previous measurements of CO clumps in two of the clouds show the CO and HI are co-moving within a few km s⁻<fnref n="1" /> in regions of high HI column density, particularly when considering just the CNM.

**Key words.** ISM: clouds – ISM: kinematics and dynamics – galaxies: dwarf – galaxies: ISM

## 1. Introduction

The interstellar medium (ISM) consists of many different components, of which one of the most abundant is neutral hydrogen (HI). The study of neutral hydrogen helps us to understand the distribution and conditions needed to facilitate star formation. Understanding the ways in which HI cools and condenses to become molecular hydrogen is a key part of the star formation process.

One factor that affects the cooling ability of the ISM is the metallicity. Understanding how metallicity affects the distribution of the neutral hydrogen, requires studying the HI in

different environments. The Magellanic Clouds provide us with an excellent, nearby laboratory in which to study HI in low-metallicity environments. The Large Magellanic Cloud (LMC) has a characteristic metallicity of 0.5 Z⊙ (Rolleston et al. 2002) and the Small Magellanic Cloud (SMC) has one of 0.2 Z⊙ (Russell & Dopita 1992). Both clouds provide low-metallicity environments, with significant HI gas reservoirs that can be well-resolved in observations. With the advent of high-resolution, high-sensitivity interferometers such as the upcoming Square Kilometre Array (SKA) (McClure-Griffiths et al. 2015), galaxies such as the Magellanic Clouds will be imaged at the same physical scales as the Milky Way had been previously, allowing for a close comparison of a broad range of Galactic environments.

<footnote n="1">¹ Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, 24 rue Lhomond, 75005 Paris Cedex 05, France</footnote>

<footnote n="2">² AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Research School of Astronomy and Astrophysics, The Australian National University, Canberra, ACT 2611, Australia</footnote>

<footnote n="4">⁴ School of Mathematical and Physical Sciences and Astrophysics and Space Technologies Research Centre, Macquarie University, 2109, NSW, Australia</footnote>

<footnote n="5">⁵ Australia Telescope National Facility, CSIRO Space and Astronomy, PO Box 76, Epping NSW 1710, Australia</footnote>

<footnote n="6">⁶ School of Physical Sciences and Nanotechnology, Yachay Tech University, Hacienda San José S/N, 100119, Urcuquí, Ecuador</footnote>

<footnote n="7">⁷ Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, 50019 Sesto Fiorentino, Italy</footnote>

<footnote n="8">⁸ School of Natural Sciences, Private Bag 37, University of Tasmania, Hobart, TAS, 7001, Australia</footnote>

<footnote n="9">⁹ Department of Physics & Astronomy, Western Kentucky University, 1906 College Heights Blvd., Bowling Green, KY 42101, USA</footnote>

<footnote n="10">¹⁰ International Centre for Radio Astronomy Research (ICRAR), Curtin University, Bentley 6102, WA, Australia</footnote>

<footnote n="11">¹¹ CSIRO Space and Astronomy, 26 Dick Perry Avenue, Kensington, 6151, WA, Australia</footnote>

<footnote n="12">¹² Department of Physics & Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA</footnote>

<footnote n="13">¹³ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="14">¹⁴ Department of Astronomy, University of Wisconsin, Madison, WI 53706-15821, USA</footnote>

<footnote n="15">¹⁵ Lennard-Jones Laboratories, Keele University, ST5 5BG, UK</footnote>

<footnote>* Corresponding author; francesbw.astronomy@outlook.com</footnote>
The SMC is an irregular galaxy, whose 3D structure is quite complex. Previous HI studies of the SMC have shown that the neutral gas is very dynamically complex, especially along lines of sight in the central areas of the galaxy. The neutral gas distribution has been shaped by the interaction of the SMC with the LMC, simulated in Diaz & Bekki (2011) and Besla et al. (2012), as well as star formation activity in the galaxy. Notably, in the HI survey conducted with a combination of Australian Telescope Compact Array (ATCA) and Murriyang (Parkes) observations (Staveley-Smith et al. 1997) just over 500 shells were catalogued, indicating that turbulence, star formation and their resulting winds, and supernovae are changing the HI landscape. These processes inject energy into the medium which can compress the neutral hydrogen, forming the dense gas needed to form the next generation of stars (Inutsuka et al. 2015; Dawson 2013). Disruption from phenomena occurring within the galaxy and the tidal forces from the larger Magellanic system interaction has made the SMC a very disturbed and dynamic system. Stellar studies have shown that the galaxy's stellar population is extended by 1–10 kpc along the line of sight (Muraveva et al. 2018). A recent study (Murray et al. 2024) comparing the HI and stellar population of the galaxy suggests that there are two main sections of the galaxy positioned at differing distances along the line of sight. This elongation, coupled with the fact that HI spectral data give only velocity information along the line of the sight, makes it particularly difficult to reliably determine distances to any discrete structures in the galaxy.

Despite the challenges in studying the HI distribution of the SMC, its low metallicity provides a valuable environment with different conditions to that of the Milky Way. Because the main cooling mechanisms of the HI are less efficient at lower metallicities (Wolfire et al. 1995), the distributions of the warm neutral medium (WNM) compared to the cold neutral medium (CNM) through the SMC should differ when compared to the HI under Galactic conditions. These different phases, WNM and CNM, are defined by their typical temperatures in the Milky Way in Wolfire et al. (1995) and Wolfire et al. (2003), where WNM (T ~ 8000 K) and CNM (T ~ 100 K). Another phase of the HI, the unstable neutral medium (UNM) is often seen in studies (McClure-Griffiths et al. 2023). This phase is HI which lies in the thermally unstable range of pressures and temperatures, and it is a transitory state of HI that will evolve into either stable phase of the HI (CNM or WNM). The exact temperatures and pressures that define all of these phases depend on the ISM conditions, as detailed in Wolfire et al. (1995) and Wolfire et al. (2003).

The CNM distribution in the SMC has been studied previously through absorption surveys, both targeted and un-targeted. Across three previous absorption studies (Dickey et al. 2000; Jameson et al. 2019; Dempsey et al. 2022) the reported mean CNM fraction for the SMC has ranged from 0.07 to 0.2. Dempsey et al. (2022) reported CNM fraction values for individual lines of sight, which show that there are a few lines of sight with a CNM fraction close to 0.5, but the majority are below 0.2. CNM is found to populate discrete filamentary structures in the ISM (McClure-Griffiths et al. 2023) which cannot be well-characterised by absorption studies alone. This is because, while absorption studies are very effective at identifying CNM, especially when the gas is very cold and easily absorbs, the detections of the CNM are limited by the background source density. Therefore absorption studies cannot provide a full map of the CNM across any specific field.

New developments in Gaussian decomposition and other HI analysis methods have helped to map the CNM fraction across

large fields. In particular, Regularized Optimization for Hyper-Spectral Analysis (ROHSA), first introduced in Marchal et al. (2019), has since been used to map the different HI phases across a number of Galactic structures (Marchal et al. 2021; Taank et al. 2022; Vujeva et al. 2023). The ROHSA algorithm takes advantage of the expected spatial coherence of the HI signal and enforces similar solutions between neighbouring sightlines. This allows the distribution of the different phases of the HI across the field to be tracked directly from the emission spectra. Of course, it is vital to have the appropriate spectral resolution in order to resolve the CNM signals that typically have linewidths <3 km s⁻¹. Another method, using Fourier transforms of the spectral axis of HI data, demonstrated in Marchal et al. (2024), estimates a lower limit for the CNM fraction with great speed across large fields of view.

Analysis of data with these tools will help refine techniques in preparation for HI surveys currently observing with the Australian Square Kilometre Array Pathfinder (ASKAP) telescope and the data that will eventually come from the SKA-MID (Square Kilometre Array MID) telescope. These telescopes provide datasets with improved angular resolution and sensitivity, as well as adequate velocity resolution to uncover the distribution of the CNM within and beyond the Milky Way. In this work we look at three discrete HI clouds within the periphery of the SMC, first identified in McClure-Griffiths et al. (2018), with data from the Pilot phase of the Galactic ASKAP HI (GASKAP-HI) survey and we analysed them using the ROHSA algorithm.

In Section 2 we outline the data we analysed for this work and the noise estimation. In Section 3 we detail the modelling of the HI data using ROHSA, describing the parameters used in the fitting process and the associated uncertainties. In Section 4 we describe the results from the fitting process and the main trends we see across the three clouds. In Section 5 we discuss the HI in relation to other tracers CO and Hα, as well as each of its phases, and we discuss formation scenarios for each of these clouds.

## 2. Data

### 2.1. GASKAP collaboration data

The data in this work comes from the GASKAP-HI survey Pilot Phase I observations of the SMC. A total of 20.9 h of integration was undertaken over two observing blocks in December 2019. The data was reduced with a joint-imaging pipeline during which it was combined with single dish observations from the Murriyang telescope, the details are contained in the data release paper (Pingel et al. 2022). After processing, the hyperspectral data (in the form of a position-position-velocity, PPV, cube) for the SMC field was imaged at an angular resolution of 30" and spectral resolution of 0.98 km s⁻¹ (smoothed from the native resolution of ≈0.3 km s⁻¹). The field of view of the datacube is ≈5 × 5 deg², which nicely encompasses the majority of the HI distribution of the SMC, including the lower column-density areas with gas that leads to the Bridge and the Counter-Arm. The velocities of the SMC datacube range from 40.0 to 253.9 km s⁻¹, capturing the full dynamic range of the SMC.

### 2.2. Noise estimation

For this dataset, the median noise level measured across the whole cube was reported as 1.1 K in Pingel et al. (2022). As we consider specific regions of the SMC as opposed to the SMC as a whole, it is necessary to have a local estimation of the noise. Thus, in this work, we take an empirical approach to obtain
<figure><img src="image_1.png" /><figcaption>Fig. 1. Column density images of all three clouds (Left: Alpha, Middle: Hook, Right: Gamma) produced by integrating each field over the velocity range indicated in each panel. The blue arrows indicate the direction of the dynamical centre of the SMC derived in Di Teodoro et al. (2019b).</figcaption></figure>

the noise values across the datacube. This method requires identifying channels in the datacube that are ‘emission-free’, from which to measure the rms noise at each pixel. The velocity range of the datacube encompasses the entirety of the expected SMC emission and extra channels either side. Through visual inspection, channels where $v < 62$ km s⁻¹ and $v > 235$ km s⁻¹ were determined to be free of SMC emission, allowing 43 of 220 channels (19.5%) to be used to calculate the rms noise level of each pixel. Taking the rms of the values across the 43 channels resulted in a median noise value of 1.6 K, with the values increasing towards the edges of datacube as the beam response decreases. The noise in the HI line increases with the strength of the emission, so we calculate the noise along each line of sight as is done in Boothroyd et al. (2011):

$$ \sigma(v) = \sigma_0 \left( \frac{T_{\text{sys}} + T_{\text{B}}(v)}{T_{\text{sys}}} \right), \qquad (1) $$

where $\sigma_0$ is the rms noise level of an individual pixel. We use the provided value of $T_{\text{sys}} = 55$ K from Pingel et al. (2022). With this added dimension, we obtain a 3D map of the noise throughout the entire datacube.

# 3. Methods

## 3.1. Field selection

The fields analysed in this work were chosen to encompass three HI clouds first identified in McClure-Griffiths et al. (2018) as outflows from massive star formation within the SMC. These clouds were chosen because their emission is very strong in the densest regions of the clouds, with column densities of the order of $10^{20}$ cm⁻². Additionally, while the distribution of CNM in the main body of the SMC can be studied in absorption, these clouds have compact structures and were not probed by any background sources in the latest absorption study (Dempsey et al. 2022), due to the intrinsic limit imposed by the background source density. Thus, this work provides an opportunity to analyse the phase

distribution just from the emission of HI clouds with high signal to noise ratios.

In Figure 1 we show the total column density for each cloud, which we name Alpha, Hook, and Gamma, respectively, over the velocity range they span. Additionally, we also indicate the direction of the SMC dynamical centre as reported in Di Teodoro et al. (2019b). Each cloud has a different morphology, with the Alpha cloud being the smallest of the three, which assuming a distance of 63 kpc (Di Teodoro et al. 2019b) measuring 240 pc across the longest diagonal (north-east to south-west), with an irregular morphology. The Hook cloud is primarily composed of the strong ridge along the west side forming a long filament that spans 610 pc. It also continues looping around the north and down the east of the field, albeit with lower column densities. The Gamma cloud has a strong core filament, and appears to be broken into multiple clumps. It covers a distance of 560 pc, similar to the Hook cloud, but with a clumpier appearance.

## 3.2. Gaussian decomposition with ROHSA

### 3.2.1. Defining a best-fit solution

We used the ROHSA algorithm developed in Marchal et al. (2019) to spectrally decompose each of the three clouds. The ROHSA tool allows us to take advantage of the expected spatial coherency of the HI signal, that the decomposition solution along a line of sight is similar to that of its neighbours. This approach helps tackle the degeneracy problem that can muddy HI decomposition efforts along complicated lines of sight with multiple components with differing spreads and centres. The degree to which this coherency is enforced is controlled by 4 hyperparameters, detailed in Section 3.2.3, and it is the fine-tuning of these hyperparameters that will have the most impact on the final solution. As part of this work, an exploration of the 4-D hyperparameter space was undertaken to determine the best solution. We define a best solution as one whose residuals are normally distributed and achieves a reduced chi-squared value close to 1 with the fewest number of components.
<figure><img src="image_2.png" /><figcaption>Fig. 2. Example spectra from the Alpha cloud field. Bottom: mean spectrum of the field before (orange) and after (blue) subtraction of the main body emission signal. Top: spectra of two separate sightlines with the same colours indicating the spectra before and after subtraction. The grey windows show the spectral range retained for the fitting described in Section 3.2.3.</figcaption></figure>

### 3.2.2. Separating cloud and main body emission

While the signal from the three target clouds is strong in the brightest areas, the emission from the main body of the SMC dominates the mean spectrum over each full field containing the clouds, shown in Figure 2. In preliminary attempts to fit the full spectra, 8–9 components were used, when only 2 components are likely necessary to model the emission from the clouds themselves. If there are not enough components to fully model the main body and cloud emission across the whole field then the weakest emission channels will be un-fit. However, adding additional components to the model can increase the degeneracy of the solution and increase the computational cost. Since the emission from the main body of the SMC is of little interest to us in this work, we opted to isolate the cloud emission to simplify the fitting efforts.

To overcome this problem, we initially took a subset of the spectrum over the velocity range of interest and proceeded to fit just this section. Although the targeted clouds are removed from the main body emission of the SMC, there are some wide (WNM) components from the main body whose tails bleed into the velocity range of interest for these clouds. If we do nothing but take a subset of the data over the velocities of interest, we risk having these tails disrupt ROHSA’s fitting efforts. If we allow ROHSA to fit this unrefined data, it will not find the true signal that this tail belongs to, as the central velocity now lies outside the velocity range it is allowed to consider, leading to one of two scenarios. Either the tail is not fit at all by the solution, increasing the residual, or; the tail is fit by a relatively thin component that will likely still leave a non-noise residual and artificially inflate the CNM fraction found by the solution. Neither scenario is preferable. To combat this, we do a preliminary fit to the entire spectrum with low amounts of regularisation to identify the bulk

<figure><img src="image_3.png" /><figcaption>Fig. 3. Change of average $\chi^2_{red}$ value as all hyperparameters are varied for the Gamma cloud field. Shaded regions show the standard deviation around the mean.</figcaption></figure>

main-body emission of the SMC. This allows us to remove the redundant emission before segmenting the spectra to the velocity range that we are interested in. Once the main body emission has been removed we discard this part of the spectrum, focusing on a subset of channels over which the feature of interest exists. In this removal stage we are not interested in individual spectral components or concerned with the degeneracy of solution, we merely aim to recreate the shape of the signal without the noise, to remove it. Now we have a reduced dataset that requires fewer components, where a stable solution is more achievable. Figure 2 shows the mean spectrum before and after subtraction of the main body emission as well as the spectra from two non-neighbouring sightlines. The grey windows in Figure 2 show the reduced spectral range used for fitting in the next section.

### 3.2.3. Fitting processed spectra

After removal of the main-body emission, the final step is to fit the features of interest. At this point, we want to optimise the hyperparameters used by ROHSA to constrain the solution. Additionally, we keep in mind from our definition of best fit solution in Section 3.2.1, that we want to use as few components as possible. The cost function of ROHSA contains a regularisation term that is set by four hyperparameters. The first three hyperparameters act to regularise the elements of each Gaussian, that is, the amplitude (A), the central velocity (μ) and the velocity dispersion (σ) within a 1 pixel radius. The fourth hyperparameter regularises the value of σ across the entire field, as it corresponds to the phase properties of that component and should not vary to a significant degree. These hyperparameters are incorporated into the cost function ($J(\theta, v)$) as such, from Marchal et al. (2019):

$$J(\theta, v) = L(\theta, v) + \frac{1}{2} \sum_{n=1}^{N} \lambda_A \|DA_n\|_2^2 + \lambda_\mu \|D\mu_n\|_2^2 + \lambda_\sigma \|D\sigma_n\|_2^2 + \lambda'_\sigma \|\sigma_n - m_n\|_2^2, \quad (2)$$

where θ represents the positional axes and v is the velocity axis, both in the hyperspectral data. N is the number of Gaussians. λ_a, λ_μ, λ_σ, and λ'_σ are the hyperparameters of ROHSA that control the spatial variance of the fitted parameters A (amplitude), μ (central velocity), and σ (velocity dispersion). m represents a mean value for μ of each N Gaussians. For further details, refer to Marchal et al. (2019).

We explored the parameter space in orders of magnitude for each hyperparameter for the Gamma cloud field. Figure 3 shows
<table><caption>**Table 1.** Fitting parameters and results for each cloud field in this work.</caption><thead><tr><th rowspan="2">Cloud</th><th rowspan="2">n</th><th rowspan="2">λ<sub>A</sub></th><th rowspan="2">λ<sub>μ</sub></th><th rowspan="2">λ<sub>σ</sub></th><th rowspan="2">λ'<sub>σ</sub></th><th rowspan="2">χ<sup>2</sup><sub>red</sub></th><th colspan="3">Component 1</th><th colspan="3">Component 2</th><th colspan="3">Component 3</th></tr><tr><th>A<br>(K)</th><th>v<br>(km s<sup>-1</sup>)</th><th>σ<sub>v</sub><br>(km s<sup>-1</sup>)</th><th>A<br>(K)</th><th>v<br>(km s<sup>-1</sup>)</th><th>σ<sub>v</sub><br>(km s<sup>-1</sup>)</th><th>A<br>(K)</th><th>v<br>(km s<sup>-1</sup>)</th><th>σ<sub>v</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Alpha</td><td>3</td><td>1</td><td>1</td><td>1</td><td>100</td><td>1.17</td><td>2.35</td><td>108</td><td>7.77</td><td>6.51</td><td>104</td><td>1.97</td><td>0.43</td><td>79.94</td><td>9.35</td></tr><tr><td>Hook</td><td>3</td><td>1</td><td>10</td><td>10</td><td>100</td><td>1.19</td><td>2.87</td><td>135</td><td>3.99</td><td>0.71</td><td>124</td><td>6.17</td><td>2.84</td><td>132</td><td>1.32</td></tr><tr><td>Gamma</td><td>3</td><td>10</td><td>10</td><td>10</td><td>10</td><td>1.12</td><td>0.99</td><td>102</td><td>6.48</td><td>1.81</td><td>113</td><td>1.29</td><td>2.65</td><td>107</td><td>1.48</td></tr></tbody></table>

**Notes.** Column 1 gives the number of components fit to the field, columns 2–5 give the values of each hyperparameter for the fit with ROHSA, column 6 gives the reduced chi-squared value for each fit. The remaining columns show for each Gaussian component that was fit to the fields, the mean values of the amplitude (A), central velocity (μ), and velocity dispersion (σ).

how the reduced chi-squared value ($\chi^2_{red}$) for the Gamma cloud field changes as a function of the hyperparameter value. Evidently, the hyperparameters that most affected the $\chi^2_{red}$ were $\lambda_A$ and $\lambda'_\sigma$. On the other hand, $\lambda_v$ and $\lambda_\sigma$ did not change the $\chi^2_{red}$ very much.

This field was explored extensively by initialising ROHSA with each hyperparameter varying over 5 orders of magnitude. This amounts to $5^4 = 625$ runs over a single field. This is very computationally expensive way to narrow down the correct hyperparameters for a single field. So for the remaining 2 fields, a more scaled-back approach was taken. We only explored a hyperparameter space one order of magnitude either side of the best fit hyperparameters that were obtained for the Gamma cloud field (indicated in Table 1).

This smaller exploration resulted in a different outcome, that the hyperparameter values did not affect the average $\chi^2_{red}$ values by an appreciable amount. We found it was difficult to ascertain, from just the average value of the $\chi^2_{red}$, what the best fit solution was. Thus, we have the stipulation in our definition that the residuals must be normally distributed. If the field is under fit the distribution would be positively skewed, and on the other hand if it was over fit, the distribution would be negatively skewed.

The number of Gaussians fit to each field were selected through a similar exploration of the parameter space. From visual inspection of the spectra in each cloud, it was estimated that at least two Gaussian components were required to capture the full signal from each cloud. From this starting point, we obtained ROHSA solutions for each cloud ranging from 2 to 4 Gaussians. We assessed the fits of differing amount of Gaussians by inspecting the map of reduced chi-squared values calculated for each field for each line of sight. In a well-fit solution, this map should resemble a grid of random noise. We found in all three fields that we were able to achieve this noise-like reduced chi-squared map with three Gaussians as a minimum. In all fields, using only two Gaussians left parts of spectra unfit and there was clear structure in the reduced chi-squared maps. We achieved very slight reductions in the mean reduced chi-squared values with four Gaussians, however the reduced chi-squared maps did not appear to change significantly. It is important to cease adding additional Gaussian components to spectra once the best-fit criteria are satisfied. It is also valuable to note that ROHSA allows Gaussian components to have an amplitude of zero if that solution minimises the cost function (Equation (2)). So, to say that a three Gaussian solution is fit to the field does not denote that all three Gaussians are actively contributing to the solution for every spectrum in the field. It should be thought of as the maximum number of Gaussians fit to any one spectrum in the field. So after inspecting the distribution of residuals and verifying that they are centred around zero we selected the

three Gaussian solutions for each of the field. These solutions are discussed in Section 4.1.

### 3.3. Error estimation

Errors are not obtained directly from the ROHSA output, so we follow the method outlined in Taank et al. (2022) in Section 4.2.2. This involves three distinct methods of estimating error and taking their quadratic sum as the total uncertainty. These methods all test the stability of the solution. The first method involves taking the solution and adding random noise onto it, the levels of which are determined by the noise calculated in Section 2.2. The ROHSA algorithm is then run on these noise realisations with the same hyperparameter values as the original solution. The second method takes one of the noise realisations and runs ROHSA with hyperparameter values that vary by up to 10% around the values used to achieve the original solution. The hyperparameter values for each of these runs were drawn from a uniform distribution within 10% of the original hyperparameter values. The third method takes the same noise realisation used in the previous method and varies the initial guess ROHSA uses to fit the top-level solution. The spread of fitted values for each component of each Gaussian fitted in the original provides a reasonable range from which to select the initial guesses, as it is sufficiently variable, without being random. The range was thus taken as the FWHM of the distribution of each Gaussian component around its mean. The initial guesses for each of these runs were drawn from a uniform distribution within the aforementioned range. For each of these three methods, 100 runs were completed, totalling 300 runs with ROHSA. To assess the uncertainty from each method, the standard deviation in each Gaussian parameter from the 100 runs was calculated. To then obtain the total uncertainty ($\sigma_{tot}$) from these three methods, we take the square root of the quadratic sum as per Equation (3),

$$ \sigma_{tot} = \sqrt{\sigma_{rn}^2 + \sigma_{hp}^2 + \sigma_{ig}^2} \quad (3) $$

where $\sigma_{rn}$ is the uncertainty from the first method with random noise, $\sigma_{hp}$ is the uncertainty from the second method with hyperparameter variation, and $\sigma_{ig}$ is the uncertainty from the third method with initial guess variation. In Taank et al. (2022), they found that the initial guess variation was the largest source of error out of the 3 for their 11 Gaussian fit. Once we calculated our errors from all 3 sources, we found the error from the noise realisations is consistently the largest source of error. The ratio of the contribution from noise realisations, hyperparameter change and initial spectrum for the Alpha field is 48:25:27, for the Hook field is 44:14:42, and for the Gamma field is 55:32:13. In all 3 fields the noise uncertainty contributes most to the overall
<figure><img src="image_4.png" /><figcaption>Fig. 4. Spatial maps of the fitted properties of the Alpha cloud from ROHSA within a $1 \times 10^{20} \text{ cm}^{-2}$ column density contour for each phase. Top: total CNM, Bottom: total WNM. Left: column density ($N_H$), Centre: central velocity ($\mu$), Right: velocity dispersion ($\sigma$). The grey ellipse in the column density maps indicates the beamsize.</figcaption></figure>

uncertainty, whereas the other two sources have differing effects on each of the fields. We find that varying the initial guess does not dominate the uncertainty in these fields as it does in Taank et al. (2022), due to the relative simplicity of our 3 Gaussian models compared to their 11 Gaussian model. Degenerate solutions become more abundant when the number of Gaussians is increased and thus solutions with fewer Gaussians will be more robust against changes to the initial guess.

# 4. Results

## 4.1. Best-fit models

After finding solutions that satisfy our previous definition of a best fit, we list the best fit parameters in Table 1. We classify the components into the different HI phases based on the mean maximum kinetic temperature ($\langle T_k \rangle$) measured from the mean Gaussian dispersion ($\langle \sigma \rangle$) for each component. The maximum kinetic temperature is the temperature of the gas if there was no turbulent broadening of the HI signal. In reality there is some contribution from turbulent broadening, so the maximum kinetic temperature is an upper limit on the true temperature of the gas. The maximum kinetic temperature is derived from the following relation from Draine (2011):

$$ \langle T_k \rangle = 121 \langle \sigma \rangle^2. \qquad (4) $$

Gaussian components with $\langle T_k \rangle < 500$ K are classified as CNM, those with $500 < \langle T_k \rangle < 5000$ K are classified as UNM, and those with $\langle T_k \rangle > 5000$ K are classified as WNM, to be consistent with the phase definitions in Heiles & Troland (2003).

For the Alpha cloud we obtained a best fit model with 3 components, with 2 of these components clearly corresponding to the cloud. The CNM component has $\langle T_k \rangle = 470$ K and the WNM component has $\langle T_k \rangle = 7305$ K. In Figure 4 the CNM component traces the strong filamentary structure of the

cloud as we see it in the integrated data, whereas the WNM component has a more extended distribution. However, the two components overlap in the velocity axis (seen in Figure 4), showing the components are related to each other in the position and velocity axes. The third unrelated component we obtain for this fit has a large average velocity dispersion from which it would be classified as WNM. On average, its velocity offset is $25-30 \text{ km s}^{-1}$ from the cloud, thus it is likely diffuse, low-level WNM emission at the extreme velocity ends of the SMC. Within the high column density areas of the cloud, this unrelated component contributes negligibly to the total column density. For this fit we achieve a column density weighted mean $\chi_{red}^2$ value of 1.17 and we show the column density weighted mean values for each component parameter in Table 1.

For the Hook cloud we obtained a best fit model again with 3 components, with just 2 of these components clearly corresponding to the cloud. In this case we do not detect significant WNM emission, but rather a cooler UNM component. The CNM component has $\langle T_k \rangle = 211$ K and the UNM component has $\langle T_k \rangle = 1926$ K. The CNM component traces the strongest part of the Hook shape along the west and into the north, but does not extend into the eastern part of the field. The UNM component envelops the CNM component and loops all the way around the field from west to east. From Figure 5 we can see that they have the same central velocity for the areas where they overlap on the sky, discussed further in Section 4.4. The 3rd component for this fit is unrelated to the Hook cloud. It does not trace any structure of the cloud, and does not follow the same velocity trend. It is likely a tracing diffuse WNM signal, as in the Alpha cloud. For this fit we achieve a column density weighted mean $\chi_{red}^2$ value of 1.19 and we show the column density weighted mean values for each component parameter in Table 1.

For the Gamma cloud we obtained a best fit model again with three components, with all three components contributing to the total HI of the cloud. There are two CNM components, with $\langle T_k \rangle = 201$ K and $\langle T_k \rangle = 265$ K and one WNM component
<figure><img src="image_5.png" /><figcaption>Fig. 5. Spatial maps of the fitted properties of the Hook cloud from ROHSA within a 5 × 10¹⁹ cm⁻² column density contour for each phase. The panels are the same as in Figure 4.</figcaption></figure>

with $\langle T_k \rangle = 5081$ K. Like in the Alpha cloud, we see that the CNM and WNM components overlap in velocity space, shown in Figure 6. This fit required 2 components to fully capture the CNM distribution as there are strong, narrow signals at offset velocities in parts of this field. Table 1 shows that they are on average separated by 6 km s⁻¹. The advantage of the regularisation conditions that ROHSA imposes is that we can trace how these components evolve individually over the spatial axes, further discussed in Section 4.4. For the purpose of mapping the velocity and velocity dispersion of the total CNM, we take the column density weighted mean of the two components in Figure 6. For this fit we achieve a weighted mean $\chi^2_{red}$ value of 1.12 and we show the column density weighted mean values for each component in Table 1.

## 4.2. Phase distribution

The CNM fraction of any particular HI structure highlights areas where the different phases dominate in column density. We calculate the column density of HI ($N_{HI}$) with the assumption the HI is optically thin (from Roberts 1975):

$$N_{HI}(x, y) = 1.823 \times 10^{18} \int T_B(x, y, v) dv, \quad (5)$$

where $N_{HI}$ is the column density, $T_B$ is the brightness temperature along the positional axes (x and y) and velocity axis (v).

We take the optically thin regime as a reasonable assumption since these are low column density areas. In the Dempsey et al. (2022) absorption study, they showed that the HI column

density correction factor for optical depth is less than 1.05 for uncorrected column densities below $10^{21}$ cm⁻².

As we model $T_B(x, y, v)$ as a sum of N Gaussians along the v axis, we can simplify Equation (5) that for any component n, the column density is:

$$N_{n,HI}(x, y) = 1.823 \times 10^{18} \sqrt{2\pi} A_n(x, y) \sigma_n(x, y), \quad (6)$$

where $A_n$ and $\sigma_n$ are the amplitude and velocity dispersion respectively of the n-th component.

Additionally we define the CNM fraction in this work as such:

$$f_{\text{CNM}} = \frac{N_{\text{(HI,CNM)}}}{N_{\text{(HI,CNM)}} + N_{\text{(HI,UNM)}} + N_{\text{(HI,WNM)}}.} \quad (7)$$

In Figure 7 in the Alpha cloud, we see that the fraction of CNM decreases from south-west to north-east, with the exception of the cold clump at the north-east end of the cloud (RA 01:09:31, Dec −71:16:47). At these velocities (100–120 km s⁻¹) the main emission from the SMC lies to the south-west of this field, so the closer to the main body of the SMC, the higher the CNM fraction. Additionally, at this position of the Alpha cloud in projection, the emission of the main body of the galaxy is seen at $v > 135$ km s⁻¹. Looking at the velocity maps in Figure 4 the south-west region deviates less from the main body than the north-east. In summation, the further we go from the main body emission, in velocity or spatially, the lower the CNM fraction typically gets. The only area that bucks this trend is the compact cold clump in the north-east. It is centred at a highly deviant velocity from the SMC and has a CNM fraction of ~0.6.
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Spatial maps of the fitted properties of the Gamma cloud from ROHSA within a $5 \times 10^{19} \text{ cm}^{-2}$ column density contour for each phase. The panels are the same as in Figure 4.</figcaption></figure>

The Hook cloud has a stronger contribution from the UNM phase than either of the other two clouds. The maximum value for the CNM fraction is $f_{\text{CNM,max}} = 0.67$ and changes as we move around the cloud. The western side of the cloud is where the CNM fraction is highest, especially towards the centre of the filament. This means the CNM is shrouded by an envelope of UNM. The main body emission peaks at $150 \text{ km s}^{-1}$ in this field and

at these velocities the main body emission lies off to the south and south-west of this field. There is no strong relation between the velocity structure and the CNM fraction, as the CNM fraction seems to decrease in the northern most part, whereas, the velocity structure follows a south to north trend, detailed further in Section 4.4. The CNM fraction seems to increase with the thickness of the filament along the western edge.
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** CNM fraction ($f_{CNM}$) of all three clouds (Left: Alpha, Middle: Hook, Right: Gamma) within a $1 \times 10^{20}$ cm$^{-2}$ total column density contour for the Alpha cloud and within a $5 \times 10^{19}$ cm$^{-2}$ total column density contour for the Hook and Gamma clouds.</figcaption></figure>

<table><caption>**Table 2.** Physical properties of each cloud, derived as explained in Section 4.3.</caption><thead><tr><th>Cloud</th><th>M<sub>CNM</sub><br>(10<sup>4</sup>M<sub>☉</sub>)</th><th>e<sub>M,CNM</sub><br>(10<sup>4</sup>M<sub>☉</sub>)</th><th>M<sub>WNM</sub><br>(10<sup>4</sup>M<sub>☉</sub>)</th><th>e<sub>M,WNM</sub><br>(10<sup>4</sup>M<sub>☉</sub>)</th><th>M<sub>tot</sub><br>(10<sup>4</sup>M<sub>☉</sub>)</th><th>e<sub>M,tot</sub><br>(10<sup>4</sup>M<sub>☉</sub>)</th><th>W&#773;<br>(pc)</th><th>n<sub>CNM</sub><br>(cm<sup>-3</sup>)</th><th>e<sub>n,CNM</sub><br>(cm<sup>-3</sup>)</th></tr></thead><tbody><tr><td>Alpha</td><td>1.20</td><td>0.19</td><td>0.21</td><td>0.03</td><td>1.40</td><td>0.19</td><td>23.1</td><td>5.14</td><td>1.38</td></tr><tr><td>Hook</td><td>1.34</td><td>0.21</td><td>3.86</td><td>0.61</td><td>5.21</td><td>0.65</td><td>33.7</td><td>1.76</td><td>0.48</td></tr><tr><td>Gamma</td><td>1.33</td><td>0.21</td><td>0.68</td><td>0.11</td><td>2.01</td><td>0.24</td><td>20.6</td><td>1.28</td><td>0.47</td></tr></tbody></table>

**Notes.** Columns 1–6 give the CNM, WNM, and total masses of each cloud and their respective uncertainties. Column 7 gives the measured width of the CNM component. Columns 8–9 give the number density of the CNM from Equation (8) and its uncertainty.

The Gamma cloud has strong contributions from the CNM and WNM, with the majority of the CNM residing at the south end of the cloud and the WNM residing at the north end. There is a definite gradient in the CNM fraction as we move up the cloud. It transitions from completely CNM to completely WNM from bottom to top as shown in Figure 7. In the southern part of this cloud, the CNM fraction approaches 1, which is the highest we record for any of the three clouds. Whereas the other clouds point to a scenario where the CNM is shrouded in a more diffuse WNM envelope, this cloud has no WNM envelope at its southern end. At this clouds velocity (100–112 km s$^{-1}$) the main body emission is south of the cloud and again at the clouds location, the main emission of the SMC starts at 125 km s$^{-1}$ and peaks at 145 km s$^{-1}$. So, similarly to the general trend seen in Alpha cloud, the Gamma cloud has a higher CNM fraction in the area that is closer to the main body of the galaxy. Additionally, as the CNM fraction decreases, the velocity increasingly deviates from the peak emission in this area.

### 4.3. H1 mass and density

In Table 2 we show the calculated masses for each cloud separated by phase and in total as well as the CNM volume density. For these structures, calculating an H1 volume density is non-trivial. Hyperspectral data gives us no information about the line of sight structure, thus we have no way of empirically measuring how extended these structures are along the line of sight. In previous studies, the H1 density has been calculated by assuming that the depth of a structure along the line of sight is the same as the width (For et al. 2016; McClure-Griffiths et al. 2006). We

adopt a similar approach. Since these structures are not spherical, we cannot fit a 2D Gaussian to the column density images of these clouds to reliably estimate the width, so we treat them as filamentary structures. To measure the width of the filamentary structures we used the radfil package from Zucker & Chen (2018). This package defines a path through a filament and fits a Gaussian or Plummer profile to the mean profile of the filament. To define these paths through the clouds we masked pixels outside the 67% contour level and passed the masked data through the radfil profile builder. The resulting mean profile is then fit with a Gaussian profile which provides a FWHM that we take as the defining width of the filament. To get a representative density, we then divide the maximum column density at a given position along the filament by the width of the filament in angular units and accounting for the distance of the SMC, as in Equation (8) (as is done in For et al. 2016; Ben Bekhti et al. 2006):

$$n(x) = \frac{N_{H1,max}(x)}{d \tan \theta(x)}, \qquad (8)$$

where $x$ is the distance along the filament, $\theta(x)$ is the FWHM at $x$, and $d$ is the distance to the cloud (63 ± 5 kpc). We used the reported distance of the SMC from Di Teodoro et al. (2019b) in this instance, even though it may not be the distance to these individual clouds. We are unable to determine the distance to the clouds themselves without stellar associations that provide this information. The uncertainty inherent in this assumption is translated into the uncertainty in the number density measurement by including the uncertainty in the SMC distance reported in Di Teodoro et al. (2019b). These values and associated
uncertainties for the measured filament widths and subsequent number densities are reported in Table 2.

Using this method, we find the minimum values for the FWHM in each filament of around 67" (20.6 pc) which is just over twice our beamsize of 30". So the thinnest parts of the filaments are barely resolved. This could indicate that the CNM structures are more compact than we are able to measure, thus our measures of density, which range from 1 to 5 cm⁻³, should be treated as lower limits. These values are on the lower end of the typical densities expected for the CNM (Wolfire et al. 1995). For comparison, in one study of CNM filaments in the Milky Way Kalberla et al. (2016) calculate a upper limit on filament thickness of 0.3 pc. In another study of HI filaments in the Riegel-Crutcher Cloud (McClure-Griffiths et al. 2006) the magnetically dominated filaments had a typical width of 0.1 pc. Using these widths would increase the density to 70–200 cm⁻³, within the high end of the expected CNM densities. To resolve structures on these scales for the SMC we would need $\le 1''$ angular resolution, which is not readily achievable with current radio interferometers.

## 4.4. Velocity structure

From the best solutions for each field we can construct velocity maps for each phase. These velocity maps are shown in the middle columns of Figures 4, 5, and 6.

In the Alpha cloud we find different velocity structures in each component. The CNM component shows an approximately $-10 \text{ km s}^{-1}$ gradient from west to east whereas the WNM component has no strong gradient from west to east. Both components are centred around the same velocity on the western side, but moving across to the eastern side end up diverging from each other by $12 \text{ km s}^{-1}$. This divergence is driven mainly by the CNM velocity gradient. The CNM velocity decreases as we move the side of the cloud furthest from the SMC main body at these velocities.

The Hook cloud has an interesting velocity structure, it has a gradient that runs south to north with decreasing velocity along both sides. The CNM and UNM follow each other very well, only offset by $0.7 \text{ km s}^{-1}$ on average where the CNM lies, along the western edge. The velocity gradient shows that the areas in which the cloud is closest to the main body emission (in the south-west direction) have the least divergent velocities from the main body. This is a trend we see in all 3 clouds, a velocity decrease as we move radially away from the SMC centre in opposition to the general velocity trend we see across the SMC. The SMC emission moves east-west across the sky as the velocity decreases from ~200 to 100 km s⁻¹.

We see a different trend in the Gamma cloud. Figure 6 shows a clear velocity gradient from $118 \text{ km s}^{-1}$ at the bottom of the cloud to $95 \text{ km s}^{-1}$ at the top of the cloud. The CNM and WNM components follow each other along this gradient, south to north along the cloud at a slight offset. They are offset from each other on average by $2.7 \text{ km s}^{-1}$. So, they both slow down relative to the SMC main body as the distance from the SMC increases. As mentioned in Section 4.2, the CNM fraction decreases along this same path, so overall in this cloud the warmer the HI, the slower the LOS velocity relative to the main body emission.

### 4.4.1. Deviation velocities

To make a measure of the deviation velocity of each cloud, a measure commonly used to categorise low, intermediate, and high velocity clouds (LVC, IVC, and HVCs) in the Milky Way,

we calculate the first moment map of the galaxy with the clouds removed. To remove the clouds, we mask every voxel within the fields analysed where the total emission from the ROHSA model is above the noise level. From this masked PPV cube, we obtain the first moment and compare the mean velocities of each cloud with the first moment velocity at their position. We define the deviation velocity ($v_{DEV}$) as in Equation (9) below.

$$v_{DEV} = v_{peak} - v_{M1} \quad (9)$$

where $v_{peak}$ is the peak velocity, the velocity at which the brightness temperature from the ROHSA model is at its maximum, and $v_{M1}$ is the first moment velocity. From Equation (9) we find that the deviation velocities of the Alpha, Hook, and Gamma clouds are $-56.4, -13.8$, and $-34.8 \text{ km s}^{-1}$ respectively.

By the definitions asserted in Wakker (2004), classifying clouds by their deviation velocity, the Alpha cloud is an IVC and the Hook and Gamma clouds are LVCs of the SMC.

### 4.4.2. Clouds in the context of the 3D morphology of the SMC

Recently, work was done in Murray et al. (2024) to investigate the distribution of the SMC HI gas along the line of sight, using distances obtained from the stellar population. They found that the stellar population of the SMC could be separated into two distinct groups, sitting at two distinct distances along the line of sight and used this information to separate the HI into ‘front’ and ‘behind’ sections. The distinction between these two components is well-determined where there are many stellar objects, which is primarily in the main body of the galaxy. The distinction becomes hazier in the outskirts of the galaxy, where there are fewer stellar objects. We compared the central velocities of the clouds in this work to the first moment maps of the front and behind components constructed in Murray et al. (2024).

From the first moment maps in Murray et al. (2024) we extracted the regions around each cloud in this work. From each subsection we calculated the median value in the front and behind components, obtaining two respective first moment velocities for each cloud region. To characterise the uncertainty on these first moment velocities, we measured the dispersion in each field. The range of first moment velocities in these fields were quite extreme, due to low column density lines of sight, so the reported uncertainty is the mean of the dispersion measured when the 10% most positive and most negative values are excluded, and the dispersion measured when lines of sight with column densities below $1.32 \times 10^{20} \text{ cm}^{-2}$ are excluded (thus above the $3\sigma$ column density sensitivity limit reported in Pingel et al. 2022).

For the Alpha field, the median velocities for the front and behind components are $166.1 \pm 5.2 \text{ km s}^{-1}$ and $120.8 \pm 22.4 \text{ km s}^{-1}$ respectively. For the Hook field, the median velocities for the front and behind components are $142.1 \pm 22.8 \text{ km s}^{-1}$ and $143.3 \pm 14.7 \text{ km s}^{-1}$ respectively. For the Gamma field, the median velocities for the front and behind components are $142.6 \pm 15.6 \text{ km s}^{-1}$ and $146.1 \pm 11.0 \text{ km s}^{-1}$ respectively. The difference between the first moment front and behind velocities in the Hook and Gamma cloud fields is relatively small, eclipsed by the uncertainties on each measurement. There is a larger distinction between the two component velocities for the Alpha cloud field, however the difference between them is still within a $2\sigma$ uncertainty level. Given that the peak velocity of the Alpha cloud is $106.4 \text{ km s}^{-1}$, it is offset from front component by $59.7 \text{ km s}^{-1}$, a $12\sigma$ level, but offset by $14.4 \text{ km s}^{-1}$ and well within the $1\sigma$ level of the behind component. This suggests that
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** CO velocity comparison ($v_{CO} - v_{HI}$) for Alpha cloud (top) and Hook cloud (bottom). Size of circle markers indicate the strength of the integrated CO flux density scaled linearly.</figcaption></figure>

the Alpha cloud belongs to the behind component of the SMC identified in Murray et al. (2024).

Working under that assumption may shed some light on the distance to the Alpha cloud. In Section 4.3 we assumed a distance to all 3 clouds of $63 \pm 5$ kpc from Di Teodoro et al. (2019b), but if the Alpha cloud belongs to the behind component of the SMC, we could also assume it sits at the mean distance of that component reported in Murray et al. (2024) of 66 kpc. If this is true, the physical scales of the Alpha cloud will increase by 5% and the HI mass will increase by 10%.

# 5. Discussion

## 5.1. CO comparison

$^{12}$CO(2 → 1) data were obtained during targeted APEX observations by Di Teodoro et al. (2019a). It covered the entirety of the Alpha cloud and a subsection of the Hook cloud field.

Within the field encompassing the Alpha cloud, there were 8 clumps of CO detected in Di Teodoro et al. (2019a). Comparison of the HI CNM component central velocities and the CO data show that four of the clumps agree in velocity within a 3 channel (~3 km s⁻¹) window with the CNM, whereas the other four are all offset by approximately 15 km s⁻¹ from the CNM,

see Figure 8. The four clumps that are consistent with the velocity structure of the CNM, are more spatially coincident with the high column density regions of CNM than the others. This divides these clumps into two distinct populations, one that likely belongs to the Alpha cloud and one that does not. Interestingly, the clumps of CO that lie off the highest column density areas of this feature we consider here can not be matched with any strong HI emission. There is very little HI emission in this region at $v < 93$ km s⁻¹ that is significant (at $3\sigma$ above the noise).

When comparing the CO data to the HI WNM component central velocities, they agree within the same velocity window as the CNM component, except on the eastern side. This is due to the divergence of the CNM from the WNM in velocity space, which we note in Section 4.4. From this, we conclude that the CO is more dynamically aligned with the CNM than the WNM.

Within the field encompassing the Hook cloud there were nine clumps detected in CO. This cloud spans three times the size of the Alpha cloud in both spatial dimensions, so the entirety of the Hook cloud was not observed in CO APEX observations from Di Teodoro et al. (2019a). The observations cover the areas of strong total HI emission in the base of the western edge of the cloud, shown in Figure 8. Seven out of nine clumps are in agreement with the CNM within 3 km s⁻¹. These same seven clumps are offset slightly more from the UNM velocities. The
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Left: fitted component central velocities for the different phases for two different paths through the Alpha cloud (top and bottom). The WNM is shown in red, the CNM in blue with the darker colours indicating higher column density. The error bars in the top panels show the $1\sigma$ uncertainties from Section 3.3. The error bars in the bottom panels show the combined uncertainties at the $3\sigma$ level. Right: path through the cloud for respective velocity relations shown in purple. 0 denotes the start of the path, and 128 and 145 denote the ends of each path.</figcaption></figure>

two clumps that strongly disagree with this trend are slightly spatially offset to the CNM filament, however one clump that agrees with the CNM velocities is also offset from the CNM filament. These two anomalies are also among the clumps with the lowest integrated flux densities. Clumps with higher integrated flux densities and more spatially coincident with the high CNM column density regions, agree best with the central velocities of the CNM component. This indicates that the strongest CO clumps are associated with the dense CNM regions of the Hook cloud. Since there is little difference in the velocities of the UNM and CNM, we could just as equally suppose that the CO is associated with the UNM envelope so that all this gas is travelling together.

Overall, in both clouds, the CO is more dynamically aligned with the CNM, significantly in the Alpha cloud and marginally in the Hook cloud, particularly in areas with high CNM column density.

## 5.2. Phase velocities

With the information obtained from the ROHSA decomposition, it is possible to trace the velocity of each phase through the cloud. We outlined the trends in the velocity gradients of the individual phases in Section 4.4, but discuss here the relationship between the two different velocity gradients in different parts of each cloud. In the Alpha cloud when plotting different paths through the cloud towards the south-west end, an offset between the two phases becomes apparent. In Figure 9 we show

two paths through the Alpha cloud, one from east to west and another from north to south. Both paths terminate at the concentration of emission in the south-west of the cloud and show a significant offset in the velocities of the two phases. This offset is most obvious in the bottom-left panel of Figure 9 where it reaches ~5 km s⁻¹. This offset is also seen in a region of high column density for both phases with relatively low uncertainties on the central velocity, providing strong evidence for a velocity offset between the phases towards this region of the cloud.

For the Hook cloud, there is very little difference in the central velocities of each phase, and as such, there is no significant consistent or increasing offset between the two phases as we move along the curved filament that defines the Hook cloud. There is a slight offset of less than 1 km s⁻¹ in the densest area of the cloud towards the southern end, but this offset is within the $3\sigma$ uncertainty window.

Finally, for the Gamma cloud the velocity difference between the phases is a little more complex. There are two CNM components that need to be averaged to obtain a representative central velocity for the total CNM. Both of these CNM components have higher uncertainties than the CNM components of the other two clouds. This is likely because these components have similar central velocities and were easily confused with each other amongst the many ROHSA runs completed for the bootstrapping uncertainty method described in Section 3.3. These larger uncertainties are compounded when taking the mean of the two values and this makes it difficult to have confidence in offsets of the
same magnitude as in the Alpha cloud. Additionally, the measure
of the offset becomes less meaningful when the column density
of either phase falls below the noise limit which poses problems
for this cloud. The CNM is most prominent at the southern end of
the cloud and the WNM at the northern end, so there is not much
area where the phases overlap at high column densities. There
are a few points in the northern end of the cloud which have sig-
nificant CNM column densities where an offset of ~4 km s⁻¹
can be measured at levels exceeding 3σ. From just these points,
it cannot be concluded that there is a consistent offset between
the phases, but it suggests that an offset is present at least in the
northern end of the cloud.

5.3. Factors affecting morphology

5.3.1. Alpha cloud

The Alpha cloud has a complex structure, the most irregular of
the three clouds considered in this work. It does not have just one
concentration of cold material, with the north-east and south-
west ends of the cloud each containing significant amounts of
CNM. This cloud was also analysed in Pingel et al. (2022), where
they looked at the mean properties of the cloud. When they took
a mean spectrum of the cloud and fit a Gaussian they obtained a
σ value of 2.93 km s⁻¹, which is higher than our obtained value
of 1.97 km s⁻¹. This discrepancy is likely due to the presence
of the WNM uncovered in this work which would broaden the
mean profile. Additionally, in this work we found that the CNM
and WNM are not consistently offset from each other in velocity,
so this would also add to the broadening of the mean spectrum.
These both explain the broader dispersion derived in Pingel et al.
(2022) compared to the CNM dispersion in this work.

The Alpha cloud also has an obvious cavity on the west-
ern side of the field. This was catalogued as an HI shell in
Staveley-Smith et al. (1997) with a heliocentric velocity of
117.7 km s⁻¹, corresponding to a local standard of rest (LSR)
velocity of 108.7 km s⁻¹. We measure a maximum velocity
of 113.4 km s⁻¹ (LSR) for the CNM component in this field
which is consistent with the shell velocity considering our veloc-
ity resolution of 0.98 km s⁻¹ and the velocity resolution in
Staveley-Smith et al. (1997) of 1.6 km s⁻¹. This cavity was also
noted in the analysis in Pingel et al. (2022) where there was no
conclusive stellar association made due to the absence of radial
velocity information from the Gaia release (Gaia Collaboration
2021). An expanding shell could explain the absence of HI,
specifically the CNM, in this region, but it does not explain the
offset in velocity between the phases seen in the south-western
side of the shell, where the column density is highest. The off-
set in velocity, coupled with the higher CNM fraction at this
end of the cloud, suggests stellar feedback from the direction the
main galaxy is acting on the cloud and has managed to move the
less dense WNM more efficiently than the CNM. This feedback
has begun to strip the cloud of its WNM envelope at this posi-
tion, leaving the CNM and CO behind, similar to the head-tail
structure seen in some HVCs and IVCs.

An additional piece of information available in this area is
the Hα data from the Magellanic Cloud Emission Line Survey
(MCELS) (Winkler et al. 2015). Figure 10 shows the MCELS
Hα where there is strong emission in the southern region of
the cloud. Around this Hα emission, there is an enhancement
of the CNM fraction, except on the southern side where there
is very little HI. The positional coincidence of this Hα emis-
sion with the HI cloud suggests they could be related, but there
is no spectroscopic1 data available for this region at sufficient

<figure><img src="image_10.png" /><figcaption>Fig. 10. MCELS Hα data with the HI column density of the Alpha cloud overlaid in contour levels from 1 × 10²⁰ to 4.2 × 10²⁰ cm⁻².</figcaption></figure>

resolution to make an association between the cloud and the Hα
emission.

Overall, the Alpha cloud seems to be shaped by multiple pro-
cesses. Star formation that is facilitating expansion, or turbulent
motions as suggested in Pingel et al. (2022), have created a cav-
ity in the centre of the cloud, while stellar winds from the main
galaxy have stripped the WNM envelope in the south-western
head of the cloud. Additionally, compact Hα emission in the
south of the cloud could indicate stellar activity may be ionising
the surrounding HI, if it is associated with the Alpha cloud.

5.3.2. Hook cloud

The Hook cloud seems to have a simpler morphology than the
Alpha cloud, with the CNM and UNM matching each other well
in position and velocity. There is no velocity gradient along the
length of the cloud as we move further away from the main
body of the galaxy. This suggests that the filament is not strongly
affected by star formation from within the SMC main body. As
it has a concentration of CNM on the western side, with a UNM
envelop that also loops around to the eastern side, it seems more
likely to be caused by an expansive force from the centre of
the field. It is not spherical, which could be due to expansion
occurring at multiple points within this area, causing uneven
elongation along the different axes.

The source of this expansive force could come from stellar
objects in the area. Martínez-Delgado et al. (2019) identified an
arm of stellar objects that extends into this field that could be
injecting energy into the diffuse medium of the SMC outskirts,
creating CNM filaments. There are multiple stellar clusters
within this field, catalogued in Bica et al. (2020), which could
provide the energy needed to create a shell. However, all of
these clusters, for which we have age measurements, are old
(t > 100 Myr) and no shell in the most recent catalogue of HI
shells in the SMC (Staveley-Smith et al. 1997) exceeds 40 Myr in
age. We can characterise the energies associated with the poten-
tial shell by using the information from the HI. We are unable
to measure an expansion velocity of this shell-like structure, as
there are no apparent front and back walls in HI spectra from the
centre of the field. However, this does not preclude the existence
of front and back walls as they may be below the noise limit.
So for the Hook cloud, we therefore assume expansion veloc-
ties of v<sub>exp</sub> = 10–20 km s<sup>-1</sup>, as this is within the range found in
Staveley-Smith et al. (1997) for the larger shells. Using the total
HI mass shown in Table 2 of 5.21 M☉, we calculate the kinetic energy of the shell to be 0.5–2.1 × 10⁵⁰ erg.

We can also define a kinetic age ($t_k$), where $t_k = R_s/v_{exp}$, using the mean of the two shell axes as the representative radius of the shell ($R_s$). The two axes of the shell measure at 360 and 680 pc, thus $R_s = 520$ pc, giving the kinetic age $t_k = 25–51$ Myr. Following the formalism in McCray & Kafatos (1987) we can also calculate a timescale and number of stars required to produce this shell if there was constant supernovae energy injection using $R_s = 97$ pc ($N_*/n_0)^{1/5} t_7^{3/5}$ and $v_{exp} = 5.7$ km s⁻¹ ($N_*/n_0)^{1/5} t_7^{-2/5}$, where $t_7$ is the time in units of 10 Myr, $N_*$ is the number of massive stars, and $n_0$ is the ambient density of the medium before the shell was formed. To estimate the ambient density, we take the total mass of the HI and divide by the volume of shell, using the minor axis of the shell as the depth along the line of sight. This gives an ambient density of 0.05 cm⁻³. With $R_s = 520$ pc, $n_0 = 0.05$ cm⁻³, and $v_{exp} = 10–20$ km s⁻¹, we get $t_7 = 1.5–3.1$ and $N_* = 7–57$. This range in timescale overlaps with the range we derived for the dynamic timescale.

These measurements are compatible with the Hook cloud being a shell formed within the last 50 Myr from a cluster of massive stars. The elongated expansion along the major axis of the shell is reminiscent of HI chimney structures, such as in Normandeau et al. (1996), Pidopryhora et al. (2007), and Dawson et al. (2008). In these cases the HI shell expanded much further away from the Galactic Plane of the Milky Way than they could into the Galactic Plane due to the density and pressure difference. It is also common to see a cap of HI that forms the top of the chimney, which could be the northern edge of the Hook in this case.

The MCELS data for this field shows low-level Hα emission along the southern edges of the Hook on both sides which was also highlighted in McClure-Griffiths et al. (2018), but as this field is at the edge of the MCELS data range, we have no information about the Hα emission for the cap of the Hook. An HI chimney-like expansion history seems the most probable explanation for the Hook cloud, which would make it one of the largest shells identified in the SMC, when comparing with the catalogue presented in Staveley-Smith et al. (1997) and expanded upon in Hatzidimitriou et al. (2005).

### 5.3.3. Gamma cloud

The Gamma cloud exhibits a head-tail structure that is commonly seen in HVCs. This could indicate that the cloud is experiencing infall, with simulations in Heitsch et al. (2022), Konz et al. (2002) and Quilis & Moore (2001) all showing that the heads of infalling clouds are typically colder than the trailing tail of the cloud. This is consistent with what we see in the Gamma cloud and is further supported by the fact that the cloud is elongated roughly along the axis that points towards the dynamical centre of the SMC, indicating the path it has taken. However, this would also be seen if a neutral cloud that already existed in the periphery was subject to ram-pressure stripping from wind from the main galaxy. The cloud would be elongated along the direction of the stellar feedback, and the cold core would remain at the base of the cloud while the WNM is stripped and pushed away to trail behind the cold core. In the outflow scenario, we are catching the cloud at a point before the unshielded CNM core has begun to dissociate. This is similar to the outflow scenario proposed in Noon et al. (2023) for smaller Milky Way clouds. The clouds analysed in Noon et al. (2023) also had varying amounts of CO present, depending on the age

of the cloud, which was not available for the Gamma cloud in this analysis.

As with the Alpha cloud, the MCELS data shows that there is significant Hα emission at the base of the Gamma cloud, shown in Figure 11, that is also spatially conincident with a stellar cluster HW32, catalogued in Maia et al. (2014). There is no radial velocity information available at sufficient resolution to confirm a dynamic association with the Gamma cloud, but the spatial coincidence is compelling. This region could be the source of the stellar activity that is ionising and pushing the WNM away from the base of HI cloud. However, the cluster would need to provide sufficient energy to elongate an HI cloud over 560pc. So a combination of feedback from this cluster and the SMC could be a possible formation scenario.

## 6. Conclusions

From this work, we have found that in the periphery of the SMC, there is a significant amount of cold material in cloud structures detached from the main body of the SMC. This result complements the results of previous absorption surveys of the SMC in building a picture of the CNM distribution for the SMC. Our main results are:

- significant maximum values of the CNM fraction ($f_{CNM}$), ranging from 0.6 to 1 across the three clouds in a low metallicity environment;

- preferential spatial and dynamical agreement of observed CO clumps with the CNM over the WNM;

- new mass measurements for each cloud showing they all have similar masses, of the order of $10^4$ M☉, in line with the results from McClure-Griffiths et al. (2018); and

- first estimates of the CNM number density providing a lower limit of the order 1 cm⁻³ for these clouds, limited by spatial resolution.

Our analysis has also shown that the Alpha and Gamma clouds are likely shaped by stellar feedback from the direction of the

<figure><img src="image_11.png" /><figcaption>Fig. 11. MCELS Hα data with the HI column density of the Gamma cloud contour levels from 7 × 10¹⁹ to 42 × 10¹⁹ cm⁻². Green marker shows the position of HW 32 a young stellar cluster of the SMC.</figcaption></figure>
main body of the SMC, as well as possible nearby stellar clusters, stripping the WNM envelope away from the shielded CNM at the bases of these clouds. However, we cannot exclude the infall scenario for the Gamma cloud. In contrast to the other two clouds, the Hook cloud is likely an HI supershell showing signs of chimney-like blowout, one of the largest shells catalogued in the SMC.

**Acknowledgements.** We would like to thank the anonymous referee for their comments that helped improve the quality of this work. We would like to acknowledge the Interstellar Institute programs I15 and I16 for fostering collaboration and discussion around the work presented in this paper. This scientific work uses data obtained from Inyarrimanha Ilgari Bundara/the Murchison Radioastronomy Observatory. We acknowledge the Wajarri Yamaji People as the Traditional Owners and native title holders of the Observatory site. CSIRO's ASKAP radio telescope is part of the Australia Telescope National Facility (https://ror.org/05qajvd42). Operation of ASKAP is funded by the Australian Government with support from the National Collaborative Research Infrastructure Strategy. ASKAP uses the resources of the Pawsey Supercomputing Research Centre. Establishment of ASKAP, Inyarrimanha Ilgari Bundara, the CSIRO Murchison Radioastronomy Observatory and the Pawsey Supercomputing Research Centre are initiatives of the Australian Government, with support from the Government of Western Australia and the Science and Industry Endowment Fund. This paper includes archived data obtained through the CSIRO ASKAP Science Data Archive, CASDA (http://data.csiro.au). EMDT was supported by the European Research Council (ERC) under grant agreement no. 101040751. This work made use of the numpy Harris et al. (2020) and matplotlib Hunter (2007) python packages. Additionally, this work made use of Astropy: http://www.astropy.org a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022).

## References

* Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

* Ben Bekhti, N., Brüns, C., Kerp, J., & Westmeier, T. 2006, A&A, 457, 917

* Besla, G., Kallivayalil, N., Hernquist, L., et al. 2012, MNRAS, 421, 2109

* Bica, E., Westera, P., Kerber, L. d. O., et al. 2020, AJ, 159, 82

* Boothroyd, A. I., Blagrave, K., Lockman, F. J., et al. 2011, A&A, 536, A81

* Dawson, J. R. 2013, PASA, 30, e025

* Dawson, J. R., Mizuno, N., Onishi, T., McClure-Griffiths, N. M., & Fukui, Y. 2008, MNRAS, 387, 31

* Dempsey, J., McClure-Griffiths, N. M., Murray, C., et al. 2022, PASA, 39, e034

* Diaz, J., & Bekki, K. 2011, MNRAS, 413, 2015

* Dickey, J. M., Mebold, U., Stanimirovic, S., & Staveley-Smith, L. 2000, ApJ, 536, 756

* Di Teodoro, E. M., McClure-Griffiths, N. M., De Breuck, C., et al. 2019a, ApJ, 885, L32

* Di Teodoro, E. M., McClure-Griffiths, N. M., Jameson, K. E., et al. 2019b, MNRAS, 483, 392

* Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton: Princeton University Press)

* For, B. Q., Staveley-Smith, L., McClure-Griffiths, N. M., Westmeier, T., & Bekki, K. 2016, MNRAS, 461, 892

* Gaia Collaboration (Luri, X., et al.) 2021, A&A, 649, A7

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

* Hatzidimitriou, D., Stanimirovic, S., Maragoudaki, F., et al. 2005, MNRAS, 360, 1171

* Heiles, C., & Troland, T. H. 2003, ApJ, 145, 329

* Heitsch, F., Marchal, A., Miville-Deschênes, M. A., Shull, J. M., & Fox, A. J. 2022, MNRAS, 509, 4515

* Hunter, J. D. 2007, CiSE, 9, 90

* Inutsuka, S.-i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49

* Jameson, K. E., McClure-Griffiths, N. M., Liu, B., et al. 2019, ApJ, 244, 7

* Kalberla, P. M. W., Kerp, J., Haud, U., et al. 2016, ApJ, 821, 117

* Konz, C., Brüns, C., & Birk, G. T. 2002, A&A, 391, 713

* Maia, F. F. S., Piatti, A. E., & Santos, J. F. C. 2014, MNRAS, 437, 2005

* Marchal, A., Miville-Deschênes, M.-A., Orieux, F., et al. 2019, A&A, 626, A101

* Marchal, A., Martin, P. G., & Gong, M. 2021, ApJ, 921, 11

* Marchal, A., Martin, P. G., Miville-Deschênes, M.-A., et al. 2024, ApJ, 961, 161

* Martínez-Delgado, D., Vivas, A. K., Grebel, E. K., et al. 2019, A&A, 631, A98

* McClure-Griffiths, N. M., Dickey, J. M., Gaensler, B. M., Green, A. J., & Haverkorn, M. 2006, ApJ, 652, 1339

* McClure-Griffiths, N. M., Stanimirovic, S., Murray, C., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 130

* McClure-Griffiths, N. M., Dénes, H., Dickey, J. M., et al. 2018, Nat. Astron., 2, 901

* McClure-Griffiths, N. M., Stanimirović, S., & Rybarczyk, D. R. 2023, ARA&A, 61, 19

* McCray, R., & Kafatos, M. 1987, ApJ, 317, 190

* Muraveva, T., Subramanian, S., Clementini, G., et al. 2018, MNRAS, 473, 3131

* Murray, C. E., Hasselquist, S., Peek, J. E. G., et al. 2024, ApJ, 962, 120

* Noon, K. A., Krumholz, M. R., Di Teodoro, E. M., et al. 2023, MNRAS, 524, 1258

* Normandeau, M., Taylor, A. R., & Dewdney, P. E. 1996, Nature, 380, 687

* Pidopryhora, Y., Lockman, F. J., & Shields, J. C. 2007, ApJ, 656, 928

* Pingel, N. M., Dempsey, J., McClure-Griffiths, N. M., et al. 2022, PASA, 39, e005

* Quilis, V., & Moore, B. 2001, ApJ, 555, L95

* Roberts, M. S. 1975, Radio Observations of Neutral Hydrogen in Galaxies (Chicago: University of Chicago Press), 309

* Rolleston, W. R. J., Trundle, C., & Dufton, P. L. 2002, A&A, 396, 53

* Russell, S. C., & Dopita, M. A. 1992, ApJ, 384, 508

* Staveley-Smith, L., Sault, R. J., Hatzidimitriou, D., Kesteven, M. J., & McConnell, D. 1997, MNRAS, 289, 225

* Taank, M., Marchal, A., Martin, P. G., & Vujeva, L. 2022, ApJ, 937, 81

* Vujeva, L., Marchal, A., Martin, P. G., & Taank, M. 2023, ApJ, 951, 120

* Wakker, B. P. 2004, in High Velocity Clouds (Berlin: Springer), 312, 25

* Winkler, P. F., Smith, R. C., Points, S. D., & MCELS Team. 2015, ASP Conf. Ser., 491, 343

* Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152

* Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278

* Zucker, C., & Chen, H. H.-H. 2018, ApJ, 864, 152
## Appendix A: Model fitting statistics

To ensure that there was no signal left unfit in addition to visually inspecting the reduced chi-squared map (outlined in Section 3.2.3, we also looked at the distribution of the residual values across the cloud subcubes. For each subcube we calculate the residuals ($R(x, y, v)$) as below:

$$R(x, y, v) = Data(x, y, v) - \sum_{n=1}^{N} A_n \exp\left(\frac{-(v - \mu_n)^2}{2\sigma_n^2}\right), \quad (A.1)$$

where $A$, $\mu$, and $\sigma$ are the fitted components (amplitude, central velocity and dispersion) of $n$ number of Gaussian components.

From this we plot the distribution of all values in the resulting residual cube and calculate the skewness ($\gamma$) (Equation A.2).

$$\gamma = \frac{\sum_{i}^{N_R} (R_i - \bar{R})^3}{(N_R - 1) \sigma_R^3}, \quad (A.2)$$

where $N_R$ is the number of residual values in the subcube, $R$ is the residual value, $\bar{R}$ is the mean residual value of the subcube, and $\sigma_R$ is the standard deviation of the residuals in the subcube.

The distributions of the residuals for all three fields are shown in Figure A.1. All of the distributions have a calculated skewness of 0–0.01 which means they are normally distributed, so we do not see any signal in the residual cubes.

<figure><img src="image_12.png" /><figcaption>Fig. A.1: Normalised distribution of the residual values of brightness temperature for all three cloud subcubes. The black dotted line indicates the expected centre at 0.</figcaption></figure>