Neutral Hydrogen (H i) 21 cm as a Probe: Investigating Spatial Variations in Interstellar Turbulent Properties

Interstellar turbulence shapes the H i distribution in the Milky Way (MW). How this affects large-scale statistical properties of H i column density across the MW remains largely unconstrained. We use the ∼13,000 deg2 GALFA-H i survey to map statistical fluctuations of H i over the ±40 km s−1 velocity range. We calculate the spatial power spectrum (SPS) of the H i column density image by running a 3° kernel and measuring the SPS slope over a range of angular scales from 16′ to 20°. Due to GALFA’s complex observing and calibration strategy, we construct detailed estimates of the noise contribution and account for GALFA beam effects on the SPS. This allows us to systematically analyze H i images that trace a wide range of interstellar environments. We find that the SPS slope varies between ∼ −2.6 at high Galactic latitudes and ∼ −3.2 close to the Galactic plane. The range of SPS slope values becomes tighter when we consider H i optical depth and line-of-sight length caused by the plane-parallel geometry of the H i disk. This relatively uniform, large-scale distribution of the SPS slope is suggestive of large-scale turbulent driving being a dominant mechanism for shaping H i structures in the MW and/or the stellar feedback turbulence being efficiently dissipated within dense molecular clouds. Only at latitudes above 60° do we find evidence for the H i SPS slope being consistently more shallow. Those directions are largely within the Local Bubble, suggesting that the recent history of this cavity, shaped by multiple supernova explosions, has modified the turbulent state of H i and/or fractions of H i phases.

Over the past two decades, numerous observational (e.g., Crovisier & Dickey 1983;Stanimirovic et al. 1999) and numerical studies (e.g., Grisdale et al. 2017Grisdale et al. , 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 re-mains 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 (τ ) 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. 2013Pingel et al. , 2018)), the δ-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. 2013Pingel et al. , 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 ∼ −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.

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 ′ × 3.8 ′ .The survey area spans approximately 13,000 square degrees, including the declination range (DEC) of −1 • 20 ′ to 38 • 02 ′ , 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 |v| ≤ 40 km s −1 .We calculate the HI column density under the assumption of HI being optically thin.The HI column density, N (HI) thin , is given by: N HI,thin = C 0 ˆTB (v)dv. (1) where C 0 = 1.823 × 10 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 .3. METHODS: SPATIAL POWER SPECTRUM (SPS) AND THE ROLLING SPS

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: where k is the spatial frequency, measured in units of wavelength and being propor-tional to the inverse of the spatial scale1 , while L = x − 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 ∼ 4 ′ .We then apply a median filter with box size of 1 • to identify pixels with column density ≥5-σ, 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 ∼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, ⟨ℜ 2 + ℑ 2 ⟩ (where ℜ and ℑ 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

√
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: where A is the power-law amplitude, and γ 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).

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: (4) where we assume a 2D Gaussian beam shape with a FWHM 4 ′ × 4 ′ size, therefore P beam (k) ∝ exp(−4π 2 σ 2 beam k 2 ), and In the case of GALFA-HI , the noise term P 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 noise (k) = C × 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×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 ≤ k 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 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 × 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.

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 • 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 • to 16 ′ per single 3degree kernel.To obtain the SPS slope γ 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 ′ 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 χ 2 for each fit.

GALFA-HI for the inner, mid and outer Galactic regions
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 × 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 • to 16 ′ .
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 ∼ −3.2 ± 0.1 relative to the high-latitude regions (G4 and G5) which have a more shallow slope, ∼ −2.9 ± 0.1.In particular, the two most extreme slopes are found for regions G1 and G6, with corresponding SPS slopes being −3.37 ± 0.05 and −2.80 ± 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 highlatitude HI and to extend the range of an-gular scales even further, in Figure 3 we select three regions, each about 40 × 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 ∼ 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 largescale end of the fitting range, the highlatitude region has a much shallower SPS The GALFA-HI LVC column density image, divided into ten equal sub-regions named as G0 to G9 where each region is 36 × 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 • and 16 ′ , which corresponds to spatial frequencies k = 10 −3 to ∼ 6 × 10 −2 arcmin −1 ).slope relative to the low-latitude regions (e.g.−2.35 ± 0.17 vs −3.73 ± 0.17).This shows that the main difference between low-and high-latitude HI SPS occurs on scales larger than ∼ 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 lowand high-latitude regions is comparable, but the HI power on large scales is a factor of ∼ 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.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 noisedominated 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 χ 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.
On a 3-degree spatial scale, the SPS slope varies between −2.6 and −3.2 (1σ 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| ≤ 10 degrees, while in blue we show data for |b| ≥ 70 degrees.The mean SPS slope for these two sub-samples is ∼ −3.0 and ∼ −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.
In Figure 7 we show the SPS slope values as a function of Galactic latitude.For Here, N ′ 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 ∼ −2.9 across

Comparison of the SPS slope with previous studies
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 sam-  ple 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 Galac-tic 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 ± 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 ± 0.14 after accounting for instrumental effects (shown as "UM" in Figure 7).

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: f = (0.47 ± 0.09) × log 10 (N (HI)/10 20 ) + (0.66 ± 0.12). (5) However, for sightlines through the Galactic plane area (|b| < 5 • ), the correction factor increased rapidly as f = (2.41 ± 0.93) × log 10 (N (HI)/10 20 ) + (2.57 ± 1.57).At low column densities below ∼ 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 • by making the SPS slope slightly more shallow.At latitudes |b| ≥ 20 • , 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 ∼ 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 • ), we hypothesize that the full correction would introduce more small-scale structure resulting in a more shallow SPS slope within 10 • 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.

Accounting for the plane parallel distribution of HI
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) × 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) × sin |b| distribution is shown in the bottom panel of the same figure.
The N (HI) × 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 ±50 km s −1 and showed that the N (HI) × sin |b| distribution is flat at 2.6 × 10 20 cm −2 up to sin |b| ∼ 0.6, and then decreases roughly by a factor of two (reaching 1.5 × 10 20 cm −2 at |b| = 85 • ).Dickey & Lockman (1990) combined data from several HI surveys over the velocity range ±250 km s −1 and found that the decrease of N (HI) × sin |b| at high latitudes follows: N (HI) = 3.84×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) × sin |b| distribution suggests that our velocity range of ±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) × sin |b| distribution goes to 0 at b ∼ 0 • as the MW HI disk has a finite size.For sin |b| ∼ 0.1-0.8 the N (HI) × sin |b| distribution is relatively uniform.For |b| > 50 • the HI column density decreases by a factor of ∼ 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)

DISCUSSION AND SUMMARY
In this study, we calculated the HI SPS of the HI column density image (integrated over local velocities ±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 ±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) × sin |b| distribution is remarkably flat (within 1-σ error bars) for |b| ≤ 60 • .Only at high Galactic latitudes, |b| > 60 • , 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 lineof-sight length caused by the plane parallel geometry of the HI disk.However, the more shallow SPS slope at |b| > 60 • , relative to lower latitudes, remains statistically significant at about 3-σ 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 ∼ 10 20 to ∼ 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 ∼ 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 ∼ 1 pc to ∼ 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 12 CO, exhibited no change in magnitude across a spatial range of ∼ 0.01 to ∼ 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 δ-variance ap-proach 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.

ACKNOWLEDGMENTS
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.
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) Once the non-HI (N off ) and HI noise (N on−HI ) images are combined, the final step in creating the noise image involves identifying and replacing any extreme outliers (hot pixels with values ≥ |5σ|).This step is necessary to eliminate extreme sharp discontinuities that can cause artificial effects on the SPS via Gibbs ringing.where P noise (k) is the SPS of the noise image, P beam (k) is the SPS of the ALFA's beam, and P 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 subimage from the full GALFA-HI noise image.To improve statistics for any give noise SPS, we parameterize P noise (k) as a function of the standard deviation value of the input noise sub-image, sd = σ(N total ) estimated for each 3-degree kernel, and write it as P noise (k) = C × 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 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 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 σ rx .Since each kernel covers a 3-degree area, σ(N total ) is used as a proxy for σ rx .In Figure 10 we plot the log of σ(N total ) as a function of the log of P 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 σ(N total ) and log(P 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 σ(N total ).As a result, we have a family of about 100 Q(k, sd) functions to represent log(P noise ). Figure 11 demonstrates the functional change of log(P noise ) as a function of σ(N total ).
Figure 12 illustrates a sample of the raw noise SPS (curves black data points) overplotted 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

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 obs (k) is shown with red data points in Figure 13.The symmetrical Gaussian beam shape, P beam (k)

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

Figure 2 .
Figure2.Top: The GALFA-HI LVC column density image, divided into ten equal sub-regions named as G0 to G9 where each region is 36 × 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 • and 16 ′ , which corresponds to spatial frequencies k = 10 −3 to ∼ 6 × 10 −2 arcmin −1 ).

Figure 5 .
Figure5.Top: The GALFA-HI rolling SPS slope, where the SPS slope is calculated for each 3 • kernel and rolled across RA and DEC by stepping every 8 ′ to cover the entire GALFA-HI region.The SPS slope was determined only for regions where the SPS power ≥ 9× noise power, and the fitting was conducted over the range of 1.5 • to 16 ′ .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 χ 2 map of the single power-law SPS fit.the entire latitude range.For |b| > 60 degrees, the slope becomes slightly more shallow.Our mean SPS slope is in excellent agreement with whatKalberla & Haud (2019) found by using the all-sky HI4PI survey.Using the entire data set they found a SPS slope of −2.943 ± 0.003 and −2.944 ± 0.005 when considering only |b| > 20 degrees (for the velocity range from −8 to 8 km s −1 ).Figure7hints 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 ∼ 0 • to b ∼ 80 • , on the southern side the SPS slope is relatively constant from b ∼ 0 • to −50 • , and then has a sharp decrease from ∼ −2.9 to ∼ −2.6.This asymme-

Figure 6 .
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| ≤ 10 • is shown in blue, whereas the off-plane histogram pertains to latitudes |b| ≥ 70 • and is shown in orange.The vertical dashed lines represent the mean value for the three samples.

Figure 8 .
Figure 8. Top panel: The HI column density integrated over ±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) × sin |b| distribution, averaged in 2-degree latitude bins, as a function of sin |b|.distribution is modified by the explosion of multiple supernovae that formed the Local Bubble.While lines of sight at |b| ≤ 50 • predominantly probe the disk HI , at least half of the path length at |b| > 50 • 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 ∼ 200 pc (Welsh et al. 2010)) and therefore has a different fraction of warm to cold HI relative to the HI at |b| < 50 • .As shown in the bottom panel of Figure 8, the SPS slope of the N (HI) ×

Figure 10 .
Figure 10.The logarithmic power for the noise SPS, P noise (k), as a function of the logarithm of σ(N 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.
is shown as grey dashed curves.The noise SPS, C × Q(k, sd) is displayed as yellow data points in the same figure.Finally, [P obs − C × Q(k, sd)]/P beam is shown as green data points.The range of k-values for fitting is restricted to 1.5-degree and 4 × FWHM for all kernels, as explained in the text.