3.1. ASKAPsoft pipeline

ASKAP observations are calibrated and imaged using a high-performance processing pipeline (Whiting Reference Whiting2020; Hotan et al. Reference Hotan2021) 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. Reference Guzman2019). 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.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.Footnote b It summarises and visualises outputs produced by the ASKAPsoft (Guzman et al. Reference Guzman2019) 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. Reference For2019; Koribalski et al. Reference Koribalski2020), 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. Reference Huynh, Dempsey, Whiting and Ophel2020). 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,Footnote 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 pipelineFootnote d that utilises a combination of Common Astronomy Software Application (CASA; McMullin et al. Reference McMullin, Waters, Schiebel, Young, Golap, Shaw, Hill and Bell2007) tasks to pre-process the measurement sets, WSClean Footnote 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 Reference Sault, Teuben and Wright1995) 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. Reference Heywood2019). We input the ASKAP primary beams measured from holography at 1.4 GHz (Hotan et al. Reference Hotan2021) 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. Reference Cornwell, Golap and Bhatnagar2005), 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 Reference Cornwell2008) 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ć Reference 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ć Reference Stanimirović2002; Cotton Reference Cotton2017), the image plane directly (Faridani et al. Reference Faridani, Bigiel, Flöer, Kerp and Stanimirović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. Reference Koda2011; Koda et al. Reference Koda, Teuben, Sawada, Plunkett and Fomalont2019). More recently, Rau et al. (Reference Rau, Naik and Braun2019) 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 (GASSFootnote f; McClure-Griffiths et al. Reference McClure-Griffiths2009) 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. Reference Kalberla2010). 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 Reference Gordon1976) 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. Reference McClure-Griffiths2009) 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 Reference Hancock, Trott and Hurley-Walker2018), 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. Reference Hotan2021). 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.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).Footnote g We measure the rms noise using the source finding software, SoFiA 2; Serra et al. Reference Serra2015; westmeier et al. Reference Westmeier2021 , 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. (Reference Johnston2007) 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. (Reference Dickey2013) (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. (Reference Dickey2013) (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. Reference Martin, Blagrave, Lockman, Pinheiro Gonçalves, Boothroyd, Joncas, Miville-Deschênes and Stephan2015; Miville-Deschênes et al. Reference Miville-Deschênes, Duc, Marleau, Cuillandre, Didelon, Gwyn and Karabal2016; Pingel et al. Reference Pingel, Lee, Burkhart and Stanimirović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. (Reference Koch, Chiang, Utomo, Chastenet, Leroy, Rosolowsky and Sandstrom2020), who used the PYMC3 package (Salvatier, Wiecki, & Fonnesbeck Reference Salvatier, Wiecki and Fonnesbeck2015), 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. (Reference Stanimirović, Staveley-Smith, Dickey, Sault and Snowden1999), 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. (Reference Martin, Blagrave, Lockman, Pinheiro Gonçalves, Boothroyd, Joncas, Miville-Deschênes and Stephan2015) demonstrated through similar SPS that the noise is effectively the combination of receiver and sky noise (Boothroyd et al. Reference Boothroyd, Blagrave, Lockman, Martin, Pinheiro Gonçalves and Srikanth2011), $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. (Reference Boothroyd, Blagrave, Lockman, Martin, Pinheiro Gonçalves and Srikanth2011) 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. Reference Hotan2021).

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.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. (Reference McClure-Griffiths2018) 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. (Reference Staveley-Smith, Sault, Hatzidimitriou, Kesteven and McConnell1997) 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. Reference Stanimirović, Staveley-Smith, Dickey, Sault and Snowden1999), 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. Reference Gordon2014). 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. Reference Winkler, Smith, Points and Team2015) and other star burst regions such as the N66/NGC 346 complex (Heydari-Malayeri & Selier Reference Heydari-Malayeri and Selier2010). 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. Reference Meixner2013). 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. Reference Rachford2009), 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. Reference Di Teodoro2019a), 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. (Reference Clark, Peek and Putman2014) 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. Reference Leroy2021). 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. Reference Bialy, Sternberg, Lee, Le Petit and Roueff2015 and references therein), traces sites of star formation. Glover & Clark (Reference Glover and Clark2012) 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. Reference Lee2012; Burkhart et al. Reference Burkhart, Lee, Murray and Stanimirović2015; Wang et al. Reference Wang2020) 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. (Reference Di Teodoro2019b), a majority of these anomalous features likely originate from strong star-formation-driven outflows (McClure-Griffiths et al. Reference McClure-Griffiths2018), 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 Reference Mac Low and Klessen2004), 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. Reference Szotkowski2019), or energy is being injected over a wide range of scales without preference at any particular scale (e.g., Norman & Ferrara Reference Norman and Ferrara1996). 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. Reference Stanimirović, Staveley-Smith, Dickey, Sault and Snowden1999; Elmegreen & Falgarone Reference Elmegreen and Falgarone1996). 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. (Reference Stanimirović, Staveley-Smith, Dickey, Sault and Snowden1999). 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. (Reference Stanimirović, Staveley-Smith, Dickey, Sault and Snowden1999), 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. (Reference Stanimirović, Staveley-Smith, Dickey, Sault and Snowden1999). 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. Reference Falgarone, Phillips and Walker1991). 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. (Reference Sun2018) 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. Reference Vázquez-Semadeni, Palau, Ballesteros-Paredes, Gómez and Zamora-Avilés2019), 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. (Reference McClure-Griffiths2018) 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. Reference Di Teodoro2019b). 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 Reference Di Teodoro and Fraternali2015), 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. Reference McClure-Griffiths2018) 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. (Reference Di Teodoro2019b). 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. (Reference Di Teodoro2019b) 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. Reference Gaia2016; Gaia Collaboration et al. Reference Gaia2021) 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 Reference Hartmann1999). 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. (Reference Daigle, Joncas and Parizeau2007) 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 Reference Pon, Johnstone and Kaufman2012; Girichidis et al. Reference Girichidis, Naab, Walch and Berlok2021), transfer mass, momentum, and energy (Vijayan et al. Reference Vijayan, Kim, Armillotta, Ostriker and Li2020; Schneider et al. Reference Schneider, Ostriker, Robertson and Thompson2020), 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. Reference Banda-Barragán, Parkin, Federrath, Crocker and Bicknell2016; Gronke et al. Reference Gronke, Oh, Ji and Norman2021). Di Teodoro et al. (Reference Di Teodoro2019b) 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 Reference Phillips, Ramos-Larios and Perez-Grana2009; Grønnow et al. Reference Grønnow, Tepper-Garca, Bland-Hawthorn and McClure-Griffiths2017; Grønnow, Tepper-Garc a, & Bland-Hawthorn Reference Grønnow, Tepper-Garca and Bland-Hawthorn2018; Banda-Barragán et al. Reference Banda-Barragán, Zertuche, Federrath, Garca Del Valle, Brüggen and Wagner2019; Banda-Barragán et al. Reference Banda-Barragán, Brüggen, Heesen, Scannapieco, Cottle, Federrath and Wagner2021).

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. Reference McClure-Griffiths2018), 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. Reference Marchal, Miville-Deschênes, Orieux, Gac, Soussen, Lesot, d’Allonnes and Salomé2019) or GaussPy+ (Lindner et al. Reference Lindner2015; Riener et al. Reference Riener, Kainulainen, Henshaw, Orkisz, Murray and Beuther2019).