# Considerations with stacking absorption spectra: cold H i gas in cirrus region of the Milky Way

Callum Lynn, Antoine Marchal, N M Mcclure-Griffiths, Marc-Antoine

Miville-Deschênes, Claire E Murray, Hiep Nguyen, James Dempsey, Enrico Di

Teodoro, Jacco Th van Loon, John M Dickey, et al.

► To cite this version:

Callum Lynn, Antoine Marchal, N M Mcclure-Griffiths, Marc-Antoine Miville-Deschênes, Claire E Murray, et al. Considerations with stacking absorption spectra: cold H i gas in cirrus region of the Milky Way. Monthly Notices of the Royal Astronomical Society, 2025, 536 (4), pp.3538-3553. 10.1093/mnras/stae2818 . hal-05409584

HAL Id: hal-05409584

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

Submitted on 10 Dec 2025

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
# Considerations with stacking absorption spectra: cold H I gas in cirrus region of the Milky Way

Callum Lynn<sup>1,★</sup> Antoine Marchal<sup>1</sup> N. M. McClure-Griffiths<sup>1</sup> Marc-Antoine Miville-Deschênes<sup>1,2</sup>

Claire E. Murray<sup>3,4</sup> Hiep Nguyen<sup>1</sup> James Dempsey<sup>1,5</sup> Enrico Di Teodoro<sup>5,6</sup>

Jacco Th. van Loon<sup>7</sup> John M. Dickey<sup>8</sup> Min-Young Lee<sup>9,10</sup> Gilles Joncas<sup>11</sup> Yik Ki Ma<sup>1</sup>

Nickolas M. Pingel<sup>12</sup> Snežana Stanimirović<sup>12</sup> Ian Kemp<sup>13</sup> Steven Gibson<sup>14</sup> and Helga Dénes<sup>15</sup>

Accepted 2024 December 16. Received 2024 December 6; in original form 2024 July 16

## ABSTRACT

We use the Milky Way neutral hydrogen (H I) absorption and emission spectra from the Galactic Australian Square Kilometre Array Pathfinder (GASKAP) Phase II Pilot survey along with toy models to investigate the effects of stacking multicomponent spectra on measurements of peak optical depth and spin temperature. Shifting spectra by the peak in emission, ‘primary’ components shifted to 0 km s⁻<fnref n="1" /> are correctly averaged. Additional components on individual sightlines are averaged with non-centred velocities, producing a broader and shallower ‘secondary’ component in the resulting stack. Peak optical depths and brightness temperatures of the secondary components from stacks are lower limits of their true average values due to the velocity offset of each component. The spin temperature however is well correlated with the truth since the velocity offset of components affects the emission and absorption spectra equally. Stacking 462 GASKAP absorption-emission spectral pairs, we detect a component with a spin temperature of 1320 ± 263 K, consistent with gas from the unstable neutral medium and higher than any previous GASKAP detection in this region. We also stack 2240 pilot survey spectra containing no Milky Way absorption, revealing a primary narrow and secondary broad component, with spin temperatures belonging to the cold neutral medium (CNM). Spatially binning and stacking the non-detections across the plane-of-sky by their distance from CNM absorption detections, the primary component’s optical depth decreases with distance from known locations of cold gas. The spin temperature however remains stable in both components, over an approximate physical plane-of-sky distance of ~ 100 pc.

**Key words:** ISM: structure – solar neighbourhood – galaxies: ISM – radio lines: ISM.

## 1 INTRODUCTION

Neutral hydrogen (H I) is the most abundant element in the Universe, making it an important tracer for large-scale galactic morphology. Because H I is the building block of stars, understanding its temper-

ature in the current epoch of the interstellar medium (ISM) is crucial to our understanding of the current state of star formation in the Milky Way.

Through a combination of observations and theory, H I has been best described with a bi-stable model with two stable temperature regimes: the warm neutral medium (WNM), T<sub>s</sub> ~ 4000–8000 K in the Solar neighbourhood, filling most of the H I volume, and the cold

<footnote n="1"><sup>1</sup>Research School of Astronomy & Astrophysics, The Australian National University, Canberra ACT 2611, Australia</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris, F-75005 Paris, France</footnote>

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

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

<footnote n="5"><sup>5</sup>CSIRO Information Management and Technology, GPO Box 1700 Canberra, ACT 2601, Australia</footnote>

<footnote n="6"><sup>6</sup>Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, via G. Sansone 1, F-50019 Sesto Fiorentino, Firenze, Italy</footnote>

<footnote n="7"><sup>7</sup>Lennard-Jones Laboratories, School of Chemical and Physical Sciences, Keele University, Keele, Staffordshire ST5 5BG, UK</footnote>

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

<footnote n="9"><sup>9</sup>Korea Astronomy and Space Science Institute, 776 Daeadeok-daero, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>

<footnote n="10"><sup>10</sup>Department of Astronomy and Space Science, University of Science and Technology, 217 Gajeong-ro, Daejeon 34113, Republic of Korea</footnote>

<footnote n="11"><sup>11</sup>Dépt. de physique, de génie physique et d'optique, Centre de recherche en astrophysique du Québec, 1045 avenue de la médecine, Université Laval, Québec, Québec, G1V 0A6, Canada</footnote>

<footnote n="12"><sup>12</sup>Department of Astronomy, University of Wisconsin–Madison, 475 N Charter St, Madison, WI 53703, USA</footnote>

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

<footnote n="14"><sup>14</sup>Department of Physics & Astronomy, Western Kentucky University, 1906 College Heights Boulevard, Bowling Green, KY 42101, USA</footnote>

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

<footnote>*E-mail: Callum.Lynn@anu.edu.au</footnote>
neutral medium (CNM), $T_s \sim 25-250$ K, that exists in dense cold clumps enveloped by the WNM. Between these two temperature regimes lies the thermally unstable neutral medium (UNM), a region thought to be produced by the condensation of the WNM via thermal instability (Field 1965; Field, Goldsmith & Habing 1969; Wolfire et al. 1995; Hennebelle & Pérault 1999, 2000; Wolfire et al. 2003; Audit & Hennebelle 2005).

Observationally, the CNM is the easiest of the three phases to detect in absorption, even through low-sensitivity observations, due to its typically high-optical depth and column density (Lazareff 1975; Dickey, Salpeter & Terzian 1977; Crovisier, Kazes & Aubry 1978; Heiles & Troland 2003; Mohan, Dwarakanath & Srinivasan 2004; Begum et al. 2010; Roy et al. 2013). However, the UNM and WNM, due to their lower opacity, require high-sensitivity observations and instruments to be observed in absorption (Carilli, Dwarakanath & Goss 1998; Dwarakanath, Carilli & Goss 2002; Patra et al. 2018). Therefore, most of their physical properties have been estimated from H1 emission alone, decomposed through methods of Gaussian decomposition (Haud 2000; Nidever, Majewski & Butler Burton 2008; Lindner et al. 2015; Kalberla & Haud 2018; Marchal et al. 2019; Riener et al. 2019; Marchal & Miville-Deschênes 2021).

It is only recently with advances in absorption survey sensitivity and field-of-view that we have the number of sightlines to be able to stack Galactic H1 absorption spectra to improve their signal-to-noise. Using absorption sources from the first half of the 21-SPONGE survey, Murray et al. (2014) stacked the spectra from 19 sightlines with Milky Way H1 absorption to improve the optical depth sensitivity to the point of detecting the Milky Way's WNM in absorption ($\sigma_r \approx 2.6 \times 10^{-4}$). The spin temperature of this component was approximately $7200_{-1200}^{+1800}$ K. With the completed 21-SPONGE survey catalogue, Murray et al. (2018) conducted a similar stacking process, detecting a spin temperature of $\sim 10^4$ K. While this was one of the first steps in using stacking techniques with H1 absorption, the 21-SPONGE survey is limited both in the small number of sightlines used as well as sources being located vast distances from each other on the plane-of-sky, yielding an average of gas that may not physically be interacting or connected.

In this work, we make use of recent observations of Galactic absorption and emission in the direction of the Magellanic system ($-45^\circ < b < -25^\circ, 270^\circ < l < 300^\circ$). The gas of the Milky Way in this direction forms two distinct filamentary structures overlapping each other (Nguyen et al. 2024; Lynn et al., in preparation). The H1 survey in this region of the sky comes from the pilot phases of the Galactic Australian Square Kilometre Array Pathfinder (GASKAP) survey, observing the region with a spatial resolution of 30 arcsec (Dickey et al. 2013). Using the absorption pipeline for GASKAP described in Dempsey et al. (2022), the absorption survey searched in an unbiased manner across the full plane-of-sky survey range, producing a large number of 462 detections of Galactic H1 absorption.

While the individual detections sightlines are used for studying the physical and spatial properties of cold gas in the cirrus gas of the Milky Way (Nguyen et al. 2024; Lynn et al., in preparation), we instead focus on utilising stacking both the detections sightlines, as well as the non-detections of Galactic H1 absorption. Through stacking these two samples independently, we probe the global properties of the gas that lies beyond the sensitivity limits of individual spectra through GASKAP's current operational sensitivity, revealing the higher temperature phases of H1 gas in this region of the sky. The large number of the non-detections also allows for stacking sources based on their spatial distribution to investigate differences across the sky, a technique that was not possible previously.

In this paper, we present analysis of the optical depth and spin temperature of H1 gas in the Milky Way towards the Magellanic system and the spatial distribution of these properties through stacking the detection and non-detection sources from the GASKAP Pilot Phase II observations (Nguyen et al. 2024). The paper is organized as follows. Section 2.1 presents the data used in this work. In Section 3, we provide a theoretical background on stacking lines-of-sights with multiple Gaussian components and the consequences for extracted parameters. Section 4 investigates stacking the detection spectra after subtracting the already detected cold gas components from Nguyen et al. (2024). Section 5 then describes the stacking analysis of our GASKAP non-detection catalogue using (i) the whole sample and (ii) spatially distinct bins based on their proximity to known locations of cold gas. A discussion and summary are provided in Sections 6 and 7.

## 2 DATA

### 2.1 Observations

We made use of observations from the GASKAP-H1 Pilot Phase II survey using the Australian Square Kilometre Array Pathfinder (ASKAP) (Dickey et al. 2013). The observations were taken from 2021 October (SBID 33047) to 2022 April (SBID 38845) with each field gaining 10 h of exposure. The Pilot Phase II observations comprise ten 25 square-degree adjacent fields of view covering the entire Large Magellanic Cloud, Small Magellanic Cloud, as well as the whole Magellanic Bridge. Put together, the fields form a large mosaic from which we have analysed the background sources seen in absorption against the Milky Way foreground gas that was mapped at local velocities.

### 2.2 Source catalogues

#### 2.2.1 Detections

Nguyen et al. (2024) utilised the GASKAP absorption pipeline developed by Dempsey et al. (2022) to process 2714 sources of H1 absorption and their spectra. Twelve of these absorption spectra were saturated and removed from the sample. The remaining 2702 sources selected had a background source flux, $S_{\text{cont}}$, greater than 15 mJy. Of these spectra, 462 were deemed to contain detections of H1 absorption from gas pertaining to the Milky Way, defined as having a point above $3\sigma$ and a neighbouring channel above $2.8\sigma$. Nguyen et al. (2024) used Gaussian decomposition to then analyse the spin temperature and optical depth properties of the cold gas in these detection sightlines.

Relevant for this work, Fig. 1 (top panels) displays one of these absorption detections (bottom) against the radio source J040848–750719 and its associated emission spectra interpolated around the source (top). It is noteworthy that there are multiple peaks seen in absorption against this source. The blue and red lines show the modelled contribution from the cold and warm gas obtained following the methodology prescribed in Heiles & Troland (2003).

The mean rms noise of the detection sample is approximately $\sigma_r = 4.2 \times 10^{-2}$ and a median of around $3.1 \times 10^{-2}$, around 15 times higher than the noise level of individual channels in more sensitive observations like 21-SPONGE (Murray et al. 2014, 2018). However, the main advantage with the GASKAP survey is the order of magnitude increase in the number of sightlines available to stack.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top panels: Brightness temperature $T_B$ in units of kelvin of the absorption source J040848–750719 from the detection catalogue (top), and its corresponding absorption profile $e^{-\tau}$ (bottom). The coloured lines show the Gaussian decomposition as documented in Nguyen et al. (2024), with the fit of the Gaussian components found in absorption (dashed blue lines) along with their total contribution in emission (solid blue), an additional component added in emission tracing warm gas not found in absorption (dashed red) and the total fit of the emission profile (green). Bottom panels: Brightness temperature $T_B$ (top) and its corresponding absorption profile (bottom) of a randomly selected line-of-sight with no detected Galactic H I absorption.</figcaption></figure>

### 2.2.2 Non-detections

In this work, we also used the remaining 2240 sources where there were no detections of H I absorption with a signal-to-noise greater than 2.8 in multiple neighbouring channels. Hereafter, we will refer to them as the non-detection sample. An example of an absorption-emission pair selected randomly from the non-detection sample is shown in Fig. 1 (bottom panels). While the sensitivity of the absorption profile is too low to detect any features, the associated emission spectrum shows prominent peaks, one of which is relatively narrow, that reveal the presence of CNM gas mixed within the warmer phases of the ISM (i.e. the UNM and the WNM).

Table 1 lists the source name, Galactic coordinates l and b, average rms noise $\sigma_{\tau}$ over all velocity channels, and velocity at the peak brightness temperature $v_{\text{peak}}$ of ten sample entries of the 2240 non-detection source catalogue. The mean rms noise of the non-detections is $\sigma_{\tau} = 0.18$ with a median rms noise of approximately 0.14. The noise level of the GASKAP absorption sightlines is several orders

<table><caption>**Table 1.** Sample of GASKAP Pilot Phase II non-detection parameters.</caption><thead><tr><th>Source</th><th>l<br>(°)</th><th>b<br>(°)</th><th>σ<sub>τ</sub><br>(×10<sup>-3</sup>)</th><th>v<sup>b</sup><sub>peak</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J055932–634614</td><td>273.2</td><td>-29.8</td><td>244.5</td><td>2</td></tr><tr><td>J044109–631421</td><td>273.9</td><td>-38.5</td><td>158.0</td><td>12</td></tr><tr><td>J053738–624208</td><td>272.0</td><td>-32.3</td><td>93.3</td><td>6</td></tr><tr><td>J060128–624704</td><td>272.0</td><td>-29.6</td><td>204.9</td><td>-1</td></tr><tr><td>J053234–631640</td><td>272.7</td><td>-32.8</td><td>149.0</td><td>5</td></tr><tr><td>J060038–674747</td><td>277.8</td><td>-29.7</td><td>154.8</td><td>-5</td></tr><tr><td>J003524–732222</td><td>304.5</td><td>-43.7</td><td>116.5</td><td>0</td></tr><tr><td>J044926–691203</td><td>280.7</td><td>-35.9</td><td>88.5</td><td>3</td></tr><tr><td>J050107–691737</td><td>280.5</td><td>-34.9</td><td>190.4</td><td>4</td></tr><tr><td>J013229–734122</td><td>299.0</td><td>-43.1</td><td>238.4</td><td>1</td></tr></tbody></table>

Notes. $^a$Mean rms noise in optical depth over all velocity channels.

$^b$Line-of-sight velocity of the peak in brightness temperature.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Scatter plot of the background source flux, $S_{\text{cont}}$, versus the optical depth noise, $\sigma_{\tau}$, for each source of the 2240 non-detections (blue) and the 462 detections (orange) of Milky Way H I absorption. All sources, regardless of detection, have a background flux higher than the cut-off of 15 mJy (Dempsey et al. 2022).</figcaption></figure>

of magnitude higher than that of 21-SPONGE (Murray et al. 2018). However, similar to the detection sample, this noise level increase is countered by having two orders of magnitude more sources than that used in the stacking analysis in Murray et al. (2014, 2018).

### 2.2.3 Flux vs sensitivity

A comparison between the two data sets of detections and non-detections is shown in Fig. 2, comparing their background source flux and the mean rms noise of their absorption spectra. Detections and non-detections are shown in orange and blue, respectively. We note that the detections typically have higher background source fluxes, where there is a higher probability of the line-of-sight (LOS) producing a detection of absorbing gas. Most of the detections are within the same parameter space as the non-detections. However when the noise level is greater than $\sigma_{\tau} \sim 3 \times 10^{-1}$, the spectra are too noisy for any gas regardless of its opacity to produce a signal beyond the $\sim 3\sigma$ signal-to-noise requirement of a detection.

Since the non-detections spectra have similar background source fluxes to most of the detections, there is nothing intrinsically different about these spectra other than the higher noise levels. Thus, while the non-detection sources did not meet the criteria for being used in regular absorption analysis (Dempsey et al. 2022; Nguyen et al.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Spatial distribution of background sources from the GASKAP-HI absorption survey overlaid on a grey-scale map of the column density of H1 gas computed in the optically thin limit from the GASS survey (McClure-Griffiths et al. 2009; Kalberla et al. 2010; Kalberla & Haud 2015). Detections of cold gas in the velocity range $-50 < v < 50$ km s$^{-1}$ are annotated with white circles. Sightlines with non-detections of cold gas in absorption are shown with red circles. The blue dashed lines indicate the arcs of the vertical and horizontal H1 filaments in this region of the Milky Way (Lynn et al., in preparation).</figcaption></figure>

2024), they should nonetheless still contain astrophysical signals buried within their noise.

HI filamentary structures in the region, defined by high H1 column density and high-brightness temperature (Lynn et al., in preparation).

## 2.3 Emission data

To pair these absorption sightlines with emission data, we made use of the H1 emission data provided by the GASKAP Pilot surveys as mentioned in Section 2.1. The emission data have a beamsize of 30 arcsec and native spectral resolution of 0.24 km s$^{-1}$. For noise reduction, the emission and absorption spectra were binned spectrally by a factor of 4 to a spectral resolution of $\sim 1$ km s$^{-1}$. As the emission sightline located at the position of H1 absorption is affected by the absorbing feature, emission spectra from nearby sightlines on the plane-of-sky were used instead, assuming that the gas within this distance is similar. These emission sightlines also have to be far away enough from each other that beam smearing does not affect their average. For this purpose, we used a method similar to that of Nguyen et al. (2024). For each absorption sightline, approximately 20 emission sightlines were selected around the absorption feature, within an annulus of inner-radius 1 beam-radius and outer-radius 2.5 beam-radii, ensuring each is at least 1 beam-radius from other sightlines. These emission sightlines are then averaged, weighted by the inverse of their plane-of-sky distance from the absorption sightline at the centre.

Fig. 3 shows the locations of the absorption sightlines, including both detections and non-detections at Milky Way velocities in the velocity range $-50 < v < 50$ km s$^{-1}$ overlaid onto the column density map of the region taken from the Galactic All-Sky survey computed in the optically thin limit (GASS, surveyed with the Parkes telescope with a beam size of 16 arcmin and 1 km s$^{-1}$ spectral resolution; McClure-Griffiths et al. 2009; Kalberla et al. 2010; Kalberla & Haud 2015). The column density map highlights the vertical and horizontal

## 3 STACKING CONSIDERATIONS

When stacking multiple spectra to produce an approximate average spectra over a region, one needs to take into consideration the effects of averaging multiple Gaussian components that are potentially misaligned in velocity. In this section, we explore, through use of the GASKAP detection sample, the biases introduced by mixing complex multiphase emission and absorption spectra through stacking.

### 3.1 Velocity alignment

When stacking spectra, the first thing we must consider is how to account for effects like Galactic rotation, or bulk motions of gas that could displace components in line-of-sight velocity. For this work, we explore shifting the emission and absorption spectra such that the peak brightness temperature is centred on 0 km s$^{-1}$. The peak in emission typically traces cold gas and so shifting to centre this ensures that at least the bulk of H1 cold gas is centred on the same velocity when stacking. With sightlines containing single components, this method of shifting perfectly aligns all components to produce a true average spectrum. We note that in reality some minor deviations are present between the peak in emission and absorption that are typically on the order of 1–2 km s$^{-1}$, as observed in the GASKAP detections (Nguyen et al. 2024). These deviations are due to a combination of the difference in noise-levels between the emission and absorption spectra as well as the emission and absorption spectra probing slightly different gas as they are not along the same sightline. These effects cause some uncertainty in the
central velocity of the peak of absorption and its alignment with the peak in emission. However, the 1–2 km s⁻¹ deviation is negligible compared to the ~10–15 km s⁻¹ variation in central velocities of cold components in this region of the sky (Nguyen et al. 2024).

Other methods in the past have utilised shifting by the first velocity moment (Murray et al. 2014, 2018) producing similar results with single components as this simplifies to centring by the peak. However, the two shifting methods deviate when we consider sightlines with multiple components. As shown in the example detection spectra of Fig. 1, the most sensitive detection sources showcase multiple components: a narrow primary component aligned well to the peak in the respective emission profile and a slightly broader secondary component offset in velocity. Therefore, velocity shifting needs to be able to accommodate for the existence of multiple clouds along a sightline. In a multiphase medium with multiple components per sightline, none of the components along a line-of-sight will be perfectly centred to 0 km s⁻¹ if shifted by the first velocity moment. When stacking these partially shifted spectra, the key features of colder and warmer phases will be heavily mixed together.

Since in a general stacking case, we are unsure of the number of components in each sightline prior to stacking, shifting by the first velocity moment will make it difficult to separate out the phases of H1 through Gaussian decomposition. If we shift instead by the peak in emission, we avoid these problems, as we standardise the individual sightlines such that the bulk of cold gas (traced by the peak in emission and henceforth referred to as the 'primary component') is centred to the same velocity. While additional components at other local maxima will be mixed together since their central velocities are not accounted for (producing what will be referred to as the 'secondary component'), shifting in this manner allows for easier separation of the phases of H1 with a major part of the cold phase aligned in velocity at 0 km s⁻¹ prior to stacking.

## 3.2 Stacked spectrum

The GASKAP detection sample provides a unique opportunity to test how stacking affects measured properties, given that the sightlines are all in the same plane-of-sky location, and sightlines contain multiple cold components (Nguyen et al. 2024).

Our stacking analysis follows the procedure described in Murray et al. (2014), where they formed a stack of 19 absorption spectra from the 21-SPONGE survey, after subtracting Gaussian fits for cold gas, to investigate properties of the WNM. For this section, however, we are using the un-subtracted detection spectra as observed, containing all the phase information of the ISM to compare the stack with the average properties of the individual sightlines.

The stacked spectrum was obtained by taking the average of the sample weighted by the noise of each absorption spectrum, with weights $\tau_{res}/\sigma_{\tau}^2$ (Treister et al. 2011) to favour spectra with lower noise. The weights simplify to $1/\sigma_{\tau}$ since $\sigma_{\tau} \propto T_B$ and for constant temperature, $\tau_{res} \propto T_B$ (Murray et al. 2014). The same weights were applied in the stack of the emission spectra. Fig. 4 shows the stacked emission and absorption spectra (top and bottom panel, respectively). The absorption stack produces a two-component profile; a deep, narrow, primary component, and a broad secondary component.

## 3.3 Gaussian decomposition

To properly Gaussian decompose the stacked spectra, we implemented the method prescribed in Heiles & Troland (2003) to model their spin temperatures individually. Two Gaussian components were

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top: Weighted-mean emission profile of the stacked GASKAP detection spectra. The original stacked spectrum is shown in black and the model of the joint Gaussian decomposition in green. The model is made by the addition of the emission produced by the narrow and broad components (shown in blue), and an additional broad component required to complete the fit (shown by the red dashed line). Middle: Residual profiles of the emission (middle top) and absorption (middle bottom) fits. Bottom: Weighted-mean absorption profile of the stacked GASKAP detections. The original stacked spectrum is shown in black. The narrow and broad components' absorption fits are shown by the blue dashed lines, with their total contribution to the absorption profile shown in the solid blue line.</figcaption></figure>

fitted to the absorption stacked spectrum. The fitting parameters were then used to model the amplitude of their emission counterparts, allowing for a ±1 km s⁻¹ deviation in their central velocities and 10 per cent deviation in their widths. An additional broader Gaussian component was added to fit the signal in emission to ensure the column density of the fit was not lower than that of the original data from GASKAP, and that the correction factor, $R = N_{HI}/N_{HI}^*$, was greater than one, where * denotes the column density computed in the optically thin regime. The additional component in emission represents gas from the WNM that is not able to be detected in absorption due to the insufficient optical depth sensitivity of our stack. The resulting decomposition is displayed in Fig. 4. It should be noted that a potential bias in adding only a single emission-only component is that it under-represents the spin temperature of the absorption feature. Adding additional components to fit the same data would lead to them being narrower in width and shallower in amplitude, allowing the absorption components to more closely fill the emission at the peak. This higher amplitude in brightness temperature would lead to a higher recorded spin temperature.

When solving the radiative transfer equation of the line, we accounted for each possible permutation of the broader emission component being in either the foreground or background relative to the colder components seen in absorption along the LOS. For the fraction of warm gas in the background (and hence absorbed by the cold gas infront of it), we used three values: 0, 0.5, and 1. Averaging these solutions to the radiative transfer yielded the mean spin temperatures. The primary and secondary component of the total stacked spectrum using all 462 detection sources have a mean spin temperature of $T_s = 68 \pm 7$ K and $159 \pm 58$ K, peak optical depth of $\tau_{peak} = 0.160 \pm 0.033$ and $0.040 \pm 0.007$, and full-width half-maximum (FWHM) of $3.4 \pm 1.4$ km s⁻¹ and $12.5 \pm 1.1$ km s⁻¹, respectively. The uncertainties are produced through a method of bootstrapping,
where new samples of 462 detection spectra were sampled from the original 462 sightlines with replacement. Repeating the stacking and fitting process of the new sample, and repeating this 10⁵ times, the standard deviations of the above measurements of $T_s$, $\tau_{peak}$, and FWHM were computed.

## 3.4 Comparison of stack versus original sample

To examine the effect of stacking spectra on the measurements properties of spin temperature, optical depth and FWHM, we can compare these values to that of the individual detection sightlines, using the properties extracted by Nguyen et al. (2024). Following the discussion in Section 3.1, absorption detection components located closest spectrally to the peak in emission in each detection sightline are considered primary components. The remaining absorption components are collected as the secondary components in each spectra. Fig. 5 shows the comparison between the observed peak $\tau$, FWHM, and spin temperatures of both stacked components and the distribution of the individual sightlines' primary and secondary components. We generate the weighted mean of these individual sightline properties using the same weighting of the inverse of the optical depth noise that produced the stacks. The uncertainty in these means were calculated by adding the uncertainty of the mean and the systematic errors in quadrature. The systematic errors were approximated as the mean of the uncertainties of individual sightline values in $\tau_{peak}$, FWHM and $T_s$. When considering the overlap of the errors in the stacked values and that of the mean individual values, we find that the primary component's optical depth and FWHM are both accurate representations of an average value. However, when it comes to the secondary component, the observed peak optical depths lie below the mean of the individual sightline peak $\tau$ measurements that produced them. Conversely, the FWHM measurements observed in the secondary component's stack lies above the mean FWHM of the individual sightline components.

The reason for this discrepancy between the secondary component's stacked values compared to the mean properties of the individual sightlines comes from the velocity offset of components. Since the components are slightly offset in velocity from each other, the width of the resulting stacked Gaussian will be wider than the average of the individual Gaussians. Due to the total flux of all the components being conserved in the average stack, if the width of the Gaussians is increased, to compensate the amplitude of the resulting stacked Gaussian must diminish. Therefore, regardless of the severity of the range of central velocities of secondary components in individual sightlines, the peaks will always be lower limits and their widths always greater than the average of the individual components themselves.

The primary components on the other hand, aligned with the peak in emission, are shifted correctly to roughly 0 km s⁻¹ with some minor deviation due to possible slight offset between absorption and emission peaks in velocity (as mentioned in Section 3.1). As such the primary component's lack of a noticeable velocity offset leads to the stacked result representing a relatively accurate average of the peak optical depth and FWHM.

When calculating the spin temperatures of the components, the primary component's stack and mean values are well aligned along with their uncertainty envelopes. The secondary component's values are also similar. While their uncertainty envelopes do not perfectly align, it should be noted that due to the potential biases in spin temperature by only including one emission component (as mentioned in Section 3.2), it is likely we are undersampling the error in the observed stacked temperature. Therefore, we conclude that both

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Distributions of the peak optical depth (top), FWHM (middle), and spin temperature (bottom) of individual absorption components found by Nguyen et al. (2024). The components are separated into primary components (blue), which are aligned with the peak in their corresponding emission spectra, and secondary components (red), containing additional components offset from the emission peak. The solid vertical lines represent the mean of these distributions in their respective colours. The dotted lines highlight the values obtained from the detection stack in Fig. 4 with the shaded regions representing the uncertainties of both the mean and observed stacked values.</figcaption></figure>

the primary and secondary components' spin temperatures are well aligned with the means of the original distributions. This is due to the fact that since the velocity offset of components is equal in both emission and absorption, the peaks of emissions and absorption are affected in a similar way. Since the spin temperature is found through the division between peaks in brightness temperature and optical depth, the velocity offset effect is cancelled out, resulting in a spin temperature that is well correlated to the true average temperature of the individual sightlines' components. An exploration into this velocity offset effect is explored in Appendix A where toy models are used instead of real detection data, yielding similar results.
## 3.5 Theoretical implications

From stacking the detection sample and comparing the result with the individual sightline properties (as well as through the use of toy models in Appendix A), we find several effects to be considered.

The uncertainty of the true peaks and widths of the stacked components means we should be cautious when calculating secondary information extracted from emission-absorption pair stacks that use these values such as the FWHM, corresponding turbulent Mach numbers, and mass fractions, as they become unreliable.

This is especially true for extragalactic observations, where individual sightlines' beamsizes pass through much larger portions of a galaxy than Galactic observations. The large amount of extragalactic ISM material mixed together at various velocities will broaden the Gaussian components in a similar manner to that of the secondary component explored in this section, with individual sightlines acting as a pseudo-stack of the region. As a result, a similar effect will occur and peak optical depth measurements of extragalactic sources should take this lower-limit effect into account when determining results.

The rest of this work will calculate and report parameters including the peak optical depth, FWHM and spin temperature. However, other quantities such as the turbulent Mach number, column density and mass fraction will not be discussed as they develop large uncertainties, as mentioned above.

# 4 STACKING MODEL-SUBTRACTED DETECTIONS

In this section, we investigate the existence of higher temperature H1 gas in the region of the GASKAP H1 Pilot survey by stacking the GASKAP detection spectra after subtracting the cold gas already detected from Nguyen et al. (2024).

## 4.1 Model subtraction

Higher temperature phases of H1 tend to have weaker signals in absorption due to lower optical depths. In order to detect these higher temperatures in absorption, we must first subtract the cold gas already detected from individual sightlines. For each of the 462 detection sightlines, Nguyen et al. (2024) catalogues the fitted model parameters of the cold gas detected in absorption and emission. The modelled Gaussian absorption and emission components, $τ_i(v)$ and $T_{B,i}$ (where i indicates the ith component) were subtracted from the respective spectra from which they were detected to produce 462 residual absorption spectra, $τ_{res}(v) = τ(v) - Στ_i(v)$, and residual emission spectra $T_{B, res}(v) = T_B(v) - ΣT_{B,i}(v)$.

## 4.2 Stacked spectrum

Since we have removed the majority of the signal through subtracting the cold gas components, we no longer have information about the velocity of the remaining underlying absorption hidden beneath the noise in individual detection spectra. We follow the same procedure as in Section 3.1, relying on the GASKAP emission data to align the spectra before performing the stack to account for Galactic rotation and large scale bulk motions in the region.

To produce the stacked model-subtracted emission and absorption spectra, we use the same weighting as in Section 3.2 following Murray et al. (2014). The resulting stacked emission and absorption spectra are shown in Fig. 6. The rms noise of the stacked spectrum is $σ_τ = 1.1 × 10^{-3}$, calculated from taking the standard deviation of 30 empty channels in the post-shifting velocity range $-50 < v < 0$.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Top: Weighted-mean emission profile of the stacked GASKAP detection spectra after subtracting the cold gas models in individual sightlines from Nguyen et al. (2024). The original stacked spectrum is shown in black and the model of the joint Gaussian decomposition in green. The model is made by the addition of the emission produced by the absorption component (shown in blue), and an additional broad component required to complete the fit (shown by the red dashed line). Middle: Residual profiles of the emission (middle top) and absorption (middle bottom) fits. Bottom: Weighted-mean absorption profile of the stacked GASKAP non-detections. The original stacked spectrum is shown in black. The absorption fit is shown by the solid blue line.</figcaption></figure>

$-20 \text{ km s}^{-1}$. The stacked detection spectra's noise level is almost 40 times smaller than the mean rms noise of the individual detection sightlines ($σ_τ = 4.2 × 10^{-2}$) and nearly 30 times smaller than their median rms noise ($σ_τ = 3.1 × 10^{-2}$). The stacked absorption profile has an equivalent width (EW) in the velocity range $-50 < v < 50 \text{ km s}^{-1}$ of EW = $\int_{-50}^{50} (1 - e^{-τ}) dv = 0.12 \text{ km s}^{-1}$ (Murray et al. 2014).

To ensure that our stack is not plagued by any outlying spectra and the detected component of the stack is not produced by only a few dominant narrow profiles in the average, we conducted two bootstrapping Monte Carlo simulation experiments (Wall & Jenkins 2003; Murray et al. 2014) to test for variations in the stack's EW. First, we repeated the stacking analysis on a new sample of 462 sources randomly chosen from our original 462 sightlines with replacement. Repeating this analysis $10^5$ times, we computed the EW for each trial. The distribution of these EW values is shown in Fig. 7. The distribution is Gaussian, with a mean value of $\langle EW \rangle = 0.12 ± 0.02 \text{ km s}^{-1}$ (68 per cent confidence), suggesting that the EW computed from the original sample is consistent with being drawn from a parent population of spectra with a comparable mean and is not influenced by outlying spectra. Secondly, we produced a similar EW distribution where in each trial we inverted (multiplied by -1) half of our randomly selected sample. The distribution, also shown in Fig. 7, is centred around $0 \text{ km s}^{-1}$, consistent with no signal. The result verifies that our stack does not produce a false signal.

We used the method described in Section 3.2 to decompose the stacked spectrum and determine the properties of the remaining gas. The stacked spectrum produces a single wide component in absorption, with a peak optical depth of $τ_{peak} = (6.6 ± 0.5) × 10^{-3}$, FWHM of $16.3 ± 1.4 \text{ km s}^{-1}$, and $T_s = 1320 ± 263 \text{ K}$. The uncertainties in these values were evaluated by repeating the same fitting procedure to each of the $10^5$ trials of the bootstrapping sample used earlier.
<figure><img src="image_8.png" /></figure>

**Figure 7.** Histograms of the EW distribution from $10^5$ trials of bootstrapping, randomly sampling 462 spectra from the model-subtracted detection sample derived in Section 4.1, with possibility of replacement. The blue histogram displays the distribution of EW values from the trials, with a mean of $\langle \text{EW} \rangle = \left\langle \int_{-50}^{50} (1 - e^{-t}) dv \right\rangle = 0.12 \pm 0.02 \text{ km s}^{-1}$ (68 per cent confidence). The orange histogram displays the same random trials except with half of the spectra inverted (multiplied by -1) before stacking to test for contamination from outlying spectra. The mean of the orange distribution is $0.00 \pm 0.02 \text{ km s}^{-1}$ (68 per cent confidence). The vertical dotted line crossing the middle of the blue histogram represents the EW of the original stack (EW= 0.12 km s$^{-1}$). The dotted line along the middle of the orange histogram represents the expected mean of 0 km s$^{-1}$ for the half-inverted stack EW. The additional dashed lines mark the 68 per cent confidence interval for the two distributions.

These uncertainties are larger than those from the variation of the fraction of background warm gas in the radiative transfer solution.

# 5 STACKING NON-DETECTIONS

Stacking the detection spectra is useful for investigating higher temperature gas in sightlines where we know cold gas already exists. However, the GASKAP survey has the unique opportunity to use 2240 non-detection spectra in the same region on the sky that likely also contains information about the multiphase ISM buried within their noise. In this section, we investigate stacking the GASKAP non-detection spectra to reveal more information about cold gas in this region of the Milky Way.

## 5.1 Velocity alignment

Since there is no information about the velocity of the underlying absorption line hidden beneath the noise in individual non-detection spectra, we rely on the peak in the GASKAP emission data to align them before performing the stack to account for Galactic rotation and large scale bulk motions in the region (as explained in Section 3.1).

## 5.2 Stacked spectrum

For the non-detection stack, we follow the same stacking analysis used in previous sections for the detection sample in Section 4.2. The top and bottom panels in Fig. 8 show the stacked emission and absorption spectra of the non-detection sample, respectively. The high sensitivity of the absorption stack reveals a two-component profile; a deep narrow primary component resembling absorption profiles of the CNM seen in the detection catalogue, and a broad secondary component.

The rms noise of the stacked absorption profile is $\sigma_{\tau} = 2.0 \times 10^{-3}$, calculated by taking the standard deviation of 30 empty channels in the post-shifting velocity range $-50 < v < -20 \text{ km s}^{-1}$. The noise

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

**Figure 8.** Top: Weighted-mean emission profile of the stacked GASKAP non-detections. The original stacked spectrum is shown in black and the model of the joint Gaussian decomposition in green. The model is made by the addition of the emission produced by the narrow and broad components (shown in blue), and an additional broad component required to complete the fit (shown by the red dashed line). Middle: Residual profiles of the emission (middle top) and absorption (middle bottom) fits. Bottom: Weighted-mean absorption profile of the stacked GASKAP non-detections. The original stacked spectrum is shown in black. The narrow and broad components' absorption fits are shown by the blue dashed lines, with their total contribution to the absorption profile shown in the solid blue line.

level has decreased from the original mean of $\sigma_{\tau} = 1.7 \times 10^{-1}$ for the non-detections by nearly two orders of magnitude. The absorption stack's noise level is also approximately 15 times lower than that of the median noise level ($\sigma_{\tau} = 3.1 \times 10^{-2}$) of the single lines-of-sight classified as detections in the GASKAP-HI absorption survey (Nguyen et al. 2024). The stacked absorption profile has an equivalent width (EW) in the velocity range $-50 < v < 50 \text{ km s}^{-1}$ of EW = $\int_{-50}^{50} (1 - e^{-t}) dv = 0.44 \text{ km s}^{-1}$ (Murray et al. 2014).

To ensure that our stack is not plagued by any outlying spectra and that the primary component of the stack is not produced by only a few dominant narrow profiles in the average, we conducted the same two bootstrapping Monte Carlo simulation experiments used in Section 4.2. We completed the stacking analysis on a new sample of 2240 sources randomly chosen from our original 2240 sightlines with replacement. The process was repeated $10^5$ times, with the EW computed in each trial. The distribution of EW values is shown in Fig. 9. The distribution is Gaussian, with a mean value of $\langle \text{EW} \rangle = 0.41 \pm 0.06 \text{ km s}^{-1}$ (68 per cent confidence), suggesting that the EW computed from the original sample is consistent with being drawn from a parent population of spectra with a comparable mean and is not influenced by outlying spectra.

We produced an additional EW distribution where in each trial we inverted (multiplied by -1) half of our randomly selected sample. The distribution, also shown in Fig. 9, is centred around 0 km s$^{-1}$, consistent with no signal. From this result, we conclude our stack does not produce a false signal.

Since there are multiple components within this stacked absorption spectrum, we followed the method described in Section 3.2 to decompose the spectra and calculate the gas' properties, tabulated in Table 2. The primary and secondary component of the total stacked spectrum using all 2240 sources have a mean spin temperature of $T_s = 98 \pm 12 \text{ K}$ and $T_s = 255 \pm 106 \text{ K}$, respectively. The uncertain-
<figure><img src="image_10.png" /></figure>

**Figure 9.** Histograms of the EW distribution from $10^5$ trials of bootstrapping, randomly sampling 2240 spectra from the total non-detection sample with possibility of replacement. The blue histogram displays the distribution of EW values from the trials, with a mean of $\langle \text{EW} \rangle = \langle \int_{-50}^{50} (1 - e^{-\tau}) dv \rangle = 0.41 \pm 0.06 \text{ km s}^{-1}$ (68 per cent confidence). The orange histogram displays the same random trials except with half of the spectra inverted (multiplied by -1) before stacking to test for contamination from outlying spectra. The mean of the orange distribution is $0.00 \pm 0.05 \text{ km s}^{-1}$ (68 per cent confidence). The vertical dotted line crossing the middle of the blue histogram represents the EW of the original stack (EW= 0.44 km s$^{-1}$). The dotted line along the middle of the orange histogram represents the expected mean of 0 km s$^{-1}$ for the half-inverted stack EW. The additional dashed lines mark the 68 per cent confidence interval for the two distributions.

ties were evaluated by repeating the same fitting procedure on the $10^5$ trials of the bootstrapping sample used above. These uncertainties are larger than those from the variation of the fraction of background warm gas in the radiative transfer solution.

## 5.3 Spatial distribution of stacked features

Given the large number of sources that have been obtained within GASKAP's relatively narrow field of view of the Magellanic system, we have the unique opportunity to search for any spatial variation in the absorption stack profile. In this work, we have chosen to bin the data by increasing source density of detections, allowing us to assess how the profile's peak optical depth and spin temperature change with distance from known locations of CNM gas. Since regions of

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

**Figure 10.** A 2-dimensional contour plot of the source density of Milky Way H1 detections found through the GASKAP Pilot II survey towards the Magellanic system (Nguyen et al. 2024). The contours are created based on a kernel density PDF using the spatial location of the detected sources (shown with white circles). The levels of the contours separate the 6 equally sized bins of non-detection sources, highlighting regions of higher absorption source density on the plane-of-sky towards the centre of the H1 filamentary structure observed in Fig. 3. The dotted blue lines indicate the H1 filaments from Fig. 3 for context.

cold gas typically show higher column densities, this also acts as a probe for binning by column density (as shown in Table 2).

### 5.3.1 Binning by source density of detections

One parameter of interest is how the profile's peak optical depth and spin temperature change with distance from known locations of CNM gas. To probe this, we first used the locations of the 462 H1 detections of Milky Way absorption in this region (Nguyen et al. 2024) to create a kernel probability distribution function (PDF) of known cold gas. We convert the density values of the PDF into source density of detections in units of sources per square-degree. We attribute to each non-detection a value corresponding to the source density grid at the non-detection's spatial location.

The non-detections were then binned into six equally sized populations based on their respective PDF values. The six regions created from this binning process are displayed via the contour map in Fig. 10, overlaid by the location of absorption detections shown in

<table><caption>**Table 2.** Fitted parameters for the total non-detection stack and spatially binned stacks.</caption><thead><tr><th>Bin</th><th>N<sup>a</sup></th><th>n<sup>b</sup><br>(deg<sup>-2</sup>)</th><th>N<sub>tot</sub><br>(×10<sup>20</sup> cm<sup>-2</sup>)</th><th>σ<sub>τ</sub><sup>c</sup><br>(×10<sup>-3</sup>)</th><th>Peak τ<br>(×10<sup>-3</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>T<sub>s</sub><br>(K)</th></tr></thead><tbody><tr><td>Total</td><td>2240</td><td>0.97±0.76</td><td>4.2</td><td>2.0</td><td>60±10</td><td>4.0±1.1</td><td>98±12</td></tr><tr><td>Non-detections</td><td></td><td></td><td></td><td></td><td>>20±10**</td><td><14.2±1.9**</td><td>255±106</td></tr><tr><td>1</td><td>374</td><td>0.08±0.06</td><td>2.5</td><td>7.0</td><td>-*</td><td>-*</td><td>-*</td></tr><tr><td>2</td><td>373</td><td>0.32±0.09</td><td>3.1</td><td>7.8</td><td>46±20</td><td>3.0±1.5</td><td>75±26</td></tr><tr><td>3</td><td>373</td><td>0.66±0.11</td><td>4.0</td><td>7.8</td><td>>17±6**</td><td><13.5±3.6**</td><td>306±117</td></tr><tr><td>4</td><td>373</td><td>1.02±0.12</td><td>4.8</td><td>7.8</td><td>62±27</td><td>4.1±1.8</td><td>77±26</td></tr><tr><td>5</td><td>373</td><td>1.48±0.15</td><td>5.5</td><td>8.1</td><td>>40±15**</td><td><10.2±4.3**</td><td>135±215</td></tr><tr><td>6</td><td>374</td><td>2.28±0.34</td><td>6.2</td><td>12.2</td><td>187±16</td><td>6.1±0.3</td><td>88±6</td></tr></tbody></table>

Notes.<sup>a</sup> Number of non-detection sources.

<sup>b</sup> Mean plane-of-sky number density of detection sources.

<sup>c</sup> Mean rms noise over 30 velocity channels from approximately -50 to -20 km s<sup>-1</sup>.

*Due to the low signal-to-noise of the data in bin 1, a fit could not be made.

**As discussed in Section 3, optical depths and FWHM of the secondary components in shifted stacks are lower limits and upper limits respectively.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Spatial distribution of the non-detection sources coloured based on their respective contour level from the contour regions shown in Fig. 10. We define the bins 1–6 starting from the outer regions of the survey, with higher bin numbers probing regions with higher source density of H1 absorption, as highlighted with the locations of the 462 Milky Way H1 absorption detections marked with white circles. The background is a grey-scale of the Milky Way H1 column density in the region of the sky from the GASS survey (McClure-Griffiths et al. 2009; Kalberla et al. 2010; Kalberla & Haud 2015).</figcaption></figure>

white used to calculate the source density map and the filamentary structures marked in Fig. 3 for context.

The spatial distribution of the individual non-detection sources after binning is displayed in Fig. 11, overlaid onto the column density map from the GASS H1 emission survey. Binning using this source density metric highlights a concentric diamond-configuration centred on the Galactic filament cross, marking it as a region densely populated by CNM detections. With each of these bins of non-

detections, we repeated the same methodology of stacking their absorption and emission spectra and fitting them using Gaussian decomposition as done in Sections 4.2 and 5.2 to extract their physical properties. We repeated this using the same Monte Carlo bootstrapping process as done previously to determine uncertainties.

The Gaussian decomposition of each bin is shown in Fig. 12. Table 2 displays the physical properties obtained including the total column density of the bins, $N_{\text{tot}} \approx 1.823 \times 10^{18} \int T_{\text{B}} \, \text{dv}$, peak optical depth, spin temperature and FWHM of each component. We note that with each bin getting closer to regions of cold gas (bin 1 to bin 6 increasing in source density), the total column density also increases, allowing us to also probe the variation of the absorption features with column density simultaneously.

The primary component's peak optical depth diminishes from approximately $\tau_{\text{peak}} = (187 \pm 16) \times 10^{-3}$ to $(46 \pm 20) \times 10^{-3}$ as the bins of non-detections move outwards, away from regions of known cold gas, before disappearing in the bin most distant from the GASKAP detections. This trend is supported further with the left plot of Fig. 13 highlighting the variation of peak $\tau$ of the primary component in each bin, as it decreases the further we move from the filament structure. The peak optical depth values from the detection sample stack in Section 3.2 is also displayed, representing the closest available bin in our GASKAP sample to cold gas. We find the detection stack value agrees with the trend of the primary components in the non-detection bins.

The secondary component values on the other hand appear to plateau in peak optical depth, with the value in bin 4 being relatively similar to that of the detection sample stack in the central region of the filament structure. While this is a lower limit of the true average (as discussed in Section 3), the plateau is still noteworthy considering

<figure><img src="image_13.png" /><figcaption>**Figure 12.** H1 emission and absorption stacked profiles of the spatially binned groups of spectra from Fig. 11 and their Gaussian decompositions using the method from Heiles & Troland (2003). The layout and colour scheme used is similar to that used in Fig. 8. We note a secondary component is detected in bins 2 and 4, whereas the primary component exists in all bins except for bin 1, increasing in amplitude in absorption from bins 2 to 6 with increased detection source density. Due to the signal-to-noise being inconsistent between bins, bin 1 could not be fitted with a statistically significant Gaussian component.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Properties of the Gaussian decompositions of the stacked profiles binned in Section 5.3 (bin 1 not shown). Left: Plot of the peak $\tau$ values for the primary components (blue circles) and secondary (red squares) components from each of the stacked absorption profiles. The primary component increases in peak $\tau$ as the stacked spectra get closer to regions of known cold gas, and increase at a faster rate relative to the secondary component, whose lower limit remains relatively stable in peak optical depth. The blue and red shaded regions represent the peak optical depth of each component in the total stack from Fig. 8 within $1\sigma$. The optical depths of the detection sample stack from Section 3.2 are shown to the right, and agree with the trends from the non-detection bins as they approach colder gas. Right: Plot of the FWHM for the primary (blue circles) and secondary (red squares) components from each of the bins’ stacked profiles. The FWHM of the primary component in each bin is consistent with that of the total stack, shown in the shaded region. The secondary component’s upper limit of the FWHM is also relatively consistent, except for the outlier in bin 4. The FWHM’s of the detection sample stack from Section 3.2 are shown to the right, and agree with values of those from the non-detection bins.</figcaption></figure>

these lower limits are linearly correlated with their true counterparts (as shown in Appendix A).

We note that due to variation of signal-to-noise between each bin, bin 1's signal-to-noise was too low to consistently Gaussian decompose the stack during the Monte Carlo process used in this work. As a result, uncertainties could not be calculated for this bin. Therefore, we treat bin 1 as a non-detection of HI absorption.

The right panel of Fig. 13 shows the FWHM of the two absorption components in each bin. We note the primary component's FWHM is consistent with that of the total stack. The secondary component's FWHM, while an upper limit, also seems consistent with the value of the total stack (except for an outlier in bin 4). With the flux of these stacked components being conserved, the FWHM upper limit is related to the lower limit of the peak optical depth. As such, this implies that the FWHM upper limit is also linearly correlated with the true average FWHM (similar to the optical depth as mentioned earlier). Therefore, while we can not confirm that the true FWHM in each bin is similar to the total stack, we can conclude that the bins' values are at least similar to each other.

Fig. 14 displays the variation of $T_s$ for the primary and secondary components with distance from cold gas. We find that both the primary and secondary components remain relatively consistent in temperature regardless of distance or changes in detection source density, and notice no overall trend. There is however a minor decline in the secondary component's spin temperature as they approach the value of the detection stack from Section 3.2. Regardless, the spin temperatures of the two components still align well with those of the total stack, with most bins' components falling within $1\sigma$ (68 per cent confidence) of the total stack's values. One might argue that the stability in spin temperature could be caused by the additional component added in emission during the fitting process accommodating any changes in the absorption properties. However, since the linewidths of the broad and narrow components (see values tabulated in Table 2), and hence their $T_{k,max} = 121\sigma_v^2$ remain relatively consistent, similar to their spin temperatures, it is likely

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Variation of the spin temperature values for the primary (blue circles) and secondary (red squares) components from each of the stacked profiles binned in Section 5.3 (bin 1 not shown). The black line denotes the cutoff between CNM and UNM used in this work of ~250 K. The primary component remains relatively consistent in spin temperature regardless of the region being probed. Both components align with their counterparts in the total stack, shown in the shaded regions with their widths indicating the $1\sigma$ standard deviation. The spin temperatures of the detection sample stack from Section 3.2 are shown to the right, and align with those from the non-detection bins.</figcaption></figure>

that this trend of constant spin temperature is physical and not a byproduct of the fitting process.

## 6 DISCUSSION

Previous absorption surveys have typically used deep observations in single LOS locations (Carilli et al. 1998; Dwarakanath et al. 2002; Patra et al. 2018) or observations in a widespread field (Heiles &
Troland 2003; Murray et al. 2018). In this work, however, we have taken advantage of the unbiased observational capabilities of ASKAP to observe the spatial distribution of the properties of cirrus H1 gas in the relatively small region of the Milky Way toward the Magellanic system.

## 6.1 Spin temperature

In Section 4, stacking the residual detection spectra reveals a broad component with a spin temperature of 1320 ± 263 K. This is representative of gas from the UNM, and is much higher than the highest temperature detected by Nguyen et al. (2024) from the 462 individual detection sightlines. Since this is a stacked average of the detection spectra, which are sampled in an unbiased manner across the spatial span of the GASKAP survey range, this component represents an average of the UNM temperature within this region of the sky. The component's spin temperature is consistent with average UNM temperatures found from previous surveys that observed relatively diffuse gas (~800 ± 570 K from Heiles & Troland 2003 and ~600 ± 480 K from Murray et al. 2018).

The fact that this broader UNM component appears simultaneously in the same region of sky as the cold components found in both Nguyen et al. (2024) as well as this work supports previous findings that CNM structures are typically surrounded by a broader warmer feature such as colder UNM (Kalberla & Haud 2018) and that they are interlinked as part of a continuous range of H1 properties, with the UNM being a transitional stage between the stable phases.

From the spin temperatures of the components in the total non-detection H1 stack in Fig. 8, we are able to discern that both the primary and secondary components with temperatures of 98 ± 12 and 255 ± 106 K originate from gas belonging to the CNM ($T_s$ <~ 250 K). We note that both values are higher than the mean spin temperature of ~50 K for the components found in single lines-of-sight (Nguyen et al. 2024). This discrepancy is likely due to the sightlines used being non-detections and as such did not contain deep absorption profiles of colder denser gas like those classified as detections. However, only a few instances of warmer gas still within the regime of CNM and under 250 K were found in the detection survey. This upper regime of cold gas has optical depths beyond the sensitivity of the current Pilot survey and as such were only found through stacking in this work to reduce the signal-to-noise. The primary component's spin temperature also aligns with the average CNM spin temperature of 80 ± 60 K from the Millenium survey (Heiles & Troland 2003) and 71 ± 52 K from the more recent 21-SPONGE survey (Murray et al. 2018) albeit slightly warmer. This indicates that while probing more diffuse gas than typical past H1 absorption surveys, the CNM spin temperatures remain relatively consistent.

The secondary component's higher spin temperature of a warmer CNM gas seems to indicate that the sensitivity reached by the stacked spectrum is high enough to detect and quantify the properties of warmer cool gas throughout the entire region of the GASKAP survey. As shown in Section A2, while the observed spin temperature of primary and secondary components is well correlated with the true average, we do not have any information as to the dispersion of the individual spin temperatures that form these averages. Due to the range of thermal pressures allowing a two-phase medium, the CNM temperature predicted can theoretically range from ~20 to 350 K in Solar neighbourhood conditions (Wolfire et al. 2003). While the upper limit of stable spin temperatures for the CNM is a function of dust extinction (McClure-Griffiths, Stanimirović & Rybarczyk

2023), we can only safely conclude that this secondary component at approximately 250 K is at the rough boundary of the cold and unstable regime of the neutral ISM. We conclude that the secondary component likely contains a mixture of warmer CNM and cooler UNM.

Unlike the detection sample stack, the non-detections when averaged did not reveal any gas with temperatures above ~1000 K. This is likely because the UNM component's optical depth found in the detection sample residual stack resides close to the ~3σ detection threshold with the non-detection stack's noise level. Since the detection and non-detection sightlines are close in location and span the same small patch-of-sky, it is likely such a component could be found with the non-detections with longer exposure time and hence a reduced noise profile.

The WNM has been observed in absorption in the past (Carilli et al. 1998; Murray et al. 2014; Patra & Roy 2024); however, our large sample of absorption signals when stacked only yield components with spin temperatures belonging to the CNM in the non-detection stack and UNM in the detection stack. This is in contrast to recent work observing WNM temperatures up to $10^4$ K in H1 absorption at optical depth sensitivities close to that of the absorption stack in this work at $\sigma_{\tau} \approx 1-2 \times 10^{-3}$ (Patra & Roy 2024). This could point to a lack of WNM in this region relative to other regions on the sky, especially with the field dominated by 2 large cold dense filamentary structures (Nguyen et al. 2024; Lynn et al., in preparation). Another possibility is that the WNM in this region resides at higher optical depths than normal. Adding to this is the optical depth sensitivity of our absorption profile stack being an order of magnitude lower than that of more sensitive observations like the 21-SPONGE stacking analysis of Murray et al. (2018) ($\sigma_{\tau} \approx 10^{-3}$ versus the 21-SPONGE stack sensitivity of $\sigma_{\tau} \approx 2 \times 10^{-4}$). Despite not detecting the WNM in absorption, we can define an upper limit on the peak optical depth of the WNM in this region of the sky of approximately $\tau_{peak, WNM} \lesssim 3.3 \times 10^{-3}$ using the ~3σ detection threshold of the fitting algorithm when applied to the residual detection spectra stack in Section 4.

## 6.2 Non-detections spatial variations

To determine whether the temperature and optical depth trends observed in Section 5.3 reflect changes in the true density and relative abundance of the H1 phases, we looked to dust. Since at high Galactic latitudes the total amount of hydrogen is in the neutral phase, dust is a good tracer of H1 at high latitudes (Low et al. 1984; Boulanger et al. 1996; Lenz, Hensley & Doré 2017). Converting dust extinction densities from three-dimensional (3D) dust maps to total hydrogen densities (dominated by H1) therefore allows us to probe any variations to 3D H1 density between each of the 6 bins defined in Figs 10 and 11.

Recently, work by Edenhofer et al. (2024) has produced a local 3D dust extinction map ranging from 69 pc to 1.2 kpc from the Sun with an angular resolution of 14 arcmin and a varying distance resolution, from 0.4 to 7 pc. We used the coordinates of the sightlines used in the non-detection absorption stacks and collected, using the dustmaps<fnref n="1" /> package, all of their respective sightlines in the 3D dust map for each of the six spatial bins used in Section 5.3.

We followed the procedure described in O'Neill et al. (2024) to convert the differential extinction measurements $A'_{ZGR_{23}}$, defined in term of unitless extinction $ZGR_{23}$ (Zhang, Green & Rix 2023), into

<footnote n="1">¹https://dustmaps.readthedocs.io/</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Weighted-mean of $n_H$ along the sightlines of each of the 6 bins described in Section 5.3. The dust densities were converted from the Edenhofer et al. (2024) 3D dust map to total hydrogen densities using a conversion factor of 1653 cm⁻³ (O'Neill et al. 2024). Note as the spatial bins sample sightlines further away from the known regions of cold gas (from bin 6 to 1), we observe a decrease in the peak density of total hydrogen and by extension a decrease in the amount of colder phases of H I as we are observing at high Galactic latitudes (Lenz et al. 2017). The vertical dashed line represents the average distance of maximum $n_H$ in the 2240 non-detection sightlines (~240 ± 45 pc), with the grey region indicating the 1σ standard deviation.</figcaption></figure>

densities of total hydrogen atoms H (molecular gas included) given by

$$n_H = C \times A'_{ZGR_{23}}, \qquad (1)$$

where $C = 1653 \text{ cm}^{-3}$. The binned $n_H$ sightlines were then averaged, weighted by the same $1/\sigma_\tau$ values as used in the stacks of Section 5.3. These weighted-mean $n_H$ spectra are shown in Fig. 15.

We find that the collections of sightlines further away from regions of known cold gas and with a lower detection source density (i.e. bin 1) contain a noticeably lower peak with lower $n_H$ densities compared to the regions close to the filament and known cold gas (i.e. bin 6) with the density spectra's peaks increasing in amplitude the higher the detection source density. This trend supports our findings in Section 5.3 where the cold gas optical depth disappears drastically as we move away from the filaments (shown in the left plot of Fig. 13). It is interesting to note that the $n_H$ spectra contain two local maxima, implying the existence of multiple physically separated layers of H I, especially in the bins closer to the location of cold gas. This result aligns well with the fact that absorption sightlines in the region contain multiple cold clouds along individual sightlines and that the region contains two distinct filament structures at different physical distances (Nguyen et al. 2024).

The fact that the optical depth of the primary cold component found in the non-detection stack increases as we move closer to the regions with higher densities of cold gas detections (matching the value of the detection stack used in Section 3.2), also highlights the unbiased nature of the GASKAP absorption survey in determining regions with higher amounts of cold gas.

The primary component's spin temperature on the other hand remains relatively stable across the region, with the value in each bin staying close to the value of the total stack as shown in Fig. 14. Since the primary component traces the peak in emission, which typically accounts for the bulk of cold gas, this result seems to suggest there is little structural variation of the CNM's temperature across this region of the sky. Assuming that the H I absorption lies

at the same physical distance as the peak in the weighted mean dust extinction spectra, the average physical distance to each H I absorption location is approximately 240 ± 48 pc (as shown in Fig. 15). Using this physical distance along with the angular extent of the entire GASKAP Pilot survey of ~25°, the temperature of cold gas within this region of the Milky Way does not vary over a plane-of-sky distance of ~100 ± 20 pc. We note that due to the width of the Local Bubble's shell in the LOS direction, this length scale likely reflects only the plane-of-sky direction and not the depth along the LOS.

The secondary component's lower limit of optical depth exhibits negligible spatial variation from each bin compared to the bulk of the cold gas in the primary components. The optical depth plateaus in the middle bins at a value close to that in the detection stack in Section 3.2 as shown in Fig. 13. The spin temperature of the secondary component also shows a similar trend as shown in Fig. 14. However, due to the lack of sample size, detecting a secondary component in only 2 of the 6 non-detection bins, we cannot make any conclusions on the spatial variation of the secondary component's properties with any confidence.

## 7 SUMMARY

We present analysis using the recent Milky Way GASKAP Phase II Pilot survey's H I absorption detection sample investigating the effect of multiple components along individual sightlines on resulting stacked spectra and reported measurements such as the spin temperature and optical depth. Taking these effects into consideration, we report the results of stacking the Milky Way H I absorption and emission spectra from both the detection and non-detections of the GASKAP survey towards the Magellanic system. The main conclusions of this work are as follows:

(i) When stacking multicomponent H I absorption and emission spectra, any secondary components not centred before stacking create a signal with a peak optical depth lower than the true average value. The FWHM of the component is conversely an upper limit to conserve the total flux of the individual stacked components.

(ii) Due to the velocity offset being equal in emission and absorption, the spin temperature from a stack of secondary components is still well correlated with average of the individual components.

(iii) Stacking the residual spectra of the 462 GASKAP absorption detection sightlines to improve the total signal-to-noise, we detected gas in H I absorption with a spin temperature of ~1320 ± 263 K belonging to the UNM.

(iv) From stacking the 2240 non-detection GASKAP spectra to find gas hidden within the noise of individual sightlines, we detected the existence of a primary narrow and secondary broad component in H I absorption.

(v) Both the primary and secondary component can be attributed to averaged CNM fields in this region of the sky, with a spin temperature of approximately 98 ± 12 and 255 ± 106 K, respectively.

(vi) When spatially binning these non-detections by their plane-of-sky distance from known regions of cold gas, the primary CNM component varies noticeably in optical depth by almost an order of magnitude, decreasing with distance from the locations of known cold gas.

(vii) The change in the primary cold component's optical depth is also reflected in the change in the mean peak extinction density traced by dust over the same lines-of-sight.

(viii) The spin temperatures of cold gas components in each bin remain constant and within $1\sigma$ of the spin temperature values
from the stack of all 2240 non-detections regardless of their spatial location, over a physical plane-of-sky distance of ~ 100 pc.

In the near-future, the GASKAP survey will observe the same fields discussed in this work with longer exposure times, yielding a higher signal-to-noise ratio and sensitivity. Assuming the same number of sources are found as a conservative estimate, the signal-to-noise of all the absorption stacks throughout this work should improve by almost an order of magnitude. From this, we will then approach a signal-to-noise rivalling that of the stacked 21-SPONGE spectra (Murray et al. 2014, 2018) and will be able to detect higher spin temperature components, possibly even belonging to the WNM.

## ACKNOWLEDGEMENTS

This scientific work uses data obtained from Inyarrimanha Ilgari Bundara/the Murchison Radio-astronomy Observatory. We acknowledge the Wajarri Yamaji People as the Traditional Owners and native title holders of the Observatory site. CSIRO's ASKAP radio telescope is part of the Australia Telescope National Facility (https://ror.org/05qajvd42). Operation of ASKAP is funded by the Australian Government with support from the National Collaborative Research Infrastructure Strategy. ASKAP uses the resources of the Pawsey Supercomputing Research Centre. Establishment of ASKAP, Inyarrimanha Ilgari Bundara, the CSIRO Murchison Radio-astronomy Observatory, and the Pawsey Supercomputing Research Centre are initiatives of the Australian Government, with support from the Government of Western Australia and the Science and Industry Endowment Fund. This research is supported by an Australian Government Research Training Program (RTP) Scholarship. This research was partially funded by the Australian Government through the Australian Research Council through the award of an Australian Laureate Fellowship (project number FL210100039) to NMc-G. The authors acknowledge Interstellar Institute's programme 'I16' and the Paris-Saclay University's Institut Pascal for hosting discussions that nourished the development of the ideas behind this work. EDT was supported by the European Research Council (ERC) under grant agreement no. 10104075.

This research made use of Astropy, a community-developed core Python package for astronomy (Astropy Collaboration et al. 2013), NumPy (van der Walt, Colbert & Varoquaux 2011), matplotlib, a Python library for the publication of quality graphics (Hunter 2007), and SciPy (Virtanen et al. 2020).

## DATA AVAILABILITY

The GASKAP-HI Pilot survey data are available on the CSIRO ASKAP Science Data Archive<fnref n="2" /> (CASDA). The non-detection data set and auxiliary data products from this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Astropy Collaboration et al., 2013, A&A, 558, A33

* Audit E., Hennebelle P., 2005, A&A, 433, 1

* Begum A., Stanimirović S., Goss W. M., Heiles C., Pavkovich A. S., Hennebelle P., 2010, ApJ, 725, 1779

* Boulanger F., Abergel A., Bernard J. P., Burton W. B., Desert F. X., Hartmann D., Lagache G., Puget J. L., 1996, A&A, 312, 256

* Carilli C. L., Dwarakanath K. S., Goss W. M., 1998, ApJL, 502, L79

* Crovisier J., Kazes I., Aubry D., 1978, A&AS, 32, 205

* Dempsey J. et al., 2022, PASA, 39, e034

* Dickey J. M. et al., 2013, PASA, 30, e003

* Dickey J. M., Salpeter E. E., Terzian Y., 1977, ApJ, 211, L77

* Dwarakanath K. S., Carilli C. L., Goss W. M., 2002, ApJ, 567, 940

* Edenhofer G., Zucker C., Frank P., Saydjari A. K., Speagle J. S., Finkbeiner D., Enßlin T. A., 2024, A&A, 685, A82

* Field G. B., 1965, ApJ, 142, 531

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

* Haud U., 2000, A&A, 364, 83

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

* Hennebelle P., Pérault M., 1999, A&A, 351, 309

* Hennebelle P., Pérault M., 2000, A&A, 359, 1124

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

* Kalberla P. M. W. et al., 2010, A&A, 521, A17

* Kalberla P. M. W., Haud U., 2015, A&A, 578, A78

* Kalberla P. M. W., Haud U., 2018, A&A, 619, A58

* Lazareff B., 1975, A&A, 42, 25

* Lenz D., Hensley B. S., Doré O., 2017, ApJ, 846, 38

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

* Low F. J. et al., 1984, ApJ, 278, L19

* Marchal A., Miville-Deschênes M.-A., 2021, ApJ, 908, 186

* Marchal A., Miville-Deschênes M.-A., Orieux F., Gac N., Soussen C., Lesot M.-J., d'Allonnes A. R., Salomé Q., 2019, A&A, 626, A101

* McClure-Griffiths N. M. et al., 2009, ApJS, 181, 398

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

* Mohan R., Dwarakanath K. S., Srinivasan G., 2004, JA&A, 25, 143

* Murray C. E. et al., 2014, ApJ, 781, L41

* Murray C. E., Stanimirović S., Goss W. M., Heiles C., Dickey J. M., Babler B., Kim C.-G., 2018, ApJS, 238, 14

* Nguyen H. et al., 2024, MNRAS, 534, 3478

* Nidever D. L., Majewski S. R., Butler Burton W., 2008, ApJ, 679, 432

* O'Neill T. J., Zucker C., Goodman A. A., Edenhofer G., 2024, ApJ, 973, 136

* Patra N. N., Kanekar N., Chengalur J. N., Roy N., 2018, MNRAS, 479, L7

* Patra N. N., Roy N., 2024, MNRAS, 529, 4037

* Riener M., Kainulainen J., Henshaw J. D., Orkisz J. H., Murray C. E., Beuther H., 2019, A&A, 628, A78

* Roy N., Kanekar N., Braun R., Chengalur J. N., 2013, MNRAS, 436, 2352

* Treister E., Schawinski K., Volonteri M., Natarajan P., Gawiser E., 2011, Nature, 474, 356

* van der Walt S., Colbert. S. C., Varoquaux G., 2011, Comput. Sci. Eng., 13, 22

* Virtanen P. et al., 2020, Nat. Methods, 17, 261

* Wall J. V., Jenkins C. R., 2003, Cambridge Observing Handbooks for Research Astronomers, vol. 3, Practical Statistics for Astronomers. Cambridge Univ. Press, Cambridge, UK

* 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

* Zhang X., Green G. M., Rix H.-W., 2023, MNRAS, 524, 1855

## APPENDIX: TOY MODEL THEORETICAL CONSIDERATIONS

In this section, we use toy models to investigate the biases introduced by mixing complex multiphase emission and absorption spectra through stacking.

### A1 Single component

We first begin our toy model by assuming the simple case of all sightlines containing a single component. Here, every component is perfectly shifted to 0 km s⁻¹ by the peak of the spectrum in emission. However, as mentioned in Section 3.1, this shifting may be inaccurate

<footnote n="2">²https://research.csiro.au/casda/</footnote>
to ~ 1–2 km s⁻¹ from the misalignment of absorption and emission peaks due to differing noise levels and sightlines used.

To test how stacking these components affects measurements such as the spin temperature, we set up a toy model with 2000 sightlines, each with an emission and absorption spectrum. In each sightline, a single component was created with a width and amplitudes in brightness temperature and optical depth typical of cold gas in the CNM. The properties for the 2000 sightlines were sampled from Gaussian distributions with the mean value uniformly chosen within 5 ≤ T<sub>B,peak</sub> ≤ 15 K, 0.02 ≤ τ<sub>peak</sub> ≤ 0.3, FWHM ~ 5 km s⁻¹ and standard deviations of σ<sub>T<sub>B</sub></sub> = 4 K, σ<sub>τ<sub>peak</sub></sub> = 0.01 and σ<sub>FWHM</sub> = 0.2. The central velocities of these components is centred around 0 km s⁻¹, with a slight deviation to account for the 1–2 km s⁻¹ deviation between emission and absorption components found in reality (a central velocity dispersion of ~ 0.5 in the toy model). Stacking these spectra, we then calculate the spin temperature, T<sub>s</sub> of the stacked component via:

$$T_s = \frac{T_{B, \text{peak}}}{1 - e^{-\tau_{\text{peak}}}} \qquad (A1)$$

Alongside this, we also compute the spin temperature of the individual sightlines and take the average to create a spin temperature that should be measured. Repeating this scenario 200 times (using a different sampled mean value of T<sub>B,peak</sub> and τ<sub>peak</sub> for each iteration), we find the spin temperature of the stack is well correlated with the average spin temperature of the individual sightlines (figure not shown).

## A2 Multiple components

As discussed in Section 3.1, the most sensitive detection sources contain multiple components: a narrow primary component aligned with the peak in the corresponding emission profile, and a broader secondary component offset in velocity.

The primary component tied to the peak in emission will be aligned well at 0 km s⁻¹ after shifting, whereas the additional peaks will merge into a broader secondary component.

To explore these effects of stacking multicomponent spectra, we create a similar toy model to that of Section A1. We create 2000 mock sightlines, each with an absorption and emission spectrum. We create primary components similar to the previous section. In addition, we add a broader component. This secondary component has a width and brightness temperature and optical depth amplitude randomly sampled within typical values of warmer CNM and colder UNM gas. The properties for the 2000 sightlines were sampled from Gaussian distributions with the mean value uniformly chosen within 3 ≤ T<sub>B,peak</sub> ≤ 15 K, 0.02 ≤ τ<sub>peak</sub> ≤ 0.03, FWHM ~ 15 km s⁻¹ and standard deviations of σ<sub>T<sub>B</sub></sub> = 4 K, σ<sub>τ<sub>peak</sub></sub> = 0.01 and σ<sub>FWHM</sub> = 0.6, respectively. The central velocity of this component is randomly sampled over a normal distribution centred on 0 km s⁻¹ but with a dispersion of ~ 15 km s⁻¹, simulating the potential of a non-zero central velocity of these secondary components after shifting.

Fig. A1 demonstrates the effects of stacking these multicomponent absorption sightlines featuring a primary component, shifted correctly to 0 km s⁻¹, and a broader secondary component whose central velocity is chosen from a normal distribution with a deviation of 15 km s⁻¹. The values for the peaks and widths of the two components were chosen based on the parameter space mentioned earlier. The primary components, being correctly centered, produce an averaged component after Gaussian decomposition equivalent to an average of the individual primary components. However, due to the velocity variation of the secondary components, the average

<figure><img src="image_17.png" /><figcaption>**Figure A1.** Comparison of the primary and secondary component from Gaussian decomposition of a toy model stacked absorption spectra with the separated stacks of the individual primary and secondary components. The total stack (black) is produced as described in Section A2. We note the primary component in the stack (cyan-dashed) is identical to the separately averaged primary component (solid blue). The secondary component (red-dashed) is shallower and has a greater FWHM than the separately averaged secondary component (solid red) where each component was centered to 0 km s⁻¹.</figcaption></figure>

secondary component found after Gaussian decomposition varies drastically in amplitude and width from the average of the individual sightlines of the secondary components if they were all centred at 0 km s⁻¹.

Repeating this experiment 100 times (with a different sampled mean value of T<sub>B,peak</sub> and τ<sub>peak</sub> for each iteration within the ranges mentioned above), we compared the measured amplitudes of the components with that of the average amplitude of the components in emission and absorption in individual sightlines. This comparison is shown in Fig. A2. Only the secondary component’s properties are displayed since the primary component’s values agree with the average of individual primary components due to the correct centring of the features prior to stacking. For the secondary components, however, the measured amplitudes are always lower. This is because the total flux is conserved in the average regardless of the Gaussian’s central velocities. Increasing the velocity spread of the Gaussians broadens the final stack, and to compensate, the averaged amplitude must diminish. Therefore, regardless of the severity of the range of central velocities of secondary components in individual sightlines, the peak of the secondary component will always be lower and its width always greater than the average of the individual sightline values.

Despite this lower limit, the stacked secondary components’ amplitudes are still linearly correlated with the average values. Since the central velocities of the secondary components are approximately equal in emission and absorption, the overall effect of the velocity dispersion on the stack is equivalent in emission and absorption. As a result, the gradient of the linear correlation of measured peak brightness temperature and optical depth with the individual sightline average is the same.

When calculating the spin temperature of the secondary component through this toy model, we find that these gradients cancel out through the division of emission and absorption amplitudes. The spin temperature of the secondary component is therefore well correlated with the average of the individual secondary components’ spin temperatures, as shown in Fig. A2.
<figure><img src="image_18.png" /></figure>

**Figure A2.** Comparison of observed peak optical depth (red) and average spin temperature (blue) of the secondary components from 100 runs of stacking 2000 two-component spectra, as described in Section A2. The peak brightness temperature is not displayed but shows a similar trend to the peak optical depth due to the velocity offsets being the same in both emission and absorption. The observed peak optical depth is always less than the average peak optical depth of the individual sightlines' secondary components. However, the observed spin temperature of the secondary component correlates well with the average spin temperature of each sightline's secondary component.

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