# GALFA-SPS

Amit Kumar Mittal  Department of Astronomy, University of Wisconsin-Madison,
475 N Charter Street, Madison, WI, 53706-1582 USA  Department Electrical and
Computer Engineering, University of Wisconsin-Madison, 1415 Engineering Dr,
Madison, WI, 53706-1582 USA  Brian L Babler  Department of Astronomy,
University of Wisconsin-Madison, 475 N Charter Street, Madison, WI, 53706-1582
USA  Snežana Stanimirović Department of Astronomy, University of Wisconsin-
Madison, 475 N Charter Street, Madison, WI, 53706-1582 USA  Nickolas Pingel
Department of Astronomy, University of Wisconsin-Madison, 475 N Charter
Street, Madison, WI, 53706-1582 USA

(July 2023)

# Neutral Hydrogen (HI ) 21 cm as a probe: Investigating Spatial Variations in
Interstellar Turbulent Properties

Amit Kumar Mittal  Department of Astronomy, University of Wisconsin-Madison,
475 N Charter Street, Madison, WI, 53706-1582 USA  Department Electrical and
Computer Engineering, University of Wisconsin-Madison, 1415 Engineering Dr,
Madison, WI, 53706-1582 USA  Brian L Babler  Department of Astronomy,
University of Wisconsin-Madison, 475 N Charter Street, Madison, WI, 53706-1582
USA  Snežana Stanimirović Department of Astronomy, University of Wisconsin-
Madison, 475 N Charter Street, Madison, WI, 53706-1582 USA  Nickolas Pingel
Department of Astronomy, University of Wisconsin-Madison, 475 N Charter
Street, Madison, WI, 53706-1582 USA

(July 2023)

###### Abstract

Interstellar turbulence shapes the HI  distribution in the Milky Way (MW). How
this affects large-scale statistical properties of HI  column density across
the MW remains largely unconstrained. We use the $\sim 13,000$ square-degree
GALFA-HI  survey to map statistical fluctuations of HI  over the $\pm 40$ km
s-1 velocity range. We calculate the spatial power spectrum (SPS) of HI
column density image by running a 3-degree kernel and measuring SPS slope over
a range of angular scales from 16′to 20$\degree$. Due to GALFA’s complex
observing and calibration strategy, we construct detailed estimates of the
noise contribution and account for GALFA beam effects on SPS. This allows us
to systematically analyze HI  images that trace a wide range of interstellar
environments.

We find that SPS slope varies between $\sim-2.6$ at high Galactic latitudes,
and $\sim-3.2$ close to Galactic plane. The range of SPS slope values becomes
tighter when we consider HI  optical depth and line-of-sight length caused by
the plane-parallel geometry of HI  disk. This relatively uniform, large-scale
distribution of SPS slope is suggestive of large-scale turbulent driving being
a dominant mechanism for shaping HI  structures in the MW and/or the stellar
feedback turbulence being efficiently dissipated within dense molecular
clouds. Only at latitudes above $60$ degrees we find evidence for HI  SPS
slope being consistently more shallow. Those directions are largely within the
Local Bubble, suggesting the recent history of this cavity, shaped by multiple
supernovae explosions, has modified the turbulent state of HI  and/or
fractions of HI  phases.

Interstellar atomic gas (833); Interstellar phases (850); Interstellar medium
(847); Milky Way Galaxy (1054); Neutral hydrogen clouds (1099)

††software:  Astropy (Astropy Collaboration et al., 2013), DS9 (Joye & Mandel,
2003), HealPy (Górski et al., 2005), KARMA (Gooch, 1996), Matplotlib (Hunter,
2007), NumPy (van der Walt et al., 2011), SciPy (Oliphant, 2007)

tablenum SIXtablenum

##  1 Introduction

Interstellar turbulence is an integral driver in structuring the interstellar
medium (ISM) (e.g., Elmegreen & Scalo, 2004; McKee & Ostriker, 2007; Lazarian,
2009). It is thought to be ubiquitous throughout the ISM and important for
many astrophysical processes, such as molecule formation (Bialy et al., 2017),
accretion within disks around planets (Klessen & Hennebelle, 2010), and even
large-scale galactic dynamics (e.g., Klessen & Hennebelle, 2010; Krumholz &
Burkhart, 2016). Many astronomical sources and processes can drive
interstellar turbulence, which then cascades and affects distributions of
interstellar gas and dust.

Over the past two decades, numerous observational (e.g., Crovisier & Dickey,
1983; Stanimirovic et al., 1999) and numerical studies (e.g., Grisdale et al.,
2017, 2018) have investigated turbulent properties in the ISM. Despite these
efforts, many questions related to the physical processes driving interstellar
turbulence remain unanswered. These include the scales on which these
processes operate and the observational signatures they imprint on the
structure of interstellar clouds over their lifetime. Additionally, the
characterization of turbulent properties remains a challenge, with open
questions regarding their manifestation across different ISM phases, such as
the cold, unstable, and warm neutral medium (CNM, UNM, and WNM). Similarly,
the influence of optical depth ($\tau$) and line-of-sight depth also remain
largely unconstrained.

Many different sources of turbulence have been suggested over the years. They
all influence turbulence properties (the slope of fluctuations, driving-scale
and modes of dominant motions) in different ways. For example, stellar
feedback sources (proto-stellar jets, stellar outflows and winds, photo-
heating, supernova explosion) operate largely on sub-cloud spatial scales (Kim
et al., 2001; de Avillez & Breitschwerdt, 2005; Joung & Mac Low, 2006;
Tamburro et al., 2009; Ostriker & Shetty, 2011; Shetty et al., 2012; Ferland
et al., 2013; Grisdale et al., 2017), while gravitational and magneto-
rotational instabilities (Wada et al., 2002; Block et al., 2010; Krumholz &
Burkhart, 2016) operate largely on cloud-size spatial scales. In addition, the
supernova-driven turbulence is more effective in producing compressive
motions, while in regions with low star formation activity solenoidal motions
are expected to be more prominent (Federrath et al., 2010).

Although the various sources of turbulence affect the density and velocity
fluctuations in three dimensions, ISM turbulence is traditionally studied
using statistical descriptors of intensity fluctuations. For instance, the
spatial power spectrum, SPS, (Crovisier & Dickey, 1983; Pingel et al., 2013,
2018), the $\delta$-variance (Stutzki et al., 1998), and different orders of
structure functions (Padoan et al., 2012; Burkhart et al., 2015) are
frequently used. The SPS, in particular, is widely applied in both Galactic
and extra-galactic studies (Dame et al., 2001; Elmegreen et al., 2001; Block
et al., 2010; Combes et al., 2012; Zamora-Avilés et al., 2012; Pingel et al.,
2013; Szotkowski et al., 2019). The properties of intensity fluctuations can
be used to derive information about underlying velocity and/or density
fluctuations, as well as to uncover the spatial scales at which energy is
injected or dissipated and the mechanisms by which this energy cascades from
large to small scales, or vice versa.

There have been many studies of the HI  turbulent properties in the MW using
the SPS. For example, the SPS method was applied on HI  images for various
regions of our Galaxy (Crovisier & Dickey, 1983; Green, 1993; Dickey et al.,
2001; Pingel et al., 2013, 2018). Similarly, HI observations of several nearby
galaxies, such as the Small Magellanic Cloud (SMC; Stanimirovic et al., 1999;
Stanimirović & Lazarian, 2001; Szotkowski et al., 2019), the Large Magellanic
Cloud (LMC; Elmegreen et al., 2001; Szotkowski et al., 2019), galaxies from
the LITTLE THINGS survey (Hunter et al., 2012; Zhang et al., 2012) have been
investigated with the SPS. The results of these studies have shown the
existence of a hierarchy of structures in the diffuse ISM, with the SPS
spectral slope varying slightly with interstellar environments. For example,
HI  integrated intensity images of several regions in the MW, which are
largely confined to the disk and likely affected by high optical depth, showed
a power-law slope of the SPS of $\sim-3$ (Dickey et al., 2001). A slightly
steeper slope of $-3.4$ was found for the HI  column density in the SMC, and
$-3.7$ for the LMC. Galaxies like the LMC and M33 showed a break in their HI
power spectra, while no departure from a single power law was observed in the
MW, the SMC (Stanimirovic et al., 1999; Szotkowski et al., 2019; Koch et al.,
2020), and the Magellanic Bridge when the integrated HI  intensity
distributions were examined. Muller et al. (2004) did notice a change in the
SPS slope across several different regions between $-2.9$ and $-3.5$.

These observational studies offer valuable insights into the distribution of
power at different spatial scales, which can be compared with theoretical and
numerical predictions to determine the characteristics of the turbulent
cascade (Block et al., 2010; Pilkington et al., 2011; Grisdale et al., 2017).
The SPS can also identify the significance of shocks (Beresnyak et al., 2005),
and the scales at which energy is injected and dissipated (Kowal & Lazarian,
2007), therefore providing crucial constraints for numerical simulations.

In this study we continue our investigation of spatial variations of turbulent
properties as a way of quantifying the importance of stellar feedback,
gravity, and ram pressure as structuring agents of the HI  distribution in the
MW. In particular, to investigate possible spatial variations of the SPS slope
we apply a consistent methodology on the HI  data from the entire Galactic
Arecibo L-Band Feed Array-HI (GALFA-HI ) survey (Peek et al., 2011), which
offer a unique compromise between a high angular resolution and a large
spatial area (resulting in a high spatial dynamic range).

This paper is organized in the following way. In Section 2 we summarize the
data sets used in this study. Section 3 explains our SPS and rolling SPS
approach for mapping spatial variability of turbulent properties across the
GALFA-HI  survey. We discuss sources of noise and present a method to
characterize the noise contribution to the SPS in the Appendix section.
Section 4 presents our SPS slope image and considers how optical depth and
plane-parallel geometry affect HI  turbulent properties. Finally, we discuss
and conclude our results in Section 5.

##  2 Data: GALFA-HI

The primary data utilized in this investigation originates from the second
data release (DR2) of the GALFA-HI  survey (Peek et al., 2018), which was
carried out using the ALFA seven-beam receiver array. The full width at half
maximum (FWHM) size of the elliptical ALFA beam is $3.3\arcmin\times
3.8\arcmin$. The survey area spans approximately 13,000 square degrees,
including the declination range (DEC) of $-1\degree 20\arcmin$ to $38\degree
02\arcmin$, and the entire range of right accession (RA). The survey covers
the velocity range of $-$650 to $+$650 km s-1, with velocity resolution of
0.18 km s-1. Further details about the GALFA-HI  survey observations and data
processing are provided in Peek et al. (2018).

The full velocity coverage includs the local MW HI  as well as intermediate
and high-velocity HI  beyond what is expected from the MW rotation curve. To
focus solely on the local HI  or the low-velocity clouds (LVC), we limit the
velocity range to $|\rm{v}|\leq$ 40 km s-1. We calculate the HI  column
density under the assumption of HI  being optically thin. The HI  column
density, $N(HI)_{\rm thin}$, is given by:

$N_{HI,\rm{thin}}=C_{0}\,\int T_{B}(v)dv.$ (1)

where $C_{0}=1.823\times 10^{\rm{18}}$ cm-2/(K km s-1) (Draine, 2011), $dv$ is
the width of a single spectral channel in km s-1, and $T_{B}$ in K is the
brightness temperature. Equation (1) is used to approximate $N(HI)$ in the
absence optical depth information. In Section 4.4 we investigate the effect of
high optical depth on the SPS slope.

Figure 1 shows the spatial coverage of the GALFA-HI  survey. The survey
includes two regions of the MW plane, inner and outer Galaxy, as well as high-
latitude HI .

Figure 1: The HI  column density calculated by integrating the HI  brightness
temperature from the GALFA-HI  survey over the velocity range $-$40 $\leq\rm
v\leq$ 40 km s-1 . The Galactic coordinate grid is overlaid.

##  3 Methods: Spatial power spectrum (SPS) and the rolling SPS

###  3.1 SPS

To calculate HI  SPS we use a methodology similar to these from earlier
studies (Crovisier & Dickey, 1983; Green, 1993; Stanimirovic et al., 1999;
Elmegreen et al., 2000). For a 2D image $I(x)$, the 2D SPS is defined as:

$P(k)=\iint\langle I(x)I(x^{\prime})\rangle e^{-i{{\bf
L}\cdot{\bf\emph{k}}}}\,dL$ (2)

where $k$ is the spatial frequency, measured in units of wavelength and being
proportional to the inverse of the spatial scale111 Spatial frequency ($k$ in
arcmin-1) varies inversely with spatial scale ($\lambda$ in arcmin) as:
$k=1/\lambda$., while ${\bf L}={\textbf{\emph{x}}}-{\textbf{\emph{x}'}}$ is
the distance between two points.

We first regrid the GALFA-HI  column density image to ensure independent
pixels, with the pixel size matching the angular size of the telescope beam
which is $\sim 4^{\prime}$. We then apply a median filter with box size of 1∘
to identify pixels with column density $\geq$5-$\sigma$, where standard
deviation was calculated for the same 1∘ box. We replace the hot/extreme
pixels with the corresponding values from the median filter image. This is
done to remove $\sim$0.005$\%$ hot pixels which would affect the Fourier
transform due to sharp discontinuities.

To calculate the SPS for an image (column density), we take the Fourier
transform and square the modulus of the transform,
$\langle\Re^{2}+\Im^{2}\rangle$ (where $\Re$ and $\Im$ are the real and
imaginary parts of the 2D Fourier transform), which we refer to as the modulus
image. We proceed by selecting a range of annuli within the modulus image,
with evenly spaced intervals determined by $0.03\log(k)$. This specific bin
size strikes a balance, as both increasing and decreasing this spacing have
drawbacks: an increase would lead to a reduction in the number of data points
in the fitting, while a decrease would result in the appearance of only a few
bins near the center of the Fourier plane. Within each annulus, we assume
azimuthal symmetry and compute the median value. Pingel et al. (2013)
demonstrated that the median provides a more accurate representation of the
average power due to occasional bright pixels in the modulus image caused by
the Gibbs phenomenon. The Gibbs phenomenon, also known as the edge effect,
arises due to sharp image boundaries and the HI  emission rarely reaching zero
at the edges of the surveyed regions. As the Fourier transform of a step
function is a $sinc$ function (Bracewell, 1986), bright pixels appear at the
center of $P(k)$ (an example is shown in Figure 2 of Pingel et al. (2018)).

We plot the median power as a function of $\log$ (spatial scale), and estimate
the uncertainties using the median absolute deviation (MAD) divided by
$\sqrt{N}$, where $N$ is the number of independent pixels in each annulus As
most previous studies of the HI  SPS in the MW have shown featureless SPS
spectra (Pingel et al., 2013; Marchal & Miville-Deschênes, 2021), we start by
fitting a single power-law function to the observed SPS:

$P_{\rm obs}(k)\propto P_{\rm model}=A\times k^{\gamma}$ (3)

where $A$ is the power-law amplitude, and $\gamma$ is the power-law SPS index
for a given SPS. When fitting the SPS we use a weighted fitting procedure
where weights for each binned data points are calculated as inversely
proportional to the square of the uncertainty for that binned data point
(Taylor, 1997).

###  3.2 Modeling of the observed SPS

Equation (3) is an ideal model which does not consider any telescope
systematic effects. However, as pointed out by several recent studies (Pingel
et al., 2022; Marchal & Miville-Deschênes, 2021; Koch et al., 2020; Kalberla &
Haud, 2019; Blagrave et al., 2017; Martin et al., 2015), an observed SPS is
affected by the telescope beam and fluctuations due to instrumental noise
(random and systematic). Therefore, we adopt the methodology proposed by Koch
et al. (2020) to account for instrumental effects:

$P_{\rm{obs}}(k)=P_{\rm{beam}}(k)\times P_{\rm{model}}+P_{\rm{noise}}(k)$ (4)

where we assume a 2D Gaussian beam shape with a FWHM
$\rm{4}\arcmin\times\rm{4}\arcmin$ size, therefore
$P_{\rm{beam}}(k)\propto\exp(-4\pi^{2}\sigma_{\rm{beam}}^{2}k^{2})$, and
$\sigma_{\rm{beam}}=\rm{FWHM}/2\sqrt{\rm{\log 2}}$. In the case of GALFA-HI ,
the noise term $P_{\rm{noise}}(k)$ is not just due to random (white) noise and
therefore can not be represented as a constant term in Equation (4). In fact,
the noise comes from several contributions, including residuals from the
telescope scanning pattern, and we show in the Appendix that the noise term
can be written as $P_{\rm{noise}}(k)=C\times Q(k,sd)$ as it because it depends
on the measured standard deviation of the total noise column density image and
brightness temperature. To estimate the noise power spectrum $Q(k,sd)$, we
first generate a noise image that takes into account both the off-HI  line
noise (modeled as white Gaussian noise) and systematic noise resulting from
scanning artifacts. From the noise image we calculate the noise power spectrum
and then parameterize it as $C\times Q(k,sd)$ to make the noise power spectrum
smoother and improve statistics. The Appendix contains a comprehensive
description and modeling of the noise power spectrum. We also show there
several examples of the noise SPS and its contribution to the HI  SPS for
different Galactic regions.

To fit the SPS in Fourier space, we only consider angular scales where $k\leq
k_{\rm{max}}$, which is defined as four times the FWHM. The purpose of this
limitation is to avoid spatial scales where the complex ALFA beam pattern
(which we model as a 2D Gaussian function) may still affect the SPS.
Conversely, we restrict the fitting process at the larger spatial scale end to
$k>k_{\rm{min}}$, which is defined as the inverse of half the size of the
image. This is done because the largest scale bin is estimated from just a few
samples in the 2D power spectrum and thus has a large uncertainty.

Prior to fitting the SPS, we compare $P_{obs}(k)$ and $P_{noise}(k)$ and
proceed with fitting a power-law only if $P_{obs}(k)>9\times P_{noise}(k)$ for
the entire range of spatial frequencies. This is equivalent to ensuring that
the signal-to-noise ratio (S/N) exceeds 3 in the image domain. As a result,
some pixels in the SPS slope image (Figure 5) are left blank, those are the
pixels that do not satisfy the noise threshold criterion.

###  3.3 Rolling SPS

We use the recently developed method – the rolling SPS – by Szotkowski et al.
(2019) which runs a moving kernel across an HI  column density image to
produce an image of the SPS slope. To ensure the best possible angular
resolution while having enough statistics to calculate the SPS we selected a
3$\degree$ kernel size. In their study, Szotkowski et al. (2019) examined the
optimal kernel size by analyzing simulated images. They found that increasing
the kernel size led to a decrease in the standard deviation of the mean value
of the SPS slope, as demonstrated in Figure A of their appendix. Additionally,
they observed that increasing the sub-image size resulted in a reduction of
usable pixels at the edge of the image, which conveyed less information about
small-scale structure variations.

In this study, we probed angular scales in the range of $1.5\degree$ to
$16\arcmin$ per single 3-degree kernel. To obtain the SPS slope $\gamma$ for
each kernel, we follow the methodology outlined in the previous section,
accounting for the effects of both the beam and noise. We then shift the
kernel by $8\arcmin$ and repeat the process until the SPS slopes are computed
across the entire GALFA-HI  field. During the SPS fitting process we also
check the goodness of the fit by evaluating the $\chi^{2}$ for each fit.

##  4 Results

###  4.1 GALFA-HI  for the inner, mid and outer Galactic regions

(a)

(b)

Figure 2: Top: The GALFA-HI  LVC column density image, divided into ten equal
sub-regions named as G0 to G9 where each region is $36\times 36$ degrees in
size. Bottom: The corresponding HI  SPS calculated for the above G0 to G9
regions. The grey area shows the range of angular scales used for fitting the
SPS (between 18$\degree$ and 16$\arcmin$, which corresponds to spatial
frequencies $k=10^{-3}$ to $\sim 6\times 10^{-2}$ arcmin-1).

To investigate spatial variations of HI  turbulent properties we start by
dividing the entire HI  column density image into ten equal regions, each
about $36\times 36$ degree in size. These regions are shown in the top panel
of Figure 2 and sample the inner MW plane HI  emission, the high latitude, and
the outer MW plane HI  emission, respectively. The bottom panels of Figure 2
show the SPS of these regions. In fitting the SPS, we used the approach
explained in Section 3. The SPS shown in Figure 2 was fitted over the angular
range of scales 18$\degree$ to $16\arcmin$.

Figure 2 shows that for all regions the HI  SPS is well represented with a
single power-law function without obvious break points. It also suggests that
the inner and outer disk regions (viewed mainly in regions G1 and G7) have a
slightly steeper SPS slope $\sim-3.2\pm 0.1$ relative to the high-latitude
regions (G4 and G5) which have a more shallow slope, $\sim-2.9\pm 0.1$. In
particular, the two most extreme slopes are found for regions G1 and G6, with
corresponding SPS slopes being $-3.37\pm 0.05$ and $-2.80\pm 0.06$
respectively. This result suggests that some spatial variations of the SPS
slope may exist across the MW.

To demonstrate the difference in the SPS slope between the disk and high-
latitude HI  and to extend the range of angular scales even further, in Figure
3 we select three regions, each about $40\times 40$ degrees large,
specifically centered on the inner Galaxy, the outer Galaxy, and at high
latitudes. The inner Galaxy region covers the longitude range 30 to 77 degrees
and for our selected velocity range includes mostly the Perseus arm (and a
small fraction of HI  in each of the near Sagittarius-Carina and the Outer
Scutum Centarus arms), therefore including the HI  mostly within few to 10 kpc
from the Sun. The outer Galaxy region covers the longitude range of 170 to 215
degrees, and includes four major spiral arms (Local, Perseus, Scutum-Centaurus
and the Norma–Outer arm which is at $\sim 10$ kpc from the Sun). While this
region includes the HI  within up to 5-10 kpc from the Sun, it also includes
some more distant HI  that reaches Galactocentric distance of up to 20 kpc, as
opposed to the inner Galaxy region which stays within the Galactocentric
radius of 10 kpc.

The corresponding SPS for the three regions is shown in Figure 4. This figure
confirms that the inner and outer Galactic plane regions (shown as red and
green lines) have similar SPS slopes that are consistently slightly steeper
than the slope of the high-latitude region (shown in yellow). If we divide the
fitting range into two, the three regions have consistent SPS slope measured
for the small-scale end of the fitting range. However, on the large-scale end
of the fitting range, the high-latitude region has a much shallower SPS slope
relative to the low-latitude regions (e.g. $-2.35\pm 0.17$ vs $-3.73\pm
0.17$). This shows that the main difference between low- and high-latitude HI
SPS occurs on scales larger than $\sim 5$ degrees. This is due to lack of
large-scale diffuse HI  emission at high latitudes. The same conclusion can be
reached by comparing top two panels of Figure 13in the Appendix, where the
noise SPS between low- and high-latitude regions is comparable, but the HI
power on large scales is a factor of $\sim 100$ higher for the low-latitude
field. It is also interesting to find that the inner and outer Galaxy regions
have consistent SPS slopes. While both regions probe HI  within a few kpc from
the Sun, the outer Galaxy region is at a higher Galactocentric radius.

(a)  Figure 3: The HI  column density images for three regions: inner Galactic
plane (left), high Galactic latitude (mid), and outer Galactic plane (right).
Each region has a size $40\times 40$ degrees. The three regions are centered
at: RA $\sim$ 19h30m (left), RA $\sim$ 12h45m (mid), and RA $\sim$ 06h00m
(right) with declination centered at $18.5$ degree for each. (a)  Figure 4:
The SPS for three HI column density images shown in Figure 3. The fitting
range, shown by the grey area, covers a range of angular scales from $\sim
20$° to 16$\arcmin$, which corresponds to spatial frequencies $k=10^{-3}$ to
$\sim 6\times 10^{-2}$ arcmin-1.

###  4.2 Rolling SPS for the entire GALFA-HI  field

The rolling SPS analysis, as explained in Section 3, is employed using a
3-degree kernel over the entire low-velocity HI  column density image. The SPS
slope is calculated for each kernel, but only for regions where the HI
brightness temperature is at least three times higher than the noise level to
avoid areas dominated by noise. As we explained in Section 3.2, the S/N check
is done in the Fourier plane, and we do not perform fitting if the S/N
threshold is not met.

Figure 5 (top) displays the resulting rolling SPS slope image, with noise-
dominated pixels being left blank (that way these pixels are not included in
further analyses). The uncertainty in the estimated slope is shown in the
middle panel of the figure, with a typical uncertainty of 0.1-0.2. No residual
structure is visible in this image. Finally the bottom panel of the figure
shows the $\chi^{2}$ image which looks relatively uniform and featureless and
we conclude that a single power-law fit provides a good fit to the SPS and
discard any need of broken power-law fitting.

(a)  Figure 5: Top: The GALFA-HI rolling SPS slope, where the SPS slope is
calculated for each $3\degree$ kernel and rolled across RA and DEC by stepping
every $8\arcmin$ to cover the entire GALFA-HI region. The SPS slope was
determined only for regions where the SPS power $\geq 9\times$ noise power,
and the fitting was conducted over the range of 1.5$\degree$ to 16$\arcmin$.
Pixels that do not satisfy this criterion are left blank (NaN). The Galactic
coordinate grid is overlaid as black dashed lines. Middle: The uncertainty in
the fitted SPS slope. Bottom: The $\chi^{2}$ map of the single power-law SPS
fit.

On a 3-degree spatial scale, the SPS slope varies between $-2.6$ and $-3.2$
(1-$\sigma$ range). This figure confirms that the Galactic plane regions have
systematically the steepest SPS slope (the image is dominated by red/orange
color), while the high Galactic latitude regions have the most shallow SPS
slope (the image is dominated by blue color). In Figure 6 we plot histograms
of pixels shown in the top of Figure 5 for the most extreme SPS slope values:
in orange we show data for $|b|\leq 10$ degrees, while in blue we show data
for $|b|\geq 70$ degrees. The mean SPS slope for these two sub-samples is
$\sim$ $-3.0$ and $\sim$ $-2.7$, respectively. As we have already seen, there
appears to be a slight difference in SPS slope values between the MW plane and
high latitude regions. However, Figure 5 also shows that pixels with low
signal/noise (S/N) are mainly found at high Galactic latitude. We excluded
extreme, low-S/N pixels while fitting the SPS. As we discussed in the
Appendix, the GALFA-HI  noise is complex and affects a range of angular
scales. If not accounted properly (see Figure 13), it would result in a more
shallow SPS slope. While some of the trend in the SPS slope being more shallow
at high latitudes could be due to the systematic decrease of S/N as we move
away from the Galactic plane, we have experimented with changing the range of
scales used for fitting. We find that the shallow SPS slope at high latitudes
persists and can not be explained as being caused by the noise SPS in low S/N
regions.

(a)  Figure 6: The histogram (shown in green) of the SPS slope image shown in
the top of Figure 5. The on-plane histogram corresponds to latitudes $|b|\leq
10\degree$ is shown in blue, whereas the off-plane histogram pertains to
latitudes $|b|\geq 70\degree$ and is shown in orange.The vertical dashed lines
represent the mean value for the three samples.

In Figure 7 we show the SPS slope values as a function of Galactic latitude.
For each 2 degree wide Galactic latitude bin the average SPS slope value is
calculated from the SPS map shown in Figure 5 (top panel). Since the rolling
SPS method uses a 3-degree wide kernel, shifted by 8 arcmins, the SPS slope
values in Figure 5 are correlated. The error bars in Figure 7 are estimated
from the standard deviation of SPS slope values within each 2 degree Galactic
latitude range and divided by $\sqrt{(N^{\prime}-1)}$. Here, $N^{\prime}$ is
an estimate of the number of independent kernels within each 2-degree Galactic
latitude zone, given by the number of pixels in Figure 5 that fall within each
2-degree latitude bin divided by the number of pixels per kernel, which is
(180/8)2.

Figure 7 shows that the SPS slope remains roughly constant at $\sim-2.9$
across the entire latitude range. For $|b|>60$ degrees, the slope becomes
slightly more shallow. Our mean SPS slope is in excellent agreement with what
Kalberla & Haud (2019) found by using the all-sky HI4PI survey. Using the
entire data set they found a SPS slope of $-2.943\pm 0.003$ and $-2.944\pm
0.005$ when considering only $|b|>20$ degrees (for the velocity range from
$-8$ to 8 km s-1). Figure 7 hints at a slight asymmetry between the northern
and southern hemispheres. While on the north side we see a gradual change in
the SPS slope from $b\sim 0^{\circ}$ to $b\sim 80^{\circ}$, on the southern
side the SPS slope is relatively constant from $b\sim 0^{\circ}$ to
$-50^{\circ}$, and then has a sharp decrease from $\sim-2.9$ to $\sim-2.6$.
This asymmetry in SPS slope between north and south latitudes likely stems
from the asymmetric HI  distribution in the MW (Burton, 1988; Levine et al.,
2006a, b; Kalberla et al., 2007; Kalberla & Kerp, 2009).

###  4.3 Comparison of the SPS slope with previous studies

(a)  Figure 7: The HI  SPS slope, averaged in 2∘ latitude bins, as a function
of Galactic latitude. Green points show results from previous studies: 1-
Pingel et al. (2013), 2- Kalberla et al. (2017), 3- Miville-Deschênes et al.
(2003), 4- Blagrave et al. (2017), 5- Martin et al. (2015), 6- Khalil et al.
(2006), 7- Dickey et al. (2001), 8- Kalberla & Haud (2019), 9- Marchal &
Miville-Deschênes (2021), 10- Chepurnov et al. (2010), 11- Pingel et al.
(2018), 12- Choudhuri & Roy (2019).

In Figure 7 we also compare our SPS slope with results from previous studies.
There have been many studies of the HI  SPS mostly focusing on smaller
regions. While most studies shown in this figure have used the Fourier
transform method to calculate the SPS, the exact details of the methods used
and the fitting range, as well as the velocity range, vary greatly, making a
comparison of this heterogeneous sample difficult. Our work is the first to
apply the same methodology uniformly over a large spatial area, about 13,000
square degrees.

We note that all these measurements were obtained from HI  images integrated
over a range of velocity channels that roughly corresponds to the local HI
(excluding intermediate and high-velocity gas) and likely represent estimates
of density fluctuations (Lazarian & Pogosyan, 2000; Yuen et al., 2021). Our
SPS slope values are generally consistent with many prior findings. Two high-
latitude HI  clouds, MBM16 and Ursa Major, stand out as having a much steeper
SPS slope of $-3.7$, as well as one region in the Galactic plane. MBM16 is
especially interesting as it is a starless cloud. Both Pingel et al. (2013)
and Miville-Deschênes et al. (2003) hypothesized that the steep SPS slope
could be caused by the lack of stellar feedback. However, it is likely that
the smallest spatial scales in these studies were affected by the beam shape
and resulted in a steep SPS. Kalberla et al. (2017) re-examined the MBM16
power spectrum and estimated the SPS slope of $-3.63\pm 0.04$ after applying
apodization, beam correction, and the noise bias correction (point 2 in Figure
7). Therefore, the MBM16 region remains as being an example of an unusually
steep HI  SPS slope. Further studies are needed to understand the origin of
this steep slope. On the other hand, for the Ursa Major field, Blagrave et al.
(2017) obtained SPS slope of $-2.68\pm 0.14$ after accounting for instrumental
effects (shown as “UM” in Figure 7).

###  4.4 Accounting for high optical depth

When calculating the SPS slope we use the column density images by assuming
that HI  is optically thin. This assumption is reasonable for high-latitude HI
, but certainly underestimates the HI  column density close to and within the
MW plane. To investigate whether the SPS slope distribution is affected by the
spatial variations of the HI  optical depth we applied the correction for high
optical depth on the entire GALFA-HI  field. As shown by Lee et al. (2015) and
Nguyen et al. (2019), there are significant regional variations of this
correction factor. For example, Nguyen et al. (2019) investigated HI  in the
vicinity of several Giant Molecular Clouds (GMCs) and found that the
correction factor ($f$) from their best fit is:

$\begin{split}f&=(0.47\pm 0.09)\times\log_{10}(N(HI)/10^{20})\\\ &+(0.66\pm
0.12).\end{split}$ (5)

However, for sightlines through the Galactic plane area ($|b|<5^{\circ}$), the
correction factor increased rapidly as $f=(2.41\pm
0.93)\times\log_{10}(N(HI)/10^{20})+(2.57\pm 1.57)$. At low column densities
below $\sim 10^{21}$ cm-2, the correction factor was close to unity.

As a first-order test of how the optical depth correction affects spatial
variations of the SPS slope, we applied the correction from Equation (5) on
the entire GALFA-HI  column density image. This affected the SPS slope only
for $|b|<20^{\circ}$ by making the SPS slope slightly more shallow. At
latitudes $|b|\geq 20^{\circ}$, the SPS slope remained unchanged. The observed
change of the SPS slope is, however, within the SPS slope uncertainties
(typical error bars shown in Figure 7 are $\sim 0.1$).

We conclude that the optical depth is unlikely to modify significantly the
observed relatively uniform large-scale distribution of the SPS slope. While
our attempted correction for high optical depth is not realistic at very low
latitudes ($|b|<10^{\circ}$), we hypothesize that the full correction would
introduce more small-scale structure resulting in a more shallow SPS slope
within $10^{\circ}$ from the plane where currently we see the steepest SPS
slope (see Figure 7). In essence, the optical depth correction would likely
further diminish spatial variations of the SPS slope seen between low and
intermediate Galactic latitudes.

###  4.5 Accounting for the plane parallel distribution of HI

(a)  Figure 8:  Top panel: The HI  column density integrated over $\pm 40$ km
s-1 velocity range multiplied by $\sin|b|$ is shown as black data points with
respect to $\sin|b|$. Bottom panel: the SPS slope of the $N(HI)\times\sin|b|$
distribution, averaged in 2-degree latitude bins, as a function of $\sin|b|$.

The HI  distribution in the MW has a rough plane-parallel geometry (Dickey &
Lockman, 1990). This results in the fact that lines of sight at different
Galactic latitudes probe HI  along different path lengths. Therefore, the HI
column density changes significantly as a function of Galactic latitude. We
investigate here how the plane-parallel geometry of the HI  disk affects the
SPS slope. We therefore multiply $N(HI)$ with $\sin|b|$ as this effectively
projects the HI  column density along the vertical axis and therefore
normalizes $N(HI)$ for the difference in path lengths. The
$N(HI)\times\sin|b|$ distribution is shown in the top panel of Figure 8 as a
function of $\sin|b|$, while the rolling SPS slope calculated for the
$N(HI)\times\sin|b|$ distribution is shown in the bottom panel of the same
figure.

The $N(HI)\times\sin|b|$ distribution (top panel of Figure 8) is in agreement
with several previous studies and shows several features (Dickey & Lockman,
1990; Knapp, 1975). Knapp (1975) used the velocity range of $\pm 50$ km s-1
and showed that the $N(HI)\times\sin|b|$ distribution is flat at $2.6\times
10^{20}$ cm-2 up to $\sin|b|\sim 0.6$, and then decreases roughly by a factor
of two (reaching $1.5\times 10^{20}$ cm-2 at $|b|=85^{\circ}$). Dickey &
Lockman (1990) combined data from several HI  surveys over the velocity range
$\pm 250$ km s-1 and found that the decrease of $N(HI)\times\sin|b|$ at high
latitudes follows: $N(HI)=3.84\times\csc|b|-2.11$ in units of $10^{20}$ cm-2.
The excellent agreement between our study and both Knapp (1975) and Dickey &
Lockman (1990) the $N(HI)\times\sin|b|$ distribution suggests that our
velocity range of $\pm 40$ km s-1 encompasses the bulk of HI  emission and
that the decline of the column density at high latitudes is not due to the
selected velocity range.

The $N(HI)\times\sin|b|$ distribution goes to 0 at $b\sim 0^{\circ}$ as the MW
HI  disk has a finite size. For $\sin|b|\sim 0.1$-0.8 the $N(HI)\times\sin|b|$
distribution is relatively uniform. For $|b|>50^{\circ}$ the HI  column
density decreases by a factor of $\sim 2$ as it is affected by the presence of
the Local Bubble which contains largely hot gas (Cox & Reynolds, 1987) and
some cooler HI  (Redfield & Linsky, 2007) but this HI  distribution is
modified by the explosion of multiple supernovae that formed the Local Bubble.
While lines of sight at $|b|\leq 50^{\circ}$ predominantly probe the disk HI ,
at least half of the path length at $|b|>50^{\circ}$ is contained inside the
Local Bubble (the scale height of the cold/warm HI  is about 300/500 pc near
the Sun (McClure-Griffiths et al., 2023), while the radius of the Local Bubble
is $\sim 200$ pc (Welsh et al., 2010)) and therefore has a different fraction
of warm to cold HI  relative to the HI  at $|b|<50^{\circ}$.

As shown in the bottom panel of Figure 8, the SPS slope of the
$N(HI)\times\sin|b|$ distribution is essentially constant at $-2.85\pm 0.02$.
This is in agreement with the mostly constant $N(HI)\times\sin|b|$
distribution caused by the plane-parallel geometry of the MW HI . Only at
$\sin|b|>0.8$ where the $N(HI)$ is confined to the Local Bubble and the
fraction of cold vs. warm HI is different, we do find a small difference,
$-2.72\pm 0.03$. While this difference is barely at a 3-sigma level, it could
suggest the different fraction of HI  phases within the Local Bubble and/or
different turbulent state of the HI  inside the Local Bubble.

##  5 Discussion and summary

In this study, we calculated the HI  SPS of the HI  column density image
(integrated over local velocities $\pm 40$ km s-1) that covers close to 13,000
square degrees. As the GALFA-HI  survey has complex scanning patterns and
noise statistics due to its receiver cluster, we developed a methodology to
quantify and parameterise the noise (and beam) contributions to the HI  SPS
across the entire survey area. Using the rolling SPS method, we mapped the SPS
slope spatially across the GALFA-HI  field by rolling a 3-degree kernel and
consistently treating beam and noise contributions to the SPS within each
kernel.

We focused only on the HI  column density, relying on the assumption that the
$\pm 40$ km s-1 velocity range depicts the bulk of the HI  mass and that the
HI column density is a reliable tracer of 3D density fluctuations. While
previous studies (e.g., Lazarian & Pogosyan, 2000) have calculated the HI  SPS
from column density, the question regarding the exact contributions of density
vs. velocity fluctuations to the observed intensity fluctuations has received
significant attention in recent years. For example, the approach based on
varying the width of velocity channels to separate density from velocity
fluctuations (Lazarian & Pogosyan, 2000) was recently updated with a more
complex treatment (Yuen et al., 2021). On the other hand, Clark et al. (2019)
and Kalberla et al. (2016) suggested that the HI  intensity fluctuations seen
at different velocities are phase dependent, with narrow velocity channels
being dominated by the CNM and causing the HI  SPS slope to be more shallow
than what is found for the integrated intensity images. In addition, Körtgen
et al. (2021) used simulations to show that the HI  SPS depends on the multi-
phase medium. While understanding of direct contributions from density,
velocity, and different HI  phases to the HI  SPS awaits further development
and tests, we emphasize that our approach provides a systematic treatment of
the HI  SPS over a range of different environments with varying S/N ratios and
therefore enables a comprehensive examination of spatial variations of the HI
SPS slope.

The key result from our study is that we find a relatively uniform SPS slope
across most of the field covered by the GALFA-HI  survey. When projecting the
HI  column density along the vertical axis, by multiplying the HI  column
density by $\sin|b|$, the SPS slope of the $N(HI)\times\sin|b|$ distribution
is remarkably flat (within 1-$\sigma$ error bars) for $|b|\leq 60^{\circ}$.
Only at high Galactic latitudes, $|b|>60^{\circ}$, we find a slightly more
shallow SPS slope. The difference between low and high Galactic latitudes is
less pronounced when we consider the high-optical-depth correction and the
difference in the line-of-sight length caused by the plane parallel geometry
of the HI  disk. However, the more shallow SPS slope at $|b|>60^{\circ}$,
relative to lower latitudes, remains statistically significant at about
3-$\sigma$ level. At high latitudes the lines-of-sight are mostly within the
local ISM where the Local Bubble represents a prominent structure. On the
other hand, at lower latitudes the lines-of-sight probe much longer path-
lengths throughout the Galactic disk. Even after scaling the path-lengths for
a direct comparison of the HI  column density, it appears that the local ISM
has distinct turbulent properties. Our results therefore suggest that HI
turbulent properties within the Local Bubble are modified by the multiple
supernovae explosions.

The relatively uniform SPS slope across the majority of the GALFA-HI field is
a somewhat surprising finding when considering that the GALFA-HI  survey
probes many different interstellar environments. For example, the HI  column
density varies from $\sim 10^{20}$ to $\sim 10^{22}$ cm-2, tracing regions
with diffuse HI  all the way to regions with a significant fraction of
molecular gas (e.g. near molecular clouds such as Perseus, Taurus and
California). Based on HI  absorption measurements, it has been shown that the
CNM fraction varies across the GALFA-HI  field between $\sim 0.2$ and 0.75
(Nguyen et al., 2019). Yet, the SPS slope remains largely around $-2.9$. This
may suggest that the phase distribution, fractions of cold vs. warm HI , have
a secondary role in explaining SPS properties.

The relatively uniform SPS slope could be caused by the dominance of large-
scale turbulent driving, on Galaxy-wide spatial scales, introduced by a
combination of Galactic rotation, gravitational instabilities, and/or infall
onto the Galaxy (Wada et al., 2002; Block et al., 2010; Krumholz & Burkhart,
2016). Numerical simulations by Yoo & Cho (2014) showed that the turbulent
energy spectrum is very sensitive to large-scale driving. Unless the energy
injection rates on small scales are much higher than the energy injection rate
on large scales, the large-scale driving will always dominate. Therefore, even
in the presence of significant stellar feedback in the plane of the MW, the
large-scale turbulent driving could be the dominant mechanism. The finding of
a relatively uniform SPS slope for the MW resembles the HI  turbulent
properties within the Small Magellanic Cloud where recently Szotkowski et al.
(2019) found a very uniform SPS slope across the entire galaxy. In addition,
Yuen et al. (2022) applied the new method for separating density and velocity
fluctuations from spectral-line observations on HI  and CO observations of
several Galactic regions and showed that the turbulent velocity cascade has a
universal spectrum over a broad range of scales (from $\sim 1$ pc to $\sim
10^{3}$ pc). They concluded that this supports the picture where turbulence is
driven on large spatial scales, while the energy injection on small spatial
scales may not influence significantly the turbulent cascade. This study found
that the velocity power spectrum slope was more uniform than the slope of the
density fields.

Furthermore, Hennebelle & Falgarone (2012) found that the kinetic energy
transfer rate, observed in the molecular cloud population traced by 12CO,
exhibited no change in magnitude across a spatial range of $\sim 0.01$ to
$\sim 100$ pc. This led them to conclude that molecular clouds are part of the
same turbulent cascade as HI . If molecular clouds of diverse star formation
and stellar feedback levels exhibit negligible changes in their internal
turbulent characteristics, it is not surprising that the HI  also experiences
limited effects from small-scale turbulent driving, which could explain the
relatively consistent HI  SPS slope. A similar conclusion was reached in
Elmegreen et al. (2022) who investigated the HI  turbulent properties in a
sample of nearby galaxies and suggested that feedback does not affect much
atomic gas, hypothesizing that most of the feedback energy goes into dense
molecular clouds where it gets radiated efficiently.

As shown in Figure 4, we do not see any difference between inner and outer
Galaxy HI  turbulent properties when considering the SPS slope. While the
outer Galaxy region includes more distant HI  the inner Galaxy region is
probing HI  largely within a few kpc from the Sun. This is different from Dib
et al. (2021) who studied HI  turbulent properties of 33 disk galaxies using
the $\delta$-variance approach and concluded that inner parts are affected
largely by supernova feedback while outer parts are mostly affected by large-
scale dynamics. We note that this study has a much lower spatial resolution
and has calculated turbulent properties by averaging data for entire galaxies.

This study is supported by the National Aeronautics and Space Administration
under Grant No. 4200766703, proposal Stanimirovic (2020). We also acknowledge
support provided by the University of Wisconsin-Madison Office of the Vice
Chancellor for Research and Graduate Education with funding from the Wisconsin
Alumni Research Foundation. We are very grateful to Daniel Rybarczyk and Trey
Wenger for useful discussions at various stages of this project. We are
grateful to Peter Kalberla for providing test data and a code for comparison
of results. We thank the referee for constructive and insightful comments and
suggestions.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Beresnyak et al. (2005) Beresnyak, A., Lazarian, A., & Cho, J. 2005, ApJ, 624, L93, doi: 10.1086/430702
  * Bialy et al. (2017) Bialy, S., Burkhart, B., & Sternberg, A. 2017, ApJ, 843, 92, doi: 10.3847/1538-4357/aa7854
  * Blagrave et al. (2017) Blagrave, K., Martin, P. G., Joncas, G., et al. 2017, ApJ, 834, 126, doi: 10.3847/1538-4357/834/2/126
  * Block et al. (2010) Block, D. L., Puerari, I., Elmegreen, B. G., & Bournaud, F. 2010, ApJ, 718, L1, doi: 10.1088/2041-8205/718/1/L1
  * Bracewell (1986) Bracewell, R. N. 1986, The Fourier transform and its applications, Vol. 31999 (McGraw-Hill New York)
  * Burkhart et al. (2015) Burkhart, B., Lazarian, A., Balsara, D., Meyer, C., & Cho, J. 2015, ApJ, 805, 118, doi: 10.1088/0004-637X/805/2/118
  * Burton (1988) Burton, W. B. 1988, in Galactic and Extragalactic Radio Astronomy, ed. K. I. Kellermann & G. L. Verschuur, 295–358
  * Chepurnov et al. (2010) Chepurnov, A., Lazarian, A., Stanimirović, S., Heiles, C., & Peek, J. E. G. 2010, ApJ, 714, 1398, doi: 10.1088/0004-637X/714/2/1398
  * Choudhuri & Roy (2019) Choudhuri, S., & Roy, N. 2019, MNRAS, 483, 3437, doi: 10.1093/mnras/sty3342
  * Clark et al. (2019) Clark, S. E., Peek, J. E. G., & Miville-Deschênes, M. A. 2019, ApJ, 874, 171, doi: 10.3847/1538-4357/ab0b3b
  * Combes et al. (2012) Combes, F., Boquien, M., Kramer, C., et al. 2012, A&A, 539, A67, doi: 10.1051/0004-6361/201118282
  * Cox & Reynolds (1987) Cox, D. P., & Reynolds, R. J. 1987, ARA&A, 25, 303, doi: 10.1146/annurev.aa.25.090187.001511
  * Crovisier & Dickey (1983) Crovisier, J., & Dickey, J. M. 1983, A&A, 122, 282
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792, doi: 10.1086/318388
  * de Avillez & Breitschwerdt (2005) de Avillez, M. A., & Breitschwerdt, D. 2005, A&A, 436, 585, doi: 10.1051/0004-6361:20042146
  * Dib et al. (2021) Dib, S., Braine, J., Gopinathan, M., et al. 2021, A&A, 655, A101, doi: 10.1051/0004-6361/202141803
  * Dickey & Lockman (1990) Dickey, J. M., & Lockman, F. J. 1990, ARA&A, 28, 215, doi: 10.1146/annurev.aa.28.090190.001243
  * Dickey et al. (2001) Dickey, J. M., McClure-Griffiths, N. M., Stanimirović, S., Gaensler, B. M., & Green, A. J. 2001, ApJ, 561, 264, doi: 10.1086/323409
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium
  * Elmegreen et al. (2000) Elmegreen, B. G., Kim, S., & Staveley-Smith, L. 2000, in American Astronomical Society Meeting Abstracts, Vol. 197, American Astronomical Society Meeting Abstracts, 112.06
  * Elmegreen et al. (2001) Elmegreen, B. G., Kim, S., & Staveley-Smith, L. 2001, ApJ, 548, 749, doi: 10.1086/319021
  * Elmegreen et al. (2022) Elmegreen, B. G., Martinez, Z., & Hunter, D. A. 2022, ApJ, 928, 143, doi: 10.3847/1538-4357/ac559c
  * Elmegreen & Scalo (2004) Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211, doi: 10.1146/annurev.astro.41.011802.094859
  * Federrath et al. (2010) Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M. M. 2010, A&A, 512, A81, doi: 10.1051/0004-6361/200912437
  * Ferland et al. (2013) Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mexicana Astron. Astrofis., 49, 137.  https://arxiv.org/abs/1302.4485
  * Gooch (1996) Gooch, R. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 80
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Green (1993) Green, D. A. 1993, MNRAS, 262, 327, doi: 10.1093/mnras/262.2.327
  * Grisdale et al. (2018) Grisdale, K., Agertz, O., Renaud, F., & Romeo, A. B. 2018, MNRAS, 479, 3167, doi: 10.1093/mnras/sty1595
  * Grisdale et al. (2017) Grisdale, K., Agertz, O., Romeo, A. B., Renaud, F., & Read, J. I. 2017, MNRAS, 466, 1093, doi: 10.1093/mnras/stw3133
  * Hennebelle & Falgarone (2012) Hennebelle, P., & Falgarone, E. 2012, A&A Rev., 20, 55, doi: 10.1007/s00159-012-0055-y
  * Hunter et al. (2012) Hunter, D. A., Ficut-Vicas, D., Ashley, T., et al. 2012, AJ, 144, 134, doi: 10.1088/0004-6256/144/5/134
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Joung & Mac Low (2006) Joung, M. K. R., & Mac Low, M.-M. 2006, ApJ, 653, 1266, doi: 10.1086/508795
  * Joye & Mandel (2003) Joye, W. A., & Mandel, E. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook, 489
  * Kalberla et al. (2007) Kalberla, P. M. W., Dedes, L., Kerp, J., & Haud, U. 2007, A&A, 469, 511, doi: 10.1051/0004-6361:20066362
  * Kalberla & Haud (2019) Kalberla, P. M. W., & Haud, U. 2019, A&A, 627, A112, doi: 10.1051/0004-6361/201834533
  * Kalberla & Kerp (2009) Kalberla, P. M. W., & Kerp, J. 2009, ARA&A, 47, 27, doi: 10.1146/annurev-astro-082708-101823
  * Kalberla et al. (2017) Kalberla, P. M. W., Kerp, J., Haud, U., & Haverkorn, M. 2017, A&A, 607, A15, doi: 10.1051/0004-6361/201629627
  * Kalberla et al. (2016) Kalberla, P. M. W., Kerp, J., Haud, U., et al. 2016, ApJ, 821, 117, doi: 10.3847/0004-637X/821/2/117
  * Khalil et al. (2006) Khalil, A., Joncas, G., Nekka, F., Kestener, P., & Arneodo, A. 2006, ApJS, 165, 512, doi: 10.1086/505144
  * Kim et al. (2001) Kim, J., Balsara, D., & Mac Low, M.-M. 2001, Journal of Korean Astronomical Society, 34, 333, doi: 10.5303/JKAS.2001.34.4.333
  * Klessen & Hennebelle (2010) Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17, doi: 10.1051/0004-6361/200913780
  * Knapp (1975) Knapp, G. R. 1975, AJ, 80, 111, doi: 10.1086/111719
  * Koch et al. (2020) Koch, E. W., Chiang, I.-D., Utomo, D., et al. 2020, MNRAS, 492, 2663, doi: 10.1093/mnras/stz3582
  * Körtgen et al. (2021) Körtgen, B., Pingel, N., & Killerby-Smith, N. 2021, MNRAS, 505, 1972, doi: 10.1093/mnras/stab1408
  * Kowal & Lazarian (2007) Kowal, G., & Lazarian, A. 2007, ApJ, 666, L69, doi: 10.1086/521788
  * Krumholz & Burkhart (2016) Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671, doi: 10.1093/mnras/stw434
  * Lazarian (2009) Lazarian, A. 2009, Space Sci. Rev., 143, 357, doi: 10.1007/s11214-008-9460-y
  * Lazarian & Pogosyan (2000) Lazarian, A., & Pogosyan, D. 2000, ApJ, 537, 720, doi: 10.1086/309040
  * Lee et al. (2015) Lee, M.-Y., Stanimirović, S., Murray, C. E., Heiles, C., & Miller, J. 2015, ApJ, 809, 56, doi: 10.1088/0004-637X/809/1/56
  * Levine et al. (2006a) Levine, E. S., Blitz, L., & Heiles, C. 2006a, ApJ, 643, 881, doi: 10.1086/503091
  * Levine et al. (2006b) —. 2006b, Science, 312, 1773, doi: 10.1126/science.1128455
  * Marchal & Miville-Deschênes (2021) Marchal, A., & Miville-Deschênes, M.-A. 2021, ApJ, 908, 186, doi: 10.3847/1538-4357/abd108
  * Martin et al. (2015) Martin, P. G., Blagrave, K. P. M., Lockman, F. J., et al. 2015, ApJ, 809, 153, doi: 10.1088/0004-637X/809/2/153
  * McClure-Griffiths et al. (2023) McClure-Griffiths, N. M., Stanimirović, S., & Rybarczyk, D. R. 2023, ARA&A, 61, 19, doi: 10.1146/annurev-astro-052920-104851
  * McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565, doi: 10.1146/annurev.astro.45.051806.110602
  * Miville-Deschênes et al. (2003) Miville-Deschênes, M. A., Joncas, G., Falgarone, E., & Boulanger, F. 2003, A&A, 411, 109, doi: 10.1051/0004-6361:20031297
  * Muller et al. (2004) Muller, E., Stanimirović, S., Rosolowsky, E., & Staveley-Smith, L. 2004, ApJ, 616, 845, doi: 10.1086/425154
  * Nguyen et al. (2019) Nguyen, H., Dawson, J. R., Lee, M.-Y., et al. 2019, ApJ, 880, 141, doi: 10.3847/1538-4357/ab2b9f
  * Oliphant (2007) Oliphant, T. E. 2007, Computing in Science and Engineering, 9, 10, doi: 10.1109/MCSE.2007.58
  * Ostriker & Shetty (2011) Ostriker, E. C., & Shetty, R. 2011, ApJ, 731, 41, doi: 10.1088/0004-637X/731/1/41
  * Padoan et al. (2012) Padoan, P., Haugbølle, T., & Nordlund, Å. 2012, ApJ, 759, L27, doi: 10.1088/2041-8205/759/2/L27
  * Peek et al. (2011) Peek, J. E. G., Heiles, C., Douglas, K. A., et al. 2011, ApJS, 194, 20, doi: 10.1088/0067-0049/194/2/20
  * Peek et al. (2018) Peek, J. E. G., Babler, B. L., Zheng, Y., et al. 2018, ApJS, 234, 2, doi: 10.3847/1538-4365/aa91d3
  * Pilkington et al. (2011) Pilkington, K., Gibson, B. K., Calura, F., et al. 2011, MNRAS, 417, 2891, doi: 10.1111/j.1365-2966.2011.19450.x
  * Pingel et al. (2018) Pingel, N. M., Lee, M.-Y., Burkhart, B., & Stanimirović, S. 2018, ApJ, 856, 136, doi: 10.3847/1538-4357/aab34b
  * Pingel et al. (2013) Pingel, N. M., Stanimirović, S., Peek, J. E. G., et al. 2013, ApJ, 779, 36, doi: 10.1088/0004-637X/779/1/36
  * Pingel et al. (2022) Pingel, N. M., Dempsey, J., McClure-Griffiths, N. M., et al. 2022, PASA, 39, e005, doi: 10.1017/pasa.2021.59
  * Redfield & Linsky (2007) Redfield, S., & Linsky, J. L. 2007, arXiv e-prints, arXiv:0709.4480, doi: 10.48550/arXiv.0709.4480
  * Shetty et al. (2012) Shetty, R., Beaumont, C. N., Burton, M. G., Kelly, B. C., & Klessen, R. S. 2012, MNRAS, 425, 720, doi: 10.1111/j.1365-2966.2012.21588.x
  * Stanimirovic (2020) Stanimirovic, S. 2020, Shaken or stirred: How stellar feedback drives interstellar turbulence?, NASA Proposal id.20-ADAP20-159
  * Stanimirović & Lazarian (2001) Stanimirović, S., & Lazarian, A. 2001, ApJ, 551, L53, doi: 10.1086/319837
  * Stanimirovic et al. (1999) Stanimirovic, S., Staveley-Smith, L., Dickey, J. M., Sault, R. J., & Snowden, S. L. 1999, MNRAS, 302, 417, doi: 10.1046/j.1365-8711.1999.02013.x
  * Stutzki et al. (1998) Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998, A&A, 336, 697
  * Szotkowski et al. (2019) Szotkowski, S., Yoder, D., Stanimirović, S., et al. 2019, ApJ, 887, 111, doi: 10.3847/1538-4357/ab53df
  * Tamburro et al. (2009) Tamburro, D., Rix, H. W., Leroy, A. K., et al. 2009, AJ, 137, 4424, doi: 10.1088/0004-6256/137/5/4424
  * Taylor (1997) Taylor, J. 1997, Introduction to Error Analysis, the Study of Uncertainties in Physical Measurements, 2nd Edition
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Wada et al. (2002) Wada, K., Meurer, G., & Norman, C. A. 2002, ApJ, 577, 197, doi: 10.1086/342151
  * Welsh et al. (2010) Welsh, B. Y., Lallement, R., Vergely, J. L., & Raimond, S. 2010, A&A, 510, A54, doi: 10.1051/0004-6361/200913202
  * Yoo & Cho (2014) Yoo, H., & Cho, J. 2014, ApJ, 780, 99, doi: 10.1088/0004-637X/780/1/99
  * Yuen et al. (2022) Yuen, K. H., Ho, K. W., Law, C. Y., Chen, A., & Lazarian, A. 2022, arXiv e-prints, arXiv:2204.13760.  https://arxiv.org/abs/2204.13760
  * Yuen et al. (2021) Yuen, K. H., Ho, K. W., & Lazarian, A. 2021, ApJ, 910, 161, doi: 10.3847/1538-4357/abe4d4
  * Zamora-Avilés et al. (2012) Zamora-Avilés, M., Vázquez-Semadeni, E., & Colín, P. 2012, ApJ, 751, 77, doi: 10.1088/0004-637X/751/1/77
  * Zhang et al. (2012) Zhang, H.-X., Hunter, D. A., & Elmegreen, B. G. 2012, ApJ, 754, 29, doi: 10.1088/0004-637X/754/1/29

##  Appendix A Noise and telescope beam contributions to the SPS

###  A.1 Noise analysis

The calculation of the HI  SPS for a given kernel must take into account both
the intricate ALFA beam and the noise in the GALFA-HI  data, which varies in
both the spectral and spatial dimensions. A detailed understanding of the
noise is particularly crucial in regions at mid and high latitudes, which
represent roughly 25$\%$ of the observed sky, as the signal to noise ratio
decreases substantially with Galactic latitude and the noise contribution can
impact the slope of the SPS.

The GALFA-HI  survey resulted from merging various HI  observations obtained
via different scanning methods. Although global calibration solutions were
derived to address the relative calibration of different regions (Peek et al.,
2011), the residuals from distinct scanning patterns persist and constitute a
form of systematic noise. To account for the noise contribution during the
fitting of the rolling HI  SPS, we generate a noise image,
$N_{\rm{total}}(x,y)$, which serves as the basis for computing the rolling
noise SPS, $P_{\rm{noise}}(k)$. The GALFA-HI  column density image encompasses
multiple noise components, including receiver (rx) noise, noise stemming from
the brightness temperature of the incident HI  emission (both associated with
the system temperature, $T_{\rm{sys}}$), as well as systematic noise
originating from calibration artifacts resulting from various scanning
techniques employed in the GALFA-HI  survey. The calibration noise consists of
both an root mean square (rms) component (stemming from the scanning and data
merging process) and a systematic offset component (arising from the zero-
point calibration). As it is impossible to completely separate the random
noise due to receiver temperature ($T_{\rm{rx}}$) from the systematic noise
associated with calibration, we choose to divide the total noise image
corresponding to the HI  column density image into two components:

$N_{\rm{total}}(x,y)=N_{\rm{off}}(x,y)+N_{\rm{on-HI}}(x,y)$ (A1)

where $N_{\rm{off}}$ is the combined contribution of random fluctuations due
to receiver temperature plus the systematic scanning calibration uncertainty
(this image can be estimated from the emission-free channels of the GALFA-HI
data cube), and $N_{\rm{on-HI}}$ is the noise contribution due to the HI
brightness temperature.

$N_{\rm{off}}(x,y):$ The noise from the ALFA receiver, characterized by
$T_{\rm{rx}}$ = 30 K (Peek et al., 2018), is independent of velocity and gives
rise to the radiometer noise, with a standard deviation of $\sigma_{\rm{rx}}$.
To obtain $\sigma_{\rm{rx}}$, we select the velocity range nearly free of HI
emission, i.e., $-520$ to $-440$ km s-1 (we refer to this subset as the “non-
HI ” data cube, which shares the same velocity width as our LVC data cube).
For each spatial pixel in this cube, we calculate the standard deviation over
the velocity axis. This yields the $\sigma_{\rm{rx}}(x,y)$ image (shown in top
3 panels of Figure 9 ), whose pixel values range from 1.4 to 2.3 K.

We obtain the HI  column density of the emission-free data cube by integrating
the brightness temperature across the HI -free velocity range of $-520$ to
$-440$ km s-1, resulting in the $N_{\rm{off}}(x,y)$ image. The fourth panel of
Figure 9 depicts this image, which has pixel values spanning from $-$7.6
$\times 10^{18}$ to $+$7.1 $\times 10^{18}$ cm-2 across 97$\%$ of the area.
Although it contains receiver noise, its spatial structure is largely
determined by the scanning and calibration artifacts.

Given the challenge of differentiating the systematic calibration noise due to
scanning and merging of different observations, we note that the non-HI  noise
image involves contributions from both the $\sigma_{\rm{rx}}(x,y)$ image
(receiver noise) and systematic calibration offsets.

$N_{\rm{on-HI}}$: The noise term due to the HI  brightness temperature is
radial velocity dependent and can be expressed by scaling the standard
deviation estimated over emission-free channels ($\sigma_{\rm{rx}}$):

$\sigma_{\rm{rx+HI}}(x,y,v)=\sigma_{\rm{rx}}(x,y)\left(1+T_{b}(x,y,v)/T_{\rm{rx}}\right)$
(A2)

where $T_{b}(x,y,v)$ is the brightness temperature at each pixel location and
each velocity channel. $\sigma_{\rm{rx+HI}}$ is essentially a noise envelope
(or spectrum) calculated at each spatial pixel. We propagate
$\sigma_{\rm{rx+HI}}(x,y,v)$ to estimate the corresponding contribution to the
HI  column density image by adding noise velocity pixels in quadrature:

$\sigma^{\prime}_{\rm{rx+HI}}(x,y)=\sqrt{\Sigma_{v}({\sigma_{\rm{rx+HI}}(x,y,v)}^{2})}$
(A3)

The above method produces an image $\sigma^{\prime}_{\rm{rx+HI}}(x,y)$, where
each spatial pixel shows the uncertainty value for the HI column density
resulting from the combined effects of the ALFA receiver, the HI  brightness
temperature, and the scanning/calibration artifacts. Typical pixel values in
this image are 1.3 to 3.3 K km s-1 .

We now calculate the uncertainty image associated purely with the HI  line,
$\sigma_{\rm HI}$, as:

$\sigma^{\prime}_{\rm{rx+HI}}=\sqrt{\sigma_{\rm{rx}}^{2}+\sigma_{\rm{HI}}^{2}}$
(A4)

and therefore:

$\sigma_{\rm{HI}}=\sqrt{\sigma_{\rm{rx+HI}}^{{}^{\prime}2}-\sigma_{\rm{rx}}^{2}}.$
(A5)

Both $\sigma^{\prime}_{\rm{rx+HI}}$ and $\sigma_{\rm{HI}}$ are shown in Figure
9. To produce the $N_{\rm{on-HI}}(x,y)$ image (shown in Figure 9), for each
spatial pixel, a random number generator is used to draw a random number from
a normal distribution with standard deviation equal to
$\sigma_{\rm{HI}}(x,y)$.

Once the non-HI  ($N_{\rm{off}}$) and HI  noise ($N_{\rm{on-HI}}$) images are
combined, the final step in creating the noise image involves identifying and
replacing any extreme outliers (hot pixels with values $\geq|5\sigma|$). This
step is necessary to eliminate extreme sharp discontinuities that can cause
artificial effects on the SPS via Gibbs ringing.

###  A.2 Constructing the noise SPS data cube from the noise image

Starting with Equation (4):

$P_{\rm{obs}}(k)=P_{\rm{beam}}(k)\times P_{\rm{model}}+P_{\rm{noise}}(k)$ (A6)

where $P_{\rm{noise}}(k)$ is the SPS of the noise image, $P_{\rm{beam}}(k)$ is
the SPS of the ALFA’s beam, and $P_{\rm{model}}$ is the SPS of the underlying
HI  distribution unaffected by instrumental effects (noise, calibration), as
explained in Section 3.2.

For the rolling SPS derivation, we derive the noise SPS for each three degree
square sub-image from the full GALFA-HI noise image. To improve statistics for
any give noise SPS, we parameterize $P_{\rm{noise}}(k)$ as a function of the
standard deviation value of the input noise sub-image, $sd=\sigma(N_{\rm
total})$ estimated for each 3-degree kernel, and write it as
$P_{\rm{noise}}(k)=C\times Q(k,sd)$. Here, $C$ is a constant which varies as a
function of spatial location and $Q(k,sd)$ is a function dependent on both $k$
and the standard deviation ($sd$) of the input noise sub-image. Below, we
explain and justify the reason for this functional form.

We first perform the rolling SPS analysis on the noise image
$N_{\rm{total}}(x,y)$ to generate what we refer to as the raw noise SPS data
cube. This cube contains the SPS of each 3-degree kernel of the the
$N_{\rm{total}}(x,y)$ image. To enhance the statistics and smooth out these
SPS functions before employing them in Equation (A6), we compute an average
noise SPS, denoted as $Q(k,sd)$. Due to the scanning and calibration
artifacts, the primary source of noise in the GALFA-HI  data is the systematic
noise (as seen in the bottom panel of Figure 9). Hence, the averaged noise SPS
is significantly influenced by $\sigma_{\rm{rx}}$. Since each kernel covers a
3-degree area, $\sigma(N_{\rm{total}})$ is used as a proxy for
$\sigma_{\rm{rx}}$. In Figure 10 we plot the $\log$ of
$\sigma(N_{\rm{total}})$ as a function of the $\log$ of $P_{\rm noise}(k)$ at
each $\log(k)$ value used in the SPS calculation. This figure illustrates that
for each $\log(k)$ value there is a linear relationship between
$\sigma(N_{\rm{total}})$ and $\log(P_{\rm noise})$. We therefore fit a linear
function for each $\log(k)$ bin, providing the means for predicting the noise
SPS of any noise sub-image by assessing $\sigma(N_{\rm{total}})$. As a result,
we have a family of about 100 $Q(k,sd)$ functions to represent $\log(P_{\rm
noise})$. Figure 11 demonstrates the functional change of $\log(P_{\rm
noise})$ as a function of $\sigma(N_{\rm{total}})$.

Figure 12 illustrates a sample of the raw noise SPS (curves black data points)
over-plotted with the parameterized $Q(k,sd)$ noise power spectrum (orange
line). As shown in Figure 9, the total noise image has complex spatial
structure and there are regions where the estimate of the local standard
deviation will not provide a good fit to the noise SPS data. There are several
such examples in Figure 12, e.g. the top left panel ($(l,b)=(354.93,83.355)$)
where the orange line is clearly offset slightly from the data. To account for
this, we estimate the constant value $C$ that is needed to shift the $Q(k,sd)$
noise power spectrum along the vertical axis and fit the data better. To
estimate $C$ for each 3-degree kernel, we calculate the mean difference
between $Q(k,sd)$ and $P_{\rm noise}(k)$ across different $k$ values. This
simple procedure results in the $C\times Q(k,sd)$ noise power spectrum that is
shown as a blue line in Figure 12.

(a)  (b)

Figure 9: (Top 3 panels) The standard deviation of the radiometer noise due to
receiver temperature ($\sigma_{\rm rx}$), the standard deviation due to the HI
— brightness temperature ($\sigma_{\rm HI}$), and the standard deviation due
to both receiver and HI brightness temperature ($\sigma_{\rm{rx+HI}}$).
(Bottom 3 panels) The noise contribution to the HI column density image:
random fluctuations due to receiver temperature plus the systematic scanning
uncertainty ($N_{\rm{off}}$), the contribution due to the brightness
temperature that is velocity dependent ($N_{\rm{on-HI}}$), and the sum of
these two contribution ($N_{\rm total}$).  (a)  Figure 10: The logarithmic
power for the noise SPS, $P_{\rm noise}(k)$, as a function of the logarithm of
$\sigma(N_{\rm total})$ for the first 20 binned $k$ values obtained for all
the 3-degree kernels are shown as grey data points. The mean values are
overplotted in black. For each $k$, a linear function $Q(sd)$ was fitted to
the data. (a)  Figure 11: The parametrized noise SPS, $Q(k,sd)$ for different
values of $\sigma(N_{\rm total})$.  (a)  Figure 12: Examples of the noise SPS
spectra for different kernels. The raw SPS data points, $P_{\rm noise}(k)$ are
shown as black data points. The parametrize noise SPS, $Q(k,sd)$ is shown in
yellow, while the final parametrized noise SPS, $C\times Q(k,sd)$, is shown in
blue. The SPS samples the range of angular scales from 3 degrees to 8 arcmin,
which corresponds to spatial frequencies $k=5\times 10^{-3}$ to $10^{-1}$
arcmin-1.

###  A.3 Final SPS after considering noise and beam effects

To demonstrate the entire process of accounting for the beam and noise
effects, in Figure 13, we present a sample of SPS curves calculated for
3-degree kernels at four distinct latitudes spanning from the Galactic plane
to the Galactic pole. The observed SPS, $P_{\rm obs}(k)$ is shown with red
data points in Figure 13. The symmetrical Gaussian beam shape, $P_{\rm
beam}(k)$ is shown as grey dashed curves. The noise SPS, $C\times Q(k,sd)$ is
displayed as yellow data points in the same figure. Finally, $[P_{\rm
obs}-C\times Q(k,sd)]/P_{\rm beam}$ is shown as green data points. The range
of $k$-values for fitting is restricted to 1.5-degree and 4 $\times$ FWHM for
all kernels, as explained in the text.

(a)  Figure 13: Examples of the HI  SPS for four different regions with
individual contributions from Equation (4). Within each panel, the grey,
yellow, red, blue, and green data points represent the beam, total noise,
observed, observed-beam, and final SPS, respectively. Furthermore, the dashed
lines in the figure indicate the range of fitting for both large and small
scales as 1.5°and 16′, which corresponds to spatial frequencies $k=10^{-2}$ to
$\sim 6\times 10^{-2}$ arcmin-1.