GASKAP-HI pilot survey science I: ASKAP zoom observations of Hi emission in the Small Magellanic Cloud

Abstract

We present the most sensitive and detailed view of the neutral hydrogen ( ${\rm H\small I}$ ) emission associated with the Small Magellanic Cloud (SMC), through the combination of data from the Australian Square Kilometre Array Pathfinder (ASKAP) and Parkes (Murriyang), as part of the Galactic Australian Square Kilometre Array Pathfinder (GASKAP) pilot survey. These GASKAP-HI pilot observations, for the first time, reveal ${\rm H\small I}$ in the SMC on similar physical scales as other important tracers of the interstellar medium, such as molecular gas and dust. The resultant image cube possesses an rms noise level of 1.1 K ( $1.6\,\mathrm{mJy\ beam}^{-1}$ ) $\mathrm{per}\ 0.98\,\mathrm{km\ s}^{-1}$ spectral channel with an angular resolution of $30^{\prime\prime}$ ( ${\sim}10\,\mathrm{pc}$ ). We discuss the calibration scheme and the custom imaging pipeline that utilises a joint deconvolution approach, efficiently distributed across a computing cluster, to accurately recover the emission extending across the entire ${\sim}25\,\mathrm{deg}^2$ field-of-view. We provide an overview of the data products and characterise several aspects including the noise properties as a function of angular resolution and the represented spatial scales by deriving the global transfer function over the full spectral range. A preliminary spatial power spectrum analysis on individual spectral channels reveals that the power law nature of the density distribution extends down to scales of 10 pc. We highlight the scientific potential of these data by comparing the properties of an outflowing high-velocity cloud with previous ASKAP+Parkes ${\rm H\small I}$ test observations.

1. Introduction

The evolution of the interstellar medium (ISM), which is driven by processes such as gas infall, star formation and subsequent stellar feedback, is fundamentally linked to the formation and overall evolution‘ of galaxies. Stellar feedback in the form of outflows from sites of massive star formation greatly influences the overall anatomy of the ISM by means of the enrichment of heavy elements (Nomoto et al. 2013), mixing (Kreckel et al. 2020), and transfer of kinetic energy (Fierlinger et al. 2016). In terms of gas infall—which ultimately fuels star formation—the accretion of diffuse gas onto the discs of galaxies is the leading explanation for how galaxies have continued to form stars while retaining a nearly constant atomic neutral hydrogen ( ${\rm H\small I}$ ) content since $z\sim 2$ (Noterdaeme et al. 2012; Madau & Dickinson 2014). Hydrogen is characterised into several phases in models of the ISM (e.g., Wolfire et al. 1995; Wolfire et al. 2003; Vázquez-Semadeni 2012), assuming solar metallicty: molecular ( $\mathrm{H}_2$ ; T $\leq 50\,\mathrm{K}$ ), cold neutral medium (CNM; $\textit{T} \sim 100\,\mathrm{K}$ ), lukewarm neutral medium (LNM; $100\,\mathrm{K} \leq T\leq 6\,000\, \mathrm{K}$ ), warm neutral medium (WNM; $T\geq 6\,000\, \mathrm{K}$ ), and ionised ( $T\geq 8\,000\, \mathrm{K}$ ). While hydrogen in all its phases dominates the composition of the ISM (Wolfire et al. 1995; Wolfire et al. 2003; Kalberla & Haud 2018), the densest and coolest regions that host and fuel site of star formation contain mostly molecular gas, including the CO-dark phase consisting of $\mathrm{H}_2$ that lies outside carbon monoxide (CO) emitting regions ( $T\sim 10\,\mathrm{K}$ ; $n\sim1\,000\,\mathrm{cm}^{-3}$ ; Gnedin, Tassis, & Kravtsov 2009).

There is an abundance of theoretical work that describes the astrophysics behind the life cycle of the multi-phase ISM (e.g., McKee & Ostriker 1977; Wolfire et al. 2003; Audit & Hennebelle 2005; Vázquez-Semadeni 2012); however, there are few observational studies that thoroughly explore the interconnected roles of dense molecular gas, dust, and the cold and warm neutral phases as traced by ${\rm H\small I}$ at equivalent spectral ( ${<}1\,\mathrm{km\ s}^{-1}$ ) and physical pc scales. This paper introduces the first images from the ${\rm H\small I}$ pilot observations of the Galactic Australian Square Kilometre Array Pathfinder (GASKAP; Dickey et al. 2013) Survey, which aims to explore the relationship between ${\rm H\small I}$ and other fundamental phases of the ISM within the Milky Way and nearby Magellanic System.

The proximity of the Large and Small Magellanic Clouds (LMC and SMC), two nearby dwarf galaxies at the respective distances of 50 and 60 kpc (Hilditch, Howarth, & Harries 2005; Pietrzyński et al. 2013), makes them a unique laboratory for the detailed study of the impact from powerful outflows on the makeup of the surrounding ISM and thus the overall regulation of star formation. The GASKAP ${\rm H\small I}$ spectral line cubes possess the sensitivity, angular resolution, and spectral resolution to study the ISM at scalesapproaching the molecular gas observed with ALMA ( ${\sim}0.5''{-}3^{\prime\prime}$ ; $0.15{-}0.90\,\mathrm{pc}$ at the distance of the SMC) and comparable to dust as revealed by space-based telescopes in the far-IR such as the Spitzer Space Telescope and the Herschel Space Observatory ( ${\sim}20''{-}40^{\prime\prime}$ ; $6{-}12\, \mathrm{pc}$ at the distance of the SMC). McClure-Griffiths et al. (2018) and Dempsey et al. (2020) identified cool ${\rm H\small I}$ outflows around the SMC in emission and absorption, respectively. Furthermore, Di Teodoro et al. (2019b) detected for the first time ${}^{12}\mathrm{CO}$ gas entrained within one of these ${\rm H\small I}$ outflows at anomalous SMC velocities. The presence of molecular gas within outflowing ${\rm H\small I}$ at velocities exceeding the escape velocity suggests that the SMC could eject its entire cool gas reservoir that fuels star formation within ${\sim}1{-}3\, \mathrm{Gyr}$ , while also enriching the surrounding circum-galactic medium and ultimately feeding this gas to the Milky Way. On the other hand, the reservoir of cold gas could be replenished by gas infall, which has been suggested for regions of the LMC from the Magellanic Bridge based on anomalous kinematic features in ${\rm H\small I}$ (Indu & Subramaniam 2015). The pilot GASKAP observations presented here are ideal for the detailed study of both inflowing and outflowing ${\rm H\small I}$ and high-velocity clouds (HVCs) around the SMC and the beginnings of the Stream and Bridge, which will ultimately inform important constraints on the future evolution of the Magellanic System and Milky Way.

Turbulenceis another key process that drives the evolution of the ISM (Elmegreen & Scalo 2004; Mac Low & Klessen 2004) on a variety of length scales. For example, the accretion of circumgalactic material (Klessen & Hennebelle 2010), galactic-scale gravitational instabilities, and magnetorotational instabilities likely inject turbulent energy on large scales ( ${>}\mathrm{kpc}$ ; Krumholz & Burkhart 2016; Ibáñez-Mejía et al 2017), while stellar feedback from outflows and supernova explosions inject energy on smaller scales ( ${\sim}10\,\mathrm{pc}$ ; Padoan et al. 2016; Grisdale et al. 2017). Similar to studies of star formation, there are a multitude of theoretical works that characterise turbulence over a range of length scales and density regimes (Kowal & Lazarian 2007; Padoan, Haugbølle, & Nordlund 2012; Agertz, Romeo, & Grisdale 2015); however, few observational studies have attempted to probe turbulent properties in the context of a multi-phase ISM over similar physical scales (Stanimirović et al. 1999; Padoan et al. 2006; Burkhart et al. 2010; Pingel et al. 2013). Specific to the SMC, Szotkowski et al. (2019) used early pilot ASKAP observations to investigate spatial variations of the ${\rm H\small I}$ power spectrum and found the spectral slopes to be relatively uniform. Kalberla & Haud (2019) investigated the turbulent properties of ${\rm H\small I}$ separated into the CNM, LNM, and WNM to reveal that each phase has distinct turbulent properties. Eden et al. (2021) analysed the two-dimensional spatial power spectrum (SPS) of maps of the dense gas mass fraction derived independently from a combination of $\mathrm{H}_2$ column density estimated from dust continuum and the ratio of CO intensities to show ${\sim}10\,\mathrm{pc}$ is the characteristic scale for the largest variations in the clump and star formation efficiencies in CO-traced clouds formed by supersonic turbulence. Pingel et al. (2018) explicitly compared the turbulent properties of several multi-wavelength tracers in the Perseus molecular cloud including dust, ${\rm H\small I}$ , and CO. Each tracer showed characteristics of a distinct turbulent environment, such as self-gravitating, supersonic medium in the dust and mostly transonic medium in the ${\rm H\small I}$ . Each of these studies were able to draw conclusions on how turbulence influences the structure and dynamics of different ISM phases. However, a true multi-phase characterisation of ISM turbulence has been limited by differences in the angular and spectral resolutions of observational data.

The GASKAP survey is designed as a high angular and spectral resolution ASKAP survey of the Milky Way Galactic Plane and the Magellanic System in ${\rm H\small I}$ and hydroxyl (OH) spectral lines at frequencies of $\lambda = 21$ and $18\, {\mathrm{cm}}$ . The ${\rm H\small I}$ component of the GASKAP survey (hereafter referred to as GASKAP-HI) will probe the properties of ${\rm H\small I}$ in our own Galaxy and the nearby Magellanic System at unprecedented spatial resolution ( ${\sim}10\,\mathrm{pc}$ at the distance of the SMC), similar to the scale of individual dark and star-forming regions in Galactic molecular clouds (e.g., Perseus, Lee et al. 2012) and spectral ( ${\sim}0.5\,\mathrm{km\ s}^{-1}$ ) resolutions. The OH component of GASKAP will serve as an important probe of the molecular content of the ISM not associated with CO emission (Allen et al. 2012; Allen, Hogg, & Engelke 2015; Busch et al. 2021). The survey’s primary objective is to follow the cycle of gas evolution from diffuse, warm Hi through cold Hi, to molecular gas as traced by OH, and finally to star formation and evolution through OH masers (Dickey et al. 2013).

Many previous surveys have looked at the ${\rm H\small I}$ in the Galaxy and Magellanic System. However, previous large-area single dish surveys of the ${\rm H\small I}$ (McClure-Griffiths et al. 2009; Kerp et al. 2011; Peek et al. 2018; Martin et al. 2015; HI4PI Collaboration et al. 2016) have all lacked sufficient spatial resolution to probe the structure and dynamics at small scales. Interferometers, on the other hand, sacrifice surface brightness sensitivity (Braun & Walterbos 1985; Stanimirović 2002) for high angular resolution in order to reveal bright, small-scale structures with narrow velocity widths. GASKAP-HI builds on the Canadian Galactic Plane Survey (Taylor et al. 2003), Southern Galactic Plane Survey (SGPS; McClure-Griffiths et al. 2005), Very Large Array Galactic Plane Survey (VGPS; Stil et al. 2006), the HI/OH Recombination-line Survey of the Inner Milky Way (THOR; Beuther et al. 2016), the Southern Galactic Centre Survey (McClure-Griffiths et al. 2012), and surveys of the Magellanic System (Kim et al. 1998; Staveley-Smith et al. 1997; Stanimirović et al. 1999) by increasing the angular resolution to arcsecond scales (a factor of two increase relative to SGPS and VGPS), increasing the spectral resolution to $ 0.5\,\mathrm{km\ s}^{-1}$ (useful to probe kinetic temperatures of ${\sim}6\,\mathrm{K}$ ), and expanding the instantaneous uv-coverage to characterise a virtually continuous range of spatial frequencies. Combining these novel GASKAP-HI data with data from the 64 m Parkes single dish telescope (Murriyang) will allow, for the first time, a view of nearby ${\rm H\small I}$ emission at similar physical scales as dust and molecular gas in the Magellanic System and supplement the aforementioned surveys of the Galactic Plane. Furthermore, coupling the detailed view of the ${\rm H\small I}$ emission structure with high spectral resolution ${\rm H\small I}$ absorption measurements will facilitate accurate measurements of the spin temperature, placing lower limits on the kinetic temperatures of the WNM and LNM. The high number densities in CNM thermalise the ${\rm H\small I}$ line in this phase such that absorption measurements will provide the true kinetic temperature.

Each antenna element of the Australian Square Kilometer Array Pathfinder (ASKAP; Hotan et al. 2021) telescope is outfitted with a Phased Array Feed (PAF) receiver and a matched set of digital equipment that forms 36 simultaneous primary beams on the sky, expanding the instantaneous field of view (FoV) from ${\sim}1$ to ${\sim}25\,\mathrm{deg}^2$ . The unique wide-field imaging capability of ASKAP presents several unique data processing challenges that the GASKAP-HI survey must address before beginning full survey operation. For example, it is paramount to remove the inherent instrumental point-spread function (PSF) from a sky brightness distribution that extends over the boundaries of multiple beams without degrading the recovered structure. It is also crucial to combine ASKAP data with data from a single dish to fill in the missing short-spacings filtered out by the fixed antenna distribution (Stanimirović et al. 1999; Stanimirović 2002).

ASKAP surveys, including GASKAP, were selected in 2009 based on scientific merit during the construction phase of the telescope. The plans and specifications of each survey have been progressively refined in the intervening years to match the telescope capabilities and increased computing power, as described in Hotan et al. (2021). In late 2019, ASKAP commenced pilot observations for each survey, allocating 100 h of observing time to each survey team to test full survey operations. In this paper, we present the first GASKAP-HI pilot observations of the SMC, obtained in December 2019. The observations are one part of a multi-faceted pilot survey, covering an extended region of the SMC, the starts of the Magellanic Stream and Bridge, as well as two fields in the Milky Way Galactic Plane. These images are also the first to be produced by our custom imaging pipeline built around WSClean, a command-line application first introduced by Offringa et al. (2014) that is specifically designed to perform efficient wide-field imaging. This software allows us to perform joint deconvolution distributed over a computing cluster, which ensures accurate reconstruction of the extended ${\rm H\small I}$ emission intrinsic to the Galaxy and Magellanic System while also managing the overall memory footprint.

This paper is structured as follows: Section 2 describes the philosophy and parameters of the GASKAP-HI pilot survey; Section 3 explains the calibration procedure, provides a thorough description of our custom imaging pipeline, including the combination with Parkes data to correct for the missing short-spacings, and discusses our quality assessment methods; Section 4 outlines the noise properties of our final image cubes and derives the global transfer function to fully characterise the represented spatial frequencies; Section 5 highlights the advantages of the increased sensitivity and angular/spectral resolution through global channel maps of the SMC and presents an analysis of a known outflowing HVC; finally, in Section 6, we summarise these results from the SMC and anticipate future results from our other pilot fields and the eventual full GASKAP-HI survey.

Table 1. Summary of GASKAP-HI pilot observations. $t_{\mathrm{int}}$ is the integration time spend on each interleave position (20.9 h total integration), $\nu_0$ is the central frequency of our measurement sets (see text) in a topocentric frame, BW is the total bandwidth, B is the native spectral resolution, and $\theta_{\mathrm{PAF}}$ is the PAF rotation, which refers to the total rotation of the footprint on the sky, including the natural $-45\,\mathrm{deg}$ rotation to align with celestial coordinates. Note that the pointing centres were cycled every 10 minutes throughout the listed UTC date and time ranges.

2. ASKAP observations and the GASKAP-HI pilot survey

GASKAP-HI makes optimum use of ASKAP’s wide field-of-view and bi-modal antenna baseline distribution with peaks near 500 and 3 000 m to simultaneously image the diffuse emission with high surface brightness sensitivity and high angular resolution for point sources, such as ${\rm H\small I}$ absorption against continuum sources. GASKAP-HI will give a detailed view of gas evolution in galaxies by targeting the three very different galaxies: the Milky Way, the LMC and SMC.

GASKAP is unique amongst the ASKAP surveys in its requirement for high spectral resolution. The standard ASKAP spectral resolution is $\Delta \nu= 1\, {\mathrm{MHz}}$ for the continuum surveys (EMU; Norris et al. 2011; POSSUM Murphy et al. 2013; VAST; Anderson et al. 2021) and $\Delta \nu =18.5\, {\mathrm{kHz}}$ , optimised for extragalactic Hi surveys like WALLABY (Koribalski et al. 2020), DINGO (Meyer et al. 2017), and FLASH (Allison et al. 2020). By contrast, GASKAP is designed to fully resolve the spectral linewidths observed in cold ${\rm H\small I}$ , OH masers, and OH absorption, requiring frequency resolution of $\Delta\nu \lesssim 5\, {\mathrm{kHz}}$ . This is achieved through a ‘zoom’ mode of the ASKAP poly-phase filterbank, through which observations can be obtained with one of several different spectral resolutions between $\Delta \nu = 0.58\, {\mathrm{kHz}}$ and the standard resolution of $\Delta \nu =18.5\, {\mathrm{kHz}}$ . GASKAP-HI plans to use a spectral resolution of $\Delta \nu = 1.15\, {\mathrm{kHz}}$ , giving a velocity resolution of $\Delta v \approx 0.3\, {\mathrm{km}\, \mathrm{s}^{-1}}$ over a total bandwidth of $18\, {\mathrm{MHz}}$ spanning 15 552 fine channels. We inclusively select the fine channels $7\,887-9\,934 \sim 400\,\mathrm{km\ s}^{-1}$ to $-100\,\mathrm{km\ s}^{-1}$ ) in the topocentric reference frame, which captures ${\rm H\small I}$ emission from the Magellanic System and Milky Way, to reduce the total amount of data that needs processing. The centre frequency of each measurement set is $1419.81\, \mathrm{MHz}$ in the topocentric reference frame.

Our observations of the SMC total 20.9 h of integration, split equally into two separate sessions. Each session is known as a schedule block (SB) and identified with IDs (SBIDs) 10941 and 10944, making it easier to identify observations in the CSIRO ASKAP Science Data Archive (CASDAa). The PAF footprint, known as closepack36, was centred on $\mathrm{J}2000\,\mathrm{RA} = 00^{\mathrm{h}}58^{\mathrm{m}}43.280^{\mathrm{s}}$ , $\mathrm{Dec} = -72^{\mathrm{d}}31^{\mathrm{m}}49.03^{\mathrm{s}}$ . This particular PAF footprint places 36 simultaneously formed beams in a $6\times 6$ hexagonal grid such that the beam centres are separated by 0.9 deg to provide more uniform sensitivity (McConnell 2016; Hotan et al. 2021). Throughout the observation, every 10 min the PAF footprint cycles through three pointing centres provided in Table 1. This process, known as interleaving, ensures uniform sensitivity across the FoV. Due to beams formed offset from the pointing centre, it is necessary to track the parallactic angle on the sky throughout the interleaving process. The ASKAP antennas achieve this through a unique third axis of rotation on the reflector itself, in addition to a simple azimuth-elevation mount. Table 1 summarises the observational setup.

3. Data reduction

3.1. ASKAPsoft pipeline

ASKAP observations are calibrated and imaged using a high-performance processing pipeline (Whiting 2020; Hotan et al. 2021) that runs on the Galaxy supercomputer at the Pawsey Supercomputing Centre. This typically includes all of the steps necessary to calibrate the data and produce images and spectral cubes of the observation, using the ASKAPsoft package (Guzman et al. 2019). For the GASKAP-HI processing, since the imaging is done in the joint-imaging pipeline described below, the ASKAP pipeline was limited to bandpass calibration, flagging, and self-calibration of the time-dependent gains.

3.1.1. Bandpass and flux scale calibration and flagging

Prior to self-calibration, the ASKAP pipeline prepares the data through calibration of the bandpass and flagging of unwanted signal. The bandpass and overall flux scale is calibrated using the primary flux calibrator PKS $\mathrm{B}1934{-}638$ . The target and flux calibrator data are then flagged to remove radio frequency interference (RFI). This employs a dynamic flagging algorithm operating in the time-domain for each channel, which flags discrepant time ranges based on the Stokes V flux.

Data is then averaged to form 1 MHz-wide channels (that is, averaging 864 channels together), which will be used to create the continuum images used in the self-calibration. Those channels that are identified to contain HI emission are then flagged for this step, so that only continuum emission contributes to the resulting images.

3.1.2. Self-calibration

The derived beamforming weights are generally updated once every day or two and require 2.5 h of telescope overhead for calibration (Hotan et al. 2021). In lieu of spending time on traditional phase calibration methods—such as regularly observing a reference source near the science target—the amount of flux within the wide ASKAP FoV at any given time ensures adequate flux to utilise an in-field self-calibration scheme to correct for time-varying phase errors. The ASKAPsoft pipeline creates continuum images for each beam from the calibrated, flagged, and averaged science visibilities to use as initial models. No modelling of the spectral index is performed due to the relative narrow bandwidth.

Due to the relatively weak continuum in this field, no baselines were excluded in the continuum imaging. For future GASKAP-HI fields on the Galactic Plane or LMC, baselines below 500 m will be excluded. An image of the field is produced with the ASKAPsoft BasisfunctionMSF solver with Hogbom clean and a shallow sky model of the field is made from that image. A total of 257 w-planes are used in this image. This model is used to calibrate the complex gains over short intervals ( ${\approx}30\,\mathrm{s}$ ) across the entire observation. To preserve the flux scale, we perform phase-only self-calibration by fixing the amplitude of the calibration parameters to 1. The resulting calibration parameters are applied to both the continuum and spectral visibility data sets, and a subsequent continuum image is made using the calibrated data.

3.2. Observation diagnostics

The size of the calibrated visibilities from a typical GASKAP-HI field, consisting of roughly 20 h of integration time, 2 048 spectral channels, and 108 total beams ( $36\,\mathrm{beams}\times 3$ interleaves) in closepack36 configuration, totals ${\sim}7\, \mathrm{TB}$ . Before expending the large effort required for spectral line imaging to create an 11 GB image cube, it is important to understand the quality of an observation. This allows us to verify that the observation is suitable for its intended science use. In particular, we need to know whether

• • the sensitivity across the field is even, so that features can be compared across the field;

• • there are enough data available from short baselines to provide good sensitivity to large-scale structure and low surface brightness features;

• • there are enough data available from long baselines to adequately sample fine structure and to detect absorption;

• • the calibration bandpasses are reasonably free of structure.

We produce a diagnostics report to describe the quality of the observation.b It summarises and visualises outputs produced by the ASKAPsoft (Guzman et al. 2019) pipeline up to the production of calibrated measurement sets, along with the continuum image and source catalogue of the field. The report provides a snapshot of the overall quality and is generally sufficient to determine whether to proceed with imaging or re-observe.

The first two sections of the report provide details of the observation and the continuum image. In the observation we need to know where the field is and how it was observed—duration, beam pattern (or PAF footprint), and the correlator mode. The footprint is also used later to determine the layout of beam based plots. Noting that the continuum data is well covered by the ASKAP continuum validation report, the continuum image report is limited to basic details.

To allow assessment of the calibration we visualise the median calibration bandpasses both by beam and by antenna. In the plots we highlight up to 10 bandpasses that either have

1. (a) a median more than one standard deviation from the median of all bandpasses, or

2. (b) an amplitude range that is more than one standard deviation from the median amplitude range of all bandpasses.

The first rule highlights antenna sensitivity issues and the second calls out structure in the calibration bandpasses. As the raw ASKAP visibilities are not retained to manage disc space, recalibration is not possible, so any problems reflected here would generally mean the dataset needs to be rejected and the field re-observed.

In the diagnostics section, we largely focus on what data have been flagged by the ASKAPSoft pipeline (see Section 3.1.1) within the measurement sets and the effect this has on the suitability of the dataset for different science cases. The flagging of the measurement sets is done in four dimensions: by baseline (or combinations of antennas), by beam, by channel, and by time. Flagging is generally used to remove bad data from the dataset, whether for an antenna that is performing incorrectly, an antenna that is in the shadow of another antenna with respect to the target source, or for RFI. We provide plots of the flagging

1. (a) by baseline, (see Figure 1), which shows if there is any bias in coverage by angular scale,

   

   Figure 1. Baseline coverage for the ASKAP-36 array showing the proportion of baselines flagged and unflagged in the observation. The flagging proportion is measured across all beams and all channels. However, any time periods where the data is entirely flagged are excluded. The blue profile shows all available baselines, while the transparent orange profile denotes unflagged baselines. A histogram of the flagged baselines (excluding the automatically flagged autocorrelations) is outlined by the dotted line. The two vertical dotted lines represent the definition for short ( $l_{\text{baseline}} < 500\,\mathrm{m}$ ) and long baselines ( $l_{\text{baseline}} > 2\,000\,\mathrm{m}$ )

2. (b) by time to show if our uv-coverage is even, and

3. (c) by percentage of each baseline flagged to highlight frequency based flagging.

For example, the histogram of the baseline coverage across all beams and all channels for the entire observation for SBID 10944 in Figure 1 shows a majority of the available baselines remain unflagged during this scheduling block. Figure 1 also demonstrates that flagged baselines are evenly distributed across the entire range of available baselines. This ensures there are no unexpected gaps in the continuous uv-coverage, which would lead to potential artefacts at particular spatial scales in the final images. The two distinct peaks between 400 and 1 200 m, and 2 and 3 km are ideal for synthesised beam sizes between 30^′′ to 60^′′ to detect Galactic and Magellanic ${\rm H\small I}$ in emission. The longest baseline of 6 km produces a beam size of ${\sim}10''$ that is necessary for absorption measurements. All diagnostic plots are sourced from the flagging summary files for each measurement set. As is done in the WALLABY diagnostics (For et al. 2019; Koribalski et al. 2020), the linked diagnostic reports also provide plots of

1. (a) the fraction flagged by beam,

2. (b) the number of antennas flagged in each beam and from these,

3. (c) the expected RMS per channel for each beam of theobservation.

To assist the assessment of quality of the observation we report five metrics on a good, uncertain, and bad scale.

1. 1. The percentage of data flagged from short baselines ( $l_{\text{baseline}} < 500\,\mathrm{m}$ ). These are crucial for large-scale emission and low surface brightness sensitivity.

2. 2. The percentage of data flagged from long baselines ( $l_{\text{baseline}} > 2\,000\,\mathrm{m}$ ). These baselines are essential for tracing fine detail and for detecting absorption.

3. 3. The percentage difference of the expected RMS across the field:

   (1) \begin{equation} \text{difference} = \frac{{\mathrm{RMS}}_{\text{exp,max}}-{\text{RMS}}_{\text{exp,min}}}{{\mathrm{RMS}}_{\text{exp,min}}}\end{equation}
   which is used to show how even the sensitivity is across the field.

4. 4. The percentage of integrations with less than 5% of channels flagged. In observations of local Hi we do not want the diffuse emission spectral lines flagged out.

5. 5. The number of unflagged antennas with a phase RMS $\geq 40^\circ$ . These are antennas where the phase solutions are severely unstable over the course of the observation, leading to poor positional accuracy.

These metrics are made available both in the report and on CASDA (Huynh et al. 2020). The metrics and their threshold values are shown in Table 2. For an observation to be acceptable for emission studies, both the short baseline and the expected RMS variance tests must be in the uncertain or good ranges. For an observation to be acceptable for absorption studies, both the long baseline and the expected RMS variance tests must be in the uncertain or good ranges.

Based on these metrics, we find that the observations presented here are suitable to use for emission line analysis. Whilst antennas (denoted by the prefix ‘ak’) ak06 and ak32 are completely flagged out, the flagging is even across baseline lengths, with flagging percentages of 16% on short baselines and 25% on long baselines (denoted by the dashed lines in Figure 1). There is increased flagging on integrations around the change of interleaves. In scheduling block 10941, the last integrations for the observation in all interleaves are heavily flagged. The calibration bandpasses exhibit structure with small amplitudes, including dips at the 1 MHz beam forming intervals. These dips are not an issue for GASKAP-HI due to our relatively narrow bandpass.

3.3. Joint-imaging pipeline

The PAF receivers equipped on each ASKAP antenna extend the typical FoV for a 12 m diameter antenna from ${\sim}1$ to ${\sim}25\,\mathrm{deg}^2$ by simultaneously forming 36 distinct primary beams. Standard full-field ASKAP images are produced by deconvolving the inherent PSF response from the sky brightness distribution of each separate beam independently and then linearly mosaicking them together. This approach works for most other ASKAP surveys, as the emission is small scale (e.g., distant galaxies with small angular size) and generally is contained within a single beam. Such techniques work less well for GASKAP-HI because the targeted diffuse ${\rm H\small I}$ emission extends up to the boundaries of the primary beam. To ensure we accurately recover the distribution of structure down to baselines corresponding to the size of the primary reflectors themselves (i.e., 12’m), fill in gaps between all the short baselines, and deconvolve the images deeper, we must utilise a joint deconvolution approach. Here, a large image is created from the ourier inversion of the sampled visibilities from all beams (i.e., the ‘dirty image’) before proceeding with deconvolution. Since the calibrated visibilities from all of the beams must be managed simultaneously in a joint deconvolution approach, the computation requirements—specifically the memory footprint—are immense relative to other imaging techniques. As of publication, ASKAPsoft 1.3.0,c the primary processing software utilised to image all other ASKAP data, is not fully optimised to run a distributed joint deconvolution scheme.

Table 2. Observation diagnostics metrics.

We instead developed a custom imaging pipelined that utilises a combination of Common Astronomy Software Application (CASA; McMullin et al. 2007) tasks to pre-process the measurement sets, WSClean e to grid the multiple ASKAP pointings onto a single grid for joint deconvolution using corresponding primary beam models for each ASKAP primary beam, and miriad (Sault, Teuben, & Wright 1995) tasks to combine ASKAP and Parkes data to correct for short-spacings. This imaging package enables efficient gridding of visibilities and stable deconvolution, as demonstrated by the high-fidelity images produced of the Galactic Centre with the MeerKAT telescope (Heywood et al. 2019). We input the ASKAP primary beams measured from holography at 1.4 GHz (Hotan et al. 2021) to ensure highly accurate gridding of the gain-calibrated visibilities and scaling of the components subtracted from the dirty image during deconvolution major and minor iterations, respectively. The pipeline takes advantage of a high-performance computing (HPC) environment, while also managing the memory footprint, by submitting simultaneous batch imaging jobs of individual channels. A majority of the processing was performed with the Avatar computing cluster located at the Australian National University Research School of Astronomy and Astrophysics. The configuration of WSClean imaging mode and resource allocation is performed by bash scripts specific to both the Portable Batch System (PBS) and Simple Linux Utility for Resource Management (SLURM) scheduling software.

Figure 2 summarises the work flow of our joint-imaging pipeline to image a single GASKAP-HI field. We download from CASDA 108 measurement set files that contain the gain-calibrated visibilities for a field observed in the closepack36 beam configuration. We then use the built-in WSClean tool chgcentre to phase-rotate the visibilities to Pointing centre 1 listed in Table 1, thus setting a common direction where all phases are zero and effectively sets the centre of the resultant image. To avoid smearing and noncoplanar baseline effects (see Cornwell et al. 2005), such as the distortion of sources sitting far away from the phase centre, the phase information in the measurement sets is updated during gridding based on the individual beam phase centres provided to WSClean via a configuration file.

3.3.1. Binning and uv-based continuum subtraction

For the purpose of increasing the signal-to-noise ratio (SNR), we first bin each measurement set by four channels, achieving a spectral resolution of $4.628\,\mathrm{kHz}$ ( ${\sim}0.98\,\mathrm{km\ s}^{-1}$ at the rest frequency of ${\rm H\small I}$ ) for a single channel.

It is also necessary to remove the contribution from continuum sources in the field to ensure we infer the correct physical properties of the ${\rm H\small I}$ gas and to aid deconvolution. We select a range of ${\rm H\small I}$ emission-free channels by plotting a mean spectrum from the visibilities contained in several central beams and baselines between core antennas, which is expected to contain the brightest emission. An example spectrum is shown in Figure 3. The continuum is removed from the uv spectra for each individual beam by fitting and subtracting a first order polynomial to both the real and imaginary visibility spectra over the binned channel range from 0 to 140 (1 418.63 to 1 419.28 MHz), 300 to 340 (1 420.01 to 1 420.20 MHz), and 400 to 500 (1 420.48 to 1 420.94 MHz), as no emission from the SMC or foreground Milky Way is observed in thesechannels.

Figure 2. A flow chart summarising the workflow of our custom pipeline used to image a typical GASKAP-HI field. In brief, the phase centres for each of the 108 measurement sets, which contain the gain-calibrated visibilities per beam ( $36\,\mathrm{beams}\times3$ interleaves), are set to single value using the chgcentre utility available through the WSClean package. These altered measurement sets can be binned and put through a uv-based continuum subtraction scheme, respectively using the CASA tasks split and uvcontsub within the scripts bin_All and uvcontsub_All, that is submitted in a distributed manner for efficient processing. Single spectral channels are then split out to ensure reasonable memory management. A series of batch jobs, each utilising 4 total CPU cores, are submitted with the 108 measurement sets per spectral channel as input to the WSClean command-line imager to produce a jointly deconvolved image. The final deconvolved images are then collated into a single cube before combination with Parkes data to fill in the missing short-spacings (see Section 3.3.3). The pipeline processes are represented as rectangles, actual input/output of these processes are represented by parallelograms, and the action of input is represented by the curved arrows.

3.3.2. Joint deconvolution

Individual channels are split out from each of the 108 measurement set files to best manage the large memory footprint required by joint deconvolution. The imaging of each channel is submitted as a separate batch job that utilises a total of 4 cores to take advantage of the built-in parallelisation of multiple compute threads. Table 3 summarises the typical configuration to process a single channel with significant emission. We perform a multiscale CLEAN deconvolution (Cornwell 2008) with: a pixel size (i.e., scale parameter) of 7^′′, multiscale scales set to 0, 8, 16, 32, 64, 128, and 256 pixels (angular range between 0^′′ to 30^′), the total number of iterations (niter) set to 25 000, and mgain set to 0.7, meaning that the peak flux must be reduced by 70% to trigger the next major iteration; typical runs go through two to three major iterations.

We set the auto-threshold parameter to 3, which sets the global stopping threshold as 3 times the measured background noise. This value is automatically computed at the end of each major iteration by measuring the median absolute deviation (MAD) of the noise. In general, these parameters allow the deconvolution to run until the peak residual reaches the global stopping threshold. We set this threshold to 5 mJy, which is ${\sim}2.5$ times the expected noise level, to avoid cleaning too deeply into the noise. The ‘smallscalebias’ parameter is set at 0.85 to bias the multiscale deconvolution towards the smaller scales and aid cleaning of the expected small-scale features, and the gridded visibilities are weighted (without excluding any baselines) with a Briggs weighting scheme with a robust parameter of 0.9 to generate a PSF with higher sensitivity but still reasonable sidelobe levels.

The final deconvolved images are spatially smoothed such that the final restoring beam is $30^{\prime\prime} \times 30^{\prime\prime}$ as a compromise between sensitivity to extended structure and final angular resolution. This is comparable to the synthesised beam size expected from our robust weighting applied to the uv-coverage. The total time to image a typical channel ranged between 4 to 8 h, depending on the complexity of the emission present in the channel and number of major iterations. We find that the processing is generally memory limited, though processing of other pilot fields will increase optimisation and determine clear limitations.

Figure 3. The mean amplitude of XX and YY correlations averaged over 10 000 s as a function of binned channel for a central ASKAP beam taken from the baseline between ak04 and07. The ${\rm H\small I}$ emission-free channel range used to fit and subtract a first order polynomial model of the continuum is denoted by the shaded regions.

3.3.3. Filling in the missing short-spacings

The size of the primary reflectors on each antenna is an inherent physical limitation of how closely individual elements of an interferometer can be placed, thus filtering out the low spatial frequencies corresponding to the dish diameter (and below) and eliminating sensitivity to large angular scales. These missing baselines are commonly referred to as the missing short-spacings. Observations of the same area of the sky with a large single dish telescope with sufficient overlap in the uv-plane can be used to fill in the baseline samples below that of the minimum spacing between antennas (Stanimirović 2002) (i.e., the shortest baseline), including the total power at u = v = 0. The total flux and structure of diffuse emission on scales larger than that sampled by the shortest baseline is then recovered in the final map.

Several established techniques are available to combine interferometric and single dish observations either in the uv-plane after deconvolution (e.g., Stanimirović 2002; Cotton 2017), the image plane directly (Faridani et al. 2018) or through approximating the single dish data as artificial visibilities to be included in the image reconstruction process applied to interferometer data (Koda et al. 2011; Koda et al. 2019). More recently, Rau et al. (2019) have developed a generic joint reconstruction algorithm called SDINT that combines aspects from several of the aforementioned approaches. Given the extremely extended nature and complexity of our targeted emission, we utilise the most commonly used combination method called feathering, where the combined image is constructed by computing the weighted sum of the single dish and interferometer data in the uv-plane.

Table 3. WSClean parameters. See https://wsclean.readthedocs.io/_/downloads/en/latest/pdf/ for comprehensive documentation on these parameters.

We use the miriad task, IMMERGE, to feather a cube extracted from the Parkes Galactic All-Sky Survey (GASSf; McClure-Griffiths et al. 2009) with our deconvolved ASKAP image. This cube was centred on the SMC and regridded to be on the same spectral and spatial grid before combination. IMMERGE combines the two data sets by Fourier transforming both the deconvolved ASKAP-only and Parkes images. The Parkes data are then scaled by the ratio of the solid angle of the two beams,

(2) \begin{equation}\alpha = \frac{\Omega_{\mathrm{ASKAP}}}{\Omega_{\mathrm{Parkes}}},\end{equation}

where we find $\alpha$ to equal $9.77 \times 10^{-4}$ , assuming the Parkes beam is a Gaussian with a full-width at half-maximum (FWHM) of 16^′ at 1.4 GHz (Kalberla et al. 2010). This factor accounts for the difference in flux density of the two data sets based solely on the differences in resolution. The ASKAP uv data are then scaled by the factor

(3) \begin{equation}\beta = 1-\mathscr{F}[B_{\mathrm{Parkes}}],\end{equation}

where the second term represents the 2D Fourier transform ( $\mathscr{F}$ ) of the Parkes beam, $B_{\mathrm{Parkes}}$ . The two scaled uv data sets are then summed and Fourier transformed back to the image plane. The scaling depicted in Equation (3) ensures the effects of the poorly sampled low spatial frequencies in the ASKAP data are smoothly removed before the well-sampled low spatial frequencies provided by the Parkes data are added. The effective beam of the resultant image is the same as the original ASKAP restoring beam. The above scaling assumes the two data sets are perfectly flux calibrated; however, technical issues introduced by stray radiation and spillover make single dish data notoriously difficult to correctly flux calibrate.

Figure 4. The SPS profiles of the spectral channel centred on $151.3\,\mathrm{km\ s}^{-1}$ made from an ASKAP-only cube (solid blue), Parkes (solid orange), and several ASKAP+Parkes cubes where the sdfactor parameter (referred to as SDF in the text) in miriad’ns IMMERGE, which applies a scale factor to the single dish data before the combination in order to correct for potential flux offsets, is varied. The orange and blue dotted vertical lines denote the maximum recoverable scale of ASKAP based on the smallest baseline distance of 22 m and restoring beam size, respectively. Note that the profile for the SDF=1.0 values lies exactly on top of the ASKAP-only profile towards the smaller angular scales.

We take an interactive approach to determine the ideal flux calibration factor, henceforth referred to as the single dish factor (SDF), by setting the factor keyword in IMMERGE to be 0.5, 1.0, and 2.0 and then divide out this same factor in the resulting cube to conserve flux. This wide range of possible SDF values is useful to investigate the extremes of overestimating and underestimating the correction. We then compute the SPS of several spectral channels. The SPS is defined as

(4) \begin{equation}P\left(k\right) = \mathscr{F}\left(k_{\mathrm{RA}}, k_{\mathrm{Dec}}\right)\times\mathscr{F}^*\left(k_{\mathrm{RA}}, k_{\mathrm{Dec}}\right),\end{equation}

where k is the magnitude ( $k=\sqrt{k_{\mathrm{RA}}^2 + k_{\mathrm{Dec}}^2}$ ) of the wavenumber ( $k_{\mathrm{RA/Dec}} = \theta_{\mathrm{RA/Dec}}^{-1}$ , where $\theta_{\mathrm{RA,Dec}}$ is the angle on the sky along RA and Dec in rad), $\mathscr{F}\left(k_{\mathrm{RA}}, k_{\mathrm{Dec}}\right)$ is the 2D Fourier transform of the image under study, and the ${}^*$ symbol represents the complex conjugate. This is a two-point correlation function that measures how power (i.e., structure) is distributed across spatial scales. In practice, the distribution of power is measured by computing the 2D Fourier transforms of the integrated intensity images and measuring the median power in progressively larger annuli.

Figure 4 shows the SPS profiles of the spectral channel centred on $151.3\,\mathrm{km\ s}^{-1}$ in the Local Standard of Rest Kinematic (LSRK; Gordon 1976) reference frame produced from combinations that vary the SDF, as well as profiles from the deconvolved ASKAP-only cube and regridded Parkes image. Clearly, there is a dearth of power at large scales in the ASKAP-only image, while there is a large decrease in power detected in the Parkes image at scales below the larger single dish beam. All combinations recover the power at large scales since the feather procedure is designed to scale the flux density on those scales to that of the single dish. Interestingly, the variation in SDF has the most notable effect on the power at small scales. When the SDF is underestimated and less than one, the power increases at small scales because the scaling to conserve flux after the combination introduces a multiplicative factor that increases the flux of the already existing small-scale structure in the ASKAP-only image. When SDF is overestimated, the power decreases at small scales because the low resolution single dish data effectively washes out the observed small-scale structure.

Figure 5. A comparison between the ASKAP-only (left) and resultant combined ASKAP+Parkes image (right) for a single spectral channel centred on $151.3\,\mathrm{km\ s}^{-1}$ in the LSRK reference frame. The optimal SDF factor of 1.0 was applied during the combination.

We find that an SDF of 1.0 produces virtually an exact agreement with the ASKAP-only profile that extends from the maximum recoverable scale of ${\sim}33'$ down to the size of the 30^′′ restoring beam. Figure 5 shows the negative bowls in the ASKAP-only data, which are caused by the filtering of spatial frequencies corresponding to structures on the largest scales, are completely filled when applying a SDF of 1.0 during the combination. A SDF of 1.0 is reassuring and ensures we do not risk misrepresenting the observed structure by over or under-weighting the power at a given scale. Most importantly, the interesting structures at small scales are not washed-out. SDF values of 0.9 and 1.1 produce noticeable offsets in the SPS profiles, strengthening our conclusion. While we do not show an uncertainty envelope on the SPS profiles for the sake of clarity, Section 4.3 demonstrates that the typical uncertainties on individual power values are on the order of 5%.

3.4. Quality assessment

The spectral line cubes produced from GASKAP-HI observations also need to be validated, both to determine their suitability for science use and also to compare different processing techniques. We aim to validate the representation of the diffuse Hi of the Galactic Plane and the Magellanic System in the following ways:

• • that large-scale emission is present in velocity ranges where it is expected and not present where not expected, based on previous surveys;

• • that the noise level in the emission-free velocity range is close to theoretical, but not below, and is even across the cube; and

• • that spectra taken at the locations of selected targets, such as against bright sources, do not show any periodic features.

We validate the spectral line data at the ASKAP-only stage rather than after they have been combined with Parkes data. This allows us to compare the ASKAP-only data with data from GASS; McClure-Griffiths et al. 2009) in our tests for expected presence and absence of emission. It also provides a more rigorous test of noise levels in the ASKAP data.

To check for the presence or absence of emission, we define two $\Delta v \approx 40\, {\mathrm{km}\, \mathrm{s}^{-1}}$ regions where Hi line emission is either present or not present in the GASS cube. We extract a slab from the ASKAP cube covering each velocity region and use BANE (part of the Aegean suite; Hancock, Trott, & Hurley-Walker 2018), to produce a large-scale emission line image of each slab. The large-scale emission line images exclude Hi point sources and smooth small-scale fluctuations based on the synthesised beam size. We assess the maximum flux density level of these images against the thresholds in Table 4.

We assess the spectral line noise by taking the standard deviation of each pixel of the offline slab (as described above) along the spectral axis. We then compare the median value of this noise map to a fiducial value, the theoretical noise, with the thresholds listed in Table 4. The theoretical noise, at the beam centre and for natural weighting, is calculated as

(5) \begin{equation}\sigma_F = \frac{{\mathrm{SEFD}}}{\epsilon_{\mathrm{c}} \sqrt{n_{\mathrm{pol}} n_{\mathrm{ant}} (n_{\mathrm{ant}}-1) \Delta_{\mathrm{f}} t_{0}}}\end{equation}

where SEFD is the System Equivalent Flux Density, $\epsilon_{\mathrm{c}}$ is the correlator efficiency, $n_{\mathrm{pol}}$ is the number of polarisations, $n_{\mathrm{ant}}$ is the number of antennas, $\Delta_{\mathrm{f}}$ is the channel width and $t_{0}$ is the duration of the observation. ASKAP has a $\mathrm{SEFD} = 1\,800\,\mathrm{Jy}$ , $\mathrm{n}_{\mathrm{pol}} = 2$ , and $\mathrm{n}_{\mathrm{ant}} = 36$ (Hotan et al. 2021). We assume that $\epsilon_{\mathrm{c}}$ is included in the SEFD as the SEFD is measured through observations, but it should be ${\sim}1$ regardless (see Table 5). A typical GASKAP-HIobservation of $t_{0} = 20\,\mathrm{h}$ with $\Delta_{\mathrm{f}} = 5\,\mathrm{kHz}$ channels gives $\sigma_F = 1.89\,\mathrm{mJy\ beam}^{-1}$ .

Table 4. Spectral Line Cube Diagnostics Metrics. These are the thresholds for the different quality tests on GASKAP-HI spectral line cubes to assess them with a rating of good, uncertain or bad.

Table 5. Noise scaling parameters.

Lastly, we validate spectra extracted both from the locations of the 15 brightest sources (as identified in the continuum catalogue) and the centre of each beam. At each position we extract a single pixel spectrum and plot these spectra for visual inspection. We assess the spectra using a partial autocorrelation function (PACF) to identify any periodic features in each spectrum, requiring a 3-sigma result to classify a feature as periodic. We rate the dataset based on the number of spectra with periodic features as shown in Table 4.

We can use these metrics on the ASKAP-only image cube for the SMC field to assess its quality. Emission is present ( $21\,\mathrm{Jy\ beam}^{-1}\,\mathrm{km\ s}^{-1}$ ) where expected in the range $119 \leq \textit{v}_{\mathrm{LSRK}} \leq 165\,\mathrm{km\ s}^{-1}$ and absent ( $0.4\,\mathrm{Jy\, beam}^{-1}\,\mathrm{km\ s}^{-1}$ ) where not expected in the range $30 \leq \textit{v}_{\mathrm{LSRK}} \leq 73\,\mathrm{km\ s}^{-1}$ . Figure 6 shows a map of the spectral noise across the entire field. The spectral noise is $2.34\,\mathrm{mJy}$ per 5 kHz, or $1.24$ times theoretical, well within the good range. The slight loss in sensitivty can be attributed to the Briggs weighting applied to the visibilities (see Section 4.2.) The noise is generally smooth across the primary area of the field, however a grid pattern is apparent in the noise at low levels. This structure is likely a numerical artefact introduced by the Fourier Transform inversion of the visibilities when gridded to a finite grid to return to the image domain during major cleaning iterations. The effect of this structure is characterised in Section 4.3. The highest noise is at the edges of the field outside the FWHM of the formed primary beams, as expected. Finally, there are no periodic features in the test spectra. With all metrics in the good range we find that the SMC cube is suitable for science use.

Figure 6. Map of the spectral line noise for the ASKAP-only SMC cube from scheduling blocks 10941 and 10944. There is high noise at the edges of the cube and a notable cross-hatch low level noise pattern in the main area of the cube. However, it shows a generally low and consistent noise level across the main area of the cube, with a median noise level of $2.34\,\mathrm{mJy\ beam}^{-1}$ , 1.24 times theoretical.

4. Data products

4.1. Image cubes

The final product of our custom imaging pipeline, after masking the cube at the 10% level of the combined primary beam response, is a $3\,901\times3\,471\times222$ pixel image cube containing jointly deconvolved ASKAP data feathered with regridded Parkes data that covers a LSRK velocity range of $40.0{-}253.9\,\mathrm{km\ s}^{-1}$ at a resolution of $0.98\,\mathrm{km\ s}^{-1}$ . This cube is publicly accessible through a unique Digital Object Identifier (DOI).g We measure the rms noise using the source finding software, SoFiA 2; Serra et al. 2015; westmeier et al. 2021 , which fits a Gaussian to the negative half of the histogram of pixel values. This method is more robust against the presence of strong, extended emission or artefacts than measuring the standard deviation of emission-free channels. The width of this Gaussian is $1.1\,\mathrm{K}$ at the restoring beam size of $30^{\prime\prime}\times 30^{\prime\prime}$ , corresponding to a $5 \sigma\, {\rm H\small I}$ column density detection limit of $4.4\times10^{19}\,\mathrm{cm}^{-2}$ over a $20\,\mathrm{km\ s}^{-1}$ linewidth. This cube is the most sensitive and highest resolution image ( $30^{\prime\prime} \times 30^{\prime\prime}$ ) of the SMC in ${\rm H\small I}$ ever achieved.

4.2. Noise as a function of angular resolution

Two key science targets of GASKAP-HI are diffuse ${\rm H\small I}$ emission distributed on angular sizes larger than the synthesised beam and cold gas as revealed through sensitive ${\rm H\small I}$ absorption towards background continuum sources. It is critical to ensure the brightness temperature rms noise, $\sigma_{\mathrm{T}}$ , is minimised and behaves as predicted at the angular resolutions relevant for key GASKAP-HI science. Johnston et al. (2007) give the survey speed (SS) of ASKAP in terms of $\sigma_{\mathrm{T}}$ as

(6) \begin{equation} SS = FBn_{\mathrm{pol}}\left(\frac{\epsilon_{\mathrm{c}}\sigma_{\mathrm{T}}}{T_{\mathrm{sys}}}\right)^2\left(\frac{f}{\epsilon_{\mathrm{s}}}\right)^2,\end{equation}

where B is the spectral resolution, F is the FoV, $n_{\mathrm{pol}}$ is the number of polarisations (2), $\epsilon_{\mathrm{c}}$ is the correlator efficiency, $T_{\mathrm{sys}}$ is the system temperature in emission-free regions, and the filling factor of the array

(7) \begin{equation} f = \frac{A\eta N\Omega\epsilon_{\mathrm{s}}}{\lambda},\end{equation}

which describes the covering factor of the array in the aperture plane. Here, A is the area of a single antenna, $\eta$ is the aperture efficiency for a single antenna, N is the total number of antennas, and the solid angle of the synthesised beam $\Omega = 1.13(\theta_{\mathrm{s}})^2$ in steradians (sr), where the size of the synthesised beam is characterised by the FWHM along the major and minor axis: $\theta_{\mathrm{s}}=\sqrt{\theta_{\mathrm{maj}}\cdot\theta_{\mathrm{min}}}$

Figure 7. The measured noise (blue points) within a subregion of our GASKAP-HI cube after applying gradually larger tapers to the baselines in the uv-plane compared with the predictions shown in Figure 4 and Table 3 from Dickey et al. (2013) (orange lines and black crosses). The dotted black line denotes the updated theoretical noise using the measured $T_{\mathrm{sys}}/\eta$ scaled to an integration time of 20.9 h.

Finally, $\epsilon_{\mathrm{s}}$ is what is known as the synthesised aperture efficiency, constrained to be $\leq 1$ and relates to the weighting of the gridded visibilities in the uv-plane before inversion to the image plane. For example, applying a Briggs weighting scheme reduces $\epsilon_{\mathrm{s}}$ , since some baselines in the uv-plane become higher weighted than others. Keeping all other parameters in Equations (6) and (7) fixed, we anticipate the rms brightness temperature to scale as $\sigma_{\mathrm{T}} \propto f^{-1} \propto \theta_{\mathrm{s}}^{-2}$ , where $\epsilon_{\mathrm{s}}$ encodes the removal of baselines from the spatial smoothing. The telescope specifications are summarised in Table 5.

Figure 7 compares the measured noise as a function of angular resolution with the predictions from Dickey et al. (2013) (orange lines and crosses) at varying angular resolutions. The predicted noise profiles are determined by solving Equation (6) for $\sigma_{\mathrm{T}}$ for various dwell times ( $t_{\mathrm{int}}$ ) by setting $SS = F/t_{\mathrm{int}}$ . We measure the noise in our cube by imaging a subregion of $4\,\mathrm{beams} \times 3$ interleaves of 20 contiguous emission-free channels with natural weighting. We then apply a Gaussian taper to the baselines in the uv-plane before the inversion to the image plane to alter the final spatial resolution. The final beam size is a fit to the main lobe of the resulting PSF after tapering. We create 30 subcubes each with a different taper size that ranges from $12\,000\lambda$ to $400 \lambda$ equally spaced in log-space. The final noise is again taken to be the width of a Gaussian fit to the negative half of the histogram of pixel value as determined by SoFiA.

The measured noise is in excellent agreement with the theoretical predictions scaled to match the current telescope specifications (dotted black line). The profiles begin to deviate from the predicted $-2$ power law at larger resolutions due to the dependence of $\epsilon_{\mathrm{s}}$ at larger uv tapering. The agreement between the expected and measured scaling of the noise is confirmation that the data products of GASKAP-HI are indeed pioneering in terms of the trade-off between rms brightness temperature with angular resolution for the study of nearby ${\rm H\small I}$ .

4.3. Represented spatial frequencies

We further characterise the range of spatial frequencies represented in our final images after multi-scale joint deconvolution by employing the SPS as a data quality assessment tool. Specifically, we investigate whether a two-dimensional Gaussian function with a FWHM of $30^{\prime\prime}$ ( ${\sim}8\,\mathrm{pc}$ at the distance of the SMC) is a good approximation for the restoring beam across all available scales in the final image.

Similar to Section 3.3.3, we calculate the two-dimensional angular power spectrum for each spectral plane of our final image cube. Many studies that apply the SPS to column density or intensity images of various observational tracers in the Milky Way show the density distribution of the ISM as a function of scale is well described with a single power law (e.g., Martin et al. 2015; Miville-Deschênes et al. 2016; Pingel et al. 2018). We therefore fit the angular power spectrum profile $P\left(k\right)$ in each velocity channel with a model

(8) \begin{equation}P\left(k\right) = B\left(k\right)Ak^{-\gamma} + QP_{\mathrm{noise}}\left(k\right),\end{equation}

where $B\left(k\right)$ is the power spectrum of our ideal 30^′′ restoring Gaussian beam to account for instrumental systematics that affect the overall shape of the SPS, $A k^{-\gamma}$ is the ISM power law term, and Q is a multiplicative factor to the noise floor. The noise template, $P_{\mathrm{noise}}\left(k\right)$ , is estimated by taking the mean SPS profile for 40 emission-free channels. After fitting for A, $\gamma$ , and Q in each velocity channel, Equation (8) can be rearranged to solve for

(9) \begin{equation}\widetilde{B}\left(k, v_{\mathrm{LSRK}}\right) = \frac{P\left(k, v_{\mathrm{LSRK}}\right)-Q\left(v_{\mathrm{LSRK}}\right)P_{\mathrm{noise}}\left(k\right)}{A\left(v_{\mathrm{LSRK}}\right) k^{-\gamma\left(v_{\mathrm{LSRK}}\right)}},\end{equation}

where A, Q, and $P\left(k\right)$ now depend on $v_{\mathrm{LSRK}}$ , the observed radial velocity in the LSRK reference frame; $\widetilde{B}\left(k, v_{\mathrm{LSRK}}\right)$ is the global transfer function, which similar to optical systems specifies how the true sky brightness distribution is represented in our image cube. For example, the absence of structure in the transfer function indicates our data is the true sky brightness distribution convolved with our 30^′′ Gaussian restoring with added radiometernoise.

Figure 8. Global transfer function of the SMC image cube. The transfer function is useful for determining whether the observed sky brightness distribution is the true sky brightness distribution convolved with our 30^′′ Gaussian restoring beam with added radiometer noise. The top panel shows the entire range of angular scales, while the bottom panel shows the response between 200^′ to the vertical dashed line at 150^′ in the top panel to highlight scales with the most significant structure in the transfer function. The pixelated regions denote mostly emission-free channels.

In order to fully characterise the uncertainty on our model parameters, we follow similar fitting procedures from Koch et al. (2020), who used the PYMC3 package (Salvatier, Wiecki, & Fonnesbeck 2015), to run a Markov Chain Monte Carlo (MCMC) sampling to fit our final models. Due to the large range of angular scales spanned by the ${\rm H\small I}$ emission, we assign uniform priors for each of our model parameters and sample the amplitude A in $log_{10}$ space:

(10) \begin{equation}log_{10}A \sim \mathcal{U}(-20, 20)\end{equation}

(11) \begin{equation}\gamma \sim \mathcal{U}(0, 5)\end{equation}

(12) \begin{equation}Q \sim \mathcal{U}(0, 10)\end{equation}

We found that the fitting routine converges quickly with these choices of priors and the final fits are not affected by reasonable changes to the range of these priors. The MAD is used as an estimate of the uncertainty in each point $P_i$ in the azimuthally averaged profile

(13) \begin{equation} \sigma_{P(k)} =\,\mathrm{median}\left(\left|P_i -\,\mathrm{median}\left(P_{\mathrm{annuli}}(k)\right)\right|\right).\end{equation}

We compute the absolute deviation from the median power for each pixel within each annulus in the spatial frequency domain and take the final uncertainty to be the median of these deviations. The MAD is more robust against bias from small sample sizes at the smallest annuli and outlying values than the standard deviation. Each value of P(k) is treated as an independent sample drawn from a normal distribution with width equal to the MAD values in each azimuthal bin. These samples are first drawn in $log_{10}$ space then converted to linear scaling to avoid negative values for the power spectrum. Finally, we restrict our fits to a range of spatial frequencies between $k < k_{\mathrm{min}}/3$ —where $k_{\mathrm{min}}$ is one half the inverse of the largest map dimension—and $3k_{\mathrm{max}}$ , where $k_{\mathrm{max}}$ corresponds to the Gaussian standard deviation of our 30^′′ ideal beam. Fitting this range of spatial frequencies avoids the fluctuations from the large-scale power and bias from small-scale statistics at the smallest annuli, while also removing influence from correlated small-scale structure below the size of the restoring beam.

Figure 8 shows the measured transfer function for our SMC image cube. The absence of structure in the transfer function at small angular scales indicates that the density distribution of ${\rm H\small I}$ in the SMC is well described by a single power law down to the resolution limit of 10 pc, consistent with the conclusions from the power spectrum analysis of Stanimirović et al. (1999), which probed ${\rm H\small I}$ in a combined ATCA and Parkes image cube down to 30 pc. We discuss the interpretation of the ${\rm H\small I}$ turbulent properties in Section 5.2

While the small-scale structure is well represented by the ISM power law, there is some structure towards the largest angular scales ( ${\gtrsim}160^{\prime}$ ), highlighted in a zoomed version in the bottom panel of Figure 8. There are several possible explanations for this discrepancy including: the presence of large-scale deconvolution residuals, imperfect combination with the Parkes data, and large-scale structure in the noise. In regards to the presence of structure in the noise, the noise template SPS profile shown in the top panel of Figure 4 does begin to slightly deviate from the flat profile expected for purely random Gaussian distributed noise at these larger angular scales and also roughly correlates with the size scales of the low-level numerical artefacts from the Fourier transforms seen in Figure 6. Consequently, the structure in the transfer function at the largest angular scales likely manifests from a combination of sources of error, including bias from small-number statistics in the azimuthal binning of the smallest annuli and the known structure in the noise. It could also trace the azimuthally averaged baseline distribution. However, the power of these artifacts in the noise is many orders of magnitude below the power of the signal at these spatial frequencies and ultimately can be accounted for by restricting the range of scales (e.g., up to ${\sim}150'$ ) fit during a power spectrum analysis.

For example, the emission SPS profile from the same representative spectral channel in the top panel of Figure 9 is exceptionally well fit by our model. Inspection of other channels with extended emission show similar convergence. By restricting the range of angular scales and allowing a multiplicative factor for the noise, we ensure the presence of artefacts is taken into account in our SPS analysis. Martin et al. (2015) demonstrated through similar SPS that the noise is effectively the combination of receiver and sky noise (Boothroyd et al. 2011), $Q = 1 + T_{0, \mathrm{b}}/T_{\mathrm{sys}}$ , where $T_{0, \mathrm{b}}$ is a characteristic brightness temperature. We slightly alter the form of this relationship to $Q = 1 + 2 T_{0, \mathrm{b}}/T_{\mathrm{sys}}$ , where the factor of 2 accounts for the differences in the definitions of the linear XX and YY polarisations between Boothroyd et al. (2011) and the ASKAPsoft calibration pipeline, i.e. (XX+YY) versus (XX+YY)/2. In the bottom panel of Figure 9, we plot the fitted Q values versus the mean brightness temperature in each spectral channel. A fit of a simple linear function to channels with significant emission returns $Q = 1.014({\pm}0.003) + 0.0182({\pm}0.0004)T_{0, \mathrm{b}}$ , which gives $T_{\mathrm{sys}} = 54.97\pm0.02\,\mathrm{K}$ . This agrees exceptionally well with the expected $T_{\mathrm{sys}}$ at 1.4 GHz (Hotan et al. 2021).

Figure 9. Top: an example fit of our single component power law to a spatial power spectrum profile from a representative emission channel with the uncertainty derived from the MAD metric denoted by the shaded region. The noise template SPS is shown in orange. Bottom: Q, the multiplicative factor applied to the template SPS noise profile, as a function of the mean $T_{\mathrm{b}}$ of each spectral channel. A linear fit to channels with significant emission returns the expected $T_{\mathrm{sys}}$ of ASKAP at 1.4 GHz.

Figure 10. The peak intensity along each line-of-sight from our combined ASKAP and Parkes image cube of the SMC.

Figure 11. ASKAP+Parkes ${\rm H\small I}$ column density image of the SMC produced from Equation (14) integrated over the LSRK velocity range $60.5{-}235.4\,\mathrm{km\ s}^{-1}$ . The solid white circles denote the locations of the 36 formed beams from interleave A. In order to demonstrate the positional offset of the PAF footprint for each interleave position, we overlay dashed circles at the locations of beam 36 in interleave B (upper) and C (lower). This interleaving scheme ensures even sensitivity, and thus a consistent noise level across the entire instantaneous $25\,\mathrm{deg}^2$ FoV. The red rectangle region contains the high velocity cloud discussed in Section 5.3.

5. Gas morphology

5.1. General properties

The peak ${\rm H\small I}$ intensity image presented in Figure 10 demonstrates the full extent of the emission captured within a single ${\sim}25\,\mathrm{deg}^2$ ASKAP footprint. While not a true representation of the projected ${\rm H\small I}$ density distribution, as we simply select the peak brightness temperature value along the spectral profile at each spatial pixel, this spectacular image reveals that an incredible wealth of small-scale structure permeates the overall smooth diffuse ${\rm H\small I}$ component of the SMC. The extensive exogalactic population of ${\rm H\small I}$ features identified in McClure-Griffiths et al. (2018) are clearly visible throughout the northern regions of this map. We further highlight the projected density distribution by computing the column density image shown in Figure 11, assuming the ${\rm H\small I}$ is optically thin:

(14) \begin{equation} N(x,y)_{\mathrm{HI}} = 1.82\times10^{18} \int T_{\mathrm{b}}\left(x,y, v\right) dv \;{\mathrm{cm}^{-2}},\end{equation}

where the brightness temperature, $T_{\mathrm{b}}\!\left(x,y\right)$ , at each spatial pixel is integrated over the spectral range $235.4{-}60.5\,\mathrm{km\ s}^{-1}$ . The boundaries of the ASKAP primary beams are overlaid to demonstrate that the size scales of emission range from a single synthesised beam element up to the extent of multiple formed primary beams. Strikingly, there are still some indications of complex filamentary structure in the column density image, highlighting that our increased angular resolution resolves a wealth of new cold narrow features even when integrating across the entire spectral range.

The single channel shown in the top two panels of Figure 12 highlights the improvement in angular resolution in our pilot data over the previous generation image made with ATCA and Parkes. It could be argued that intermediate scale features appear better recovered in the previous generation ATCA+Parkes image. This is a combined effect of including emission over a slightly larger spectral range ( $1.65\,\mathrm{km\ s}^{-1}$ ), spatial smoothing to a lower angular resolution, and the difference in the uv-weighting applied during imaging. Staveley-Smith et al. (1997) applied a robustness parameter of 0.0, giving slightly less weight to sparsely sampled (i.e., long) baselines, bringing out structure on these intermediate scales.

Figure 12. A visual comparison between a single spectral channel centred on $151.3\,\mathrm{km\ s}^{-1}$ in the LSRK reference frame from the previous generation ATCA+Parkes of the SMC in ${\rm H\small I}$ (top; Stanimirović et al. 1999), the new GASKAP-HI (ASKAP+Parkes) ${\rm H\small I}$ image (middle), and the $250\,\unicode{x03BC}\mathrm{m }$ dust image from the Herschel space observatory (bottom; Gordon et al. 2014). The white rectangle outlines the Bar region of the SMC and the white arrow points to resolved ${\rm H\small I}$ plumes in the Wing region. The extent and detail of the small-scale ${\rm H\small I}$ features in this individual spectral channel of the GASKAP-HI image are similar to those seen in the dust continuum emission.

The abundance of small-scale ${\rm H\small I}$ features in the Bar of the SMC, extending from the north to the south-west and outlined by the white rectangle, trace the stellar feedback effects from the strong ongoing star formation in this region as seen in $\mathrm{H}\alpha$ (Winkler et al. 2015) and other star burst regions such as the N66/NGC 346 complex (Heydari-Malayeri & Selier 2010). The detail in small-scale features in a single spectral channel from our new GASKAP-HI data approaches those traced by the far-infrared $250\,\unicode{x03BC}\mathrm{m}$ dust image produced as part of the HERschel Inventory of The Agents of Galaxy Evolution (HERITAGE) project (Meixner et al. 2013). A comparison of the properties of the ${\rm H\small I}$ from a single spectral channel with dust emission is not straightforward given the continuum nature of the re-radiated emission arising from the cool dust. Nevertheless, assuming the gas and dust are well mixed due to the correlation between colour excess and total hydrogen column density (Rachford et al. 2009), future studies will utilise the improved angular resolution to provide insights into spatial variation of the dust-to-gas ratio across the entirety ofthe SMC.

Towards the so-called Wing of the SMC, which extends towards the LMC, the superior spectral and spatial resolution of these data reveal a complex network of discrete linear plumes (see the arrow in the middle panel of Figure 12) with projected lengths of ${\sim}0.5\,\mathrm{kpc}$ . A similar plume stemming just east of the northern portion of the Bar extends almost 1 kpc. The spectral resolution of our new image cube ensures these plumes and other interesting features, such as the supershell centred on $\mathrm{RA} = 00^{\mathrm{h}}50^{\mathrm{m}}$ , $\mathrm{Dec}=-73^{\mathrm{d}}10^{\mathrm{m}}$ , are fully resolved over 5–20 channels. Section 5.3 demonstrates how our improved insight into the kinematic properties of these discrete features helps with constraining important physical properties of the SMC.

The three-colour images of selected spectral channel maps in Figure 13, which stack emission from adjacent channels, highlight the complexity and narrow spectral properties of the newly revealed small-scale structure. At more redshifted velocities relative to the systemic velocity of $148\,\mathrm{km\ s}^{-1}$ LSRK (Di Teodoro et al. 2019a), the ${\rm H\small I}$ morphology is relatively smooth at large scales but gradually transitions to a filamentary form on arcminute scales as the velocity decreases. The ${\rm H\small I}$ distribution becomes increasingly intricate near the systemic velocity, taking the form of plumes, filaments, knots, and shells. The discrete filamentary features are analogous to the slender ${\rm H\small I}$ fibres shown by Clark et al. (2014) to be aligned with the interstellar magnetic field of the Milky Way, while the small-scale knots and arcs trace structure on 10 pc scales, similar to the typical size of star-forming regions. The multi-colour nature of these features demonstrates a sharp velocity gradient over just several adjacent spectral channels with linewidths comparable to $\mathrm{CO}(2 \rightarrow 1$ ) emission observed in individual molecular clouds in other nearby external galaxies as reported by the PHANGS-ALMA survey (Leroy et al. 2021). In our own Galaxy, the atomic-to-molecular ( ${\rm H\small I}$ -to- $\mathrm{H}_2$ ) phase transition in the dense regions of Giant Molecular Clouds (GMCs), arising from the combined effect of dust attenuation of far-ultraviolet radiation and $\mathrm{H}_2$ self-shielding (Bialy et al. 2015 and references therein), traces sites of star formation. Glover & Clark (2012) argue that the transition from the WNM to CNM phase is an important step towards star formation. Studying the gas transition from ${\rm H\small I}$ to $\mathrm{H}_2$ -dominated regions that fuel star formation within the Bar and comparing with observational studies of Milky Way GMCs (e.g., Lee et al. 2012; Burkhart et al. 2015; Wang et al. 2020) will be an invaluable step towards a complete theory of star formation.

Figure 13. Selected channel maps of the combined ASKAP+Parkes ${\rm H\small I}$ image cube. Red, green, and blue are assigned sequentially to adjacent velocity channels, each displayed with an arcsin colourmap scaling over a brightness temperature range of $T_{\mathrm{b}} = 1$ and $70\,\mathrm{K}$ . The colourscale of all panels is the same. The LSRK velocity range is listed in each panel.

Towards blueshifted velocities (e.g., $131\,\mathrm{km\ s}^{-1}$ ), a prominent population of compact, spatially anomalous ${\rm H\small I}$ features surround the SMC. Considering the general head-tail morphology of associated molecular gas by Di Teodoro et al. (2019b), a majority of these anomalous features likely originate from strong star-formation-driven outflows (McClure-Griffiths et al. 2018), though the uncertain geometry of the SMC conspires against a definitive measure of the escaping gas. Nevertheless, these novel observations will inform the comparison between kinematic and geometric properties of infall/outflow models.

5.2. Hi turbulent properties

The extension of a single component power law down to physical scales three times smaller than previously measured demonstrates processes such as star formation, which injects energy on several hundred pc scales (Mac Low & Klessen 2004), do not produce any special feature or deviation in the large-scale SPS of the SMC. Either stellar feedback is not the dominant source of energy injection at these scales and large-scale driving is more dominant (at least within the SMC; Szotkowski et al. 2019), or energy is being injected over a wide range of scales without preference at any particular scale (e.g., Norman & Ferrara 1996). This is also congruous with the notion that the power law form of the power spectrum arises from the fractal geometry of the ISM; that is, the ISM displays hierarchical structuring over the entire range of spatial scales (Stanimirović et al. 1999; Elmegreen & Falgarone 1996). Interestingly, the average spectral index measured for all slopes over velocity range of $90{-}190\,\mathrm{km\ s}^{-1}$ is $\gamma = -2.85\pm 0.07$ , which is slightly shallower than $\gamma = -3.04\pm 0.04$ measured by Stanimirović et al. (1999). Shallower slopes across the entire velocity range indicate more power on smaller scales, possibly from more abundant cold ${\rm H\small I}$ that is resolved by our higher angular resolution GASKAP-HI observations. On the other hand, the differences in the average slopes could also be driven from a statistical bias (e.g., small number statistics for the large-scale power). Following the application of the discussion from Falconer (1997) and Stanimirović et al. (1999), we exploit the Fourier Transform relationship between the autocorrelation function and power spectrum to connect $\gamma$ to the fractal function with the relation

(15) \begin{equation} -\gamma-2=4-2s,\end{equation}

where s is the box-counting dimension. This value quantifies the fractal nature of a subset of three-dimensional Euclidean space, $\mathrm{R}^3$ , and provides an upper limit for the fractal dimension (also known as the Hausdorff dimension). For example, a perfectly smooth square has a fractal dimension of 2, while an object with a non-integer fractal dimension indicates self-similarity over a range of scales and thus the presence of hierarchical structure. We determine the s of the combined ASKAP and Parkes images of ${\rm H\small I}$ emission to be $1.58\pm 0.04$ , in agreement with the fractal analysis in Stanimirović et al. (1999). The fractal dimension observed in molecular clouds is $1.36\pm 0.02$ from scales ranging from 100 to 0.02 pc (Falgarone et al. 1991). This statistically significant difference indicates separate mechanisms drive the dynamics in the WNM and molecular clouds.

Analogously, a clustering analysis applied to the upper main-sequence stars in the SMC shown by Sun et al. (2018) shows young stellar structures, much like the gas, are also organised in a hierarchical fashion. The consistent description for both the gas and stars indicates similar mechanisms, such as turbulence and/or hierarchical gravitational fragmentation of the denser gas that forms stars (Vázquez-Semadeni et al. 2019), influence the properties of the ${\rm H\small I}$ across spatial scales extending from kpc down to the scales of individual molecular clouds on sub-parsec scales. We note, however, that a comprehensive analysis of the fractal structure, especially at individual spectral channels, should account for the effects of self-absorption and difficulty in disentangling the relative contributions of the CNM and WNM to the overall ${\rm H\small I}$ emission.

5.3. Properties of an anomalous Hi feature

The test observations presented in McClure-Griffiths et al. (2018) revealed a considerable population of spatially and kinematically anomalous ${\rm H\small I}$ features that extend upwards of $2\,\mathrm{kpc}$ from the main body of the SMC. The increase in spectral resolution and overall sensitivity of our pilot survey observations facilitate more robust constraints of the physical properties of these Magellanic HVCs—such as their mass, kinematics, and morphology. These in turn trace fundamental properties of the circumgalactic environment around the SMC and Milky Way halo like the total amount of outflowing gas and halo pressure by relating head-tail morphologies with ram pressure and/or instabilities.

We illustrate the power of our new pilot data to trace anomalous HVCs in the Magellanic System by focusing on an outflowing HVC centred on $\mathrm{RA} = 01^{\mathrm{h}}08^{\mathrm{m}}45^{\mathrm{s}}$ , $\mathrm{Dec} = -71^{\mathrm{d}}20^{\mathrm{m}}55^{\mathrm{s}}$ that entrains several compact ${}^{12}\mathrm{CO}(2\rightarrow1$ ) clouds (Di Teodoro et al. 2019b). The left and right panels of Figure 14 compare the integrated intensity images of this HVC from the pilot and test observations, respectively. We utilise the SEARCH task in BBarolo, a software package for fitting 3D tilted-ring models to emission line cubes (Di Teodoro & Fraternali 2015), in order to exclude noisy pixels from the integration of the signal. We set the SNR threshold for genuine signal to any pixel with SNR $> 5$ over the span of two or more adjacent channels, applied within the LSRK velocity range of $95{-}119\,\mathrm{km\ s}^{-1}$ . The remaining task parameters use default values. The images from the pilot and test observations show similar large-scale morphology. However, the small-scale structures in the pilot data displays more pronounced density gradients in individual clumps. This is due to the twofold increase in sensitivity, fourfold increase in spectral resolution (facilitating more narrow line features to be resolved), and more complete sampling of higher spatial frequencies from the additional antennas at longer baselines. For example, let us focus on the northeast condensation with a peak at $\mathrm{RA} = 01^{\mathrm{h}}09^{\mathrm{m}}35^{\mathrm{s}}$ , $\mathrm{Dec} = -71^{\mathrm{d}}17^{\mathrm{m}}15^{\mathrm{s}}$ . The mean column densities are roughly similar at $1.3\times10^{20}$ and $8.2\times10^{19}\,\mathrm{cm}^{-2}$ in the pilot and test images, respectively, giving a volume density $n \sim 1.6\,\mathrm{cm}^{-3}$ and $1.0\,\mathrm{cm}^{-3}$ assuming spherical symmetry, an angular extent of 90^′′, and a distance of 60 kpc.

Figure 14. Combined ASKAP+Parkes Hi masked integrated intensity images of an HVC centred on $\mathrm{RA} = 01^{\mathrm{h}}08^{\mathrm{m}}45^{\mathrm{s}}$ , $\mathrm{Dec} = -71^{\mathrm{d}}20^{\mathrm{m}}55^{\mathrm{s}}$ produced by the source finding algorithm in BBarolo. The projected physical extent of this region is 300 pc at a distance of 60 kpc. The image made from our improved pilot survey is on the left and the image from the test observations (McClure-Griffiths et al. 2018) is on the right. In both cases, the detected emission was integrated over $90 \leq \textit{v}_{\mathrm{LSRK}} \leq 118\,\mathrm{km\ s}^{-1}$ . The contour levels increase in steps of $29.5\,\mathrm{K\ km\ s}^{-1}$ and range from $5\,\mathrm{K\ km\ s}^{-1}$ to $300\,\mathrm{km\ s}^{-1}$ , inclusive. The red circle with a radius corresponding to 13 pc denotes the approximate size of a shell-like feature identified in position-velocity space. The orange contours in the left panel show the CO integrated intensities at 0.02 and $0.04\,\mathrm{K\ km\ s}^{-1}$ from Di Teodoro et al. (2019b). The size of the restoring beam for the ASKAP data sets is denoted by the circle in the lower left corner.

The full potential for constraining the physical properties of these HVCs is realised through the increased spectral resolution. For example, we break up a $2.5'\times2.5'$ region centred on the peak intensity of this condensation to the northeast for the pilot and test data into a $4\times 4$ grid to independently sample the emission across 16 restoring beam elements and compute the integrated spectrum for each sample. We then shift each integrated spectrum to the same peak velocity and average the profile. Fitting a single Gaussian function to these averaged profiles returns $v_{\mathrm{FWHM}}$ of $6.9\pm 0.5$ and $13 \pm 2\,\mathrm{km\ s}^{-1}$ for the pilot and test data, respectively; the differences arising from the differences in spectral resolution. These $v_{\mathrm{FWHM}}$ provides an upper limit for the kinetic temperature of the gas, $T_{\mathrm{max}, \mathrm{k}} = 21.866 \times v_{\mathrm{FWHM}}^2$ of $1\,000\pm200$ and $4\,000\pm1\,000\, \mathrm{K}$ . We can use the $T_{\mathrm{max}, \mathrm{k}}$ and n to obtain an estimate of the pressure of this condensation through the approximation of an ideal gas, $P/k = nT_{\mathrm{max}, \mathrm{k}}$ , where P is the pressure and k is the Boltzmann constant. The resulting upper limit pressures for the pilot and test data are respectively $P/k = 1\,600\pm300\,\mathrm{K\ cm}^{-3}$ and $P/k = 4\,000\pm1\,000\,\mathrm{K\ cm}^{-3}$ . The respective decrease of nearly 60% and 70% between the upper limit pressures and final uncertainties demonstrates that the superior spectral and angular resolution of the pilot data is crucial for constraining important physical properties. This also demonstrates how increased angular and spectral resolutions work together to improve estimates of intrinsic line widths.

Figure 15. Combined ASKAP+Parkes Hi channel maps of an HVC centred on $\mathrm{RA} = 01^{\mathrm{h}}08^{\mathrm{m}}45^{\mathrm{s}}$ , $\mathrm{Dec} = -71^{\mathrm{d}}20^{\mathrm{m}}55^{\mathrm{s}}$ . The black contours denote emission from the pilot cube, while the white contours trace emission from the lower spectral resolution test cube. Both sets of contours range between 5 and 30 K and increase in increments of 6.25 K. The increased spectral resolution of the pilot cube resolves several notable features including a Hi hole between 110.3 and $105.4\,\mathrm{km\ s}^{-1}$ . The approximate projected size of the shell at these velocities is denoted by a black circle, assuming uniform expansion and spherical symmetry. The restoring beam of the pilot data is shown in the lower left corner of each panel.

The CO integrated intensity contours from the observations of Di Teodoro et al. (2019b) overlap well both spatially and spectrally (see their Figure 3) with the southern edge of the most pronounced ${\rm H\small I}$ density gradient in this cloud centred on $\mathrm{RA} = 01^{\mathrm{h}}08^{\mathrm{m}}33^{\mathrm{s}}$ , $\mathrm{Dec} = -71^{\mathrm{d}}24^{\mathrm{m}}21^{\mathrm{s}}$ . The emission distribution in individual spectral channels in Figure 15 reveal a shell-like structure appearing first at a velocity of $110.3\,\mathrm{km\ s}^{-1}$ before completely breaking up at $103.4\,\mathrm{km\ s}^{-1}$ along the northern edge of this density gradient.

Figure 16. Combined ASKAP+Parkes Hi position-velocity distribution for three slices along RA. The top panel represents a slice 30^′′ above the centre of the circle in the left panel of Figure 14 that intercepts the Hi shell at an angle of $50^{\circ}$ from a horizontal slice through the centre located at $\mathrm{RA} = 01^{\mathrm{h}}08^{\mathrm{m}}38^{\mathrm{s}}$ , $\mathrm{Dec} = -71^{\mathrm{d}}23^{\mathrm{m}}00^{\mathrm{s}}$ . The middle panel represents a slice directly through the centre, and bottom panel shows a slice 30^′′ below the centre that intercepts the shell at a ${-}50^{\circ}$ angle below a horizontal slice through the centre. Assuming the geometry outlined for a uniformly expanding shell with $v_{\mathrm{exp}} = 2.0\,\mathrm{km\ s}^{-1}$ (much smaller than the typical thermal line widths expected for the WNM), we show the predicted position-velocity structure as a red ellipse. The emission contours are at levels of 3, 5, 10, and 20 K.

The distribution of emission in position-velocity space of several RA cuts in Figure 16 hints at a characteristic ellipse structure expected from an expanding shell, indicating spherical symmetry in 3D space. We determine the centre of the apparent shell-like feature to be $\mathrm{RA} = 01^{\mathrm{h}}08^{\mathrm{m}}38^{\mathrm{s}}$ , $\mathrm{Dec} = -71^{\mathrm{d}}23^{\mathrm{m}}00^{\mathrm{s}}$ and measure $v_{\mathrm{min}} = 106 \pm 1\,\mathrm{km\ s}^{-1}$ and $v_{\mathrm{max}} = 110 \pm 1\,\mathrm{km\ s}^{-1}$ and thus a systemic velocity $v_{\mathrm{sys}} = 108\pm 1\,\mathrm{km\ s}^{-1}$ . Assuming a uniformly expanding shell, the projected sizes in position-velocity space are respectively

(16) \begin{equation} v_{\mathrm{proj}} = v_{\mathrm{exp}}\cos{\theta_{\mathrm{cut}}}+v_{\mathrm{sys}}\end{equation}

and

(17) \begin{equation} R_{\mathrm{pos}} = R_{\mathrm{ang}}\cos{\theta_{\mathrm{cut}}},\end{equation}

where $\theta_{\mathrm{cut}}$ represents the angle at which a horizontal cut is above and below another horizontal cut through the centre of the feature in position-position space and $R_{\mathrm{ang}}$ is the angular projected radius at the centre and is determined to be 13 pc at the distance of adopted distance of 60 kpc. The line-of-sight (LOS) geometry is estimated using

(18) \begin{equation} \phi = \cos^{-1}\left(\frac{v_{\mathrm{LSRK}} - v_{\mathrm{sys}}}{v_{\mathrm{exp}}}\right)\end{equation}

(19) \begin{equation} R_{\mathrm{proj}}=R_{\mathrm{ang}}\sin{\phi},\end{equation}

Figure 17. Combined ASKAP+Parkes Hi maps of the intensity-weighted velocity field (left) and velocity dispersion (right) of an HVC centred on $\mathrm{RA} = 01^{\mathrm{h}}08^{\mathrm{m}}45^{\mathrm{s}}$ , $\mathrm{Dec} = -71^{\mathrm{d}}20^{\mathrm{m}}55^{\mathrm{s}}$ produced by the source finding algorithm in BBarolo. The velocity field contours range inclusively between the levels of 92.5 and $115\,\mathrm{km\ s}^{-1}$ in steps of $2.5\,\mathrm{km\ s}^{-1}$ . The contours of velocity dispersion begin at a level of $1\,\mathrm{km\ s}^{-1}$ and increase up to $4\,\mathrm{km\ s}^{-1}$ in steps of $0.75\,\mathrm{km\ s}^{-1}$ . The gradient in the velocity field across the entirety of the cloud and complexity of the velocity dispersion contours—especially towards the northwest spur—are qualitatively consistent with this HVC being a wind-swept cloud from a Galactic wind emanating from the SMC. The arrow denotes the direction towards the nearest prominent star formation region, NGC 371/395. The size of the restoring beam is shown in the lower left corner of both panels.

where $v_{\mathrm{LSRK}}$ is again the observed radial velocity in the LSRK reference frame. Using these $v_{\mathrm{min}}$ and $v_{\mathrm{max}}$ and measuring the velocity of the emission at the location of a horizontal cut that intercepts the shell at an angles of $\theta_{\mathrm{cut}}=50^{\circ}$ and $\theta_{\mathrm{cut}}=-50^{\circ}$ above and below the centre, we estimate the expansion velocity $v_{\mathrm{exp}}$ to be $2\pm 1\,\mathrm{km\ s}^{-1}$ . Using Equations (16)–(19), we show the projected size in the relevant channel maps and in position-velocity space.

A $v_{\mathrm{exp}}$ of $2 \pm 1\,\mathrm{km\ s}^{-1}$ is roughly five times smaller than typical thermal line widths for the WNM, making it difficult to discern likely origin scenarios based on input energy estimates for possible progenitors. A search for bright (G $< 20$ ), blue (colour $-0.5 <\,\mathrm{BP}-\mathrm{RP} < 0$ ) stars with parallaxes corresponding to distances larger than 55 kpc within the Gaia Early Data Release 3 catalogue (Gaia Collaboration et al. 2016; Gaia Collaboration et al. 2021) shows dozens of stars in the field. However, we cannot point to an isolated population as the potential progenitor of this expanding shell without further filtering based on radial velocity measurements, which are currently unavailable. Moreover, the measured mass within the aperture encompassing the feature is ${\sim}2\,240\,\mathrm{M}_{\odot}$ , implying a corresponding thermal energy of $8.9\times10^{47}$ erg. This is much lower than the input from a single core-collapse supernova explosion (Hartmann 1999). The absence of an isolated stellar association, small $v_{\mathrm{exp}}$ relative to the thermal linewidth of the WNM, and correspondingly low thermal energy point to a purely turbulent origin. Daigle et al. (2007) note the presence of structures at similar scales as this feature in simulated ${\rm H\small I}$ image cubes with pure turbulent motions. Regardless of the true origin, this investigation demonstrates that our GASKAP-HI data provide a novel view into the kinematic properties of the diverse population of exogalactic clouds surrounding the SMC.

The presence of galactic-scale external winds can explain the global observed morphological and kinematic properties of this cloud. Galactic winds shape the structures of the ISM, altering the kinematic and morphological properties of the gas. Hot winds interacting with cold clumps of dense gas can excite molecular species (Pon, Johnstone, & Kaufman 2012; Girichidis et al. 2021), transfer mass, momentum, and energy (Vijayan et al. 2020; Schneider et al. 2020), and strip the outer surface layer, forming elongated, filamentary structures that trail the main body of the cloud (e.g., Banda-Barragán et al. 2016; Gronke et al. 2021). Di Teodoro et al. (2019b) suggests a scenario in which this particular complex is launched from NGC 371/395, the nearest prominent star formation region located south-west. The intensity-weighted velocity field shown in the left panel of Figure 17 shows a clear gradient across the entirety of the cloud. Such velocity gradients in discrete clouds are often interpreted as rotation. However, the closed isovelocity contours in the western side of the cloud that overlap with the expanding bubble, absence of a clear major axis, and disturbed structure in the velocity dispersion map instead indicate these motions are more likely due to shear. The population of discrete cloudlets to the northwest in the individual spectral channels and moment maps is also consistent with the expectation of the formation of head-tail structures from the stripping of outer gas layers seen in simulations of wind–cloud interactions and observations (Phillips, Ramos-Larios, & Perez-Grana 2009; Grønnow et al. 2017; Grønnow, Tepper-Garc a, & Bland-Hawthorn 2018; Banda-Barragán et al. 2019; Banda-Barragán et al. 2021).

The integrated spectrum shown in Figure 18, measured over the field of view shown in Figures 14 and 15, further demonstrates the advantage of the increased spectral resolution in deriving physical properties of HVCs around the SMC. The pilot flux profile shows the multi-component kinematic properties of this cloud, which is only hinted at in the coarser spectrum. Inspection of the individual channel maps for velocities $105{-}97\,\mathrm{km\ s}^{-1}$ shows the flux offset between the two profiles occurs because a significant amount of emission happens to sit at the boundary of a coarser spectral channel in the test data, thus artificially extending the velocity range of the emission distribution. Again assuming a distance of 60 kpc (McClure-Griffiths et al. 2018), we achieve a more reliable measure of the ${\rm H\small I}$ mass to be $15\,100\pm300\, \mathrm{M}_{\odot}$ , as opposed to an overestimated value of $19\,600\pm100$ $\mathrm{M}_{\odot}$ . This discrepancy is likely related to improvements made in the calibration and processing pipeline between November 2017 and December 2019.

Figure 18. Combined ASKAP+Parkes Hi integrated spectrum of the HVC measured over the area shown in Figures 14 and 15.

These data enable us to derive a diverse range of physical properties for this single HVC. Future detailed cataloging and study of the HVC population around the SMC will inform mass outflow/inflow models to help determine whether such features are outflowing or accreting, aid in constraining the overall dynamics of the SMC, and trace the pressure equilibrium between cold and warm ${\rm H\small I}$ phases through Gaussian decomposition algorithms, such as the Regularized Optimization for Hyper-Spectral Analysis (ROHSA) package (Marchal et al. 2019) or GaussPy+ (Lindner et al. 2015; Riener et al. 2019).

6. Summary

Utilising a combination of ASKAP and Parkes GASS ${\rm H\small I}$ data to obtain a cube covering all angular scales, we presented the most sensitive and detailed ${\rm H\small I}$ images of the SMC ever achieved. These images were produced by our custom imaging pipeline, which is designed to optimise a joint deconvolution approach in a cluster environment in order to recover an accurate representation of the extended sky brightness distribution. In addition to a full description of the calibration procedures, imaging pipeline, and quality assessments, we demonstrated that:

• • the SPS is a powerful tool in determining the correct scaling factor applied to the total power observations when combining interferometric and single dish data. For example, an overestimated scaling factor results in a decrease in power at small scales because the combination effectively washes out the small-scale structure present in an interferometer-only image.

• • The final SMC image cube possesses a rms noise level of $1.1\, \mathrm{K}$ per $0.98\,\mathrm{km\ s}^{-1}$ spectral channel and an angular resolution of 30^′′ (8 pc at the distance of the SMC). The noise properties of the final data are generally consistent across the FoV until the sensitivity begins to drop at the edge of the PAF footprint. A low-level grid pattern is apparent, likely stemming from numerical artifacts from the Fourier Transforms that invert the gridded visibilities to the image domain during major iterations. The noise scales as predicted by Dickey et al. (2013) with the resolution of the restoring beam, indicating these GASKAP-HI data are a breakthrough in terms of the compromise between rms brightness temperature and angular resolution.

• • Our power spectrum analysis shows the distribution of ${\rm H\small I}$ emission in the SMC is well described by a single component power law, consistent with the interpretation that the ISM is fractal in nature and that stellar feedback is probably not highlighted at any particular scale. Stochastic processes, such as turbulence, therefore influence the structure of the ISM from scales ranging from kpc down to sub-parsecs.

• • A comparison of the basic properties of an outflowing ${\rm H\small I}$ complex derived from test data with ${\sim}4 \times$ worse spectral resolution relative to these new pilot observations demonstrates the power of our pilot data to constrain important ISM properties. The increase in spectral resolution allows us to resolve an expanding ${\rm H\small I}$ hole. The low thermal energy and expansion velocity suggests a purely turbulent origin for this feature.

These new ${\rm H\small I}$ images of the SMC show spectacular details on scales from 5 kpc extending down to 10 pc and reveal previously unresolved narrow line features. A separate paper will present the abundance of absorption detections in this same field. These pilot data facilitate the study of how ${\rm H\small I}$ contributes to the conditions necessary for star formation, probe the hierarchical turbulent spectrum down to physical scales three times smaller than previously measured, and truly characterise the multi-phase make up of the ISM. GASKAP-HI provides an unequivocal view of the nearby ${\rm H\small I}$ associated with the detailed structures observed in the dust and molecular gas. Additionally, the lessons learned from imaging extremely extended emission for the GASKAP-HI survey will inform the data processing of future ultra wide-field surveys to be undertaken by the Square Kilometer Array that will face identical imaging challenges. The increase in brightness temperature sensitivity, coupled with the unparalleled angular and spectral resolution, ensures the GASKAP-HI survey will provide touchstone data products for the study of ${\rm H\small I}$ in the Milky Way and Magellanic System for the next decade and beyond.