Kinematics of the Atomic ISM in M33 on 80 pc scales

Eric Koch, Erik Rosolowsky, Felix Lockman, Amanda Kepley, Adam Leroy, Andreas Schruba, J. Braine, Julianne Dalcanton, Megan Johnson, Snežana

Stanimirović

► To cite this version:

Eric Koch, Erik Rosolowsky, Felix Lockman, Amanda Kepley, Adam Leroy, et al. Kinematics of the Atomic ISM in M33 on 80 pc scales. Monthly Notices of the Royal Astronomical Society, 2018, 479 (2), pp.2505-2533. 10.1093/mnras/sty1674. hal-01829863

HAL Id: hal-01829863

https://hal.science/hal-01829863v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Kinematics of the atomic ISM in M33 on 80 pc scales

Eric W. Koch,<fnref n="1" />* Erik W. Rosolowsky,<fnref n="1" /> Felix J. Lockman,<fnref n="2" /> Amanda A. Kepley,<fnref n="3" />
Adam Leroy,<fnref n="4" /> Andreas Schruba,<fnref n="5" /> Jonathan Braine,<fnref n="6" /> Julianne Dalcanton,<fnref n="7" /> Megan
C. Johnson<fnref n="8" /> and Snežana Stanimirović<fnref n="9" />

Accepted 2018 June 20. Received 2018 June 18; in original form 2018 March 22

## ABSTRACT

We present new L-band (1–2 GHz) observations of the nearby spiral galaxy M33 with 80 pc resolution obtained with the Karl G. Jansky Very Large Array. The H I observations, combined with H I measurements from the Green Bank Telescope, improve the spectral resolution and sensitivity (2.8 K rms noise in a 0.2 km s⁻<fnref n="1" /> channel) compared to previous observations. We find individual profiles are usually non-Gaussian, harbouring line wings, multiple components, and asymmetries. Given this spectral complexity, we quantify the motions in the atomic interstellar medium through moment analysis of the spectra and fit to aligned, stacked profiles. The measured value of the H I line width depends strongly on the method used, with the velocity stacked profiles aligned to the peak velocity giving the minimum value of σ = 7 km s⁻<fnref n="1" /> and all other methods giving higher values (σ ~ 10 km s⁻<fnref n="1" />). All measurements of the line width show a shallow radial trend, with σ decreasing by ~2 km s⁻<fnref n="1" /> from Rgal = 0 to Rgal = 8 kpc. We consider a number of energy sources that might maintain the line width against turbulent dissipation, but no single source is adequate. We find excess emission relative to a Gaussian in the stacked profile line wings, ranging from 9 per cent to 26 per cent depending on how the spectra are aligned. By splitting the line wings into symmetric and asymmetric components, we find that the lagging rotational disc accounts for one-third of the line wing flux. We also find emission far from the rotation-axis of the galaxy in multiple discrete H I clouds, including a filament with a projected length of ~8 kpc.

**Key words:** ISM: kinematics and dynamics–galaxies: individual: M33–galaxies: ISM–radio lines: ISM.

## 1 INTRODUCTION

H I 21-cm emission is an ideal tracer of kinematics in the atomic interstellar medium (ISM) because of its ubiquitous distribution that extends beyond the optical extent of galactic discs. On large scales, H I emission can be used to study disc kinematics and rotation, while small-scale variations trace the kinematic and turbulent nature of the ISM.

Modelling the kinematics of the atomic ISM remains a significant challenge. The H I detected within the Milky Way exhibits kinematic features related to the evolution of clouds on small scales. Our limited perspective within the Galaxy makes it difficult to place these small-scale motions in the context of large-scale motions and the galactic potential. However, since the kinematics on both scales are linked, a consistent understanding of ISM kinematics requires sampling scales from molecular clouds to the entire disc (Dobbs et al. 2014). The external view of nearby galaxies can overcome these limitations in Galactic observations. However, these observations often lack spatial resolution and high sensitivity, yielding studies of

<footnote n="1">¹Department of Physics, University of Alberta, 4-183 CCIS, Edmonton, AB T6G 2E1, Canada</footnote>

<footnote n="2">²Green Bank Observatory, PO Box 2, Green Bank, WV 24944, USA</footnote>

<footnote n="3">³National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903-2475, USA</footnote>

<footnote n="4">⁴Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, OH 43210, USA</footnote>

<footnote n="5">⁵Max-Planck-Institut für extraterrestrische Physik, Giessenbachstraße 1, D-85748 Garching, Germany</footnote>

<footnote n="6">⁶Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>

<footnote n="7">⁷Department of Astronomy, Box 351580, University of Washington, Seattle, WA 98195, USA</footnote>

<footnote n="8">⁸United States Naval Observatory, 3450 Massachusetts Ave NW, Washington, DC 20392, USA</footnote>

<footnote n="9">⁹Department of Astronomy, University of Wisconsin, 475 N Charter St., Madison, WI 53706, USA</footnote>

<footnote>* E-mail: koch.eric.w@gmail.com (EWK); rosolowsky@ualberta.ca (EWR)</footnote>
21-cm emission that blend the small-scale motions seen in Milky Way studies.

This trade-off of resolution and scale between Galactic and extra-
galactic observations is also important for discerning the multiple
thermal states in the atomic medium. While superseded by later
work, the Field, Goldsmith & Habing (1969) model predicts the
presence of a warm neutral medium (WNM, *n* ~ 1 cm⁻³, *T* ~
1000 K) and cool neutral medium (CNM, *n* ~ 10 cm⁻³, *T* ~ 100 K)
over a range of pressure. Subsequent studies of the atomic medium
qualitatively support the two-phase picture (Wolfire et al. 1995;
Wolfire et al. 2003). However, in observational studies, it remains
difficult to separate these two gas phases. The optically thin WNM
is typically measured in H₁ emission, while the optically thick CNM
can only be unambiguously separated through absorption studies.

Decomposing H<sub>I</sub> spectra into thermal components has best been studied in the Milky Way through absorption towards bright extra-galactic sources (Heiles & Troland 2003; Murray et al. 2015), or by tracing H<sub>I</sub> self-absorption (HISA; Gibson et al. 2005). These studies find that the atomic ISM is found in both the WNM and CNM states, though a non-negligible fraction of components are in an unstable intermediate state (Heiles & Troland 2003; Murray et al. 2015). Extragalactic absorption studies towards ~50 lines of sight in Local Group galaxies suggest that the fraction of the atomic medium in the CNM varies between local systems (Dickey, Hanson & Helou 1990; Dickey et al. 1994, 2000).

Most extragalactic studies focus on interpreting the H<sub>I</sub> emission line profiles at the telescope resolution (typically 200–1000 pc), though different modelling approaches are used. Fitting individual spectra with a Gaussian is the most straight-forward approach (Boulanger & Viallefond 1992; Mogotsi et al. 2016), however non-Gaussian line features or multiple components require a more sophisticated model, such as a two-Gaussian model (Young & Lo 1996; Warren et al. 2012). A simpler approach is to only estimate the line width of a spectrum using the second moment (Tamburro et al. 2009), though this approach has been found to overestimate the line width in many cases (Mogotsi et al. 2016).

To study kinematics in in faint regions, several studies stack
velocity-aligned spectra to increase the signal-to-noise (S/N) in the
data. Nearly all studies that utilize this method find a common line
shape of a central Gaussian peak with enhanced line wings (Young
& Lo 1996; Braun 1997; Petric & Rupen 2007; Ianjamasimanana
et al. 2012; Stilp et al. 2013a). To characterize this shape, the stacked
spectra are either fit with a two-Gaussian model (Young & Lo 1996;
Ianjamasimanana et al. 2012) or as a Gaussian peak with a non-
parametric measure of the enhanced line wings (Stilp et al. 2013a,
b). The H1 line widths from a single Gaussian range from 6 to 10
km s⁻¹, while the two-Gaussian models have a narrow component
width of 3–6 km s⁻¹, and 9–25 km s⁻¹ in the wider component.

The results of stacked profile studies, and the models chosen to
explain the profile shape, have led to diverging physical interpreta-
tions of the atomic medium. Studies using a two-Gaussian model
argue that the narrow and broad components are naturally explained
by the CNM and WNM, respectively (Young & Lo 1996; Ianjamasi-
manana et al. 2012). On the other hand, Stilp et al. (2013a) argues
that the central Gaussian peak represents emission from a turbulent
mixture of CNM and WNM, and the enhanced line wings result
from stellar feedback. Alternatively, Braun (1997) proposes that the
brightest emission arises from a 'high-brightness network' (HBN)
of narrow filamentary structure across the disc. The HBN then arises
from a dense and optically thick CNM component, where the cen-
tral peak becomes flattened in individual spectra due to the higher
optical depth (Braun et al. 2009; Braun 2012).

The use of different methods and their disparate interpretations makes it difficult to create a clear connection with Milky Way studies. Yet, with sufficient sensitivity and resolution, extragalactic H₁ observations should recover the small-scale complexity observed in the Milky Way. New observations of nearby galaxies can provide this connection between galactic and extragalactic approaches. Observations of the Local Group provide the best means for determining this connection as they allow for scales of individual molecular clouds to be resolved, similar to studies within the Milky Way and Magellanic Clouds (e.g. Fukui et al. 2009; Wong et al. 2009). In particular, M33 is an ideal target due to its proximity (840 kpc, Freedman et al. 2001), moderate disc inclination (55°), and relatively small angular size. To that end, we have conducted new observations of M33 with the NSF's Karl G. Jansky Very Large Array (VLA) that focus on high sensitivity coupled with a high spectral resolution.

Recent work on M33 has focused on deep, lower resolution H<sub>I</sub> observations seeking the origins of H<sub>I</sub> halo gas and high-velocity clouds (HVCs) in the M31 group (Lockman, Free & Shields 2012; Keenan et al. 2016), or the structure of the dark matter halo (Corbelli et al. 2014; Kam et al. 2017). There has been comparatively less attention on the nature of the atomic gas in M33’s star-forming disc. The first complete H<sub>I</sub> map with the resolution to discern small-scale structure was presented by Deul & van der Hulst (1987), which unveiled the filamentary structure of M33’s H<sub>I</sub> morphology. More recently, Gratier et al. (2010) use archival H<sub>I</sub> VLA observations from Thilker, Braun & Walterbos (2002) to compare the H<sub>I</sub> properties around giant molecular clouds (GMCs). Druard et al. (2014) further perform a stacking analysis of the H<sub>I</sub> and CO to study the relationship of atomic and molecular gas. Imara, Bigiel & Blitz (2011) use an independent reduction of these same data to examine the H<sub>I</sub> environments around 45 GMCs in M33.

In this paper, we present new L-band (1–2 GHz) observa-
tions from the VLA of M33 with a focus on high velocity-
resolution observations of the 21-cm H1 line. Taken in the VLA’s
C-configuration, the beam size at the H1 line is 18–20 arcsec or a
physical size of ~80 pc. The VLA’s correlator allows for simultane-
ous observations of the 21-cm H1 line, four OH transitions, several
hydrogen radio recombination lines (RRLs), and the polarized radio
continuum. Here we focus on the H1 observations to analyse the H1
profile shapes on 80 pc scales and to determine the structure and
kinematics of the H1 disc. We detect no RRL emission (Section 3.1)
and find a single OH maser (Koch et al. 2018). The polarized radio
continuum will be presented in a future paper.

In Section 2 we present the VLA observations and the Green
Bank Telescope (GBT) H1 data used to provide short-spacing in-
formation (Lockman et al. 2012). Aspects of the H1 imaging, signal
masking, and upper limits on the hydrogen RRL emission are given
in Section 3. We explore the atomic gas properties of M33 in Sec-
tion 4 and extra-planar H1 structures in Section 5. In Section 6,
we critically evaluate the meaning of 21-cm H1 line profiles gener-
ated in extragalactic observations. Appendix A provides a detailed
description of the imaging and short-spacing combination.

2 OBSERVATIONS

2.1 VLA

Using the VLA in C-configuration, we observed a 13-point mosaic
covering the disc of M33 out to a radius of 12 kpc (Project 14B-
088). The data were taken in 12 tracks, split equally before and after
transit, for a total of 52 h. We used a hexagonal grid of pointings at
<table><caption>**Table 1.** The VLA SPW setup for spectral lines used in the observations.</caption><thead><tr><th>Line</th><th>Rest freq. (GHz)</th><th>Channel width (km s<sup>-1</sup>)</th><th>Channels</th></tr></thead><tbody><tr><td>H<sub>I</sub></td><td>1.42</td><td>0.206</td><td>8192</td></tr><tr><td>OH(1612)</td><td>1.612</td><td>1.45</td><td>512</td></tr><tr><td>OH(1665)</td><td>1.665</td><td>1.41</td><td>512</td></tr><tr><td>OH(1667)</td><td>1.667</td><td>1.40</td><td>512</td></tr><tr><td>OH(1720)</td><td>1.720</td><td>1.36</td><td>512</td></tr><tr><td>H(172)<sup>α</sup></td><td>1.28</td><td>1.26</td><td>512</td></tr><tr><td>H(166)<sup>α</sup></td><td>1.42</td><td>1.58</td><td>512</td></tr><tr><td>H(164)<sup>α</sup></td><td>1.48</td><td>1.83</td><td>512</td></tr><tr><td>H(158)<sup>α</sup></td><td>1.65</td><td>1.64</td><td>512</td></tr><tr><td>H(153)<sup>α</sup></td><td>1.82</td><td>1.42</td><td>512</td></tr><tr><td>H(152)<sup>α</sup></td><td>1.85</td><td>1.29</td><td>512</td></tr></tbody></table>

the H<sub>I</sub> frequency to ensure uniform sensitivity across the mosaic. To remedy a pointing error that omitted one of the mosaic centres, our final observation used a single ~5 h pointing on that location, yielding nearly equal integration times and noise properties across the mosaic. Due to its close proximity to M33, we use 3C 48 as the flux, delay, and gain calibrator. 3C 138 is observed as the polarization leakage calibrator for the continuum data.

Using the full capabilities of the VLA correlator, our spectral setup covers the entire 1–2 GHz L-band polarized continuum, the H<sub>I</sub> line, four OH lines, and six hydrogen RRLs, enabling a variety of science cases to be explored. The setup of the line spectral windows (SPWs) are shown in Table 1. The H<sub>I</sub> SPW has a high spectral resolution of 977 Hz, corresponding to a velocity of 206 m s<sup>-1</sup>. This high-spectral resolution is required to detect HISA, based on prior Milky Way observations (e.g. Gibson et al. 2005).

## 2.2 GBT

We reprocessed the high spectral resolution H<sub>I</sub> data presented in Lockman et al. (2012, project AGBT09A_17) to provide the short and zero spacing information on H<sub>I</sub> emission. The angular resolution of the GBT in the 21-cm line is 9.1 arcmin and the spectral resolution is 0.16 km s<sup>-1</sup>. Spectra were calibrated and corrected for stray radiation as described in Boothroyd et al. (2011), and a first-order polynomial was fitted to emission-free regions of the spectra to remove any residual instrumental baselines.

The data were originally collected as four separate maps with minimal spatial overlap, so we gridded the data into a cube using a Gaussian kernel, rather than the preferred Gaussian-Bessel kernel. This choice eliminates edge effects in the data at the centre of the galaxy, but lowers the effective angular resolution to 9.8 arcmin (see Section A1.1), which is a linear resolution of ~2.3 kpc at the distance of M33. We used the *gbtpipe* package<fnref n="1" /> to build the data cube, which performs spectral pre-processing to eliminate bad scans and a kernel based gridding approach developed for on-the-fly data (Mangum, Emerson & Greisen 2007). The noise in the final cube varies with position ranging from 50 to 120 mK in a 0.16 km s<sup>-1</sup> channel giving a median column density error of 3.1 × 10<sup>17</sup>cm<sup>-2</sup> for a 20 km s<sup>-1</sup> full width at half-maximum (FWHM) line. The shortest baselines in our VLA data provide information on scales of 16.5 arcmin, providing sufficient overlap in uv-space with the GBT beam to calculate comparison statistics (Appendix A1).

## 3 IMAGING AND CALIBRATION

We calibrated the data using a modified version of the VLA pipeline (version 1.3.0) with CASA 4.2.2.<fnref n="2" /> Following observatory recommendations, we used a modified version of the pipeline to better handle line SPWs, namely not using Hanning smoothing or automated radio frequency interference flagging on the line SPWs to maintain the spectral resolution and avoid flagging of narrow emission features (e.g. OH maser emission). Using 3C 48 as the sole calibrator (excluding polarization) gives a high S/N for all calibration scans, and the automated calibration solutions found by the pipeline are excellent. After the initial pipeline run, we manually flagged the data before running the pipeline once more. Most of the line SPWs required little manual flagging beyond the pipeline solutions, though one RRL SPW could not be recovered due to RFI.

Here we present upper limits for RRL emission and detail the imaging and masking approaches adopted for the H<sub>I</sub> data. A single OH(1665) maser is detected and is presented separately (Koch et al. 2018). We find no detections in the 1612, 1667, or 1720 MHz OH lines.

### 3.1 Radio recombination lines

We detect no RRL emission from the six observed transitions. One RRL SPW, the H(172)<sup>α</sup> line, was dominated by RFI and was unrecoverable. Prior to imaging, we subtracted a constant continuum background level using the CASA UVCONTSUB task by fitting the velocity ranges found to be emission-free in the H<sub>I</sub> imaging (Section 3.2). RRL emission is expected to be very faint and based on a prior RRL detection towards NGC 604 (Araya, Baan & Hofner 2004, see below), we do not expect to detect any signal in a single line. To lower the noise level, we average the individual line cubes together.<fnref n="3" /> To use this process, we generated data cubes for each line at a common spectral resolution (~10 km s<sup>-1</sup>) and corrected these by the primary beam pattern of the mosaic. Each line was then convolved to a common spatial resolution, set by the lowest resolution line, H(166)<sup>α</sup>. Because the mosaic pointings were set based on the primary beam size at the H<sub>I</sub> frequency, the primary beam pattern and sensitivity differs between the RRLs. However, based on the peak H<sub>α</sub> flux from the giant H<sub>II</sub> region NGC 604, we place limits where the RRL intensity is expected to be brightest. Because of the pointing issue mentioned in Section 2.1, the region around NGC 604 also has the best primary beam coverage for all RRLs observed. Using these regions, we set a 3<sup>σ</sup> upper limit on the RRL intensity at 3.0 mJy within a 60 arcsec region that encompasses the optical extent of the H<sub>II</sub> region. This upper limit in agreement with observations by Araya et al. (2004) with Arecibo in C band, where they detected the H110<sup>α</sup> line in NGC 604 with a peak line flux of 1.36 ± 0.19 mJy in a 58 arcsec beam.

### 3.2 H<sub>I</sub> imaging

We subtract a constant continuum component from the H<sub>I</sub> data using the CASA task UVCONTSUB. The continuum level is fitted based on identified emission-free channels in both the VLA and GBT data. The VLA data are inspected in the uv-plane to determine emission-free channels. We then use the GBT data to determine the velocities

<footnote n="1">¹v0.1.2, https://github.com/low-sky/gbtpipe</footnote>

<footnote n="2">²https://science.nrao.edu/facilities/vla/data-processing/pipeline/scripted-pipeline</footnote>

<footnote n="3">³https://casaguides.nrao.edu/index.php?title=Stacking_Multiple_Spectra_LLines_at_Same_Position</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Total emission profiles of the VLA-only and the combined VLA+GBT data before and after masking. The curve for the unmasked combined data set is indistinguishable from an unmasked GBT-only spectrum. The thick dashed vertical line indicates the systemic velocity of M33 at −179.2 km s⁻¹. The flux density in the masked VLA map is well recovered by the VLA data in the Northern half of M33, but the Southern half lacks a significant portion of the emission without short-spacing information.</figcaption></figure>

affected by Milky Way H I emission. The large bandwidth covered by the H I line gives an ample number of H I-free lines, and the fitted background level is well-constrained.

Examining the extent of Milky Way H I emission in the GBT data is imperative for the red-shifted H I emission at the southern-most tip of M33, as they become spatially coincident near $v_{LSRK} \sim -71$ km s⁻¹. This is indicated by the shaded region in Fig. 1. This overlap region is clearly shown in the H I Arecibo data from Putman et al. (2009, see their fig. 3). Based on this, we begin imaging the VLA data at a velocity of −73.0 km s⁻¹ to best avoid Milky Way emission. We estimate that excluding this spectral region removes ~50 Jy km s⁻¹ of flux from M33, assuming that the missing emission is similar to the blue-shifted side. This flux corresponds to 0.5 per cent of the total emission from M33.

The resulting data cube is large with 1178 velocity channels, each of which has a grid size of 2560² pixels. This presents significant computational barriers for imaging; our imaging process is described in Appendix A. Prior to imaging, we determine optimal CLEAN settings by running numerous combinations on a single channel, identifying the imaging parameters that yield the lowest peak residual without the algorithm diverging. Only natural weighting is tested as we prioritize maximizing sensitivity to extended structure. A CLEAN mask is defined for each channel based on the $3\sigma$ limit in the GBT data, which covers nearly all of the emission in the VLA data. Some channels with emission near the mask edge used an expanded mask to ensure all emission was included.

We use the multi-scale CLEAN algorithm (Cornwell 2008) for deconvolution with six scales, ranging from a point response to a quarter the largest recoverable scale (~970 arcsec). Changing the range and the specific scales did not have a significant effect on the resulting image. Based on our single channel tests, we deconvolve each channel until reaching 3.8 mJy beam⁻¹ (7.1 K). This limit is 2.5σ times the noise level in the final cube.

We then use the GBT data to provide short-spacing information for the deconvolved VLA cube by feathering the VLA and GBT data cubes together with the UVCOMBINE package.<fnref n="4" /> A detailed explanation

of the combination is provided in Appendix A1. The resulting cube fully recovers the total emission in the GBT data.

The final VLA H I cubes, with and without short-spacing data, have a $1\sigma$ sensitivity of 2.8 K per 0.2 km s⁻¹ channel and a beam size of 19 arcsec × 17 arcsec. This corresponds to a hydrogen column density of $1.0 \times 10^{18}$ cm⁻², assuming optically thin emission. The mosaic is cut-off at a primary beam coverage of 0.5 to avoid noise-dominated regions near the edge.

### 3.3 Signal masking

We define a multistep process for creating a reliable signal mask enabled by the high spectral resolution of our data. This masking is critical for our analysis to ensure that moment-based estimations are not influenced by noise. When using the VLA-only data, the masking also removes the influence of negative bowls in the data from missing short-spacing information.

Assuming a typical Gaussian line width of $\sigma = 6$ km s⁻¹ for the broad H I emission across the disc (see Section 4.4), each spectral feature will have an FWHM that spans about 70 channels. For the purpose of characterizing emission features alone (not HISA), we can substantially smooth in the spectral dimension to highlight low surface brightness features.

To determine the mask, we first spectrally smooth the data using a median filter with a width of 31 channels (~6 km s⁻¹). The noise in each smoothed spectrum is found using the median absolute deviation (MAD), iteratively rejecting points beyond $2\sigma_{MAD}$ until convergence. We then search for valid spectral components by requiring that each component has a maximum intensity above $5\sigma_{MAD}$ and has 30 consecutive channels above $2\sigma_{MAD}$. To ensure low surface-brightness line wings are included, the edges of the valid components are extended until the intensity reaches $1\sigma_{MAD}$. This procedure is performed on all spectra within the cube.

Next, we consider the spatial connectivity. To remove spurious spectral components found in the first step, we use the morphological opening and closing operators with a top-hat kernel equivalent to the FWHM of the beam (Shih 2009). The opening operator will erode the mask edges, removing all features smaller than the beam kernel. The closing operator then dilates the remaining components in the mask to restore their original area. The combination of these two masking steps yields a robust signal mask [$M(x, y, v)$], consistent in both the spectral and spatial dimensions.

This procedure is essential for high spectral-resolution data. Throughout the map, we find many lines of sight with either multiple, resolved components, or significantly skewed profiles. Not recovering the line wings in the masking procedure may hide these asymmetries at low-surface brightness. This also ensures negative bowls in the VLA-only cube are rejected from the signal mask and will not bias measurements of the asymmetry of a spectrum.

In Fig. 1, we compare the total GBT emission measured within the VLA mosaic to the emission retained after applying the masking procedure. For the VLA-only data, 60 per cent of the total H I emission is contained within the mask. The combined VLA+GBT data contains 87 per cent of the total emission in the mask. The missing flux density in the masked version arises from low surface-brightness features that do not satisfy the masking criteria. We discuss these low surface-brightness features in Section 5. The ratio of flux density retained in the mask is roughly constant across all channels for the combined VLA+GBT data, indicating that only low surface brightness H I is excluded by the masking procedure.

<footnote n="4">⁴https://github.com/radio-astro-tools/uvcombine</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** H I column density map of the VLA + GBT mosaic, assuming optically thin emission. The map is masked using the technique described in Section 3.3, and shown using an arcsinh stretch. The mosaic recovers emission well out to a galactic radius of 10 kpc.</figcaption></figure>

# 4 PROPERTIES OF THE ATOMIC MEDIUM

In this section, we examine the properties of the H I emission derived from the new data in the context of previous H I studies (Corbelli et al. 2014; Druard et al. 2014; Kam et al. 2017).

We calculate an atomic mass of $9 \pm 2 \times 10^8$ M$_{\odot}$ for the VLA data and $1.3 \pm 0.3 \times 10^9$ M$_{\odot}$ for the combined data, including a factor of 1.4 for He and heavier elements. When including regions outside of the VLA mosaic, the GBT data give a total atomic mass of $1.8 \pm 0.2 \times 10^9$ M$_{\odot}$, consistent with other independent H I mass measurements (e.g. Putman et al. 2009).

Fig. 2 shows the column density map of the VLA data, highlighting M33’s flocculent structure. The inner region of the disc ($R_{\text{gal}} < 2$ kpc; see Fig. 3 for context) is dominated by small-scale H I shells and lacks large-scale (~kpc) structures. The ‘mid-disc’ ($R_{\text{gal}} = 2–4$ kpc) is dominated by spiral arms prominent in the optical (e.g. fig. 2 in Corbelli et al. 2014). The H I shows similar spiral arm structure in the Northern half, while the Southern arm is dominated by one of the brightest H I clumps in the galaxy (Rosolowsky et al. 2007). The H I in the outer disc ($R_{\text{gal}} > 4$ kpc) is nearly ubiquitous, with bright spiral-arm segments throughout. The outer edge of our signal-masked H I map is near the radius where the warp in M33’s disc becomes significant (Corbelli et al. 2014).

Fig. 3 shows four spectra extracted from different positions in the H I data cube. Each of the spectra shows a non-Gaussian profile

with either multiple velocity components, extended line wings, or both. To describe these features, we use moment-based descriptions of the line profiles and velocity-aligned spectral stacking. Given the lines are clearly non-Gaussian, direct Gaussian fits to the lines do not produce a complete description of the emission.

We note that all uses of the line width in this paper are defined as the Gaussian standard deviation, not the FWHM.

## 4.1 Rotation curve and disc parameters

Using the masked cube from Section 3.3, we calculate the velocity at the peak intensity, shown in the left-hand panel of Fig. 4, to derive a rotation model. We define this velocity at peak intensity as $v_{\text{peak}} \equiv \text{argmax}_{v} I(v)$, which is commonly referred to as the ‘peak velocity’. To find $v_{\text{peak}}$, we smooth the data with a $1 \text{ km s}^{-1}$ (~5 channel) Gaussian kernel to minimize noise before identifying the velocity of the peak. In Appendix B, we discuss the difference in using $v_{\text{peak}}$ versus the centroid (first moment) velocity, $v_{\text{cent}}$, in deriving the rotation curve. Briefly, $v_{\text{cent}}$ becomes biased by asymmetric line wings (Section 4.3); whereas $v_{\text{peak}}$ is not biased by the line shape and should provide a more accurate representation of the rotation velocity.

We fit the velocity surface using DISKFIT (Spekkens & Sellwood 2007; Sellwood & Spekkens 2015) with a circular velocity model.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Example H₁ spectra compared between the archival VLA (grey) and the new VLA + GBT H₁ data (black). The spectra are centred at the velocity of peak intensity, calculated by smoothing with a 2 km s⁻¹ Gaussian kernel. The location of the spectra are indicated by the arrows on the column density map in the left-hand panel. The new data highlight the complex line shapes near the peaks and demonstrate that signal seen in the extended line wings is real. The red-dashed box on the left-hand panel shows the position where a HVC is found (Section 5.3) at a velocity of −30 km s⁻¹ from the galaxy’s rotation. The cyan contours on the map indicate galactocentric radii of 2, 4, 6, and 8 kpc, respectively, in the plane of the galaxy.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left: H₁ *v*<sub>peak</sub> surface used for fitting the rotation curve. The masked regions from the signal mask (Section 3.3) are shown in white. Centre: Residual velocities from subtracting the DISKFIT rotation model from the peak velocity surface. Right: The line width map derived from the second moment.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Rotation curve fit by DISKFIT from the VLA + GBT $v_{\text{peak}}$ surface (blue solid). The rotation curves from Corbelli et al. (2014, green-dashed) and Kam et al. (2017, red-dotted) are also shown. Our model is consistent with the Kam et al. (2017) model, but has a lower velocity than the Corbelli et al. (2014) curve in the inner 6 kpc.</figcaption></figure>

<table><caption>**Table 2.** Galactic disc parameters and fit statistics from DISKFIT for the VLA + GBT peak velocity surface. Errors are the $1\sigma$ intervals based on the DISKFIT bootstrapping. The $n$, $v_{\text{max}}$, and $R_{\text{max}}$ parameters are the fit parameters to the Brandt (1960) model (equation 1).</caption><tr><td>Points used</td><td>39 970</td></tr><tr><td>Degrees of freedom (d.o.f.) for error</td><td>39 663</td></tr><tr><td>χ<sup>2</sup>/DOF</td><td>1.44</td></tr><tr><td>Iterations</td><td>4</td></tr><tr><td>Centre R.A.</td><td>23.4607 ± 0.0042°</td></tr><tr><td>Centre Dec.</td><td>30.6583 ± 0.0032°</td></tr><tr><td>Position angle</td><td>201.12 ± 0.47°</td></tr><tr><td>Inclination</td><td>55.08 ± 1.56°</td></tr><tr><td>v<sub>sys</sub> (LSRK)</td><td>-179.18 ± 0.76 km s<sup>-1</sup></td></tr><tr><td>n</td><td>0.56 ± 0.04</td></tr><tr><td>v<sub>max</sub></td><td>110.0 ± 1.5 km s<sup>-1</sup></td></tr><tr><td>r<sub>max</sub></td><td>12.0 ± 1.3 kpc</td></tr></table>

DISKFIT derives a global rotation model by simultaneously fitting the whole velocity surface over a given set of radial bins. This provides better constraints than minimizing a set of individual rings, but requires that the disc be characterized by a common set of disc parameters (e.g. position angle, inclination). The rotation curve is shown in Fig. 5, and the fit parameters and statistics are given in Table 2. The rotation velocities in each bin are given in Table C1 in the Appendix. The parameter uncertainties are calculated from 200 bootstrap iterations in DISKFIT (Sellwood & Sánchez 2010). The $\chi^2$ calculation in DISKFIT includes an ISM line width parameter for the expected model dispersion, which we set to $8 \text{ km s}^{-1}$ (see Section 4.2). However, since this is a constant factor applied to every position in the fit, it will not affect the resulting disc parameters (Kuzio de Naray et al. 2012). Our rotation model does not include a warp component since a warp only becomes prominent at galactocentric radii beyond ~8 kpc (Corbelli et al. 2014) and we do not fit for a bar component. We use radial bins with a constant width of 100 pc in the plane of the galaxy, just larger than the VLA beam size, out to a galactocentric radius of 7.5 kpc. Because the radial bins are larger than the beam, we do not include a correction for beam smearing. Beyond this radius, the uncertainty in the derived moments increases substantially as the integrated intensity decreases. The fitted position of the galactic centre is 7 arcsec from

the 2MASS position (Skrutskie et al. 2006), within the $1\sigma$ parameter uncertainties.

Fig. 5 shows the circular rotation velocities, $v_{\text{rot}}$, from our rotation model, along with the recent H I models from Corbelli et al. (2014) and Kam et al. (2017). There is excellent agreement with the Kam et al. (2017) results, though our model has rotation velocities below the Corbelli et al. (2014) rotation curve up to 6 kpc. The difference may be due to how Corbelli et al. (2014) apply finite disc corrections (see their Appendix B), but since both the data and modelling approaches differ for the three rotation curves, it is difficult to be sure of the cause.

The residual velocity surface from the difference between $v_{\text{peak}}$ and $v_{\text{rot}}$ is shown in Fig. 4. The residual velocity surface shows only small-scale variations, suggesting that additional velocity components in the rotation model are not needed. However, the outer edges of the mosaic show what is likely the beginning of the warp. Large deviations (~20 km s$^{-1}$) both above and below the model values are evident along the edges, similar to the deviations found by Kam et al. (2017).

To remove bin-by-bin variations in the fitted rotation curve, we follow the prescription used by Meidt et al. (2008) by fitting a Brandt (1960) rotation curve:

$$v_{\text{rot}}(r) = \frac{v_{\text{max}}(r/r_{\text{max}})}{[1/3 + 2/3(r/r_{\text{max}})^n]^{3/2n}} \quad (1)$$

The fit parameters for the Brandt model are given in Table 2. Using the analytical approximation, we create a smooth version of the rotation surface for use as a model. This smooth rotation model surface is used for creating rotation-subtracted versions of the H I cube.

## 4.2 Surface density profiles

Using the galactic disc parameters in Table 2, we create a radial profile of the H I surface density, corrected for the disc inclination, in 100 pc radial bins out to 10 kpc. Assuming optically thin emission, we use a mass conversion factor of $0.019 \text{ M}_\odot \text{ pc}^{-2} (\text{K km s}^{-1})^{-1}$ to find the H I surface density, which is the standard H I column density conversion factor with a factor of 1.4 for the mass of heavier elements. Fig. 6 shows the radial surface density profile of the combined VLA and GBT data. The profiles presented here use smaller bin widths than those in Druard et al. (2014), Corbelli et al. (2014), and Kam et al. (2017), but are consistent with all of these previous works.

The profile has an average value consistent with a surface density of ~8 $\text{M}_\odot \text{ pc}^{-2}$ out to 7 kpc, where the surface density begins to taper off. The radius where the tapering begins is comparable to where the stellar surface density equals the gas surface density, based on fig. 10 in Corbelli et al. (2014). Within this radius, the average stellar surface density exceeds the gas component.

When averaged over galactocentric radii of < 7 kpc, the atomic gas surface density is nearly constant, consistent with the surface density profile from Kam et al. (2017). However, averaging over smaller radii highlights the deviations driven by the large-scale disc structure. When averaging with the VLA-only data, the average value is ~6 $\text{M}_\odot \text{ pc}^{-2}$ for the inner 7 kpc, roughly consistent with the fraction of emission recovered without the GBT data.

The profile shows two significant peaks located around 2 and 4.5 kpc; these occur at radii that enclose the prominent spiral arm structure. Within ~2 kpc, the disc lacks spiral structure and is dominated by wind-driven shells (Rosolowsky et al. 2007), making it difficult to connect these large-scale averages to the morphology
of the emission. The outer disc also lacks significant coherent large-scale structure, with the emission dominated by multiple spiral arm fragments.

To further examine these variations in the surface density profile, we create radial profiles split into the Northern and Southern halves. The aforementioned peaks are far more prominent in the Northern half than the Southern. This is even more dramatic for the dip at 3 kpc, which is entirely driven by variations in the Northern half. These differences are driven by the asymmetric main spiral arms: the Northern arm, in H1 emission, is more prominent and distinct than the southern arm.

The surface densities between the two halves are approximately equal between a radius of 1.5–2 kpc and beyond 5.5 kpc. The region from 1.5 to 2 kpc contains the beginning of the main spiral arms. The H1 arms begin at ~1.8 kpc are offset by nearly 180°, and are nearly symmetrical, despite tracing different extents at larger radii. The other matching region at 5.5 kpc is beyond the extent of the main spiral arms where the disc becomes dominated by flocculent structure. These matching regions over all angles suggest that this middle region of the disc is dominated by an arm driving mechanism and distinct from the inner and outer regions in the disc. However, the spatial scale of these variations may be different due to the large area averaged over at larger radii.

The outer region of the disc shows less variation between the halves compared to the inner and middle disc. While this region contains variations in the surface density due to the spiral arm fragments, the area averaged over is also larger and will tend to smooth out variations more than the regions at smaller radii.

## 4.3 Exploring spectral complexities

With the high spectral resolution in our data, there are a sufficient number of channels to use higher order spectral moments for describing line profile shapes. Fig. 3 demonstrates that simple spectral models (one or two Gaussians) are inappropriate to represent the data at this resolution. We quantify these features by calculating maps of the skewness and kurtosis of the velocity profiles. These higher order moments require a larger number of channels than the line width (second moment) to reduce uncertainties, which has limited their use in previous extragalactic H1 studies with coarser resolution. Instead, previous studies have relied on adding equivalent parameters to Gauss–Hermite polynomials (e.g. Young et al. 2003; Oh et al. 2011). However, this analytic form is limited in how well it can explain extended line wings, making the moments-based description better suited for our data.

The skewness is calculated from the third moment of each spectrum, normalized by the line width cubed:

$$S(x, y) = \frac{\sum_{v \in \mathcal{M}(x, y, v)} [v - M_1(x, y)]^2 I(x, y, v) \delta v}{M_0 \sigma_v^3}, \quad (2)$$

where $M_0$ is the integrated line intensity (zeroth moment), and $M_1$ is the line centroid (first moment). The line width here is defined from the second moment: $\sigma_v^2 = M_2(x, y)$. A positive skewness indicates a spectrum with a red-shifted tail, while a spectrum with negative skewness has a blue-shifted tail.

Kurtosis is similarly calculated from the fourth moment. We report the *kurtosis excess* by subtracting the value of 3 expected for a normal distribution:

$$K(x, y) = \frac{\sum_{v \in \mathcal{M}(x, y, v)} [v - M_1(x, y)]^4 I(x, y, v) \delta v}{M_0 \sigma_v^4} - 3. \quad (3)$$

<figure><img src="image_8.png" /><figcaption>**Figure 6.** The combined GBT and VLA atomic gas surface density profile (black) split into the Northern (red dot-dashed) and Southern (blue dashed) halves in 100 pc bins and corrected for the disc inclination. The surface density has an average value of $8 M_{\odot} \text{ pc}^{-2}$ in the inner 7 kpc, with substantial variation from galactic structure. The error bars are shown for the Northern and Southern halves; the errors for the entire disc (not shown) are $\sqrt{2}$ smaller. These uncertainties are the standard deviations in each radial bin, corrected by the number of beams (i.e. the independent samples) within the annulus (e.g. Druard et al. 2014). The shaded region indicates the inner 0.5 kpc where beam smearing and a small number of samples in each bin gives large uncertainties.</figcaption></figure>

A spectrum with a positive kurtosis will have strong tails and have a peak sharper than a Gaussian.

The skewness maps are shown in Fig. 7 for the VLA-only and combined data. The addition of the total power component from the GBT data has a significant effect on the skewness. The VLA-only map shows small-scale structure across much of the disc, indicating variations of the spectral shape down to the resolution of the data. This map also shows some larger scale variations, particularly for the outer disc in the Northern half. With the inclusion of the GBT data, these large-scale variations become far more prominent. The skewness structure is split into a positively skewed Northern half and a negatively skewed Southern half. This separation is highlighted in Fig. 8, where the radial profiles of skewness are shown for the Northern and Southern halves. The spectral shapes are skewed to the red-shifted side in the Northern half and to the blue-shifted side in the Southern half. This is consistent with the position–velocity slices in Section 5 that show a lagging rotational H1 component across the disc.

The kurtosis maps are less interesting due to the influence of the line wings, and this drives a strong correlation to the peak temperature map. Fig. 9 shows the close resemblance between the peak temperature and kurtosis maps for the combined data. The excess line wings bias the kurtosis, making it a measure of the spectral shape relative to the line wings rather than the bright H1 component. A higher peak temperature then leads to a larger positive kurtosis. We discuss this relation further in Section 4.3.1. Since kurtosis is velocity-weighted by the fourth power, the assumption that line shapes are well-described by single-peaked near-Gaussian breaks down and it provides less-useful empirical relations compared to the skewness.

The skewness profiles in Fig. 8 are nearly mirror images of each other in the Northern and Southern halves. By averaging over the entire disc, this distinction between the halves is lost, and the average values are near to a skewness of zero over most radii. This hides interesting variations seen separately in the halves, particularly near
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The skewness maps of the VLA (left) and combined data (right). Positive skewness indicates a red-shifted tail, while a negative value is a line profile with a blue-shifted tail. The low surface-brightness component, primarily added to the line wings, from the GBT data drives the large-scale skewness properties.</figcaption></figure>

3 kpc, where the average skewness in both halves is close to zero. In the surface density profiles (Fig. 6), this radius corresponds to a dip in surface density, though it is more prominent in the Northern half. This region contains a large portion of the Northern spiral arm, including some of the more active star-forming regions (e.g. NGC 604). Apart from this, it is unclear why only this region has a skewness close to zero and equal between the halves.

Though this asymmetry dominates both skewness maps, there are ‘limbs’ in the outer disc that do not follow the same trend. The western side of the Northern half in Fig. 7 has a predominantly blue-shifted component along most of the map edge, while the eastern side of the Southern half is mostly red-shifted. These features are nearly symmetric and are roughly pointed towards the position angle of the warped disc’s semimajor axis (~165°; Corbelli et al. 2014; Kam et al. 2017). The skewness ‘limbs’ may indicate where the warped disc becomes more prominent than the main disc. Indeed, the northern limb has larger line widths relative to the surrounding structure in Fig. 4. These enhanced line widths are also consistent with the line width map from Kam et al. (2017), which includes the structure at larger radii. Note that this does not contradict the lack of a warped disc component in the rotation model (Appendix C) since the peak velocities will not be sensitive to extended tails.

The kurtosis profiles in Fig. 8 do not show variation between the halves and are consistent within the uncertainties with the whole disc average at most radii. The positive kurtosis values indicate that the typical H₁ spectrum in the inner 6 kpc has strong tails that do not vary significantly in shape. The negative kurtosis beyond 7 kpc occurs where the brightest H₁ emission tapers off in the surface density profile (Fig. 6). A negative kurtosis value in these regions indicates a lack of narrow and bright H₁ components.

We note that the profile values in the inner ~0.5 kpc are more uncertain due to the low number of spectra in these bins and may be the cause of the negative skewness slope out to 1 kpc. If the spectra were preferentially red-shifted at the centre, the PV-slices in Section 5.1 would show this asymmetry.

### 4.3.1 Correlations with skewness and kurtosis

In Table 3, we show the Spearman correlation coefficients calculated by comparing either the kurtosis or skewness with other H₁ properties. Most lines of sight have H₁ detected along them, giving a population of ~1.2 million spectra to compare properties.<fnref n="5" /> Since the distribution of these parameter values needs not be Gaussian, we perform a bootstrap to calculate the uncertainty of the correlation and the *p*-value. In each bootstrap iteration, we randomly choose 30 per cent of the sample and re-calculate the correlations; the standard deviation of the resulting distribution from all iterations gives the uncertainty. The *p*-value calculation uses a randomly shuffled version of one of the parameters to calculate the correlation with the other parameter. Both steps are run for 1000 iterations, and all correlation values are found to be significant with small uncertainties.

Skewness is strongly correlated with $v_{cent}$ and $v_{peak}$ for the combined VLA and GBT data. This is well-explained by the north–south asymmetry in skewness evident in Figs 7 and 8.

The kurtosis shows strong correlation with the peak temperature, zeroth moment, and the spectral width of the signal mask. We test whether these correlations are driven only by the shape of the signal mask by recalculating the correlations using points with similar mask widths (e.g. the correlation between kurtosis and peak temperature only for points with a mask width of 300–320 pixels). We find that the correlations with peak temperature and the zeroth moment persist when controlling for the mask width, indicating the correlations are not solely dependent on the signal mask.

The correlations with kurtosis are instead the result of the line width bias due to extended line wings. As discussed in Appendix E, the presence of line wings in the profile leads to an overestimate of the line width. This results from having a larger number of

<footnote n="5">⁵We do not correct for small-scale correlations over the beam size since the map extent is significantly larger. This will not affect correlations on large scales.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Radial profiles of the skewness (top) and kurtosis (bottom) of the VLA + GBT data averaged over the entire disc (black), and the Northern (blue) and Southern (red) halves. The shaded region indicates the inner 0.5 kpc where beam smearing and a small number of samples in each bin gives large uncertainties. The skewness profile demonstrates the differing skewness properties between the halves of the galaxy. The excess kurtosis indicates that the average H I spectrum has heavier tails relative to a Gaussian, though the bias in the second moment requires caveats for interpreting the kurtosis (Section 4.3.1).</figcaption></figure>

channels associated with the line wings than the bright compo-
nent (e.g. spectrum d in Fig. 3), and the kurtosis suffers from the
same issue. The measured kurtosis is then with respect to the shape
of the line wings, rather than the shape of the bright component.
The positive kurtosis reflects the strong line wings present in most
spectra relative to a Gaussian shape. Since peak temperature and
zeroth moment are measures of how ‘peaked’ that bright compo-
nent is, they are positively correlated with the kurtosis. This high-
lights where the naive moments approach to examining complex
line shapes breaks down and that the results should be treated with
caution.

Neither skewness nor kurtosis shows a strong correlation with
the line width. This is encouraging since both depend on the line
width, with increasing power, and a strong correlation may imply
the variations in skewness and kurtosis are only driven by changes
in the line width. Our findings above also show that correlations
with kurtosis are not driven by the mask shape. These results show
that our measurement of significant non-zero values in the higher
order moments is robust and that single-Gaussian descriptions of
the H I line becomes inadequate at 80 pc resolution. However, since
moment-based descriptions are subject to biases in the centroid,

we further explore H I line profile shapes using stacked profiles
below.

4.4 Stacking spectra

Stacking spectra is a common procedure to enhance the signal-to-noise in regions of faint emission. By aligning spectra to a common central velocity, combining the spectra will coherently add emission components, while the noise will add incoherently. This technique has often been used for studying typical H I line properties in nearby galaxies.

We create stacked profiles using three different definitions of the line centre: the rotational velocity ($v_{\text{rot}}$; Section 4.1), the centroid velocity ($v_{\text{cent}}$), and the peak velocity ($v_{\text{peak}}$). From the resulting stacked profiles, we can assess how the reference velocity affects the stacked profile shape. A cursory test with Gauss-Hermite polynomials, as used in several previous works (e.g. Petric & Rupen 2007; Ianjamasimanana et al. 2012; Stilp et al. 2013a), gave poor results due to spectra with multiple bright H I components and those with extended line wings. The typical S/N in the data used in these studies prohibited directly estimating $v_{\text{peak}}$, which is well-defined in our higher S/N data (see also Braun 1997).

The different definitions of the line centre will highlight differ-
ent ISM properties. Previous studies have typically stacked spec-
tra based on $v_{\text{peak}}$, which will minimize the stacked profile width
and recover the typical spectrum properties at the resolution of
the data (Caldú-Primo & Schruba 2016). This is ideal for attempt-
ing to recover the thermal and turbulent properties of the medium.
Stacking with $v_{\text{rot}}$ yields stacked profiles with widths sensitive to
large-scale turbulence and non-circular motions. For example, vari-
ations from the circular rotation velocity are clear from the resid-
ual velocity map in Fig. 4. Using $v_{\text{cent}}$ as the common velocity
gives, for an optically thin tracer, the mass-weighted velocity av-
erage. As discussed in Section 4.3, this can make $v_{\text{cent}}$ biased if
the line profiles are asymmetrical. This bias in $v_{\text{cent}}$ makes a phys-
ical interpretation less clear. None the less, we explore these three
stacking methods since they are commonly used in spectral-line
studies.

We shift the spectra using Fourier transforms before co-adding them to create stacked profiles. Fourier shifting allows shifting by fractions of the channel width without interpolating, minimizing artificial broadening due to a finite-sampled grid. We then add the aligned spectra to create the stacked profile. The uncertainty in each channel of the stacked profile is the sum-in-quadrature of the uncertainty from each PPV-pixel. Since the noise is approximately constant in our data, this reduces to the noise of one PPV-pixel multiplied by the square root of the number of beam elements in the stacked region. We include all of the spectra in the stacked profiles where $v_{\text{cent}}$ and $v_{\text{peak}}$ can be calculated after applying the signal mask (i.e. the signal mask contains at least one component along the line of sight). This excludes only the small masked regions in Fig. 4.

To model the stacked spectra, we adopt the half-width-half-max (HWHM) scaling method from Stilp et al. (2013a). We favour this approach rather than other models due to the large variation in the shapes of the stacked profiles we find using different definitions of the line centre. In particular, we found that the two-Gaussian model used in other works (e.g. Young & Lo 1996; Ianjamasimanana et al. 2012) cannot account for how peaked the $v_{\text{peak}}$ stacked profiles are and the ratios between the two components are highly sensitive to covariances between parameters. We provide a more thorough comparison in Appendix D.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The peak H1 temperature (left) and the kurtosis (right) maps of the combined data. A positive kurtosis indicates a shape with strong tails and more peaked relative to a Gaussian. There is a strong correlation between the two maps that results from the influence of the line wings on the kurtosis (Section 4.3.1).</figcaption></figure>

<table><caption>**Table 3.** Correlation between skewness and kurtosis with other H1 properties. The uncertainties on each correlation are ±0.002, using bootstrap resampling of one-third of the data to recompute the correlation 1000 times. We also perform a permutation test with 1000 iterations. All correlations have a p-value of <0.001. The mask width is the number of pixels in the spectral dimension contained within the signal mask.</caption><thead><tr><th></th><th>Skewness</th><th>Kurtosis</th></tr></thead><tbody><tr><td>Mask width</td><td>0.019</td><td>0.408</td></tr><tr><td>Zeroth moment</td><td>0.065</td><td>0.682</td></tr><tr><td>Peak temperature</td><td>0.088</td><td>0.800</td></tr><tr><td>Centroid</td><td>-0.463</td><td>-0.155</td></tr><tr><td>Peak velocity</td><td>-0.516</td><td>-0.162</td></tr><tr><td>Line width</td><td>-0.033</td><td>-0.050</td></tr></tbody></table>

The HWM method assumes that the stacked profile within the FWHM of the peak intensity can be modelled as a Gaussian line profile. Since these stacked profiles are created over large areas, the S/N is extremely high and the FWHM is a reliable measure of the width. Additional properties describing the stacked profile shape are defined based on the residual between the stacked profile and the assumed Gaussian profile for the peak. We use the definitions for four of these parameters from Stilp et al. (2013a). The effective Gaussian width, $\sigma_{\text{HWM}}$, is estimated by scaling to a Gaussian shape within the FWHM of the stacked profile. This assumed Gaussian shape is centred at the velocity of the maximum intensity, $v_{\text{centre}}$. The line wings are defined as the excess relative to this Gaussian peak beyond the FWHM points. The fraction of excess flux is given as $f_{\text{wing}}$:

$$f_{\text{wings}} = \frac{\sum_{|v|>\text{HWM}} [S(v) - G(v)]}{\sum_{v} S(v)}, \quad (4)$$

where $S(v)$ is the stacked profile, $G(v)$ is the Gaussian model within the FWHM and $\text{HWM} = \text{FWHM}/2$. The sum in the denominator is over the entire profile. The line wing excess can be used to define

an effective width of the wings equivalent to the second moment:

$$\sigma_{\text{wing}}^2 = \frac{\sum_{|v|>\text{HWM}} [S(v) - G(v)] v^2}{\sum_{|v|>\text{HWM}} [S(v) - G(v)]}. \quad (5)$$

Note that the definition of the second moment relies on a Gaussian line profile. Since the residual $S(v) - G(v)$ is not Gaussian, this may not have a clear connection to a Gaussian line width. We include it here to compare with the values presented in Stilp et al. (2013a).

We define two additional parameters to describe the profile shapes. While higher order moments present one avenue for characterizing deviations from a Gaussian shape, we favour developing empirical measures that do not weight the statistic by deviations from the line centroid. These metrics provide a more uniform description of the line profile over its entire extent. First, we define the asymmetry of a profile as the difference between the total flux at velocities smaller than the reference velocity and the total flux at velocities greater than the reference velocity. The difference is then normalized by the total flux over all velocities:

$$a = \frac{\sum_{v>\text{centre}} S(v) - \sum_{v<\text{centre}} S(v)}{\sum S(v)}. \quad (6)$$

The sum in the denominator is over the whole profile. Stilp et al. (2013a) define a different asymmetry parameter that, while similar, returns an absolute value. We choose this alternate definition since it retains information on the side of the profile containing more flux, which is determined by the asymmetry in the line wings. In this manner, $a$ is analogous to the skewness, without weighting by the velocity offset and line width.

The final parameter, $\kappa$, provides the fractional difference between the central peak and the Gaussian model within the FWHM, similar to $f_{\text{wings}}$:

$$\kappa = \frac{\sum_{|v|<\text{FWHM}} [S(v) - G(v)]}{\sum_{|v|<\text{FWHM}} G(v)}. \quad (7)$$
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Stacked spectra over the entire disc using $v_{\text{rot}}$ (blue solid), $v_{\text{cent}}$ (green dashed), and $v_{\text{peak}}$ (red dot-dashed) as the definitions of the line centre. The peak intensities are set to unity to emphasize the line shapes and the velocity axis is centred at the systemic velocity. Stacking based on $v_{\text{rot}}$ and $v_{\text{cent}}$ give similar profiles with near Gaussian peaks within the FWHM. However, using $v_{\text{peak}}$ as the reference velocity yields a sharper peak shape and more prominent line wings relative to a Gaussian.</figcaption></figure>

This is another parametrization for the kurtosis: a negative value indicates the profile is more peaked than the Gaussian, while a positive value is flatter. We note, however, that the kurtosis used in Section 4.3 will primarily be sensitive to the line wing structure, whereas $\kappa$ describes the shape of the peak since it is only measured within the HWHM relative to a Gaussian of equivalent width.

Fig. 10 shows the stacked spectra over the entire disc with different reference velocities, and Table 4 gives the fitted model parameters. The parameter uncertainties are calculated in one of two ways. For $\sigma$ and $v_{\text{centre}}$, we assume that the uncertainty is set by the channel size; thus the uncertainty is half the channel width ($\sim 0.1$ km s$^{-1}$). For the rest of the parameters, we perform 100 bootstrap iterations where (1) the data are re-sampled by adding Gaussian noise within the flux uncertainty for each channel, and (2) allowing the estimated peak velocity and Gaussian width to vary by their uncertainties. This propagates the uncertainty in the flux and the assumed Gaussian peak to the other parameters. Due to the high S/N in the stacked profiles, most of the uncertainty comes from the latter source.

Fig. 10 shows how the use of different reference velocities significantly changes the stacked profile shapes. Stacking based on $v_{\text{rot}}$ and $v_{\text{cent}}$ gives the widest profiles and both have similar shapes. The shapes of their peaks are close to Gaussian ($\kappa \sim 0$; Table 4) and the widths are $\sim 10$ km s$^{-1}$ in both. The widths of these profiles are dominated by the residual velocity in the $v_{\text{rot}}$ model (Fig. 4) and the bias in $v_{\text{cent}}$ due to extended line wings evident in the skewness maps (Fig. 7). Stacking using $v_{\text{cent}}$ or $v_{\text{peak}}$ gives profiles with non-zero $a$ parameters (equation 6). This may be highlighting an asymmetry in the low surface-brightness component of the line wings not evident in the position-velocity (PV) slices in Section 5.

The stacked profile from using $v_{\text{peak}}$ is clearly non-Gaussian in its central shape ($\kappa = -0.059$). The profile is both narrower and more peaked relative to the other profiles and is similar to the shape found in other studies (e.g. Petric & Rupen 2007). The narrow and peaked profile contains a larger fraction of flux in the line wings relative to a Gaussian ($f_{\text{wings}} = 0.26$). The difference in this stacked profile compared to those from $v_{\text{cent}}$ and $v_{\text{rot}}$ is the source of broadening in the width, which we discuss further in Section 6.3. If the broadening source has a Gaussian distribution of velocities around the peak, like

the residuals from $v_{\text{rot}}$, the central-limit theorem gives that the shape will be close to Gaussian within the FWHM. When using $v_{\text{peak}}$, the dominant source of broadening is the uncertainty on $v_{\text{peak}}$ ($\sim 0.1$ km s$^{-1}$), which is small compared to typical H I line widths.

The width and shape of the $v_{\text{peak}}$ profile is relatively unchanged without the total power component. In Fig. 11, we show the $v_{\text{peak}}$ stacked profiles with and without the GBT data. Including the GBT data affects the line wings, which are far more prominent in the combined VLA + GBT data set. These differences highlight that most of the line wing emission is driven by a low surface brightness, extended component. This combination of a steep peak within the FWHM and heavy tails makes these shapes difficult to model. They are clearly non-Gaussian, but neither of these profiles has tails strong enough to be well-modelled by a Lorentzian.

The assumed Gaussian shape for sharply peaked profiles violates some of the assumptions for the HWHM model. The fractional difference in the peak shape, defined with $\kappa$, indicates that the assumed Gaussian shape overestimates the emission in the peak by 6 per cent for the $v_{\text{peak}}$ stacked profile. Though the two-Gaussian model also cannot account for the peak shape (Appendix D), this highlights the difficulty in defining a simple model for the stacked profile shapes. We discuss the factors that lead to this difficulty in Section 6.3.

### 4.4.1 Radial variations in stacked profiles

We next examine how the stacked profile parameters vary with galactocentric radius, focussing on the $v_{\text{peak}}$ stacked profiles. We create stacked profiles in galactocentric rings with 100 pc widths, as used for the surface density profiles (Section 4.2). Given the opposing line wing directions found in the skewness map (Section 4.3), we also examine stacked profiles over the Northern and Southern halves separately.

Fig. 12 shows the three profile parameters describing asymmetries ($a$) and shape ($\kappa$ and $f_{\text{wings}}$) for the $v_{\text{peak}}$ stacked profiles; we explore the line width variations in Section 4.5. The asymmetry parameter $a$ shows a clear difference between the two halves, matching the variations identified in the skewness analysis (Section 4.3). The Southern half of the galaxy has an excess of emission at negative velocities while the Northern half has an excess at positive velocities. This is qualitatively the same as is shown in the skewness radial profile in Fig. 8, though we find less variation with radius in $a$. The decreased radial variation in $a$ results from the $\sigma^3$ weighting in the skewness, which makes it sensitive to small line width variations, and from the difference in the order of operations, since stacked profiles are already an averaged quantity.

The shape within the FWHM of the stacked profiles, parametrized with $\kappa$, is largely consistent between the two halves. However, $\kappa$ varies by $\sim 33$ per cent with radius in the inner 6 kpc. There are two dips at 2.5 and 3.5 kpc where $\kappa$ increases, indicating a larger discrepancy from a Gaussian. There are similar features in the kurtosis profile in Fig. 8, however, those variations suggest that the typical spectrum becomes closer to a Gaussian profile. This apparent discrepancy may result from an increase in the fraction of multi-component spectra in these regions. These radii contain portions of the optical spiral arms, particularly the region around NGC 604 in the northern arm and the bright clump in the southern arm (see spectrum c in Fig. 3). The kurtosis is driven by the shape with respect to the line wings (Section 4.3.1) and multiple components will give a smaller excess kurtosis since the bright component is effectively wider. We see the opposite effect in $\kappa$ because the spectra
<table><caption>**Table 4.** HWHM fit values to the stacked profiles over the entire disc. These profiles are shown in Fig. 10. The uncertainties are propagated based on the uncertainty in each channel of the profiles, and assuming an uncertainty of half the spectral resolution (~0.1 km s⁻¹) for the peak velocity and FWHM line width. Note that σ<sub>wing</sub> may not be treated as an equivalent Gaussian width due to the non-Gaussian shape of the residuals; see equation (5).</caption><thead><tr><th></th><th>Rotation vel. (v<sub>rot</sub>)</th><th>Method of stacking spectra<br>Centroid vel. (v<sub>cent</sub>)</th><th>Peak vel. (v<sub>peak</sub>)</th></tr></thead><tbody><tr><td>σ<sub>HWHM</sub> (km s<sup>-1</sup>)</td><td>10.3 ± 0.1</td><td>9.7 ± 0.1</td><td>6.6 ± 0.1</td></tr><tr><td>v<sub>centre</sub> (km s<sup>-1</sup>)</td><td>-0.5 ± 0.1</td><td>-0.5 ± 0.1</td><td>0.0 ± 0.1</td></tr><tr><td>f<sub>wing</sub></td><td>0.11 ± 0.01</td><td>0.09 ± 0.01</td><td>0.26 ± 0.01</td></tr><tr><td>σ<sub>wing</sub> (km s<sup>-1</sup>)</td><td>28.9<sup>+0.5</sup><sub>-0.4</sub></td><td>26.3<sup>+1.5</sup><sub>-0.4</sub></td><td>23.4<sup>+0.3</sup><sub>-0.2</sub></td></tr><tr><td>a</td><td>-0.009<sup>+0.011</sup><sub>-0.010</sub></td><td>0.027<sup>+0.004</sup><sub>-0.003</sub></td><td>0.021<sup>+0.005</sup><sub>-0.014</sub></td></tr><tr><td>κ</td><td>-0.013<sup>+0.004</sup><sub>-0.004</sub></td><td>-0.011<sup>+0.004</sup><sub>-0.004</sub></td><td>-0.059<sup>+0.005</sup><sub>-0.004</sub></td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Stacked profiles of the VLA (top) and VLA + GBT data (bottom) when centred using v<sub>peak</sub>. The dashed green lines are the FWHM line widths and the shaded grey regions show the difference between the equivalent Gaussian profile and the data. The GBT data mostly add emission to the line wings and do not affect the shape within the FWHM.</figcaption></figure>

are added coherently and multiple components will tend to enhance the line wings.

The fraction of flux in the line wings is on average f<sub>wings</sub> ~ 0.22 in the inner 6 kpc, but there are substantial (~50 per cent) variations. There are peaks at radii of 2.5 and 3.5 kpc, consistent with variations in κ being driven more by lines of sight having multiple spectral components. However, peaks in f<sub>wings</sub> result from an increase in one half of the galaxy, whereas the variations in κ are consistent between both halves. The variations in f<sub>wings</sub> cannot then be entirely driven by an increase in multicomponent spectra. Another process, possibly outflows from star formation in the spiral arms, may lead to the increase in the line wing fraction.

None of the parameters show a monotonic trend with galactocentric radius in the inner 6 kpc, suggesting that the variations are due to local processes (i.e. outflows, local galactic structure).

A trend is clear, however, in all of the parameters beyond 6 kpc. The line profiles, over both halves, become more asymmetric (a < 0), have sharper and narrower peaks (κ and σ decrease), and have increased line wing fractions. A similar trend is seen for the kurtosis profile (Fig. 8), though the kurtosis decreases for the reasons discussed above. At ~6.5 kpc, the surface density profile strongly decreases to the edge of the map (Fig. 6). We attribute these trends

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Top: Asymmetry, a, of v<sub>peak</sub> stacked profiles in 100 pc radial bins. The parameter values are shown for stacked profiles over the whole disc (black), and the Northern (blue dashed) and Southern (red dot-dashed) halves. Middle: The peak shape κ. Bottom: The line wing fraction f<sub>wings</sub>. The shaded region in both plots highlights the inner 0.5 kpc, where the lack of data points makes the estimates more uncertain. The uncertainties in this region are underestimated since they do not account for the additional velocity uncertainty that results from spatial smearing where the rotation curve is steep (Stilp et al. 2013a). The top two panels are qualitatively similar to the skewness and kurtosis profiles from Fig. 8: the asymmetry varies between the halves of the galaxy, but the peak shape is consistent between both.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** The $v_{peak}$ stacked profiles over the whole galaxy (solid black), the Northern half (blue dashed), and the Southern half (red dot-dashed). The line wings asymmetry between the halves are highlighted with the shaded regions, where the colours correspond to the velocity the profiles are skewed towards.</figcaption></figure>

to the lack of a bright H I component when averaged on large scales. The change in asymmetry is due to the beginning of M33's warped disc, which starts to dominate the velocity surface at ~7 kpc, consistent with increases in the line width at the edge of the map (Fig. 4). The variations we find beginning at 6 kpc suggests the warped disc has a significant effect on the spectral shapes before it dominates the large-scale kinematics.

### 4.4.2 Asymmetry between the Northern and Southern halves

The skewness maps in Fig. 7 indicate a large-scale asymmetry in profile shapes between the Northern and Southern halves of the galaxy. By comparing the line wing excess in each half of the galaxy, we can estimate what fraction of $f_{wings}$ is symmetric and asymmetric in velocity.

The $v_{peak}$ stacked profiles in each half and over the entire disc are compared in Fig. 13. The peaks of all three profiles are essentially identical in shape and width, but the line wings differ in each half. We can re-write the line wing fraction (equation 4) as

$$f_{\text{wing}} = \frac{\sum_{|v|>\text{HWHM}} [S_{\text{N}}(v) + S_{\text{S}}(v) - G(v)]}{\sum_{v} S_{\text{N}}(v) + S_{\text{S}}(v)}, \quad (8)$$

where $S_{\text{N}}(v)$ and $S_{\text{S}}(v)$ are the stacked profiles of the Northern and Southern halves, respectively. The asymmetric line wing component is the total excess on opposite sides of the peak, located at $v = 0$ km s$^{-1}$ in this case. Based on the total profiles, the Southern half will have an excess towards negative velocities, while the Northern half will have an excess at positive velocities, both pointing towards the systemic velocity. Thus we define the asymmetric line wing component as

$$f_{\text{asym}} = \frac{\sum_{v<-\text{HWHM}} [S_{\text{S}}(v) - S_{\text{N}}(v)] + \sum_{v>\text{HWHM}} [S_{\text{N}}(v) - S_{\text{S}}(v)]}{\sum_{v} S_{\text{N}}(v) + S_{\text{S}}(v)}. \quad (9)$$

Note that the asymmetric component is not dependent on the assumed Gaussian peak $G(v)$, though the limits on the sums are set by the HWHM. The symmetric line wing component is the difference

<table><caption>**Table 5.** The line wing fractions from the peak velocity stacked profiles split into the Northern and Southern halves. The total line wing fraction is the same as given in Table 4.</caption><thead><tr><th></th><th>Line wing fraction</th></tr></thead><tbody><tr><td>$f_{\text{wing}}$</td><td>0.26 ± 0.01</td></tr><tr><td>$f_{\text{symm}}$</td><td>0.18 ± 0.01</td></tr><tr><td>$f_{\text{asym}}$</td><td>0.08 ± 0.01</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** The total (blue solid), symmetric (green dashed), and asymmetric (red dot-dashed) line wing fractions determined from equations (8) to (10) for $v_{peak}$ stacked profiles in 100 pc radial bins. The inner 0.5 kpc is shown in grey, where the small area and beam smearing cause deviations in the profile properties. The asymmetric component accounts for ~30 per cent of the line wing fraction and is roughly constant throughout the disc. The variations in the total line wing fraction are mostly driven by variation in the symmetric component.</figcaption></figure>

between the total and asymmetric components:

$$f_{\text{sym}} = f_{\text{wing}} - f_{\text{asym}}. \quad (10)$$

Because the asymmetric component is not defined as the absolute value of the difference, it is possible to have a negative value if the line wing excesses are opposite of the definition assumed based on Fig. 13.

Table 5 shows the asymmetric and symmetric components of the total profiles from Fig. 13. Relative to the total line wing fraction, we find that the ~70 per cent arises from a symmetric component while 30 per cent stems from an asymmetric one. In both halves, the asymmetric component points towards the systemic velocity of M33 (since $f_{asym} > 0$) and could trace a rotationally lagging disc component. We discuss this further in Section 5.

Fig. 14 shows the symmetric and asymmetric line wing fractions calculated from profiles in 100 pc radial bins. The asymmetric component is largely constant throughout the entire disc. The symmetric line wing component contains most of the variations previously seen in Fig. 12 for the total line wing fraction. In particular, the increase in $f_{wings}$ beyond 6 kpc is entirely due to the symmetric component.

### 4.4.3 Stacking profile with peak H I temperature

We also explore the properties of stacked profiles binned by their peak H I temperature. We created 20 stacked profiles in bins with a width of 5 per cent of the peak H I temperature distribution. This provides an equal number of spectra (48 900) in each stacked profile.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** The H I line width of profiles stacked based on their peak H I temperature. Each stacked profile covers a range of 5 per cent in the peak temperature distribution – the symbols in the plot indicate the middle of each bin.</figcaption></figure>

As in the previous section, we consider the profiles stacked with respect to $v_{\text{rot}}$, $v_{\text{cent}}$, and $v_{\text{peak}}$.

Fig. 15 shows the line widths for each stacking method in 5-percentile bins. Similar to radially binning the data (Section 4.5), we find an offset between the widths when using different definitions of the line centre, with the peak velocity stacking providing the smallest widths. Each stacking method has a line width that decreases with increasing peak temperature, though the change in the widths of the $v_{\text{peak}}$ stacked profiles is only ~1 km s⁻¹.

Despite the small changes we find in the $v_{\text{peak}}$ stacked profiles over the range in peak temperature, the shape parameters change substantially. The line wing fraction $f_{\text{wing}}$ monotonically decreases with increasing peak temperature, ranging from 0.17 to 0.32. The shape within the FWHM of the profile, described by $\kappa$, also increases with the peak temperature, with values ranging from −0.15 to −0.02. Note that this is the opposite behaviour from the radially binned profile parameters (Fig. 12). The combination of the trends in these two parameters and the near-constant line width demonstrates a shift in the profile shape that is near-Gaussian at higher peak temperatures and significantly non-Gaussian to low-peak temperatures.

The profile shapes at low-peak temperatures (<20 K) show two distinct components: a narrow sharply peaked component on top of a broad near-Gaussian component. Since, by area, most of the lines of sight with peak temperatures of <20 K are located at $R_{\text{gal}} > 7$ kpc, where the average surface density begins decreasing (Fig. 6), we interpret the increased importance of the broad component as arising from an increase in the proportion of emission from M33's warped disc component. This is consistent with our analysis of the skewness maps (Section 4.3), where the reversal in the sign of the skewness near the edges of the map in Fig. 7 is roughly aligned with the position angle of the warped disc component (Corbelli et al. 2014; Kam et al. 2017). The increase in $f_{\text{wings}}$ from the radially binned stacked profiles in Fig. 12 is also consistent with most spectra with low-peak temperature (<20 K) arising at $R_{\text{gal}} > 7$ kpc.

We also find no clear evidence for H I spectra with flattened peaks in the highest peak temperature bins as would result from optically thick H I emission (Braun et al. 2009; Braun 2012). If the stacked spectra had flattened tops, $\kappa$ should decrease with increasing peak temperature, opposite of the consistent increase we find.

<figure><img src="image_18.png" /><figcaption>**Figure 16.** H I line width profiles in 100 pc bins for the azimuthally averaged line width map (solid line; Fig. 4) and the stacked profiles from $v_{\text{rot}}$ (green diamonds), $v_{\text{cent}}$ (red circles), and $v_{\text{peak}}$ (purple triangles). The shaded grey region highlights the inner 0.5 kpc where the lack of samples and beam smearing cause large uncertainties in the line widths. The different methods give a large range in the H I line width. However, none of the methods shows a strong radial decrease.</figcaption></figure>

## 4.5 Comparing estimates of line width

Fig. 16 shows the line width as a function of galactocentric radius as derived from the azimuthally averaged second-moment line width map (Fig. 4) and the stacked spectra aligned with $v_{\text{rot}}$, $v_{\text{cent}}$, and $v_{\text{peak}}$. The enhanced line widths in the inner 0.5 kpc result from the small areas averaged over and from beam smearing where the rotation curve is steep.

Stacking by $v_{\text{peak}}$ gives the narrowest profiles, consistent with the whole-disc stacked spectra in Fig. 10. The typical width is ~7 km s⁻¹ in the inner 6 kpc with a shallow gradient; beyond 6 kpc, the line widths decrease to ~6 km s⁻¹. Our line widths are consistent with the stacked profile widths found by Druard et al. (2014) using the archival H I VLA observations presented in Gratier et al. (2010).

The $v_{\text{rot}}$ and $v_{\text{cent}}$ stacked profiles also have similar widths (~10 km s⁻¹) to the whole-disc stacked profiles. The widths are consistent with each other throughout most of the disc, suggesting that the line wing bias in $v_{\text{cent}}$ in the Northern and Southern halves of the disc are similar to the dispersion of velocity residuals from the $v_{\text{rot}}$ model. At radii from 2 to 4 kpc, where the spiral arms dominate the galactic structure, the $v_{\text{rot}}$ stacked profiles have larger widths than the $v_{\text{cent}}$ stacked profiles. This indicates enhanced motion in the spiral arms due to multicomponent spectra or coherent flows across the spiral arms, which are not included in the rotation model. Apart from the spiral arm region, there is no trend between $v_{\text{rot}}$ or $v_{\text{cent}}$ stacked line widths with radius.

We find that the azimuthally averaged line widths from the second moment have a typical value of ~12 km s⁻¹ and are larger than all of the stacked profile widths. Since stacking based on $v_{\text{rot}}$ yields the typical dispersion to non-circular motion and the $v_{\text{cent}}$ is biased by asymmetric line wings, it is worrying that the second moment estimates are substantially larger. Previous studies have shown that the second moment can overestimate the line width relative to stacked profiles (Ianjamasimanana et al. 2012) and Gaussian fits to individual spectra (Mogotsi et al. 2016), particularly in the inner discs of nearby galaxies. The discrepancy we find shows that, when extended line wings are detected throughout the disc, the second moment line widths overestimate the line width. In Appendix E, we find that the second moment line widths
increase by 50 per cent between the VLA and VLA+GBT data sets.

# 5 EXTRA-PLANAR VELOCITY COMPONENTS

Many nearby galaxies have velocity components offset from galactic rotation, either in the form of discrete HⅠ clouds or extended line wings related to the main disc structure (Sancisi et al. 2008). Both features relate to neutral gas accretion on to the main discs and may play an important role in galaxy evolution. We explore the preferential directions of line wings, as shown in the skewness maps (Fig. 7), in Section 5.1, and discrete HⅠ sources offset in velocity from the main disc in Section 5.2.

## 5.1 Anomalous velocity component

Fig. 17 shows a position–velocity (PV) slice along the major axis of the disc from the North to South. Similar to Kam et al. (2017), we find several features either unassociated with the M33’s main disc or lagging relative to the rotation velocity. The bright emission from the main disc, shown in grey-scale, is well-matched to the rotation model (green triangles). The contours highlight the extent of the low-surface brightness features. In particular, the 1 K contour, shown in blue, shows the asymmetry in the line wings directed towards the systemic velocity, consistent with what we find from individual spectra (Section 4.3). The extents of the blue contours in Fig. 17 suggests that the line wing structure is similar on large scales within the disc, implying that the line wings do not entirely arise from localized outflows from stellar feedback.

The blue contours in Fig. 17 highlight asymmetrical features that point towards the systemic velocity. This indicates a rotationally lagging disc component that is consistent with gas on high-latitude orbits in a thick disc. The extensive excesses near the centre of the disc may be related to the halo structure. Kam et al. (2017) find extensive ‘beard’-like structures – velocity elongations in the PV-slice – in this region, which are more evident at lower surface brightnesses.

The excess intensity at lagging velocities extends across the disc, as is also indicated by the skewness of the line wings (Fig. 7). Similar large-scale skewed profiles are found by Ianjamasimanana et al. (2012) for a sample of nearby galaxies in the THINGS survey (see their fig. 17). Their measure of skewness for individual profiles comes from Gauss–Hermite fitting, where higher order skewness and kurtosis terms can be included directly in the fit. They find a clear bias in the skewness between the approaching and receding halves for a number of galaxies in their sample.

Similar to the asymmetrical ‘limbs’ of opposing skewness along the edge of the map in Fig. 7, Ianjamasimanana et al. (2012) find several other systems with similar features towards the outer regions of galactic discs. In Section 4.3, we suggest that these limbs are roughly aligned with the position angle of M33’s warped disc component (e.g. Corbelli et al. 2014) and such features may arise from similar warped or lagging components in other nearby galaxies. For example, the direction of the asymmetrical regions found by Ianjamasimanana et al. (2012) in NGC 2403 is consistent with the position angle of its lagging rotational disc (Fraternali et al. 2001, 2002a).

Lagging disc components have been found in a number of nearby galaxies, including NGC 2403 (Fraternali et al. 2002a), NGC 4559 (Barbieri et al. 2005; Vargas et al. 2017), NGC 891 (Barnabè et al. 2006), and NGC 925 (Sancisi et al. 2008; Heald et al. 2011). As mentioned above, NGC 2403 provides an interesting comparison

to M33 because of its flocculent structure. Fraternali et al. (2002a) suggest the lagging component has an origin internal to NGC 2403, rather than being primordial, due to the coherent structure and connection to the main disc. NGC 2403 also appears to have been undisturbed by dynamical interactions (e.g. Williams et al. 2013). Fraternali et al. (2002a) suggest that the lagging component results from a galactic fountain mechanism, whereby hot ionized gas ejected into the halo from the fountain cools into atomic clouds that accumulate and fall back to the disc. Observations of ionized gas tracers confirm the presence of the lagging component (Fraternali et al. 2002b; Fraternali, Oosterloo & Sancisi 2004), consistent with the galactic fountain model. For a sample of nearby galaxies, Sancisi et al. (2008) also find that the extra-planar gas likely results from a galactic fountain, with a small component accreted from the intergalactic medium. Barnabè et al. (2006) find that the lagging component in NGC 891 is well-modelled by a baroclinic pressure model dependent on both the density and temperature of the medium. The model treats the medium as a mixture of a hot homogeneous component surrounding discrete cold HⅠ clouds. From this, they reproduce the observed velocity gradient of the lagging component in NGC 891.

A lagging disc component from the galactic fountain mechanism may be related to cold HⅠ ‘low-velocity clouds’ (LVCs), which have been observed in the Milky Way. Stanimirović et al. (2006) find a population of small (~ few pc) LVCs with velocities of <30 km s⁻¹ relative to the Milky Way rotation velocity (see also Lockman et al. 2002; Stil et al. 2006). Adopting their typical column densities of clouds (>2 × 10¹⁹ cm⁻²), we compare this to the average column density of material at v ~ 30 kms⁻¹ off M33 projected rotation speed (<2 × 10¹⁸ cm¹⁸). If the off velocity gas were organized into such clouds, it would have a reasonable filling factor of f < 0.1 and be broadly distributed across the disc. Our observations of this gas are consistent with an unresolved population of small clouds at high altitudes above the disc. The material could be raining back down as per a galactic fountain model, or the material could be directly accreting from the hot medium.

Only 1/3 of the line wing emission (~8 per cent of the total atomic gas mass) is found in the asymmetric wing components, which would be identified as the ‘lagging’ gas. The remaining 2/3 of the line wing emission in the symmetric component may require other explanations, especially since the wings are significantly wider than the thermal width of the warm neutral gas. One possible explanation is neutral gas accretion from the halo. Recently, Zheng et al. (2017a, b) favour a model where $M \sim 10^8 M_\odot$ of ionized gas in M33’s halo is accreting on to the disc at ~100 kms s⁻¹. They suggest that these features result from the infall of material produced by a galactic fountain or past interactions between M33 and M31. A key observation in this analysis is the split between the symmetric and asymmetric components of the line wings (Section 4.4.2). Our symmetric line wings are consistent with this accretion flow.

## 5.2 Discrete HⅠ clouds

Substantial HⅠ structure has been found within a few degrees of M33 (Grossi et al. 2008; Putman et al. 2009; Keenan et al. 2016), including discrete features in an HⅠ bridge that stretches towards M31 (Lockman et al. 2012). These structures are thought to comprise a population of HVCs or halo gas surrounding M33, similar to the HVC population in the Milky Way (Sancisi et al. 2008). Studies of these structures require excellent surface-brightness sensitivity and have mostly been limited to single-dish surveys
<figure><img src="image_19.png" /><figcaption>**Figure 17.** A position–velocity slice spatially averaged over a 200 arcsec region centred along the major axis. The *x*-axis is the offset angle measured from the rotation centre of the disc (Table 2). The grey-scale, indicated by the colourbar, is in K, while the contours are at 1, 2, and 3 K, respectively. The green triangles are the velocities from the rotation model (Appendix C). The horizontal dot–dashed line is at the systemic velocity, and the vertical dot–dashed line indicates the centre of M33. The region near the centre shows significant emission off the galactic rotation curve. The lowest contour highlights two features near the systemic velocity that are independent with emission from the main disc down to the noise level. Across the slice, the contours show faint emission that is preferentially skewed towards the systemic velocity (Table 2).</figcaption></figure>

at low resolution. Kam et al. (2017) provides the first-resolved study of these features nearest to the main disc, including elon- gated velocity structures connected to M33’s disc and discrete H I clouds offset in velocity. The data used by Kam et al. (2017) have a resolution of 120 arcsec, much better than previous deep H I studies.

We search for similar discrete H I structures within the spatial limits of our map. Though most of the H I emission in the aforemen- tioned surveys is outside this area, we find significant extra-planar structure within our data. Fig. 18 shows the extra-planar emission with integrated intensity maps in the Northern and Southern halves after subtracting the circular rotation model from the data cube. While most of the previously known H I clouds from single-dish studies fall outside of the area mapped our survey area, we find two structures that are consistent with the catalogue from Keenan et al. (2016).

We find a total atomic mass of $1.3 \pm 0.5 \times 10^7 M_{\odot}$ in the velocity-integrated regions shown in Fig. 18, which accounts for about 20 per cent of the total HVC and disc-halo mass in M33 (Lockman 2017), and $\sim 1$ per cent of the atomic mass of the main and lagging discs. This mass is a factor of $\sim 3$ less than the estimated HVC mass in the Milky Way (Putman, Peek & Joung 2012).

The left-hand panel of Fig. 18 shows the extent of the blue-shifted components and the right-hand panel shows the red-shifted com- ponents. The atomic mass over these regions is $5 \pm 2 \times 10^6$ and $8 \pm 3 \times 10^6 M_{\odot}$ for the blue- and red-shifted components, respec- tively. In both cases, the extra-planar H I emission has significant small-scale structure, with features that were unresolved in previ- ous studies of extra-planar gas in M33 (Westmeier, Braun & Thilker 2005; Grossi et al. 2008; Putman et al. 2009; Keenan et al. 2016). At 80 pc resolution, the structures are filamentary with ‘blob’-like concentrations. This is particularly prominent for the blue-shifted components, which appear highly elongated.

The red-shifted integrated intensity shows four structures, one of which is partially cut-off by the edge of the map. The structure nearest to the centre of M33 is a portion of the lagging rotation component highlighted in Fig. 17, shown by the deviation in the green contour at an offset of $-0.05^\circ$. The last two structures near the top of the disc ($1^{\prime}34^{\prime}00^\prime s +31^{\prime}05^\prime s$) have been catalogued by Grossi et al. (2008) as AA1 and Keenan et al. (2016) as AGESM33-1. These structures appear blended at the resolutions of both studies,

but are shown here to have multiple faint filamentary structures surrounding the brighter centre. Keenan et al. (2016) report a mass of $1.18 \times 10^6 M_{\odot}$ for the blended structure, about 1/8 of the total mass we find over this velocity range.

A second structure identified as an HVC by Keenan et al. (2016), AGESM33-22, is within the spatial extent of our data but is too faint to be detected in our data at 80 pc resolution. The HVC, located at a velocity of $-338 \pm 15 \text{ km s}^{-1}$, has a column density of $7.2 \times 10^{17} \text{ cm}^{-2}$ at a resolution of 3.5 arcmin, well below the sensitivity of our full resolution data.

The structure of the blue-shifted component is dominated by a single filamentary feature that stretches across most of the disc. Thus its mass is that of the blue-shifted region: $5 \pm 2 \times 10^6 M_{\odot}$. The portion nearest to M33’s centre is comprised of two connected fila- mentary structures, forming a loop, that then extends further south. This feature has not been discussed in previous studies, though it is noted in fig. 7 of Sancisi et al. (2008) using the Thilker et al. (2002) archival VLA data. It can also be seen in fig. 3 of Putman et al. (2009) in their Arecibo data at velocities of $-159$ to $-123 \text{ km s}^{-1}$. The Arecibo data indicate that the two components of the filament are indeed connected. The filament then has a projected length of $\sim 8$ kpc. Similarly, the velocity of the northern tip of the filamentary structure approaches the velocity of the lagging rota- tional component and the extended structure between 0 and $0.2^\circ$ in Fig. 17. Filamentary H I structures with similar lengths have been found in NGC 2403 and NGC 891 (Sancisi et al. 2008). Fraternali et al. (2001, 2002a) find a mass of $1 \times 10^7 M_{\odot}$ for the filament in NGC 2403, about twice the total mass we find in the blue-shifted off-axis emission in M33. Structures like this H I filament in M33 are strong evidence for accretion on to the main disc (Sancisi et al. 2008).

## 5.3 An H I cloud impacting the main disc?

An interesting H I cloud near the northeast edge of the map overlaps in velocity with the main disc, possibly indicating an interaction between an infalling cloud and the disc. We show the spatial struc- ture of this cloud in the top panel of Fig. 19. The morphology of the cloud is similar to the larger scale clouds in Section 5.2, with a central blob-like structure with connecting filamentary structure. The prominent filamentary feature in the figure reaches towards the
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Integrated intensity maps of components displaced in velocity from M33's main disc in the Southern half (left) and the Northern half (right). The maps are integrated over the rotation-subtracted data cube; a velocity of $v = 0$ km s$^{-1}$ is set at the systemic velocity of M33. The emission from the main disc is nearly all confined to velocities of -60 to +50 km s$^{-1}$ (see Section 4.4). The left shows the emission structure blue-shifted from M33 (Southern half), while the right shows the red-shifted components (Northern half).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Top: An integrated intensity map highlighting the structure of the H I cloud. The red dashed box in Fig. 3 shows this position relative to the rest of the galaxy. The velocity range is chosen to minimize emission from the main disc. Bottom: The average spectrum within the red-dashed circle in the top panel. The velocity axis is defined with respect to the $v_{\text{rot}}$ model from Section 4.1. We fit two Gaussian components to the spectrum, shown as the green-dashed and red dot-dashed components.</figcaption></figure>

northern edge of the map and overlaps in velocity with the main blob. The blob is about 600 pc in diameter.

The bottom panel of Fig. 19 shows the average spectrum over the cloud, where the velocity axis is defined with respect to the circular rotation velocity at $v = 0$. There are clearly two spectral components and possibly a third faint component around 60 km s$^{-1}$. We model the two bright components with Gaussians. The widths of the fitted components are $15.2 \pm 0.5$ km s$^{-1}$ for the H I cloud and $9.5 \pm 0.2$ km s$^{-1}$ for the disc.

Putman et al. (2009) show that this cloud overlaps in velocity with a portion of M33’s warped disc, labelled as the ‘northern warp’ in their fig. 3, outside of our survey area. What we detect as a separate H I cloud may be part of a larger structure related to the warped disc.

Treating the H I cloud as a discrete structure, we use the cloud’s spectral component from Fig. 19 with the total emission over the cloud to derive a mass of $6.5 \pm 0.5 \times 10^6$ M$_{\odot}$. This is similar to the total mass found in the blue- and red-shifted H I clouds in Section 5.2. With this cloud included, the total mass of off-rotation components becomes $\sim 2 \times 10^7$ M$_{\odot}$ and there are likely other fainter off-rotation H I clouds blended with the lagging rotational disc that will also increase this total mass. A full search of these features will require a more thorough spectral decomposition of the data.

# 6 DISCUSSION

Measuring the kinematic properties of the atomic ISM is vital for evaluating and comparing to theoretical models. We compare two common techniques for estimating kinematic properties in extra-galactic observations: the second moment of the line shape and spectral-stacking. Our results demonstrate that these common techniques are inherently limited and can yield vastly different results. For example, the method used to estimate the line width can yield values that differ by a factor of two (Section 4.5). In this discussion, we compare our results to those found in nearby galaxies and present a simple Gaussian Mixture Model for describing stacked profile shapes. This model recovers the qualitative profile shape found in previous studies, which we use to compare with previous model interpretations of stacked H I profiles.
## 6.1 Previous work and interpretations

Previous work on extragalactic H₁ spectra focuses on the H₁ line width, measured using the second moment (Tamburro et al. 2009), fitting Gaussian models to individual spectra (Petric & Rupen 2007; Warren et al. 2012; Mogotsi et al. 2016), or stacking spectra over large spatial regions (Young & Lo 1996; Braun 1997; Ianjamasimanana et al. 2012; Stilp et al. 2013a). The first method is attractive since second moments are simple to compute, though it assumes the profile is close to a Gaussian. Tamburro et al. (2009) use line widths from the second moment for the THINGS survey and find typical values of $10 \pm 2 \text{ km s}^{-1}$ across the sample. Since these line widths are larger than the thermal line width expected for the WNM ($\sim 8 \text{ km s}^{-1}$), the profiles are interpreted as being broadened due to turbulent motions.

Modelling individual spectra by fitting Gaussian components has the advantage of measuring variations down to the resolution of the data, and being able to account for multiple components. Young et al. (2003) and Warren et al. (2012) find that most H₁ spectra can be modelled by up to two Gaussian components. Both studies find a range of line widths in both components, suggesting emission components arise from both CNM and WNM phases. The narrow components have widths of $<6 \text{ km s}^{-1}$ and are attributed to the CNM, with the line width broadened substantially from the thermal width by turbulence. However, Warren et al. (2012) do not find a clear connection between the location of narrow components and star-forming regions. Using a similar data set as Tamburro et al. (2009), Mogotsi et al. (2016) fit a single Gaussian model to all H₁ spectra where CO is detected. They find line widths of $10 \pm 2 \text{ km s}^{-1}$, consistent with estimates from the second moment.

Alternatively, non-Gaussian line profiles have been modelled using Gauss-Hermite polynomials, which can account for variations in the skewness and kurtosis. While this method often provides a good analytical description, the connections between the model parameters and the underlying physical properties of the medium are unclear (Young et al. 2003). Young & Lo (1996), Young et al. (2003), Ianjamasimanana et al. (2012), and Stilp et al. (2013a) use this model to estimate the velocity at the peak intensity in their data, which are then shifted to a common centre and stacked.

The H₁ stacking method has the advantage of vastly increasing the S/N in the data, allowing low surface brightness emission to be studied at the expense of removing spatial variations. Nearly all studies find the same basic stacked profile shape: a near-Gaussian peak with enhanced line wings. This shape is not well-described by physically motivated profile shapes. The enhanced line wings lead to poor single Gaussian fits, yet they are not prominent enough to match a Lorentzian profile. This has led to different modelling approaches with differing interpretations of the results. Young & Lo (1996) model the stacked profiles as a wide and narrow Gaussian component. These two components have widths of $9 \text{ km s}^{-1}$ and $3 \text{ km s}^{-1}$, respectively. Since these are similar, though slightly larger, than the expected thermal line widths for WNM and CNM phases, they attribute the model as separating out the two different phases, with some broadening from turbulence. A larger sample of 34 galaxies from the THINGS survey was analysed by Ianjamasimanana et al. (2012) using the same two-Gaussian model; the sample overlaps with the independent analyses by Tamburro et al. (2009) and Mogotsi et al. (2016). They find similar results to Young & Lo (1996), with narrow components arising from the CNM with Gaussian widths of $3.4–8.6 \text{ km s}^{-1}$, and wide WNM components with widths from $10.1–24.3 \text{ km s}^{-1}$. These ranges imply that both the CNM and WNM must be highly turbulent

since the line widths are much larger than the expected thermal line widths of these phases. We note that, for comparison with our results, the data from the THINGS survey has a coarser spectral resolution compared to our data, with channel sizes from 1.3 to $5.2 \text{ km s}^{-1}$.

Petric & Rupen (2007) propose an alternate interpretation of the enhanced line wings in stacked profiles. They argue that the flux ratio of the narrow component, if exclusively tracing the CNM, to the wide component should steeply decline with galactocentric radius, matching the observed star formation activity. Instead, they find that the ratio stays roughly constant to large radii. They attribute the central component to a mixture of atomic gas, with enhanced line wings that describe bulk motions, possibly from stellar feedback or infalling halo gas. Stilp et al. (2013a) test this interpretation by modelling the stacked profiles as a single Gaussian peak with line wings measured by their fractional increase relative to the Gaussian. For a sample of 24 dwarf galaxies, they find line widths ranging from $\sim 6$ to $10 \text{ km s}^{-1}$ and $f_{\text{wings}} \sim 0.05–0.2$.

An alternative model to the previous studies is introduced by Braun et al. (2009) and Braun (2012) where spectra with flattened tops are attributed to a cool H₁ layer sandwiched between two warmer layers (Braun 1997). It is difficult to make comparisons to this model without accounting for multicomponent spectra when calculating the line width. We defer further discussion of this model to future work.

Most methods have also been used to study radial trends in the H₁ properties. Tamburro et al. (2009) examine radial variations in the H₁ line width, finding that the line width tends to decrease rapidly in $R_{\text{gal}}$ to $R_{25}$,<fnref n="6" /> similar to the radial decline in the star formation rate surface density. The slope beyond this radius in some galaxies flattens significantly. Stilp et al. (2013b) and Ianjamasimanana et al. (2015) examine radial trends by creating stacked profiles in radial bins. They find the same qualitative result despite using different models: the line width tends to decrease with radius. This trend is found in most nearby galaxies, though with three exceptions: NGC 2403, NGC 2976, and NGC 628 (Ianjamasimanana et al. 2015; Mogotsi et al. 2016). These systems have a shallow radial decrease in the line width throughout the disc. The first two galaxies are flocculent spirals similar to M33. NGC 628 is classified as a grand spiral galaxy, though its mass and star formation rate is within a factor of 2 compared to M33 (Walter et al. 2008).

These studies cover the H₁ properties of most nearby galaxies, spanning a range of galaxy properties and types. Young & Lo (1996), Young et al. (2003), and Stilp et al. (2013a) restrict their sample to dwarf galaxies, where the low-metallicity properties and lack of large-scale galactic structure may lead to differences in the ISM properties compared to the canonical view of the ISM. Tamburro et al. (2009), Ianjamasimanana et al. (2012), and Mogotsi et al. (2016) include a mix of galaxy types in their sample. Ianjamasimanana et al. (2012) find no difference in the H₁ line widths of stacked profiles between dwarf and spiral galaxies. Furthermore, Caldu-Primo et al. (2013) analyse the radial line width trends in the H₁ and CO for a similar set of nearby galaxies and find no common radial variation. These results are consistent with all stacked profiles in these studies having the same qualitative shape (Petric & Rupen 2007), suggesting that the physical mechanism governing typical line profiles is similar regardless of galaxy type.

<footnote n="6">⁶Defined as the semimajor axis of the 25 mag arcsec⁻² isophote in B band.</footnote>
Druard et al. (2014) present an analysis of H I stacked profiles in M33, using the archival VLA H I data presented in Gratier et al. (2010). In this work, they emphasize comparisons of the H I to stacked profiles of CO(2–1). In 1 kpc radial bins, they find an average H I stacked profile width of ~6.5 km s⁻¹ by fitting a single Gaussian to the stacked profiles. The line width decreases slowly with radius – a decrease of just ~2 km s⁻¹ between the bins at R<sub>gal</sub> = 0–1 kpc and the 7–8 kpc bin – similar to the three galaxies noted above.

## 6.2 Key results in M33

Our results are broadly consistent with previous studies. In particular, we find a line width of 6.7 ± 0.2 km s⁻¹ from the v<sub>peak</sub> stacked profile, consistent with the M33 study from Druard et al. (2014). The line widths of the stacked profiles over radial bins, ranging largely between 5.5 and 7.5 km s⁻¹, are also consistent with Druard et al. (2014). This is encouraging since: (1) the stacked profiles of the archival H I data should have a similarly high S/N as our data, with the only major difference being the lower spectral resolution of the archival data (1.2 km s⁻¹), and (2) a different methodology is used in parametrizing the line shapes. The stacked profiles between the two studies also have a consistent shape, with a kurtosis excess in the central peak and prominent line wings. The kurtosis excess is similar to that found by Braun (1997) and more extreme than other recent studies (e.g. Ianjamasimanana et al. 2012; Stilp et al. 2013a).

The line widths from the v<sub>peak</sub> stacked profiles are near the lower limit of the range of line widths found by Ianjamasimanana et al. (2012) and Stilp et al. (2013a) for nearby galaxies (see Appendix D). However, since these profiles have non-Gaussian shapes, we note that the equivalent Gaussian width should be treated as an upper limit on the line width. While the HWHM-scaling method from Stilp et al. (2013a) is appealing because it does not require fitting an analytical model, it still has limitations on how well the profile shape can be described. We discuss the source of these difficulties in Section 6.3.

There are some discrepancies between the range of values we find for the other HWHM model parameters compared to the set of dwarf galaxies in Stilp et al. (2013a, b) :

(i) We find consistently larger f<sub>wing</sub> values, ranging from 0.25 to 0.35 in radial bins, than the 0.1–0.25 range in the Stilp et al. (2013b) sample (Fig. 12). This range is, however, consistent if we only consider the symmetric line wing fraction (Fig. 14). Thus the symmetric line wing component, possibly arising from inflowing or outflowing material (Section 5), is similar with those measured in nearby dwarfs and the difference is due to the presence of a lagging rotational disc in M33.

(ii) Stilp et al. (2013b) find a large range of asymmetries when stacking spectra in radial bins.<fnref n="7" /> We only find large asymmetries in M33 when stacking over radial bins in the Northern or Southern half, respectively (Fig. 12). The shape within the HWHM of the stacked profiles in Stilp et al. (2013b) appear close to symmetrical (see their fig. 4), and most of the asymmetry occurs in the line wings.

The estimated line widths from the second moment have an average of ~12 km s⁻¹, consistent with the upper limit of the range

found by Tamburro et al. (2009). These widths are 50 per cent larger than largest widths of the stacked profiles (Fig. 16). In Appendix E, we demonstrate that the second moment is sensitive to the line shape, and can drastically alter the second moment line width when extended line wings are present. Thus, we discourage the use of the second moment for estimating line widths. Similar differences between other line width methods and the second moment are noted in several studies (Ianjamasimanana et al. 2012; Stilp et al. 2013a; Mogotsi et al. 2016).

## 6.3 Mixture models for stacked profiles

As we explain in Section 6.1, the shape of stacked profiles has led to different modelling approaches and different interpretations of their meanings. To illustrate why stacked profiles have shapes that are difficult to model using simple analytic forms, we invoke a Gaussian mixture model as a framework for understanding the line wings and widths. We consider the case where a stacked profile S(v) is an average over an ensemble of Gaussian line profiles with purely thermal line widths:

$$S(v) = \sum_{i} w_{i} \exp \left[ -\frac{(v - v_{i})^{2}}{2\sigma_{T,i}^{2}} \right], \qquad (11)$$

where w<sub>i</sub> is the weight of each spectrum, v<sub>i</sub> is the central velocity, and σ<sub>T</sub> is the thermal line width for temperature T, i.e. σ<sub>T</sub> = √kT/μm<sub>H</sub>. The weights are normalized such that Σ<sub>i</sub>w<sub>i</sub> = 1. We can describe the shape of the stacked profile S(v) based on the moments of the distribution. We use the derivation from Wang & Taaffe (2015) to describe the first four moments. The centroid (first moment) of the stacked profile is v<sub>0</sub> = Σ<sub>i</sub>w<sub>i</sub>v<sub>i</sub>, which depends only on the distribution of the line centre v<sub>i</sub>. The variance (second moment) of the stacked spectrum is then

$$\sigma^2 = \sum_i w_i (\sigma_{T,i}^2 + v_i^2) - v_0^2 \qquad (12)$$

and the skewness (third moment) and kurtosis (fourth moment) are respectively

$$\text{skew} = \frac{1}{\sigma^3} \sum_i w_i (v_i - v_0) [3\sigma_{T,i}^2 + (v_i - v_0)^2] \qquad (13)$$

$$\text{kurt} = \frac{1}{\sigma^4} \sum_i w_i [3\sigma_{T,i}^4 + 6(v_i - v_0)^2\sigma_{T,i}^2 + (v_i - v_0)^4] - 3. \quad (14)$$

The critical result of these moments is that the stacked profile will depend on the distributions of the line centre and width of the individual spectra that were averaged over. For example, the shape of the stacked profile could be reproduced by adopting a reference probability density function for w<sub>i</sub> = f(T, v<sub>i</sub>) encoding the amount of material found at a given temperature and velocity offset from the reference velocity.

We investigate how the common features in stacked profiles found in different studies can be attributed to the dependence on the distributions of the line centre and width of individual spectra that comprise the stacked profile.

(i) Velocity offsets contribute to the line width. – Equation (12) has an equal dependence on the line width of the spectra and the distribution of their central velocities. The differences between the v<sub>rot</sub>, v<sub>cent</sub>, and v<sub>peak</sub> stacked profile widths is due to this dependence (Fig. 16). Stacking based on v<sub>peak</sub> minimizes the variation of the central velocity to within the uncertainty of v<sub>peak</sub> (0.1 km s⁻¹). It

<footnote n="7">⁷Note that the asymmetry parameter in their study returns the fraction of asymmetry, bounded between 0 and 1: a<sub>full</sub> = $\frac{\sum_{i} \sqrt{(S(v) - S(-v))^2}}{\sum_{i} S(v)}$.</footnote>
represents the optimal recovery of the mixture of the line widths. The $v_{cent}$ and $v_{rot}$ stacked profiles are instead dominated by the distribution of the central velocities, not the distribution of the line widths. The purpose of the different approaches is discussed further in Section 6.6.

(ii) Skewness depends on the distribution of line centres. – From equation (13), an asymmetric stacked profile can only occur if the distribution of line centres is also asymmetric since the sign of the skewness is determined by the ($v_i - v_0$) term. Stacked profiles using $v_{rot}$ or $v_{cent}$ are both close to symmetric about $v = 0$, but $v_{cent}$ can be biased by the asymmetric line wings. To get asymmetric stacked profiles from the former two line centres, an additional component (i.e. the lagging disc) must be present in the line profiles.

(iii) Stacked profiles of the atomic medium almost always have a line wing excess. – The kurtosis from equation (14) has a strong dependence on the distribution of line centres and the line widths. Since the terms all have even powers, non-Gaussian features in either distribution will always increase the kurtosis, causing the mixture model to have an excess in the tails. For the atomic medium, the distribution of temperatures, which only affect $\sigma_{T,i}$ will be non-Gaussian since it will contain a mixture of CNM and WNM temperatures, with some fraction in the unstable intermediate regime. Thus, without requiring a separate physical process to produce line wings, stacked profiles of the atomic medium will always have an excess in the line wings. The only case where the kurtosis will decrease is if the majority of the spectra have multiple line components. This case is not explicitly handled in this toy model, but would tend to broaden the profile and flatten the peak, resulting in a kurtosis deficit.

Without requiring additional physical processes, this mixture model can qualitatively reproduce properties of observed stacked profiles based only on a set of Gaussian components. This accounts for the near universal stacked line shape noted in Petric & Rupen (2007) and Stilp et al. (2013a) when stacking based on the peak velocity.

By adopting simple distributions for the line centres and widths, we test whether the mixture model of thermal-Gaussian components can produce simulated stacked profiles with the quantitative properties of observed stacked profiles. We find that reasonable distributions of H I temperatures creates narrower stacked profiles with smaller line wing fractions than the observed stacked profiles. We simulate a stacked profile from the mixture model by sampling $10^4$ spectra from a uniform temperature distribution ranging from $10^2$ to $10^4$ K. We also sample the central velocities from a Gaussian distribution with width $\sigma = 1 \text{ km s}^{-1}$ to compare with the $v_{peak}$ stacking in our data. The resulting stacked profile is similar in shape to the $v_{peak}$ stacked profiles, with a sharp peak and enhanced line wings. Using the same model used for the stacked profiles, the mixture model profile has a width of $\sigma = 4.8 \text{ km s}^{-1}$, $\kappa = -0.03$, and $f_{wings} = 0.09$. When using a wider central velocity distribution to compare with the $v_{rot}$ stacked profiles, we find that the mixture model also has a narrower width and smaller $f_{wings}$.

We create another simulated stacked profile, but alter the temperature distribution to range from $10^3$ to $10^4$ K. This increases the width of the profile to $\sigma = 5.2 \text{ km s}^{-1}$, but makes the shape closer to a Gaussian than the observed profiles, with $\kappa = -0.01$ and $f_{wings} = 0.06$. This suggests that a cold atomic component in the stacked profiles is responsible for the high kurtosis excess in observed stacked profiles.

The mixture model demonstrates that the qualitative stacked profile shape can be reproduced from a mixture of Gaussian compo-

nents and a range of plausible thermal temperatures. This highlights the difficulty in using simple analytical models for describing stacked profiles. We find that observed stacked profiles are wider and have more prominent line wings than the simplified mixture model presented here (Fig. 3). These differences must arise from a combination of the CNM/intermediate/WNM fractions, the lagging rotational disc, multicomponent spectra, and turbulence. We conclude that stacked profiles reflect this full set of physical effects in M33, rather than being solely the distribution of CNM/WNM temperatures. In particular, there must be significant contributions from a multimodal distribution of line centres (e.g. from multiple components) to explain the profiles. Because our data can resolve multiple velocity components both spatially and spectrally, we will revisit this analysis in future work through the multicomponent fitting to our data, mirroring the approaches used in Galactic studies (e.g. Lindner et al. 2015; Henshaw et al. 2016).

## 6.4 Interpreting line wings

The various models used for H I stacked profiles are driven by different interpretations of the line wings. The two-Gaussian models used by Ianjamasimanana et al. (2012) argue that the line wings trace a turbulent WNM component, while Stilp et al. (2013a) assign the wings as the product of feedback. Our findings suggest that line widths result from a mixture of different physical processes, making it difficult to uniquely identify the source of the line wing excess without a detailed knowledge of the processes that can broaden the line profiles. The mixture model presented in Section 6.3 shows that the excess line wings naturally result from combining a set of Gaussian components. Since both turbulent and thermal motion produce Gaussian profiles, the effect of each component cannot be determined from the line wing fraction without additional information.

Both in individual spectra and stacked profiles (over the approaching or receding halves of the galaxy), we find that line wings are not symmetric and are consistently skewed towards the systemic velocity (Fig. 8). Using the $v_{peak}$ stacked profiles in the Northern and Southern halves, we find that about 1/3 of the line wing fraction results from this asymmetric component (Fig. 14). Based on our comparison of off-rotation and extra-planar structure in Section 5, we identify this component as arising from the lagging rotational disc. Fig. 14 shows that the lagging rotational component has roughly an equal contribution to the line wing fraction throughout the disc.

The symmetric component of the line wing fraction, however, shows substantial variation with galactocentric radius, though there is no consistent trend. Instead, the variations appear to result from local processes within the disc. For example, there is an increase in the mid-disc ($R_{gal} = 2–4 \text{ kpc}$) where the spiral arms dominate, possibly driven by enhanced recent star formation in the region and bulk non-circular motions driven by the arm mechanism.

The radial trend is inconsistent with the symmetric line wing component being driven solely by stellar feedback. Line wings resulting from stellar feedback should correlate well with tracers of the star formation rate, which drops exponentially with radius in M33 (Heyer et al. 2004; Boquien et al. 2015), as does the H$_2$ surface density (Druard et al. 2014). Stilp et al. (2013a) find a correlation between the SFR and line wing fraction for stacked profiles over the full galaxy of their sample of dwarf galaxies. However, the correlation is tenuous when the stacking is performed in radial bins (Stilp et al. 2013b). Many of the galaxies in their sample show relatively constant $f_{wings}$, or otherwise peak far from
the centre of the disc, as we find in M33. The lack of a radial trend may also be inconsistent with the turbulent WNM interpretation, since turbulent motion is expected to decrease with radius. We discuss this in context with the line width of the stacked profiles in Section 6.5.

The increase in $f_{\text{wings}}$ beyond 6 kpc, driven entirely by the symmetric component, emphasizes that multiple processes must contribute to the line wings. The increase in this region is consistent with the behaviour in skewness and kurtosis analyses in Section 4.3, and we attribute the variations in Section 4.4.1 to M33's warped disc component. Corbelli et al. (2014) and Kam et al. (2017) both find evidence for the kinematic effects of the warp at similar radii.

## 6.5 A shallow radial gradient in the line width

As mentioned in Section 6.1, the central regions of many nearby galaxies have a steep radial gradient in the velocity dispersion, with notable exceptions including the flocculent spirals NGC 2403 and NGC 2976. From Fig. 16, we find that all of the line width measurements have a shallow radial trend. The line widths from the $v_{\text{peak}}$ stacked profiles decrease by just ~1.5 km s⁻¹ between a radius of 2 and 8 kpc. This variation is consistent with the radial decrease found in M33 by Druard et al. (2014) using 1 kpc bins.

The measured line widths from the $v_{\text{peak}}$ stacked profiles are narrower than the expected thermal line width of ~8 km s⁻¹ for the WNM (for a temperature of ~8500 K; Wolfire et al. 1995). And as found using the Gaussian Mixture model, the excess kurtosis universally found in stacked profiles requires a contribution of narrow Gaussian components with equivalent thermal widths of $T < 1000$ K. The line widths of the stacked profiles should then be some combination of varying thermal line widths, ranging from CNM to WNM temperatures, with turbulent broadening, which plausibly varies as well.

The line width traces the distribution of thermal and kinetic energy in the atomic medium. Cooling and dissipative mechanisms will remove this kinetic energy, so maintaining these line widths will require energy injection. Previous work has used the energy balance in the ISM of other galaxies to determine sources for the energy injection (e.g. Tamburro et al. 2009; Stilp et al. 2013b). All 11 galaxies examined in Tamburro et al. (2009) have a steep line width gradient, which tend to flatten towards the outer disc. They find that feedback via supernova is not sufficient to maintain the observed line widths in the outer regions, but magneto-rotation instabilities (MRI) plausibly could.

In particular, M33 has an increasing (not flat) rotation curve (Fig. 5), which contrasts with models used in the previous work. Similarly the star formation rate declines exponentially with radius (Heyer et al. 2004; Corbelli et al. 2014), providing insight into the rate of energy injection by core collapse supernova. We estimate the volumetric energy dissipation rate by assuming motions are turbulent with an outer scale set by the disc scale height of the atomic ISM in the galaxy (H, taken to be 100 pc; Tamburro et al. 2009) and dissipate on crossing times ($\tau_c = H\sigma^{-1}$) for the system:

$$ \dot{u}_{\text{diss}} = \frac{3 \Sigma \sigma^2}{2 H \tau_c} = \frac{3 \Sigma \sigma^3}{2 H^2}, \qquad (15) $$

where $\Sigma$ is the surface density of the atomic ISM. The rate of MRI energy injection has been estimated as

$$ \dot{u}_{\text{MRI}} = 0.6 \frac{B^2}{8\pi} R \frac{d\Omega}{dR} \qquad (16) $$

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Radial profiles of energy injection for core-collapse supernova (SN; solid red line) and MRI (dashed blue line) and estimates of the turbulent dissipation rates in the atomic ISM. Neither mechanism, by itself, matches the range and constant profile of the observed dissipation rates.</figcaption></figure>

following Mac Low & Klessen (2004). Here $B$ is the characteristic magnetic field strength, taken to be a constant $B = 8$ µG with radius based on modelling of the non-thermal radio continuum by Tabatabaei et al. (2008). We use our rotation curve model to evaluate orbital frequency $\Omega$ at galactocentric radius $R$. We take the energy injection from core-collapse supernova from Mac Low & Klessen (2004):

$$ \dot{u}_{\text{SN}} = \epsilon_{\text{SN}} E_{\text{SN}} f_{\text{SN}} \frac{\dot{\Sigma}_{*}}{\langle m \rangle H}. \qquad (17) $$

Here $\epsilon_{\text{SN}} \sim 0.1$ is the efficiency with which supernovae inject energy into the ISM, $E_{\text{SN}} = 10^{51}$ erg is the energy of a supernova, $\dot{\Sigma}_*$ is the star formation rate surface density, $\langle m \rangle = 10^2 M_{\odot}$ is the average mass of stars that forms per core collapse supernova (Tamburro et al. 2009), and $H$ is the scale height of star formation, taken to be 100 pc again. We use the star formation rate maps from Boquien et al. (2015), specifically the combined maps from the *Galex* FUV and *Spitzer*/MIPS 24 µm data. We generate a radially averaged star formation rate profile as the median of data in annular bins with width of 50 pc.

In Fig. 20, we plot the radial profiles of the turbulent energy dissipation adopting a range of line widths from 6 to 12 km s⁻¹ and a constant surface density of 8 $M_{\odot}$ pc⁻², which describe our data well. Neither injection mechanism can explain the line widths in isolation. Because of the rising rotation curve and comparatively low shear, MRI driving rates are significantly below the dissipation rate estimates, but the MRI mechanism shows a flat profile consistent with the observed shallow trends in line width and surface density. The star formation rates fall off exponentially and cannot explain the flat profile of the line widths and surface density but it does broadly agree with the range of required energy injection rates. Appealing to Type Ia supernovae as an additional source of turbulent motion will not change the profile substantially since the old stellar disc in M33 also declines exponentially with a scale length of 1.4 kpc (Regan & Vogel 1994).

This model adopts a constant atomic gas scale height throughout the observed region, but the scale height likely increases with galactocentric radius. While the radial variation of the scale height of gas in M33 has not been measured, other disc galaxies show that the atomic gas scale height increases linearly by a factor of two
over the optical radius of the galaxy, $R_{25}$ (e.g. Yim et al. 2014). For M33, $R_{25} = 9$ kpc (de Vaucouleurs et al. 1991), implying that $\dot{u}_{diss}$ would decrease by a factor of <4 over the plotted region. The supernova driving would also show a steeper decline, dropping by an additional factor of two. However, the thickness of the star-forming disc is relatively constant (Yim et al. 2014) so the scale height used in $\dot{u}_{SN}$ may remain constant with radius. Even with this assumption, a flaring H I disc does not clearly match the factor of $10^2$ decline in supernova driving in this simple model. Using higher resolution H I data, it will become possible to map spatial variations in the scale height through changes in the turbulent properties of the gas (Padoan et al. 2001) and investigate these effects in more detail.

Without a more detailed study, we thus arrive at a similar conclusion as Stilp et al. (2013b), who examine a number of energy sources for driving turbulence in their sample of dwarf galaxies and find that a single source cannot maintain the observed line widths.

## 6.6 Interpretation of different stacking models

The properties of stacked profiles depend on how the spectra are aligned. The appropriate choice of which definition of the line centre to use – the rotation velocity ($v_{rot}$), centroid velocity ($v_{cent}$), or velocity at peak intensity ($v_{peak}$) – depends on the goal of the analysis.

Stacking based on $v_{rot}$ maximizes contributions from non-circular motions in the disc. In M33, the width of the velocity residuals is significant when compared to the typical WNM thermal width. The line width of this stacked profile will measure the average kinetic energy in the ISM in the galactic potential.

The properties of the atomic ISM components at the spatial resolution of the data are best recovered by stacking based on $v_{peak}$ since this minimizes the variance in the line centre distribution for the bright components. Line wings are also the most significant in this method but trace distributions of temperature in the ISM as well as multiple velocity components along the line of sight (Section 6.3). Assuming the statistical uncertainty on $v_{peak}$ is small, the resulting stacked profile will recover the average line shape at the resolution of the data. We explicitly test this in a forthcoming paper (Koch et al., in preparation) and our preliminary results show that individual H I components have a median line width of ~6 km s⁻¹ based on a Gaussian fit, which matches the 6–7 km s⁻¹ found for the $v_{peak}$ stacked profiles (Table 4). However, in our upcoming work, we find that individual components have widths ranging from 4 to 10 km s⁻¹. While the $v_{peak}$ profiles appear to recover the typical spectrum’s properties, it necessarily removes spatial variations in the widths.

The line width of the $v_{cent}$ stacked profile can be interpreted as tracing the kinematics of the internal motions of the neutral medium, since, for an optically thin tracer, the centroid is the mass-weighted average velocity. The dispersion around this velocity represents the energy of the gas-to-gas relative motions, as opposed to the gas-to-galactic potential motions traced by the rotation-stacked profile. Because $v_{cent}$ is biased by the line shape, however, caution must be used if the line shapes have a preferential skewness as this will artificially broaden the velocity distribution. Creating a map of the skewness (Fig. 7) or checking the asymmetry of stacked spectra over limited regions of the galaxy (Fig. 12) can be used to diagnose how substantial the bias will be.

Thus, the line width in each of the stacked profile types result in a different but physically meaningful measure of typical atomic ISM

properties. Differences in the interpretations of these line widths present in the literature is partly attributable to real variations, but also arise from the method used for modelling. Most of the literature aims to characterize components of the ISM, however we caution that these interpretations are necessarily uncertain, even with good velocity and spatial resolution. We demonstrate this uncertainty in the source of stacked profile properties with the Gaussian Mix-ture Model in Section 6.3. Reasonable distributions of ISM prop-erties – for example the temperature – can produce the qualitative properties of observed stacked profiles. However, a realistic mix-ture model must also include distributions for other physical pro-cesses, such as turbulent broadening, outflows from stellar feedback, bulk non-circular motion and lagging rotational motion, amongst others.

For the M33 data, the best recourse is to perform a multicomponent Gaussian decomposition on a per line-of-sight basis. Such methods can realistically only work for good-quality observations of nearby galaxies. These detailed analyses are able to separate the effects of different processes by retaining spatial relationships. We recommend that studies seeking to separate physical processes either model individual spectra (e.g. Warren et al. 2012) or average spectra over the size scale of the targetted processes. For more distant targets where the scale of most processes will be unresolved, or for poor S/N data, stacking methods are still useful. Though the stacked profile will combine many physical processes together, it can be used as a powerful measure of similarity between different galaxies. For example, the consistent line widths found by Ianjamasimanana et al. (2012) between dwarf and spiral galaxies suggests that some atomic ISM properties are unchanged by galaxy type.

## 7 SUMMARY

We present new L-band observations of M33 taken by the VLA in C-configuration. The new data set yields a spatial resolution of ~18 arcsec, tracing ~80 pc scales, and significantly improves the sensitivity and spectral resolution compared to previous observations with an rms brightness temperature of 2.8 K in a 0.2 km s⁻¹ channel.

(i) The new H I VLA data recover 72 per cent of the total H I emission in M33 compared to estimates from single-dish data within the VLA survey region. After combining the VLA data with GBT observations, the total emission matches the emission from the GBT data alone. This gives a total atomic mass of $1.3 \pm 0.2 \times 10^9$ M$_\odot$. The combined data cube is fully sampled down to 80 pc scales.

(ii) We fit a circular rotation model to constrain the kinematic parameters on 100 pc scales. The rotation curve is well represented by a Brandt (1960) model with $v_{max} = 110.0 \pm 1.5$ km s⁻¹ at $r_{max} = 12.0 \pm 1.3$ kpc. In general, we find good agreement with previously published rotation curves for the inner 8 kpc of the disc.

(iii) For galactocentric radii $R < 7$ kpc, the azimuthally averaged atomic gas profile has a nearly constant average surface density of $\Sigma_{H I} = 8$ M$_\odot$ pc⁻², though the observations highlight local variations of 25 per cent with radius. These variations are seen in both limbs (north and south) of the galaxy.

(iv) The H I line profiles are consistently non-Gaussian. We parametrize this by calculating higher order moments – skewness and kurtosis – for the line profiles. We find that H I profiles are asymmetrically skewed towards the systemic velocity, resulting in the Northern and Southern halves having oppositely signed skewness. By examining PV-slices, we find evidence for a lag-
ging rotational component in M33, consistent with Kam et al. (2017). We find a typical excess in kurtosis of ~0.2, indicating that the typical spectrum has an excess in the line wings relative to a Gaussian.

(v) We stack spectra over the entire galaxy, and in 100 pc radial bins. By modelling the profiles as Gaussian peaks with enhanced line wings, we find that stacking based on the velocity at the peak intensity ($v_{peak}$) of spectra gives the smallest widths of ~7 km s⁻¹. This is consistent with arising from a combination of cool and warm atomic gas with some turbulent broadening. Stacking based on the centroid ($v_{cent}$) and rotation ($v_{rot}$) velocities give large line widths of ~10 km s⁻¹, which provides an estimate of internal motions in the atomic medium. However, the centroid is biased by the asymmetric line wings, which will tend to broaden the stacked profile. Line width estimates from the second moment are larger than all of the stacked profile widths and are sensitive to extended line wings.

(vi) All estimates of the line width show a shallow decrease with radius, dropping ~2 km s⁻¹ over the inner 8 kpc, though with significant (20 per cent) fluctuations. This is atypical for most nearby galaxies, which have steeper slopes at small galactocentric radii. We find that simple estimates of the volumetric energy dissipation rates from core-collapse supernova and MRI cannot explain the radial trend in the line width. The rising rotation curve and low shear in M33 leads to MRI energy injection rates that are significantly lower than the estimated range from the observed line widths. More careful measurements of the scale height of atomic gas in M33 are needed to make these conclusions robust.

(vii) The fraction of excess emission in the stacked profile line widths ranges from 9 per cent to 26 per cent depending on the choice of line centre. There is no clear radial trend in the inner 6 kpc, but local variations of 50 per cent are significant. We split the line wing fraction into asymmetric and symmetric parts based on the stacked profiles in the Northern and Southern halves. We find the asymmetric part, from the lagging disc, accounts for 1/3 of the line wing excess and is nearly constant with radius. The radial variations in the line wings are almost entirely driven by the symmetric part. Beyond 6 kpc, the symmetric part steadily increases, which we attribute to the start of M33's warped disc.

(viii) We present a Gaussian Mixture Model to explain stacked profile shapes. Using only a set of Gaussian components with thermal line widths, we can qualitatively reproduce the stacked profile shape found in many previous studies. In particular, the highly peaked centre of the profile can only be reproduced by including line widths expected for the CNM. However, the model line widths are smaller than the observed range, suggesting that a combination of the CNM/intermediate/WNM fractions, turbulence, a lagging rotational disc, and multicomponent spectra act to broaden the observed stacked profiles. This large possible combination of physical processes cannot be extracted from the stacked profile without additional information.

(ix) We identify discrete high-velocity gas structures on the blue- and red-shifted sides of the disc. These structures have a total mass of 1.3 ± 0.5 × 10⁷ M☉ in this high-latitude component, about 1 per cent of the total atomic mass in the system. Some of these structures were previously identified with single-dish studies. Here, we can resolve these structures and their extent across the main disc. All high-velocity structures appear as a clump surrounded by fainter filamentary structure. We find two features of particular interest: a long filament in the Southern half of the disc with a project length of 8 kpc, and a cloud 600 pc in diameter overlapping in velocity with the main disc, indicating a possible interaction point. These

structures highlight the complexities of M33's halo, and possibly a connection with the warped disc component.

(x) We find no detections of RRL emission, based on stacking six RRL lines. We set a 3σ upper limit of 3.0 mJy in a 60 arcsec region towards NGC 604.

In conclusion, the high spectral resolution observations of M33 show the galaxy harbours a kinematically rich atomic medium. With complete spatial sampling down to 80 pc scales, our analysis shows that several common analysis paths used in extragalactic observations have limitations in describing these spectra. We highlight that overcoming these limitations requires decomposing the spectra into their individual components to gain a full understanding of the atomic ISM kinematics. This provides new opportunities for exploring the data using tools from the Milky Way community.

The spectral-line VLA data and derived analysis products are available at https://doi.org/10.11570/18.0003, including derived moment maps and azimuthally averaged properties used throughout this paper. Scripts to reproduce the reduction, imaging, and analysis are available at https://github.com/e-koch/VLA_Lband<fnref n="8" />.

## ACKNOWLEDGEMENTS

We are grateful for helpful discussions with Doug Johnstone in developing this paper, as well as Adam Ginsburg and Hauyu Baobab Liu for discussions on imaging. We also thank the anonymous referee for their feedback and comments that improved the paper. EWK is supported by a Postgraduate Scholarship from the Natural Sciences and Engineering Research Council of Canada (NSERC). EWK and EWR are supported by a Discovery Grant from NSERC (RGPIN-2012-355247; RGPIN-2017-03987). This research was enabled in part by support provided by WestGrid (www.westgrid.ca), Compute Canada (www.computecanada.ca), and CANFAR (www.canfar.net). The National Radio Astronomy Observatory and the Green Bank Observatory are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

**Code bibliography:** CASA (version 4.4 and 4.7; McMullin et al. 2007) – ASTROPY (Astropy Collaboration 2013) – RADIO-ASTRO-TOOLS (spectral-cube, radio-beam, uvcombine; radio-astro-tools.github.io) – MATPLOTLIB (Hunter 2007) – SEABORN (Waskom et al. 2017) – DISKFIT (Spekkens & Sellwood 2007; Sellwood & Spekkens 2015) – IMAGE REGISTRATION (http://image-registration.readthedocs.io)

## REFERENCES

Araya E., Baan W. A., Hofner P., 2004, ApJS, 154, 541

Astropy Collaboration, 2013, A&A, 558, A33

Barbieri C. V., Fraternali F., Oosterloo T., Bertin G., Boomsma R., Sancisi R., 2005, A&A, 439, 947

Barnabè M., Ciotti L., Fraternali F., Sancisi R., 2006, A&A, 446, 61

Boothroyd A. I., Blagrave K., Lockman F. J., Martin P. G., Pinheiro Gonçalves D., Srikanth S., 2011, A&A, 536, A81

Boquien M. et al., 2015, A&A, 578, A8

Boulanger F., Viallefond F., 1992, A&A, 266, 37

Brandt J. C., 1960, ApJ, 131, 293

Braun R., 1997, ApJ, 484, 637

Braun R., 2012, ApJ, 749, 87

Braun R., Thilker D. A., Walterbos R. A. M., Corbelli E., 2009, ApJ, 695, 937

Caldú-Primo A., Schruba A., 2016, AJ, 151, 34

<footnote n="8">⁸v1.1; https://doi.org/10.5281/zenodo.1300005</footnote>
Caldu-Primo A., Schruba A., Walter F., Leroy A., Sandstrom K., de Blok
W. J. G., Ianjamasimanana R., Mogotsi K. M., 2013, AJ, 146, 150

Corbelli E., Thilker D., Zibetti S., Giovanardi C., Salucci P., 2014, A&A,
572, A23

Cornwell T. J., 2008, *IEEE J. Select. Top. Signal Process.*, 2, 793

de Blok W. J. G., Walter F., Brinks E., Trachternach C., Oh S. H., Kennicutt
R. C. J., 2008, AJ, 136

de Vaucouleurs G., de Vaucouleurs A., Corwin H. G., Jr, Buta R. J., Paturel
G., Fouqué P., 1991, Third Reference Catalogue of Bright Galaxies.
Volume I: Explanations and references. Volume II: Data for galaxies
between 0<sup>h</sup> and 12<sup>h</sup>. Volume III: Data for galaxies between 12<sup>h</sup> and 24<sup>h</sup>.

Deul E. R., van der Hulst J. M., 1987, A&AS, 67, 509

Dickey J. M., Hanson M. M., Helou G., 1990, ApJ, 352, 522

Dickey J. M., Mebold U., Marx M., Amy S., Haynes R. F., Wilson W., 1994,
A&A, 289, 357

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

Dobbs C. L. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P.,
Henning T., eds, Protostars and Planets VI, Univ. Arizona Press, Tucson,
AZ, p. 3

Druard C. et al., 2014, A&A, 567, A118

Field G. B., Goldsmith D. W., Habing H. J., 1969, ApJ, 155, L149

Fraternali F., Oosterloo T., Sancisi R., van Moorsel G., 2001, ApJ, 562, L47

Fraternali F., van Moorsel G., Sancisi R., Oosterloo T., 2002a, AJ, 123, 3124

Fraternali F., Cappi M., Sancisi R., Oosterloo T., 2002b, ApJ, 578, 109

Fraternali F., Oosterloo T., Sancisi R., 2004, A&A, 424, 485

Freedman W. L. et al., 2001, ApJ, 553, 47

Fukui Y. et al., 2009, ApJ, 705, 144

Gibson S. J., Taylor A. R., Higgs L. A., Brunt C. M., Dewdney P. E., 2005,
ApJ, 626, 195

Gratier P. et al., 2010, A&A, 522, A3

Grossi M., Giovanardi C., Corbelli E., Giovanelli R., Haynes M. P., Martin
A. M., Saintonge A., Dowell J. D., 2008, A&A, 487, 161

Heald G. et al., 2011, A&A, 526, A118

Heiles C., Troland T. H., 2003, ApJ, 586, 1067

Henshaw J. D. et al., 2016, MNRAS, 457, 2675

Heyer M. H., Corbelli E., Schneider S. E., Young J. S., 2004, ApJ, 602, 723

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Ianjamasimanana R., de Blok W. J. G., Walter F., Heald G. H., 2012, AJ,
144, 96

Ianjamasimanana R., de Blok W. J. G., Walter F., Heald G. H., Caldu-Primo
A., Jarrett T. H., 2015, AJ, 150, 47

Imara N., Bigiel F., Blitz L., 2011, ApJ, 732, 79

Kam S. Z., Carignan C., Chemin L., Foster T., Elson E., Jarrett T. H., 2017,
AJ, 154, 41

Keenan O. C., Davies J. I., Taylor R., Minchin R. F., 2016, MNRAS, 456,
951

Koch E., Rosolowsky E., Johnson M. C., Kepley A. A., Leroy A., 2018,
Res. Notes AAS, 2, 24

Kuzio de Naray R., Arsenault C. A., Spekkens K., Sellwood J. A., McDonald
M., Simon J. D., Teuben P., 2012, MNRAS, 427, 2523

Lindner R. R. et al., 2015, AJ, 149, 138

Lockman F. J., 2017, in Fox A., Davé R., eds, Astrophysics and Space
Science Library, Vol. 430, Gas Accretion onto Galaxies, Springer, Cham.
p. 49

Lockman F. J., Murphy E. M., Petty-Powell S., Urick V. J., 2002, ApJS,
140, 331

Lockman F. J., Free N. L., Shields J. C., 2012, AJ, 144, 52

Mac Low M.-M., Klessen R. S., 2004, Rev. Modern Phys., 76, 125

Mangum J. G., Emerson D. T., Greisen E. W., 2007, A&A, 474, 679

McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw
R. A., Hill F., Bell D. J., eds, ASP Conf. Ser. Vol. 376, Astronomical Data
Analysis Software and Systems XVI. Astron. Soc. Pac., San Francisco,
p. 127

Meidt S. E., Rand R. J., Merrifield M. R., Shetty R., Vogel S. N., 2008, ApJ,
688, 224

Mogotsi K. M., de Blok W. J. G., Caldu-Primo A., Walter F., Ianjamasi-
manana R., Leroy A. K., 2016, AJ, 151, 15

Murray C. E. et al., 2015, ApJ, 804, 89

Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt R. C., Jr, 2011,
AJ, 141, 193

Padoan P., Kim S., Goodman A., Staveley-Smith L., 2001, ApJ, 555, L33

Petric A. O., Rupen M. P., 2007, AJ, 134, 1952

Putman M. E. et al., 2009, ApJ, 703, 1486

Putman M. E., Peek J. E. G., Joung M. R., 2012, ARA&A, 50, 491

Regan M. W., Vogel S. N., 1994, ApJ, 434, 536

Rosolowsky E., Keto E., Matsushita S., Willner S. P., 2007, ApJ, 661, 830

Sancisi R., Fraternali F., Oosterloo T., van der Hulst T., 2008, A&AR, 15,
189

Sellwood J. A., Sánchez R. Z., 2010, MNRAS, 404, 1733

Sellwood J. A., Spekkens K., 2015, preprint (arXiv:1509.07120)

Shih F. Y., 2009, Image Processing and Mathematical Morphology. Funda-
mentals and Applications. , CRC Press, LLC

Skrutskie M. F. et al., 2006, AJ, 131, 1163

Spekkens K., Sellwood J. A., 2007, ApJ, 664, 204

Stanimirović S. et al., 2006, ApJ, 653, 1210

Stanimirovic S., 1999, PhD thesis, Univ. Western Sydney Nepean, Australia

Stil J. M. et al., 2006, AJ, 132, 1158

Stilp A. M., Dalcanton J. J., Warren S. R., Skillman E., Ott J., Koribalski
B., 2013a, ApJ, 765, 136

Stilp A. M., Dalcanton J. J., Skillman E., Warren S. R., Ott J., Koribalski
B., 2013b, ApJ, 773, 88

Tabatabaei F. S., Krause M., Fletcher A., Beck R., 2008, A&A, 490, 1005

Tamburro D., Rix H. W., Leroy A. K., Mac Low M. M., Walter F., Kennicutt
R. C., Brinks E., de Blok W. J. G., 2009, AJ, 137, 4424

Thilker D. A., Braun R., Walterbos R. A. M., 2002, in Taylor A. R., Lan-
decker T. L., Willis A. G., eds, ASP Conf. Ser. Vol. 276, Seeing Through
the Dust: The Detection of HI and the Exploration of the ISM in Galax-
ies. Astron. Soc. Pac., San Francisco, p. 370

Vargas C. J. et al., 2017, ApJ, 839, 118

Walter F., Brinks E., de Blok W. J. G., Bigiel F., Kennicutt R. C. J., Thornley
M. D., Leroy A., 2008, AJ, 136, 2563

Wang J., Taaffe M., 2015, INFORMS J. Comput., 27, 193

Warren S. R. et al., 2012, ApJ, 757, 84

Waskom M. et al., 2017, mwaskom/seaborn: v0.8.1 (September 2017),
https://doi.org/10.5281/zenodo.883859

Westmeier D., Braun R., Thilker D., 2005, A&A, 436, 101

Wilcox R., 2010, Fundamentals of Modern Statistical Methods, 2nd edn.
Springer, New York, NY

Williams B. F., Dalcanton J. J., Stilp A., Dolphin A., Skillman E. D.,
Radburn-Smith D., 2013, ApJ, 765, 120

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

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

Wong T. et al., 2009, ApJ, 696, 370

Yim K., Wong T., Xue R., Rand R. J., Rosolowsky E., van der Hulst J. M.,
Benjamin R., Murphy E. J., 2014, AJ, 148, 127

Young L. M., Lo K. Y., 1996, ApJ, 462, 203

Young L. M., van Zee L., Lo K. Y., Dohm-Palmer R. C., Beierle M. E.,
2003, ApJ, 592, 111

Zheng Y., Peek J. E. G., Werk J. K., Putman M. E., 2017a, ApJ, 834, 179

Zheng Y., Werk J. K., Peek J. E. G., Putman M. E., 2017b, ApJ, 840, 65

APPENDIX A: IMAGING APPROACH

At a spectral resolution of 0.2 km s⁻¹ and a spatial grid size of 2560² pixels needed to cover the entire mosaic, imaging and deconvolution requires significant computational time and power. The size of the resulting cube is 1178 × 2560 × 2560 pixels, giving a size of ~29 GB. Rather than imaging the cube as a whole, we split the data into individual velocity channels, image and deconvolve each channel separately, then recombine the imaged channels into a final data cube. When working on a cluster, this approach allows channels to be simultaneously imaged, providing a significant
speed-up in the time required to image the entire cube. Furthermore, since the size of the channel measurement set (MS) is much smaller than the original MS, slow I/O operations are relatively minimized. Using CASA 4.4, we found that imaging a single channel from the complete MS was ~5x slower than when imaging from a split channel MS.

This approach to imaging large data cubes has a significant bottleneck during the splitting stage: the time required to split all 1178 channels from the complete MS was nearly 2 weeks on the Jasper cluster,<fnref n="9" /> which uses a lustre-based file system. In hindsight, we note that this is not an optimized choice of operation; a binary-split algorithm that progressively splits an MS would likely achieve a significant speed-up in the operation, particularly since this can be naturally parallelized. Furthermore, the storage requirements are more than double the size of the original MS. The new CASA task TCLEAN allows for the MS to be opened in a read-only mode, largely mitigating the need to split off individual channels. This approach will be used in the future.

## A1 Image combination

We explore the effects of combining the VLA data with Arecibo (Putman et al. 2009) and Green Bank Telescope (GBT; Lockman et al. 2012) observations to provide short-spacing information. Both are well suited to be combined with our VLA data as they have similar spectral resolution and at least a factor of two spatial overlap in the uv-plane. We note that older GBT data were combined with the archival VLA observations presented in Braun (2012), but these data have a spectral resolution of 1.42 km s⁻¹ and are not well-suited for combination with the 0.2 km s⁻¹ resolution VLA data presented here.

The Arecibo data from Putman et al. (2009) have a spatial resolution of 3.4 arcmin and a spectral resolution 0.4 km s⁻¹. This provides a significant spatial overlap with the VLA data, but also requires up-sampling in the spectral domain, thereby increasing the channel-to-channel correlation. We first spatially register the data using the cross-correlation method in IMAGE_REGISTRATION<fnref n="10" /> and find that no correction is needed. We use the FEATHER_simple<fnref n="11" /> task from the UVCOMBINE package to combine the data. This method of combining the data assumes the single-dish beam is well-approximated by an isotropic Gaussian kernel. However, the Arecibo data contain significant side-lobe structure, which leads to enhanced negative bowling in the feathered image on scales of ~1 arcmin. Further, the power spectrum of the Arecibo data, clipped to match the region covered by the VLA mosaic, shows a drop in power on the scale of the entire image, indicating a large-scale ripple. This ripple is far less significant in the larger region around M33 used by Putman et al. (2009). For these reasons, we did not pursue further image combination with the Arecibo data.

The GBT data from Lockman et al. (2012) have a nominal spatial resolution of 91 arcmin and a spectral resolution 0.16 km s⁻¹. Unlike the Arecibo data, no spectral up-sampling is required to match the VLA data. These data were taken in four 2° × 2° regions centred on M33; aspects of the GBT gridding are presented in Section 2.2. The effective resolution of the data is 9.8 arcmin (see below), which gives sufficient uv-overlap for combining with

the VLA data. The GBT beam has minimal side-lobe structure after calibration and the combined images with the VLA data do not show the enhanced bowling we encountered with the Arecibo data. Before feathering, we test if the data are spatially registered and find that the GBT data require a 3 arcsec shift in declination. Visually, this shift appears to provide a better combination with the data compared to when no shift is applied. However, we note that 3 arcsec is within the pointing error for the GBT data.

### A1.1 Combination tests

Following chapter 3 of Stanimirovic (1999), we run two tests on the uv-amplitudes where the spatial coverage of the VLA and GBT data overlap to (1) ensure a bias is not added from using an incorrect single-dish beam model, and (2) derive a relative calibration factor to obtain a consistent flux-density scale. We define the uv overlap region as all points between the 9.8 arcmin beam size of the GBT data and 16.5 arcmin – for the shortest baseline of 44 m in the VLA data. This gives ~200 overlap points per channel at the grid size of the VLA data (3 arcsec × 3 arcsec). To get the GBT amplitudes used below, we deconvolve its Fourier transform by dividing by the Fourier transform of the GBT beam. The GBT amplitudes in the overlap region are then multiplied by the ratio between the beam areas to account for the difference in resolution (Stanimirovic 1999).

In section 3.2.1 of Stanimirovic (1999), an approximate relation between the scale factor and the uv-distance k is given:

$$f_{cal} = \left[ 1 + \frac{\Delta\theta(2\theta_0 + \Delta\theta)}{4\ln2} k^2 \right], \quad (A1)$$

where θ₀ is the true FWHM single-dish beam size and Δθ is the deviation from the true beam size and assumed to be small. This predicts a linear relation between fcal and k² when an incorrect beam size is assumed. The correct beam size should have a slope of zero. We perform this fit for each velocity channel in the data using a robust Theil-Sen estimator to find the slope. A robust fitting method is required since the combined effects of the noise and different emission structure in the channels will naturally produce significant outliers. The Theil-Sen method computes the slope from the median of all slopes between each pair of data points, as is defined in the scipy implementation.<fnref n="12" /> The slopes for each channel are shown in Fig. A1 using a GBT beam size of 9.8 arcmin. This beam size maximizes the number of channels with slopes consistent with zero. The variations in the slope with different channels – notably the first 100 and last 200 channels – are driven by the lack of emission structure and correspondingly have larger uncertainties. The channels dominated by signal in the uv-overlap range are from 200 to 800, and since this region has slopes largely consistent with zero, we adopt a GBT beam size of 9.8 arcmin for feathering.

With the correct GBT beam size established, we now address the scaling factor between the VLA and GBT amplitudes. We follow a similar procedure to the one described in Section 3.2 in Stanimirovic (1999). The emission structure in both the VLA and GBT smoothly tapers near the map edges and no additional edge tapering is required. We first fit a line between the GBT and VLA amplitudes, similar to the method used by IMMERGE in MIRIAD. The issue

<footnote n="9">⁹https://www.westgrid.ca/resources_services</footnote>

<footnote n="10">¹⁰http://image-registration.readthedocs.io</footnote>

<footnote n="11">¹¹This matches CASA’s feather implementation.</footnote>

<footnote n="12">¹²https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.theilslop.html#scipy.stats.theilslopes</footnote>
<figure><img src="image_23.png" /><figcaption>**Figure A1.** Derived slopes of the *uv*-distance versus the ratio between the VLA and GBT amplitudes in the overlap region. A slope consistent with zero suggests the beam model is correct; an incorrect beam model will incorrectly bias the amplitudes, giving a non-zero slope.</figcaption></figure>

with the linear fitting approach is consistently dealing with outliers. We adopt the Theil-Sen method rather than the L1-minimization used in *IMMERGE* since the former tends to be insensitive to outliers in both the *x* and *y* directions, while the latter is only robust against outliers in the *y* direction (Wilcox 2010). We fit a relation to the amplitudes in groups of five velocity channels to increase the number of points. Despite the robust nature of the Theil-Sen fit, we did not find consistent results across the channel groups. The issue is the severity of the outliers and their fraction relative to the entire data set: if there are too many extreme outliers, the Theil-Sen fit still has some sensitivity to the outliers. Increasing the number of channels simultaneously fit did not show improvement, and fitting the entire set of amplitudes is prohibitive for the Theil-Sen method since it computes the slope for every pair of data points.

Using a distribution-based method, we are able to provide more reliable constraints on the scaling factor. Motivated by Section 3.2 in Stanimirovic (1999), we first examine the distributions of the GBT amplitudes, the VLA amplitudes, and their ratio across all channels. The GBT and VLA amplitudes reasonably follow a log-normal distribution.<fnref n="13" /> Since the ratio of two normal random variables follows a Cauchy (or Lorentzian) distribution, the log of the amplitude ratios can be fit to this form. The scaling factor is then the location of the distribution's peak. Fig. A2 shows the distribution of the amplitudes across all channels (blue) and the best-fit Cauchy distribution (green). This approach has the significant advantage that outliers are included in the expected model (the distribution tails) and do not require special treatment. The Cauchy distribution is also preferable to the Rayleigh distribution used in Stanimirovic (1999), which gives more weight to the outliers in the right-tail. Adopting a maximum-likelihood approach, we find a scaling factor of 1.02 ± 0.06 using the ratios from all channels. As this is consistent with 1, we do not apply a scaling factor to the GBT data before feathering. We also perform the fitting on groups of five channels and find scaling factors consistent with one between channels 200 to 800, where most of the emission is contained.

<figure><img src="image_24.png" /><figcaption>**Figure A2.** The log of the VLA-to-GBT amplitude ratios in the *uv*-overlap region across all channels. The histogram of the ratios is normalized to one and the green line is the best-fitting Cauchy distribution. The peak of the distribution is at 1.02 ± 0.06, consistent with a scale factor of one.</figcaption></figure>

## APPENDIX B: CHOICE OF VELOCITY SURFACE FOR THE ROTATIONAL MODEL

Galactic rotation curves, with a few exceptions, are typically derived from velocity surfaces. A key question explored in several publications – notably de Blok et al. (2008) – is: which velocity surface optimally traces galactic rotation? Commonly used methods include the centroid ($v_{cent}$), the velocity at the peak intensity ($v_{peak}$), and Gauss–Hermite fitting. We create velocity surfaces of the former two methods for the VLA-only and combined H I cubes and fit a rotation curve using *DISKFIT* with the same parameter settings. We did not fit a rotation model for the Gauss–Hermite velocity surface for two reasons. First, using the VLA-only cube, the difference between the Gauss–Hermite centre and $v_{peak}$ was much smaller than the difference to the $v_{cent}$ surface and will yield similar results. The second reason is that the Gauss–Hermite fitting did not perform well with the enhanced line wings for the combined data.

We find that the rotation models using the VLA-only $v_{peak}$ and $v_{cent}$ falls within the uncertainty of the combined VLA and GBT $v_{cent}$ surface, whose rotation curve is presented in Section 4.1. However, the $v_{cent}$ surface from the combined data has a shallower rotation curve. The reason for this discrepancy, and the agreement of the other rotation models, is due to the bias from the line wings. Without the component from the GBT, the VLA-only $v_{cent}$ is not as biased by line wings and largely matches the $v_{peak}$.

## APPENDIX C: ROTATION VELOCITIES

We provide a table of the calculated rotational velocities from *DISKFIT* presented in Section 4.1.

## APPENDIX D: MODELLING SUPERPROFILES

We demonstrate an alternative modelling approach for the super-profiles: fitting two stacked Gaussian profiles. Fig. D1 shows the $v_{peak}$ stacked profiles presented in Section 4.4 with a fitted two-Gaussian model and the model residuals. The model parameters are given below. Note that we do not list uncertainties due to the large covariance between parameters, which we discuss below.

These two stacked spectra drove our decision to use the HWM scaling from Stilp et al. (2013a), since the central peak is non-Gaussian and the resulting two-Gaussian fit does not provide a

<footnote n="13">¹³Note that this does not imply a log-normal distribution for the true signal. The log-normal shape results from the mixing between signal, noise, and other telescope effects.</footnote>
<table><caption>**Table C1.** Circular rotation velocities derived by DISKFIT (Section 4.1). The rotation model is fitted to the peak velocity ($v_{peak}$) surface of the VLA + GBT data.</caption><thead><tr><th>Radius (arcsec)</th><th>Circ. velocity (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>9</td><td>1.58 ± 6.24</td></tr><tr><td>27</td><td>11.88 ± 5.26</td></tr><tr><td>45</td><td>22.71 ± 5.25</td></tr><tr><td>63</td><td>30.53 ± 5.12</td></tr><tr><td>81</td><td>29.95 ± 4.47</td></tr><tr><td>99</td><td>34.20 ± 4.03</td></tr><tr><td>117</td><td>36.76 ± 3.99</td></tr><tr><td>135</td><td>44.61 ± 3.79</td></tr><tr><td>153</td><td>47.80 ± 3.49</td></tr><tr><td>171</td><td>51.67 ± 3.08</td></tr><tr><td>189</td><td>52.73 ± 2.96</td></tr><tr><td>207</td><td>53.24 ± 3.06</td></tr><tr><td>225</td><td>55.30 ± 3.22</td></tr><tr><td>243</td><td>57.58 ± 3.30</td></tr><tr><td>261</td><td>58.70 ± 3.42</td></tr><tr><td>279</td><td>59.47 ± 3.50</td></tr><tr><td>297</td><td>60.30 ± 3.59</td></tr><tr><td>315</td><td>64.32 ± 3.78</td></tr><tr><td>333</td><td>67.66 ± 3.75</td></tr><tr><td>351</td><td>69.07 ± 4.03</td></tr><tr><td>369</td><td>72.05 ± 4.49</td></tr><tr><td>387</td><td>75.89 ± 4.82</td></tr><tr><td>405</td><td>74.09 ± 4.91</td></tr><tr><td>423</td><td>73.44 ± 3.29</td></tr><tr><td>441</td><td>77.06 ± 2.76</td></tr><tr><td>459</td><td>77.06 ± 2.80</td></tr><tr><td>477</td><td>77.82 ± 2.60</td></tr><tr><td>495</td><td>79.01 ± 2.48</td></tr><tr><td>513</td><td>80.96 ± 2.59</td></tr><tr><td>531</td><td>80.34 ± 2.50</td></tr><tr><td>549</td><td>80.35 ± 2.65</td></tr><tr><td>567</td><td>82.90 ± 2.66</td></tr><tr><td>585</td><td>86.23 ± 2.66</td></tr><tr><td>603</td><td>86.92 ± 2.59</td></tr><tr><td>621</td><td>86.13 ± 2.63</td></tr><tr><td>639</td><td>87.43 ± 2.57</td></tr><tr><td>657</td><td>87.89 ± 2.63</td></tr><tr><td>675</td><td>87.82 ± 2.35</td></tr><tr><td>693</td><td>91.08 ± 2.44</td></tr><tr><td>711</td><td>90.37 ± 2.25</td></tr><tr><td>729</td><td>88.22 ± 2.09</td></tr><tr><td>747</td><td>90.67 ± 2.36</td></tr><tr><td>765</td><td>92.32 ± 2.48</td></tr><tr><td>783</td><td>93.31 ± 2.26</td></tr><tr><td>801</td><td>94.15 ± 2.15</td></tr><tr><td>819</td><td>93.40 ± 2.37</td></tr><tr><td>837</td><td>94.02 ± 2.44</td></tr><tr><td>855</td><td>94.91 ± 2.50</td></tr><tr><td>873</td><td>95.57 ± 2.48</td></tr><tr><td>891</td><td>95.15 ± 2.22</td></tr><tr><td>909</td><td>93.19 ± 2.46</td></tr><tr><td>927</td><td>94.52 ± 2.63</td></tr><tr><td>945</td><td>95.00 ± 2.77</td></tr><tr><td>963</td><td>95.72 ± 2.90</td></tr><tr><td>981</td><td>96.89 ± 2.81</td></tr><tr><td>999</td><td>98.58 ± 2.73</td></tr><tr><td>1017</td><td>98.09 ± 2.61</td></tr><tr><td>1035</td><td>99.87 ± 2.53</td></tr><tr><td>1053</td><td>99.10 ± 2.41</td></tr><tr><td>1071</td><td>99.01 ± 2.40</td></tr><tr><td>1089</td><td>97.58 ± 2.21</td></tr><tr><td>1107</td><td>98.60 ± 2.26</td></tr><tr><td>1125</td><td>99.61 ± 2.22</td></tr></tbody></table>

<table><caption>**Table C1 – continued**</caption><thead><tr><th>Radius (arcsec)</th><th>Circ. velocity (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1143</td><td>99.61 ± 2.49</td></tr><tr><td>1161</td><td>102.02 ± 2.40</td></tr><tr><td>1179</td><td>101.84 ± 2.52</td></tr><tr><td>1197</td><td>102.76 ± 2.38</td></tr><tr><td>1215</td><td>102.66 ± 2.62</td></tr><tr><td>1233</td><td>102.57 ± 2.68</td></tr><tr><td>1251</td><td>104.08 ± 2.58</td></tr><tr><td>1269</td><td>103.24 ± 2.39</td></tr><tr><td>1287</td><td>103.17 ± 2.42</td></tr><tr><td>1305</td><td>103.33 ± 2.50</td></tr><tr><td>1323</td><td>103.41 ± 2.64</td></tr><tr><td>1341</td><td>103.91 ± 2.56</td></tr><tr><td>1359</td><td>102.91 ± 2.63</td></tr><tr><td>1377</td><td>105.50 ± 2.47</td></tr><tr><td>1395</td><td>104.05 ± 2.47</td></tr><tr><td>1413</td><td>106.03 ± 2.47</td></tr><tr><td>1431</td><td>104.89 ± 2.57</td></tr><tr><td>1449</td><td>106.68 ± 2.48</td></tr><tr><td>1467</td><td>105.53 ± 2.36</td></tr><tr><td>1485</td><td>105.80 ± 2.33</td></tr><tr><td>1503</td><td>104.03 ± 2.23</td></tr><tr><td>1521</td><td>105.58 ± 2.23</td></tr><tr><td>1539</td><td>105.69 ± 2.26</td></tr><tr><td>1557</td><td>106.32 ± 2.17</td></tr><tr><td>1575</td><td>105.97 ± 2.32</td></tr><tr><td>1593</td><td>105.58 ± 2.20</td></tr><tr><td>1611</td><td>106.18 ± 2.26</td></tr><tr><td>1629</td><td>106.08 ± 2.28</td></tr><tr><td>1647</td><td>106.24 ± 2.36</td></tr><tr><td>1665</td><td>106.39 ± 2.37</td></tr><tr><td>1683</td><td>106.47 ± 2.43</td></tr><tr><td>1701</td><td>107.19 ± 2.44</td></tr><tr><td>1719</td><td>105.88 ± 2.58</td></tr><tr><td>1737</td><td>106.58 ± 2.52</td></tr><tr><td>1755</td><td>106.09 ± 2.56</td></tr><tr><td>1773</td><td>106.30 ± 2.35</td></tr><tr><td>1791</td><td>107.23 ± 2.61</td></tr><tr><td>1809</td><td>104.01 ± 2.82</td></tr><tr><td>1827</td><td>104.50 ± 2.59</td></tr></tbody></table>

better description for the data. Of particular concern is that the model components cannot account for the central peak. The sum of the amplitudes in both fits is ~0.95. The narrow and wide amplitudes are 0.51 and 0.55 for the narrow component, and 0.45 and 0.4 for the wide component, for the VLA and combined profiles, respectively. The difference in the components is about the same as the missing peak intensity, making comparisons of their ratios uncertain.

We find widths of 4.1 and 10.5 km s⁻¹ for the VLA-only profile, and widths of 4.7 and 14.4 km s⁻¹ for the combined profile. The HWHM widths of 5.9 and 6.6 km s⁻¹ for these two profiles are ~50 per cent larger than the narrow component. This is consistent with comparing the common galaxies examined with these two methods in Stilp et al. (2013a) and Ianjamasimanana et al. (2012).

Since Gaussian components do not form an orthogonal set, there are large covariances between the parameters. Stilp et al. (2013a) explore this issue in depth and demonstrate that it is this covariance that leads to the overall good fits presented in Ianjamasimanana et al. (2012). Stilp et al. (2013a) show that the individual components from the two-Gaussian fit do not account for the overall profile shape. They demonstrate this by scaling all profiles to a common width, based on the different models, and find that the shapes from
<figure><img src="image_25.png" /><figcaption>**Figure D1.** Two Gaussian fits to the $v_{\text{peak}}$ stacked profiles shown in Fig. 11. The left-hand panel is a stack of the VLA-only data and the right uses the combined data. The model residuals are shown in the bottom panels.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure E1.** The azimuthally averaged line width estimated from the second moment in 100 pc bins. The combined VLA+GBT data (green dot dashed) are 30 per cent larger than the VLA-only values (blue solid). From the stacking analysis (Fig. 11), however, we find that adding the GBT component tends to only affect the shape of the line wings. The large increase from the inclusion of the GBT data highlights that estimates from the second moment are very sensitive to the line wing structure.</figcaption></figure>

the broad and narrow components have significantly more scatter compared to the HWHM estimate.

While the two-Gaussian fits provide an adequate representation of the profiles, there are model dependencies that impact the connection to the underlying physical parameters, to the extent that the superprofiles can give. Thus we use the HWHM modelling throughout this work.

## APPENDIX E: ISSUES WITH THE SECOND MOMENT FOR ESTIMATING THE LINE WIDTH IN LOCAL GROUP GALAXIES

In Fig. E1, we show that there is a stark difference between the averaged second moment line widths with and without the total power component added – the values increase by ~30 per cent. This highlights the impact of extended line wings when using the second moment to estimate the line width, making it difficult to directly connect with the underlying physical conditions. The large discrepancy seen here is a product of two factors. The first is the limit of small angular scales in the VLA data for a relatively nearby system. Fig. 1 shows that the combined VLA and GBT data recover significantly more of the large-scale disc structure in the Southern half than the VLA data only. This issue is less important for more distant systems (e.g. Walter et al. 2008). The second factor is the extensive extra-planar component in the warped disc of M33. As we suggest from the skewness maps (Section 4.3), the warped disc component influences the line shape near the edges of our map, increasing the line widths estimated from the second moment. The extent of this factor depends on the galaxy’s environment. However, these average line widths are consistent with the range found in other H1 studies of nearby galaxies (Tamburro et al. 2009; Mogotsi et al. 2016), showing that this effect is not strong enough to make M33 an outlier, or that line widths from the second moment in other galaxies are similarly affected.

This paper has been typeset from a TeX/LATEX file prepared by the author.