© J. D. Soler et al. 2022

1 Introduction

The diffuse neutral atomic hydrogen (H I) is the matrix within which star-forming clouds form and reside and the medium that takes a large portion of the energy injected by the massive stars during their lifecycles (see, Dickey & Lockman 1990; Kalberla & Kerp 2009; Ballesteros-Paredes et al. 2020, for a review). Therefore, the investigation of the H I distribution and dynamics is crucial to understanding the cycle of energy and matter in the interstellar medium (ISM, Ferrière 2001; Klessen & Glover 2016). In this paper, we present a study of the H I emission toward the Galactic plane, mainly focusing on the orientation of its filamentary structure and the Galactic dynamics that it reveals.

The H I emission at the 21 cm wavelength has been instrumental for the study of the Milky Way. This radio line is not limited by the dust absorption that renders optical methods ineffective at distances greater than a few kiloparsecs (van de Hulst 1953). Face-on maps of the outer galaxy have been constructed by using the Doppler shift of the emission line in combination with the Galactic rotation curve (van de Hulst et al. 1954; Westerhout 1957; Oort et al. 1958). The concentration of the H I gas in spiral arm features has been evident since the first face-on projections, although the detailed reconstruction of the density structures is complicated by the departures from circular motion, the complexity of the velocity field, and the overlapping of multiple components (Kerr 1969).

Much of H I is observed to be either a warm neutral medium (WNM) with T ≈ 10⁴ K or a cold neutral medium (CNM) with T ≈ 10² K (Heiles & Troland 2003). In the classical two-phase picture, both phases exist in pressure equilibrium with stable phases for temperatures of T ≈ 300 K (CNM) and T ≈ 5000 K (WNM) embedded in a hot ionized phase (McKee & Ostriker 1977; Wolfire et al. 2003). A portion of the CNM is sampled by the absorption of background H I emission by cold foreground H I, which is generically known as H I self-absorption (HISA; Heeschen 1955; Gibson et al. 2000; Wang et al. 2020b). However, most of our understanding of the multiphase structure of H I comes from the study of absorption toward continuum sources (see, for example, Strasser et al. 2007; Stanimirović et al. 2014; Murray et al. 2018). Comparisons between 21-cm H I emission and absorption measurements indicate that, in the vicinity of the Sun, the WNM has roughly the same column density as the CNM (Falgarone & Lequeux 1973; Liszt 1983). Observations of the absorption against strong continuum sources indicate that the CNM to WNM mixture may be constant throughout the Galaxy (Dickey et al. 2009, Dickey et al. 2021).

The first all-sky H I observations, provided by the 36′-resolution Leiden/Argentine/Bonn (LAB) survey (Kalberla et al. 2005), displayed the Milky Way as a nonaxisymmetric spiral system of four arms that are traced out to at least 25 kpc from the Galactic center (Levine et al. 2006a; Nakanishi & Sofue 2016; Koo et al. 2017). These four arms are designated Sagittarius-Carina (Sgr-Car), Scutum-Centaurus (Sct-Cen), Perseus, and Outer (or Norma-Cygnus). Their existence, however, is in tension with the star counts toward the tangent points, which favor the existence of only two spiral arms (Drimmel 2000; Benjamin 2009; Churchwell et al. 2009).

Further studies of the H I emission have also identified additional spiral arm structures. Using the LAB survey observations, Dame & Thaddeus (2011) identified a presumed spiral arm structure lying beyond the Outer Arm in the first Galactic quadrant (QI) at approximately 15 kpc from the Galactic center. This arm appears to be the continuation of the Scutum-Centaurus Arm in the outer Galaxy, as a symmetric counterpart of the nearby Perseus Arm, for which it is designated as the Outer Scutum-Centaurus (OSC) arm. Using the 2′-resolution H I observations in the interferometric Southern Galactic Plane Survey (SGPS, McClure-Griffiths et al. 2005), McClure-Griffiths et al. (2004) identified a possible distant spiral arm in the fourth quadrant (QIV) of the Milky Way. This distinct emission feature can be traced for over 70° as the most extreme positive velocity feature in the longitude-velocity diagram. Assuming circular motions around the Galactic center, this structure is located at Galactocentric radii (R_gal) between 18 and 24 kpc.

The H I emission also reveals the bending of the Galactic plane noticeable at R_gal > 9 kpc and up to 40 kpc (Henderson et al. 1982; Levine et al. 2006b; Kalberla et al. 2007). This “Galactic warp” is also observed in the dust emission, in the distribution of H II regions, in the CO emission, and the kinematics of stars beyond the solar neighborhood (see, for example, Drimmel & Spergel 2001; Sun et al. 2020; Cheng et al. 2020; Paladini et al. 2004). Furthermore, the H I emission also reveals that the average disk thickness shows a pronounced flaring that can be approximated by an exponential function in the range 5 < R_gal < 35 kpc (Kalberla & Dedes 2008).

In addition to the spiral-arm structure, the warp, and the flare, the distribution of the H I emission also carries the imprint of the energy input by main sequence stars, supernovae (SNe), and other energetic process in the ISM. The first photographic representations of the H I emission toward the Galactic plane revealed a series of shells, arcs, and filaments (Weaver & Williams 1973; Heiles 1979). The H I shells are thought to be formed by the combined effects of stellar winds and SNe sweeping up gas and dust (Cox & Smith 1974; Castor et al. 1975; McKee & Ostriker 1977). The vertical H I filaments, or “worms” seemingly “crawling away from the Galactic plane”, have been interpreted as the walls of bubbles which have broken through the thin gaseous disk (Heiles 1984; Koo et al. 1992).

Dedicated studies of the filamentary structures in the H I emission, based on single-dish observations at 4′ and 30′ resolutions, have shown that they have a preferential orientation parallel to the interstellar magnetic field in the local ISM (Clark et al. 2014; Kalberla et al. 2016). This alignment may be the product of magnetically-induced velocity anisotropies, the collapse of material along field lines, shocks, or anisotropic density distributions (see, for example, Lazarian & Pogosyan 2000; Heitsch et al. 2001; Chen & Ostriker 2015; Inoue & Inutsuka 2016; Soler & Hennebelle 2017; Mocz & Burkhart 2018). More recently, Soler et al. (2020) presented the analysis of the 40″-resolution H I interferometric observations in The Hl/OH/recombination-line (THOR) survey (Beuther et al. 2016; Wang et al. 2020a), revealing that the majority of the filamentary structures in the H I emission are aligned with the Galactic plane, but that there are portions of the first Galactic quadrant dominated by vertical H I structures that are most likely resulting from the combined effect of SN feedback and the interstellar magnetic fields.

This paper presents the extension of the Soler et al. (2020) analysis to the whole Galactic disk, hence extending in Galactic longitude and Galactic latitude. The filament orientation is a characteristic of the H I emission that is less explicitly dependent on the angular resolution than other properties, such as the filament width and the length. If the filaments are understood as features in the intensity field and not as physical objects, this analysis characterizes the anisotropy in the emission distribution. This characteristic is particularly relevant given the natural axis of symmetry in the Galactic plane. The gravitational potential of the Galaxy and the stretching due to the circular motions and the spiral arms could naturally produce horizontal structures. Departures from this natural anisotropy may indicate energetic processes that lift the gas in the disk or infall. Thus, rather than defining the filaments as objects, this study aims at studying the general anisotropy of the atomic gas distribution across the Galactic disk.

To cover the entire Galactic plane, we used the state-of-the-art H I emission observations that we introduce in Sect. 2. We describe the methods for identifying the filamentary structures and quantifying their preferential orientation in Sect. 3. We present the main results of our analysis in Sect. 4. We provide a discussion of our results in Sect. 5 and present our general conclusions in Sect. 6. We reserve additional analyses for a set of appendices as follows. Appendix A presents details on the Hessian matrix filament identification method, detailing the use of circular statistics and the selection of the analysis parameters. We present the method for de-projecting the radial velocities into Galactocentric radii and the calculation of the disk’s warp and flaring in Appendix B. Appendix C presents the comparison between the filamentary structures found in the HI4PI and GALFA-H I surveys. Finally, Appendix D presents the results of the filament analysis toward two high-velocity cloud complexes.

2 Data

We used the publicly-available observations of the emission by atomic hydrogen at 21cm wavelength in the H I 4π (HI4PI) survey (HI4PI Collaboration 2016). The HI4PI survey is based on data from the Effelsberg-Bonn H I Survey (EBHIS, Kerp et al. 2011) and the Galactic All-Sky Survey (GASS, McClure-Griffiths et al. 2009; Kalberla et al. 2010). This data comprises 21-cm neutral atomic hydrogen data over the whole sky in a radial velocity range −600 <υ_LSR < 600 km s⁻¹ for declination δ > 0° and −470 < υ_LSR <470 km s⁻¹ for δ < 0°, as observed with the 100-m radio telescope at Effelsberg and the 64-m radio telescope at Parkes. These observations have full spatial sampling, thus overcoming the major issue of the LAB survey, and provide a final data product with an angular resolution of 16.′2 and sensitivity of 43 mK per 1.29-km s⁻¹ velocity channel.

We used the data distributed in FITS-format binary tables containing lists of spectra sampled on a Hierarchical Equal Area isoLatitude Pixelization (HEALPix¹) grid with N_side = 1024, which corresponds to a pixel angular size of 3′.44 (Górski et al. 2005; Zonca et al 2019) We arranged these spectra in all-sky HEALPix maps corresponding to each of the 1.29-km s⁻¹ velocity channels. Using the tools in the Python healpy package, we produced a Cartesian projection for each one of these velocity channels to obtain maps covering the range −180° < l < 180° and −10° < b < 10° with a pixel size δl = δb = 3′.44. We combined these velocity-channel maps into a single position-positionvelocity (PPV) cube, which we segmented for the analysis of particular l and b ranges using the spectral-cube² package.

We also used the publicly-available observations from the Galactic Arecibo L-Band Feed Array H I survey (GALFA-H I) data release 2 (DR2, Peek et al. 2018). This survey covers the H I emission from −650 to 650 km s⁻¹, with 0.184 km s⁻¹ channel spacing, 4′ angular resolution, and 150 mK root mean square (rms) noise per 1 km s⁻¹ velocity channel. The GALFA-H I DR2 observations cover the entirety of the sky available from the William E. Gordon 305 m antenna at Arecibo, from δ = −1°17′ to +37°57′ across the whole right ascension (α) range. We chose these data to extend our analysis to two portions of the Galactic plane at higher angular and spectral resolution. In the range |b| ≤ 10°, GALFA-H I DR2 provides observations in the ranges 36°.5 < l < 70°.5 and 175°.3 < l < 209°.4.

We obtained the GALFA-H I DR2 distributed in FITS-format cubes in the “narrow” velocity range (υ ≤ 188 km s⁻¹). through the GALFA collaboration website³. We used the Python spectral-cube and reproject⁴ packages to arrange these cubes into two mosaics covering the aforementioned l and b ranges with a pixel size δl = δb = 1′.0 and spectral resolution of 0 184 km s⁻¹

3 Methods

3.1 Analysis of filamentary structure

We selected the filamentary structures using the Hessian matrix technique as implemented in Soler et al. (2021) and coded in the hessiananalysis routine of the magnetar software package⁵ (Soler 2020) This method is less computationally demanding than other filament identification methods, such as the rolling Hough transform (RHT, Clark et al. 2014) or FILFINDER (Koch & Rosolowsky 2015), and yields similar results for the filament orientation, as detailed in Appendix C of Soler et al. (2020). The main parameters used in the analysis are summarized in Table 1 and were applied as follows.

The HI4PI data was divided in 3° × 20° × 1.29-km s⁻¹ velocity channel maps centered on Galactic latitude b = 0° and covering the Galactic longitude in the range 4° < l < 356° and the radial velocity range −200 < υ_LSR< 200 km s⁻¹. Throughout this paper we refer to these velocity channels as “tiles”. The rectangular shape of the tiles aims to maximize the b coverage and maintain the fine sampling in l, covering both the Galactic warp and flaring. This selection does not imply any loss of generality in our results, as we demonstrate for different tile sizes in Appendix A.

For each tile centered on the position l₀ and b₀ and υ_LSR = υ₀, I(l₀, b₀, υ₀), we estimated the derivatives with respect to the local coordinates (x, y) and built the Hessian matrix, $$\text{H}\left(x,y\right)\equiv \left[\begin{array}{cc}{H}_{xx}& H_{xy}\\ H_{yx}& H_{yy}\end{array}\right],$$(1)

where H_xx ≡ ∂²I/∂x², H_xy ≡ ∂²I/∂x∂y, H_yx ≡ ∂²I/∂y∂x, H_yy ≡ ∂²I/∂y², and x and y are related to the Galactic coordinates (l, b) as x ≡ l cos b and y ≡ b, so that the x-axis is parallel to the Galactic plane. We obtained the partial spatial derivatives by convolving I(l₀, b₀, υ₀) with the second derivatives of a two-dimensional Gaussian function. In practice, we used the gaussian_filter function in the open-source software package SciPy (Virtanen et al. 2020). We report the results obtained with a 18′ FWHM derivative kernel, but present the results of different kernel sizes in Appendix A.

The two eigenvalues (λ_±) of the Hessian matrix were found by solving the characteristic equation, $${\lambda }_{\pm }=\frac{\left(H_{xx}+H_{yy}\right)\pm \sqrt{{\left(H_{xx}-H_{yy}\right)}^2+4H_{xy}{H}_{yx}}}{2}.$$(2)

Both eigenvalues define the local curvature of I(l₀, b₀,υ₀). In particular, the minimum eigenvalue (λ₋) highlights filamentary structures or ridges, as detailed in Planck Collaboration Int. XXXII (2016). The eigenvector corresponding to λ₋ defines the orientation of intensity ridges with respect to the Galactic plane, which is characterized by the angle $$\theta =\frac{1}{2}\text{arctan}\left[\frac{H_{xy}+H_{yx}}{H_{xx}-H_{yy}}\right].$$(3)

We estimated θ for each of the pixels in I(l₀, b₀, υ₀). This angle, however, is only meaningful in regions of the map that are rated as filamentary according to selection criteria based on the values of λ₋ and on the noise properties of the data.

3.2 Data selection for statistical analysis

In each tile, we selected the filamentary structures based on the criterion ${\lambda }_{-}<{\lambda }_{-}^C<0$, where ${\lambda }_{-}^C$ is a curvature threshold that we determined from the data. Additionally, we chose regions where I(l₀, b₀, υ₀) > 7σ_I. We picked σ_I = 0.43 mK to match the brightness temperature noise level for a 1.29 km s⁻¹ velocity channel reported in HI4PI Collaboration (2016), which was estimated using emission line-free spectral channels to compute the standard deviation and iteratively clip outliers following the procedure described in Winkel et al. (2016). We empirically tested different signal-to-noise ratio (S/N) thresholds for a fixed ${\lambda }_{-}^C$ and found no significant differences in our results for I(l₀, b₀, υ₀) > 3σ_I, 5σ_I, and up to 9σ_I, as shown in Appendix A.2. Most of the differences introduced by the intensity threshold were found in a few tiles beyond the terminal velocities, which are dominated by fainter emission and are dismissed with the highest S/N cutoffs. We found that the distribution of V across tiles is more sensitive to ${\lambda }_{-}^C$ than to the intensity S/N, thus, we chose 5σ_I as a reference value to charaterize the curvature threshold.

Following the method introduced in Planck Collaboration Int. XXXII (2016), we estimated ${\lambda }_{-}^C$ in noise-dominated velocity channels. For each tile group centered on a particular Galactic longitude, we estimated λ₋ in five velocity channels with low S/N and determined the minimum value of λ₋ in each of them. We used the median of these five λ₋ values as the threshold value, ${\lambda }_{-}^C$. We employed the median to reduce the effect of outliers, but in general the values of λ₋ in the noise-dominated channels are similar and this selection does not imply any loss of generality. In the filament orientation analysis, we exclusively considered regions of each velocity channel map where ${\lambda }_{-}<{\lambda }_{-}^C$, which corresponds to the selection of filamentary structures with curvatures in I(l₀, b₀, υ₀) larger than those present in the noise-dominated channels. We report the results of applying a common ${\lambda }_{-}^C$ estimated from the values found in all of the tile groups, but also evaluated the selection of a different ${\lambda }_{-}^C$ for each tile group, which we present in Appendix A.2. The outcome of the two approaches is very similar and the selection of one instead of the other does not have a significant effect on the results of our analysis.

3.3 Circular statistics

Once the filamentary structures were selected, we used the angles derived from Eq. (3) to study their orientation with respect to the Galactic plane. For a systematic evaluation of the preferential orientation, we applied the projected Rayleigh statistic (V) (see, for example, Batschelet 1981), which is a test to determine whether the distribution of angles is not uniform and peaked at a particular angle. This test is equivalent to the modified Rayleigh test for uniformity proposed by Durand & Greenwood (1958) for the specific directions of interest θ = 0° and 90° (Jow et al. 2018), such that V > 0 or V < 0 correspond to preferential orientations parallel or perpendicular to the Galactic plane, respectively. It is defined as $$V=\frac{{\displaystyle \sum {}_{ij}{w}_{ij}\text{cos}\left(2{\theta }_{ij}\right)}}{\sqrt{{\displaystyle \sum {}_{ij}{w}_{ij}}/2}},$$(4)

where the indices i and j run over the pixel locations in the two spatial dimensions (l, b) for a given velocity channel and w_ij is the statistical weight of each angle θ_ij. In our application, we accounted for the spatial correlations introduced by the telescope beam by choosing w_ij = (δl/Δ)², where δl is the pixel size and Δ is the diameter of the derivative kernel that we selected to calculate the gradients. We present the results of the filament orientation trends in terms of the mean orientation angle, 〈θ〉, and the Rayleigh test for uniformity, Z, in Appendix A.1.

In this paper, we report the results obtained with an 18.′0 FWHM two-dimensional Gaussian derivative kernel. We chose this kernel size to be close to the 16.′2 FWHM resolution of the HI4PI survey, which sets the angular scale of the independent measurements with the Hessian matrix. The selection of a larger kernel reduces the number of independent gradients and thus, the significance of the orientation trends quantified by V, as illustrated in Appendix A.2. We chose a kernel slightly larger than the angular resolution of the observations motivated by an initial choice that mitigated the effect of noise in the analysis of the interferometric observations in the THOR survey (Soler et al. 2020). This early selection was carried forward in this work although its results are essentially identical to those found using 16.′2 and 17.′0 FWHM derivative kernels.

The values of V indicate a significant preferred orientation if there is sufficient clustering around the orientations θ = 0° and 90°. The null hypothesis implied in V is that the angle distribution is uniform. In the particular case of independent and uniformly distributed angles, and for a large number of samples, values of V ≈ 1.64 and 2.57 correspond to the rejection of the null hypothesis with a probability of 5% and 0.5%, respectively (Batschelet 1972). Thus, a value of V ≈ 2.87 is roughly equivalent to a 3σ confidence interval.

The statistical significance represented by V is determined numerically by taking random samples from a von Mises probability distribution, which is a circular analog of a Gaussian distribution (Mardia & Jupp 2009). It is defined as $$ℱ\left(\theta |\mu ,\kappa \right)=\frac{e^{\kappa \text{cos}\left(\theta \right)-\mu }}{2\pi I_0\left(\kappa \right)},$$(5)

where I₀(κ) is the modified Bessel function of order 0, µ is the mean of the distribution, and 1/κ² characterizes the dispersion. Following a procedure like the one described in Jow et al. (2018), the significance is determined by drawing N samples from a von Mises distribution and computing V for a given null hypothesis rejection probability, for example, 5% and 0.5%. In a similar way to the chi-square test probabilities, V and its corresponding null hypothesis rejection probability are reported in the classical circular statistics literature as tables of “critical values”, for example in Batschelet (1972), or computed in the circular statistics packages, such as circstats in astropy (Astropy Collaboration 2018).

4 Results

The general results of the H I filament orientation analysis are presented in Fig. 1. The longitude-velocity (lυ) diagrams show the values of the mean intensity, 〈I〉, and the filament orientation quantified by the projected Rayleigh statistic, V, for each of the 3° × 20° × 1.29 km s⁻¹ tiles around the Galactic plane. For the sake of reference, we have included lines that indicate the position of the main features that have been traditionally assigned to the Sgr-Car, Sct-Cen, Perseus, and Outer spiral arms, as derived from the trigonometric parallaxes and proper motions for masers associated with young high-mass stars measured in the Bar and Spiral Structure Legacy (BeSSeL) survey (Reid et al. 2016). We also included the OSC arm and the distant arm reported in McClure-Griffiths et al. (2004).

The 〈I〉 lυ-diagram, shown on the left-hand side panel of Fig. 1, corresponds to the average of the intensity over the area of each tile. It is presented as a reference for the interpretation of the V lυ diagram and not for the particular study of the spiral arm structure. However, the 〈I〉 lυ diagram clearly shows the tracks assigned to the Perseus and Outer arms in the first and second Galactic quadrants. Other features, such as the Sct-Cen and the Sgr-Car arms, are less conspicuous.

The right-hand side panel of Fig. 1 suggests that there is no apparent correlation between the locations of the spiral arms in the lυ-diagrams and the preferential filament orientations characterized by V. The most prominent features in the V lυ-diagram, shown on the right-hand side panel of Fig. 1, are two groups of filaments parallel to the Galactic plane at the minimum and maximum radial velocities toward the first two (QI and QII) and the last two (QIII and QIV) Galactic quadrants, respectively. These bands, characterized by V > 5, extend over a few tens of kilometers per second and appear on radial velocities beyond the 〈I〉 peaks and the locations of the OSC and the McClure-Griffiths et al. (2004) spiral arms. The high-V band in QIII and QIV appears to extend into QII in the range 0 <υ_LSR < 30 km s⁻¹ down to l ≈ 110°. Both high-V bands are distinguishable around l ≈180°, where they appear to be separated by a gap of around 20 km s⁻¹. The extension in l of these two bands depends on the b coverage of the analysis; that is, they appear truncated or are less prominent when considering ranges |b| ≤ 1.5, 3.0, or 5.0, as shown in Appendix A.

Around 10% of the H I emission over the whole sky is coming from intermediate- and high-velocity clouds (IVCs and HVCs, respectively), which are parcels of gas at velocities incompatible with a simple model of differential Galactic rotation (see Wakker & van Woerden 1997; Putman et al. 2012, for reviews). Most of these structures are located at |b| > 10°. The majority of the HVCs in the range |b| < 10° are too small to dominate the trends reported in this work, for example, the objects in the GALFA-H I Compact Cloud Catalog (Saul et al. 2012). However, there are two exceptions: (i) the complex H (Hulsbosch 1975), located around 110° < l < 150°, −15° < b < 15° and −220 < υ_LSR < −80 km s⁻¹ ; (ii) the anticenter (AC) complexes (Wakker & van Woerden 1991) around 155° < l < 190°, −55° < b < −5° and −340 < υ_LSR < −135 km s⁻¹. Figure 1 indicates that both of these regions are beyond the high-V bands and are composed of mainly horizontal H I filaments. We present the results of the filament analysis for these two HVC complexes in Appendix D.

Figure 2 indicates that there is no evident correlation between the mean intensity in the tiles, 〈I〉, and the orientation of the filamentary structure, characterized by V. Around 83% of the tiles show values V > 0 and roughly 23% present V > 2.87. Only 0.34% of the tiles display V < −2.87, which indicates that filamentary structures perpendicular to the Galactic plane are rare. The tiles with the highest 〈I〉 generally indicate a preference for filaments parallel to the Galactic plane, as shown by V > 0 in the range 〈I〉 ≳ 60 K. However, most of the tiles with a high significance in preferential orientation, |V| > 2.87, are located in the range 〈I〉 ≲ 60 K, thus suggesting that the overlapping of multiple components in the tiles with the highest 〈I〉 randomizes the orientation of the filamentary structures.

Figure 2 also shows the filament coverage, which we define as the number of pixels identified as filaments according to the criteria described in Sect. 3.2 over the total number of pixels in each tile. We found that the tiles with the highest 〈I〉 show a larger coverage of filamentary structures, which may result from the selection based on the intensity S/N. Still, both high and low V values appear in tiles with similar filament coverage, which indicates that the preferential orientations are not the product of the number of filaments identified in a tile. Tiles with low 〈I〉 and filament coverage are more likely to be dominated by a single intensity structure, which explains why they display the most significant values of V.

Figure 3 presents an example of the filamentary structure in the emission toward the high-V bands reported in Fig. 1. The emission maps show that the filament orientation trends are not dominated by a monolithic object but rather by a collection of filaments where the orientation parallel to the Galactic plane seems to be the most frequent. The most prominent horizontal structures in the emission maps appear over a background of fainter and seemingly randomly oriented filamentary structures. In general, the width of the filamentary structures matches the kernel size used to calculate the Hessian matrix. This fact makes the general definition of the length of the filaments problematic since the apparent overlap of the most extended structures may not be a product of spatial coherence but an artifact of the available angular resolution.

Fig. 1 Longitude-velocity (lυ) diagrams of the mean intensity (〈I〉, left) and the filament orientation quantified by the projected Rayleigh statistic (V, right). Each of the pixel elements in the diagrams corresponds to a 3° × 20° × 1.29 km s−1 velocity channel map centered on b = 0°, which we call tile throughout out this paper. These results correspond to the Hessian analysis performed using an 18′ FHWM derivative kernel. Values of V > 0 (red) or V < 0 (blue) indicate a preferential orientation of the filaments parallel, θ = 0°, or perpendicular, θ = 90°, to the Galactic plane. The 3σ statistical significance for these two orientations corresponds to V > 2.87 or V < −2.87, respectively. The overlaid curves correspond to the main spiral arms features presented in Reid et al. (2016), the outer Scutum-Centaurus arm (OSC, Dame & Thaddeus 2011), and the extended outer arm (M-G2004, McClure-Griffiths et al. 2004).

Fig. 2 Mean H I emission intensity (〈I〉) and projected Rayleigh statistic (V, Eq. (4)) for the 3° × 20° × 1.29 km s−1 tiles. The colors indicate the tile percentage covered by the filamentary structures selected following the criteria described in Sect. 3.2. The horizontal dotted lines show V = ±2.87, which indicate the 3σ significance in the preferential orientation parallel, θ = 0°, or perpendicular, θ = 90°, to the Galactic plane, respectively. The vertical dashed line indicates the intensity threshold.

Fig. 3 Examples of the H I emission and filamentary structures toward the l and υLSR ranges with predominantly horizontal filamentary structures. The top image in each panel corresponds to the H I emission intensity in three velocity-channel maps presented in red, green, and blue colors. The bottom image corresponds to the second eigenvalue of the Hessian matrix, λ−, which highlights the filamentary structures in each velocity-channel map. The white disk in the lower right corner corresponds to the size of the derivative kernel used to calculate the Hessian matrix.

4.1 Distribution of H I filamentary structures in the Galactic plane

We analyzed the distribution of H I filamentary structures in the Galactic plane by assuming circular motions around the Galactic center and a standard rotation curve, as detailed in Appendix B. First, we focused on the behavior of the filament orientations with heliocentric radius to study the combined effect of the angular resolution and distance. Second, we compared the filament orientations at different heliocentric kinematic distances. Finally, we studied the variations with Galactocentric radius and produced a face-on view of the H I filament orientation in the Galactic plane.

4.1.1 Filament orientations and heliocentric distance

Figure 4 shows the distribution of V and the far kinematic heliocentric distance, d_far, derived for each tile. We chose to report d_far to focus on the effect of distance on V, but similar conclusions can be drawn from the near kinematic distances, d_near, as shown in Appendix B. We found an increase in V with heliocentric distance that is particularly prominent toward QIII and QIV. However, we did not find that the tiles corresponding to the largest distances present notably large V, thus suggesting that the increase in V is a product of heliocentric distance alone.

The mean V values in 1-kpc d_far bins, shown by the black line in the bottom panel of Fig. 4, reveal a general preference for V > 0. However, these mean V values are below the 2.87 limit corresponding to a 3σ significance of a preferential orientation parallel to the Galactic plane. Most of the scatter in the distribution of V for a particular d_far seems to come from the azimuthal variations, as can be inferred from the comparison between the general trend and that toward each Galactic quadrant. Toward QII and QIII, we found a systematic increase in V up to d ≈ 15 kpc followed by a systematic decrease at larger distances. This may be the product of the decrease in the emission at the edge of the Galactic disk. Toward QI and QIV, we found a decrease in V up to d_far ≈ 10 kpc, or d_near ≈ 5 kpc, as illustrated in Appendix B.

To further establish if the prevalence of V > 0 at large distances is an effect introduced by the distance, we repeated the analysis using larger derivative kernels in Appendix A.2.3. Using derivative kernels with diameters two and three times larger, which is equivalent to locating the sources two or three times further away, results in a systematic reduction of V rather than a bias toward V > 0. We interpret that result as a strong indication that the relative orientation trends are not exclusively the product of the heliocentric distance. However, we also conducted a test using a derivative kernel fixed to a particular physical scale, which we describe in the following section.

Fig. 4 Kinematic heliocentric distance (d) and projected Rayleigh statistic (V, Eq. (4)) for the 3° × 20° × 1.29 km s−1 tiles. We report the largest (far) solution to the heliocentric distance equation, as detailed in Appendix B. From top to bottom, the panels correspond to each Galactic quadrant and to the whole Galactic longitude range. The colors indicate the Galactic longitude of the tiles. The horizontal dashed lines indicate V = ±2.87, which correspond to the 3σ significance in the preferential orientation around θ = 0° or θ = 90°, respectively. The solid line in the bottom panel corresponds to the mean V values in 1-kpc dfar bins.

Fig. 5 Mean intensity (〈I〉, top), projected Rayleigh statistic (V, middle), and derivative kernel size (Δ, bottom) for 3° × 20° × 1.29 km s−1 toward the indicated positions. The two lines in the middle panel indicate the results obtained with a derivative kernel with a fixed angular size and a kernel with a fixed physical size, as indicated on the legend. The kernel with a fixed physical size is exclusively used in the radial velocity range corresponding to kinematic distances where the angular size of the kernel is below 1°.

4.1.2 Filament orientations near and far

By fixing the size of the derivative kernel in the Hessian analysis, we are effectively selecting a physical size scale that changes with the distance. For example, at a distance of approximately 4.5 kpc the 18′ FWHM derivative kernel would select filamentary structures of minimum widths around 25 pc. Those structures would be blended into 100-pc-width structures sampled at 18 kpc with the same beam.

We evaluated the effect of distance by considering a derivative kernel with a fixed physical size. We ran this test toward the positions l ≈ 90° and l ≈ 270°, where the mapping between υ_LSR and d has the largest dynamic range. We limited the kernel angular size to Δ < 1° to avoid artifacts produced by the limited sampling of the 3° × 20° × 1.29 km s⁻¹ tiles. This limit excludes the radial velocity range where most of the emission lies but covers the range in which the tiles with the highest V are found. Figure 5 shows the comparison of the results obtained with the 18′ kernel and a kernel with varying angular size corresponding to Δ = 160 pc at different distances.

We found that the increasing kernel sizes produced lower values of V. This indicates that rather than producing a bias toward higher V, as would be the case if multiple structures with a random orientation were blended into a single horizontal filament, the smoothing results in a loss of the preferential orientation parallel to Galactic plane. This trend toward lower values of V was also found when using derivative kernels with larger angular sizes, as detailed in Appendix A. These results strongly suggest that the high-V bands in the lυ-diagrams in Fig. 1 cannot be explained by an effect exclusively related to the heliocentric distance.

Figure 3 presents the emission toward the range in l and υ_lsr corresponding to the high-V bands in Fig. 1. Rather than being concentrated, the emission appears as a series of horizontal filaments distributed in a broad range of Galactic latitudes. This distribution is the product of the Galactic flaring and indicates that the smoothing by the derivative kernel is not blending the emission into a single monolithic structure at large distances. Therefore, we conclude that the high values of V are not the result of an observational effect but a property of the H I structure.

Fig. 6 Face-on view of the mean H I emission intensity (〈I〉; left) and the filament orientation quantified by the projected Rayleigh statistic (V; right). Values of V > 2.87 or V < −2.87 indicate a statistical significance above 3σ in the preferential orientation around θ = 0° or θ = 90°, respectively. The yellow circular marker indicates the position of the Sun. The overlaid curves correspond to the spiral arm models identified in Reid et al. (2019), shown in black color, and their extrapolation to larger Galactocentric radii, shown in green.

4.1.3 Filament orientations and Galactocentric distance

We computed a de-projection of the 〈I〉 and V lυ-diagrams into a face-on view of the Milky Way, shown in Fig. 6. This reconstruction does not aim to produce a complete 3D density distribution but rather identify the Galactocentric radii and azimuth to which the features in the V lυ-diagrams belong. Following preceding works, such as Levine et al. (2006b) and Kalberla et al. (2007), we excluded from the analysis the Galactic longitude ranges l < 15°, l > 345° and 165° < l < 195°, as these regions have radial velocities that are too small with respect to their random velocities to establish reliable distances. We also excluded from this reconstruction points with R_gal < 10 kpc, to avoid the complications introduced by the near-far distance ambiguity and the artifacts produced by the emission around υ_LSR ≈ 0 km s⁻¹.

The face-on reconstruction of the mean intensity, shown on the left-hand-side panel of Fig. 6, displays some of the spiral arm features identified in the LAB survey around 90° ≲ ϕ ≲150° and 200° ≲ ϕ ≲ 250° around 10 < R_gal < 15 kpc (Levine et al. 2006a; Nakanishi & Sofue 2016; Koo et al. 2017). The identification of these arms is not the main focus of our analysis, but they appear highlighted by the emission concentration. The arms in QI and QIV are less conspicuous and not clearly identifiable without additional processing, such as, the median filtering applied in Levine et al. (2006a) or the spectral peak identification used in Koo et al. (2017).

The face-on reconstruction of the filament orientation trends, presented on the right-hand-side panel of Fig. 6, shows a clustering of high V at R_gal ≳ 15 kpc on both sides of the Galaxy. This collection of high-V structures, which corresponds to the bands in the lυ-diagrams in Fig. 1, does not appear simply as a ring but seems to have a pitch angle. This pitch angle impression is produced by the displacement of the peak V values from R_gal ≈ 18 kpc at ϕ ≈ 210° to R_gal ≈ 23 kpc at ϕ ≈ 270°. The pitch angle is less noticeable in the other side of the Galaxy, where there is a prominent exception to the general V > 0 trend around R_gal ≈ 17 kpc and l ≈ 95°. This feature, which corresponds to the negative V tiles at d ≳ 25 kpc in Fig. 7, was also identified in the Galactic warp and flare studies reported in Levine et al. (2006a).

Figure 6 also shows the asymmetry in the H I between the first two and the last two quadrants. Toward QI and QII there is a drop in the amount of emission at R_gal ≈ 25 kpc, which is noticeable in the lack of filaments and relative orientation trends beyond that radius in the range 45 ≲ ϕ ≲ 130°. However, there is enough emission to suggest a drop in V at R_gal ≈ 25 kpc with respect to its peak value, which we explored further by computing radial profiles.

Figure 7 shows the distribution of V with respect to the R_gal derived from the radial velocity of each tile and the average radial profile computed in 1-kpc-width bins in the range 2 ≤ R_gal ≤ 30 kpc. We included the observations at R_gal < 10 kpc since the R_gal estimates are not affected by the heliocentric distance ambiguity. The general trend is a progressive increase of V with growing R_gal until it reaches its maximum values around R_gal ≳ 20 kpc. This tendency is equivalent to the continuous transition from H I filaments slightly perpendicular or having no preferred orientation to mostly parallel to the Galactic plane with increasing distance from the Galactic center. The significance of the relative orientation in the tiles dominated by mostly perpendicular filaments increases when considering broader tiles in Galactic longitude, as illustrated in Appendix A.2.1. However, the broader tile size reduces the sampling in Galactic longitude and the resolution in the face-on reconstructions presented in Fig. 6.

Figure 7 also shows the central locations of the H I shells in the catalogs produced with the LAB survey (Ehlerová & Palouš 2013) and the SGPS observations (McClure-Griffiths et al. 2002) in the range |b| < 10°. Most of the H I shells are found within the range R_gal < 10 kpc, and their number tends to decrease with increasing Galactocentric radius. For R_gal > 15 kpc, where the values of V reveal a very significant prevalence of H I filaments parallel to the Galactic plane, there are very few identified H I shells. The likely relationship between the values of V and the density of H I shells with R_i suggest a connection between the H I filament orientation and SN feedback, which we discuss further in Sect. 5.

In the vicinity of the solar circle, R_gal = 8.15 kpc (Reid et al. 2019), we found both negative and positive values of V that show a significant departure from the general trend. Some of these fluctuations are due to H I emission in the solar neighborhood, whose proximity makes the local atomic gas dynamics more prominent than the large scale trend. The peak in negative V around R_gal ≈ 10 kpc is most likely associated with the large H I shells identified in the range 250° ≲ l ≲ 360° in radial velocities that roughly correspond to that position in the Galactic plane (McClure-Griffiths et al. 2002). Additional fluctuations in V are also introduced by noncircular motions, which are particularly dominant around the anticenter, thus explaining why the largest V are observed around l ≈ 160°. It is likely that the emission around that Galactic longitude is produced further away and that the distance estimates based on the line kinematics are inaccurate (see, for example, Burton 1971; Peek et al. 2022).

At R_gal ≳ 16 kpc the dispersion in V increases and displays relatively large values of V and an excursion into V < 0 for some of the tiles. The points corresponding to V < 0 around R_gal ≈ 20 kpc can be traced into an individual feature around ϕ ≈ 100° that is also found as an outlier in the Galactic warp (Levine et al. 2006b). Most of the tiles with high values of V at R_gal > 20 kpc are found in QIV, where the high V band in the lυ-diagram is the broadest, but for most of the lines of sight V decreases beyond R_gal > 20 kpc.

Fig. 7 Galactocentric radius (Rgal) and projected Rayleigh statistic (V, Eq. (4)) for the 3°× 20° × 1.29 km s−1 tiles. From top to bottom, the panels correspond to each Galactic quadrant and to the whole Galactic longitude range. The colors indicate the Galactic longitude of the tiles. The horizontal dashed lines indicate V = ±2.87, which indicate the 3cc significance in the preferential orientation around θ = 0° or θ = 90°, respectively. The vertical dotted line indicates the radius of the solar orbit. The solid line in the bottom panel corresponds to the mean V values in 1-kpc Rgal bins. The vertical solid lines mark the central locations of the H I shells identified in the LAB survey (Ehlerová & Palouš 2013) and the SGPS (McClure-Griffiths et al. 2002) in the range |b| < 10°, shown in gray and green colors, respectively.

4.2 Scale height of the H I filaments

In addition to the analysis of the filament orientation distribution with R_gal, we considered two other H I emission properties that also change with the distance from the Galactic center: the midplane height and the dispersion around the midplane. We computed the Galactic height z for each of the tiles and obtained its first and second intensity-weighted moments, 〈z〉 and σ_z following the procedure detailed in Appendix B. We estimated 〈z〉 and σ_z both for all of the H I emission in each tile and exclusively for the emission in the structures classified as filaments; hence we obtained two separate estimates of the warp and the flare.

Figure 8 shows the two face-on reconstructions of 〈z〉. The midplane height map obtained using all of the H I emission reproduces some of the main features obtained in the previous analysis of lower resolution data (Levine et al. 2006b; Kalberla et al. 2007). First, the warp of the Milky Way, which corresponds to the variation in midplane height from 〈z〉 ≈ 4 kpc toward QI and QII to 〈z〉 ≈ −1 kpc at R_gal ≲ 20 kpc toward QIII and QIV. Second, the finding of 〈z〉 > 0 and 〈z〉 < 0 for adjacent tiles toward QIII and QIV at R_gal ≳ 20 kpc, which has been identified as Galactic disk “scalloping” (Kulkarni et al. 1982; Levine et al. 2006b). Finally, a sharp dip around ϕ ≈ 90° and R_gal ≳ 20 kpc.

The 〈z〉 found with the H I filaments shows features that are very similar to the midplane heights obtained using the bulk of the H I emission. However, the contrast in the warp, the scalloping, and the sharp dip is much higher. This result is most likely a consequence of filtering the diffuse emission implicit in the filament identification.

The 〈z〉 radial profiles, presented in Fig. 9, are very similar for all the H I emission and for the H I filaments in the range R_gal ≲ 20 kpc. For larger R_gal, the H I filaments reveal that the Galactic warp extends up to approximately 4 kpc for tiles in the first Galactic quadrant. The H I filaments also show higher values of 〈z〉 at R_gal ≳ 20 kpc toward QIII and QIV, which strongly suggest that specific filamentary structures produce the apparent disk scalloping.

Figure 10 shows the two face-on reconstructions of σ_z. The dispersion of the H I emission around the midplane shows the flaring of the Galactic disk at R_gal > 15 kpc, both in all the emission and just in the filaments. However, the radial profiles shown in Fig. 11 indicate that the average trend is for the filamentary structures to have a lower dispersion than the bulk of the H I emission by about a factor of two at R_gal > 15 kpc. This general trend is followed in most of the l range, particularly around 250° < l < 350°, but there are clear exceptions where the σ_z values traced by the filaments drop abruptly. This can be attributed to the lower emission in radial velocities corresponding to larger distances toward QI, which makes difficult the filament identification at R_gal ≳ 25 kpc in that azimuth range.

Fig. 8 Face-on view of the H I midplane height (〈z〉) reconstructed from all the emission in the range |b| < 10° (left) and the emission in filamentary structures (right).

Fig. 9 Midplane height (〈z〉) and Galactocentric radius (Rgal) for the 3° × 20° × 1.29 km s−1 tiles. The left-hand side panel corresponds to all the H I emission in the |b| < 10° range. The right-hand side panel corresponds to the emission in the structures classified as filaments. The colors indicate the Galactic longitude, l. The vertical dotted line indicates the radius of the solar orbit around the Galactic center.

4.3 Global properties of the H I filaments

The results of the Hessian analysis provide a description of the H I emission morphology across velocity channels, which does not imply that the highlighted structures are filaments or sheets in three-dimensional space. Our finding is that the high-aspect-ratio features in the H I emission have a preferential orientation with respect to the Galactic plane. Without any loss of generality, we can also estimate what fraction of the H I column density is in these high-aspect-ratio features.

We calculated the H I column density in each tile using the assumption of optically thin emission (see, for example, Kulkarni & Heiles 1987), such that, $$N_{\text{H}}=1.83\times {10}^{18}\left(\frac{T_{\text{b}}}{\text{K}}\right)\left(\frac{\text{Δ}\upsilon }{\text{km}{\text{s}}^{-1}}\right){\text{cm}}^{-2}$$(6)

where T_b is the brightness temperature and Δυ is the width of the velocity channel. This assumption reduces our calculation to the ratio between the sum of T_b in the filamentary structures and the total T_b in each tile. We note, however, that the optically-thin assumption underestimates the amount of H I in the Galactic plane (Bihr et al. 2015; Wang et al. 2020a). The optical depth mainly affects the inferred CNM column density, since typically τ ≪ 1 for the WNM (see, for example, Heiles & Troland 2003). Insofar as the ratio of CNM to WNM in the filaments may be different than in the total gas, our results for $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ are only indicative of the global trend.

Figure 12 shows the distribution of $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ as a function of distance from the Galactic center. The $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ radial profile is tightly distributed around 80% for R_gal ≲ 18 kpc. For larger R_gal, the mean value of $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ decreases sharply, although we found a large dispersion around this general trend. The general decrease in $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ for R_gal ≳ 18 kpc is also found when fixing the I S/N threshold using larger derivative kernels, as illustrated in Appendix A.2.3. Thus, it is most likely that this trend represents a transition toward less structured H I emission rather than a drop due to the decrease in the emission.

The azimuthal variations in the $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ radial profile are most likely a product of the asymmetry in the Galactic disk. In the range l ≲ 90° we found a sharp decrease in $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ that corresponds to the drop in emission toward that particular direction. In contrast, we obtained $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ around 70% up to R_gal ≈ 30 kpc for some tiles around l ≈ 150°.

The sharp drop in the mean values of $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ with Galactocentric radius roughly coincides with the transition from mostly horizontal to randomly oriented H I filaments at R_gal ≈ 18 kpc, illustrated in Fig. 7. Likely, the large dispersion in V found at R_gal ≳ 18 kpc is the product of a reduction in the number of filamentary structures, as sampled by $N_{\text{H}}^{\text{fil}}$, and the azimuthal variations in the emission.

Fig. 10 Face-on view of the dispersion around the midplane height (σz) reconstructed from all the emission in the range |b| < 10° (left) and the emission in filamentary structures (right).

Fig. 11 Same as Fig. 9, but for the dispersion around the midplane height (σz).

Fig. 12 Radial profiles of the ratio between the column density in the filamentary structures $N_{\text{H}}^{\text{fil}}$ and the total column density (NH). From top to bottom, the panels correspond to each Galactic quadrant and to the whole Galactic longitude range. The colors indicate the Galactic longitude of the tiles. The vertical dotted line indicates the radius of the solar orbit. The solid line in the bottom panel indicates the mean values of $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ in 1-kpc Rgal bins.

4.4 Comparison with GALFA-H I

We have considered the orientation of the H I filaments found with the Hessian method fixing a particular angular scale through the selection of the derivative kernel size and in velocity channels set by the spectral resolution of the HI4PI survey. There is no survey covering the same b range across the whole Galactic plane. Still, the 4′ and 0.184 km s⁻¹ resolution GALFA-H I provides an insight into the effects of the angular and velocity resolutions, which we discuss in what follows.

4.4.1 Effect of the angular resolution

We performed the Hessian analysis in the two portions of the first and the third Galactic quadrants covered in the 4′-resolution GALFA-H I data. To independently study the effects of the difference in angular and spectral resolution, we projected the GALFA-H I observations into the same spectral axis of the HI4PI data. We performed the Hessian analysis using an 5′ FWHM derivative kernel and applied the selection criteria presented in Table 2. Further details on the implementation of the Hessian matrix method in the GALFA-H I data and the results obtained with a larger kernel are presented in Appendix C.

The results of the GALFA-H I data analysis, presented in Fig. 13, indicate at first glance that the global trends in relative orientation are roughly the same in the higher resolution data. The most substantial difference between the analysis results in the two datasets is the increase in the significance in the GALFA-H I, which is further illustrated in Fig. 15. The global trend is for V to increase by approximately a factor of two, which is most likely the combination of two effects. First, the relative orientation trends observed in HI4PI are generally the same at the scales sampled by GALFA-H I. If the primarily horizontal filaments responsible for the highest values of V in HI4PI resolved into randomly oriented structures, we would find that the high-V tiles would produce low V in GALFA-H I, which is not what we observed in most of the tiles. Second, the increase in V is related to the higher angular resolution, because a smaller beam produces more independent gradient measurements for the same tile area. This does not imply that there is a higher number of filaments in the higher resolution but that the statistical trend is still found when including more independent orientation samples. We present examples of the filaments identified tiles in the two data sets in Fig. 14.

The higher values of V in the higher resolution data indicate that the distribution of filament orientations still peaks at 0° and the prominence of that peak is more elevated. These results do not exclude cases in which the higher resolution observations reveal no preferential orientation, as in some tiles with V > 2.87 in HI4PI and V ≈ 0 in GALFA-H I. Still, the significance in V is generally higher in the GALFA-H I data by roughly a factor of two.

The increase in V with the higher resolution observations produces a more significant portion of the lυ-diagram, as illustrated in Fig. 13. This general increase in V does not imply a departure from our observation in the HI4PI data that the filaments tend to be more parallel to the Galactic plane with increasing R_gal. It means that the change in filament orientation is smoother and more significant at higher angular resolution. This trend is confirmed by the analysis of the THOR survey (Soler et al. 2020), where the increase in angular resolution by a factor of 24 with respect to HI4PI reveals a significant preference for H I filaments parallel to the Galactic plane at υ_lsr < 0 km s⁻¹ in the range 17° < l < 67° and |b| <1°25.

Fig. 13 Longitude-velocity (lυ) diagrams of the filament orientation quantified by the projected Rayleigh statistic (V, right) in the two portions of the Galactic plane covered in the range |b| ≤ 10° by both HI4PI and GALFA-H I. TheGALFA-H I observations were analyzed at 1.29 and 0.184 km s−1 resolution using a 5′ FHWM derivative kernel. The curves correspond to the same spiral arm features introduced in Fig. 1.

Fig. 14 Examples of the H I emission and filamentary structures identified in the HI4PI and the GALFA-H I data, shown in the top and bottom panels respectively. The disk in the lower right corner of each panel indicates the size of the derivative kernels used for the calculation of the Hessian matrix, 18′ and 5′ for the HI4PI and the GALFA-H I observations, respectively.

Fig. 15 Projected Rayleigh statistic (V) calculated for the filamentary structure in the HI4PI and GALFA-H I observations. Each point corresponds to a 3° × 20° × 1.29 km s−1 tile. The gray dashed lines indicate V = 2.87 and V = −2.87, which correspond to a 3σ significance in the preferential orientation around θ = 0° or θ = 90°, respectively.

4.4.2 Effect of radial velocity resolution

We also performed the Hessian analysis in the GALFA-H I observations at their native 0.184 km s⁻¹ resolution using an 5′ FWHM derivative kernel. The results, also presented in Fig. 15, show no global change in the results of our analysis when considering the higher radial velocity resolution in the GALFA-H I data. This is not unexpected given that emission from the narrowest H I components, corresponding to the CNM, is at least a few kilometers per second wide and the emission structures are likely to be coherent across several narrow velocity channels.

We illustrate the effect of the angular and radial velocity resolution for one example tile in Fig. 16. At the same 1.29km s⁻¹ resolution and S/N threshold, the results of the HI4PI analysis at 18′ FWHM and GALFA-H I at 5′ FWHM indicate a slight change in the orientation trends. This change is not unexpected, as the orientation at the large scale does not uniquely determine the orientation of the structures at the small scale. Most of the differences are found in the velocity channels where the bulk of the emission is located, υ_lsr ≈ 0 km s⁻¹. At the largest υ_lsr, which also correspond to the largest d and R_gal, the increase in angular resolution does not introduce any significant changes in the relative orientation trends, thus supporting that the large V values found there are not a product of the distance to the emission sources. For illustration, we also present the mean orientation angles 〈θ〉, which further confirms the dominant trend of filaments aligned with the Galactic plane at both angular resolutions.

The results of the GALFA-H I analysis at 0.184km s⁻¹ resolution are very similar to those at 1.29km s⁻¹ resolution. The orientation trends are coherent across a few kilometers per second in both cases. The most prominent exception is the loss of signal around at υ_LSR ≈ 80 km s⁻¹ and υ_LSR > 150 km s⁻¹, which is the result of the S/N decrease in the narrow velocity channels. In general, we do not find significant differences in the analysis of the narrow velocity channels that would critically change the filament orientation trends obtained at 1.29km s⁻¹ resolution.

5 Discussion

5.1 Atomic filament orientations as a tracer of Galactic dynamics

The trends reported in Fig. 1 indicate that most of the filaments in the H I emission are parallel to the Galactic plane. This preferential orientation appears naturally from the combined effect of the gravitational potential and Galactic rotation, as confirmed by numerical simulations of Galactic disks (see, for example, Smith et al. 2020; Tress et al. 2020). The gravitational potential constrains the gas in the disk. The Galactic rotation stretches the parcels of gas by shearing them apart. However, these effects cannot explain by themselves the slight prevalence of vertical H I filaments toward the inner Galaxy and the clear predominance of horizontal ones toward the outer Galaxy.

One observational effect that may influence the observed H I filament orientation change with Galactocentric radius is the pileup of emission from a broad range of distances into the same channel in velocity space. The portion of line of sight covered by a fixed-width velocity channel increases with Galactocentric radius. For example, at l = 90° and under the assumption of circular motions, a 1.29-km s⁻¹ velocity channel corresponds to ΔR_gal ≈ 200 pc and ΔR_gal ≈ 500 pc for υ_LSR = 100 and 150 km s⁻¹, respectively. The velocity pile up can be a plausible explanation for the decrease in V at R_gal ≳ 20 kpc reported in Fig. 7, by mixing portions of the line of sight where H I structures are mostly parallel to the Galactic plane and randomly oriented into the same velocity channel. However, it does not offer a straight forward explanation for the increase of V with R_gal for 10 ≲ R_gal ≲ 20 kpc or the high-V bands in Fig. 6.

The analysis of the GALFA-H I observations indicates that the regions dominated by filaments parallel to the Galactic plane, V > 2.87, do not resolve into structures with random orientations. Instead, the higher angular and velocity resolution observations revealed that this preferential orientation is even more significant, as illustrated in Fig. 13. The increase in angular resolution does not result in a similar increase in the orientation significance for tiles slightly dominated by filaments perpendicular to the Galactic plane, V < 0. We interpret this as indication that vertical H I filaments are less frequent than horizontal ones.

The experiments with derivative kernels fixed to a physical length, such as the one presented in Fig. 5, indicate that the progressive transition from V < 0 to V > 0 with increasing R_gal is not due to a mixture of scales introduced by the fixed angular size of the derivative kernel. Rather than producing larger V, the introduction of derivative kernels with larger angular sizes results in the reduction of V, as further illustrated in Appendix A.2.3. Hence, we consider that the progressive transition from V < 0 to V > 0 with increasing R_gal reported in Fig. 7 has an origin in the physical conditions of the Galactic disk.

5.1.1 Vertical H I filaments

Figure 17 shows an example of the emission and the filaments toward one of the inner Galaxy regions with prevalence of vertical structures, or V < 0. The majority of the structures that would be identified as HISA are constrained to relatively low Galactic longitudes, where the background emission is brightest, and are unlikely to be the dominant object producing the global filament orientation. This is also true for other regions in the inner Galaxy, as expected from the distribution and angular size in HISA catalogs (see, for example, Gibson et al. 2005). However, it does not dismiss the importance of a dedicated study of the HISA structures orientation.

Most of the structures contributing to V < 0 are worms and arcs above and below the bright filamentary structure that corresponds to the midplane. These structures extend into high Galactic latitudes, some even beyond the studied Galactic latitude range, |b| < 10°. However, the higher-resolution but lower b-coverage H I observations in the THOR survey indicate that predominantly vertical filamentary structures are also found in the |b| < 1°0 range (Soler et al. 2020).

One physical explanation for the observed vertical H I filaments is the Galactic fountain, in which hot gas is produced by the cumulative effect of supernovae in the disk; thermal expansion lifts it above the Galactic plane, where it becomes unstable; cool clouds condense and fall back ballistically (Shapiro & Field 1976; Fraternali & Binney 2006). However, direct SN feedback events can also cause vertical structures. Strong clustered SN feedback causes large scale bubbles, which are consecutively elongated along the vertical direction due to bouyancy and Parker instability (Parker 1966). Cosmic rays (CRs) add an additional source of slowly cooling energy, which builds up a vertical pressure gradient that helps lifting the gas out of the midplane (Breitschwerdt et al. 1993; Girichidis et al. 2016; Ruszkowski et al. 2017).

Numerical simulations indicate that CR-lifted gas is primarily warm rather than hot, which will increase the amount of H I compared to hot only SN-heated gas (see, for example, Girichidis et al. 2018). CRs mainly move along the magnetic field lines, which further amplifies the effects of Parker loops and the related instabilities (Hanasz & Lesch 2003). Consequently, the gas along the buoyantly rising magnetic field lines also shows a preferentially vertical orientation. Overall, strong SN feedback in the presence of a vertical stratification of gas and a dynamically relevant magnetic field will favor the formation of vertical structures. Additionally, spiral shock waves may lift and stretch the atomic gas in the vertical direction (Suchkov & Shchekinov 1974; Kim et al. 2010) and Parker instability that can create buoyant loops of magnetic field at a kiloparsec scale (Rodrigues et al. 2016; Heintz et al. 2020).

The progressive change from mostly vertical H I filaments to mostly horizontal with increasing R_gal, reported in Fig. 7, indicates that the dominant process behind the H I filament orientation also has a radial dependence. Spiral shocks trace the position of the spiral arm structure, thus they should carry an azimuthal dependence and produce differences between arm and inter-arm regions that are not noticeable in our results. The axisymmetric variations in filament orientation strongly suggest the prevailing role of SN feedback.

Supernovae are very effective in lifting the gas off from the Galactic plane, as seen in numerous numerical experiments (see, for example, Joung & Mac Low 2006; Hennebelle & Iffrig 2014; Gatto et al. 2015; Kim & Ostriker 2015). The radial variations in the clustering and frequency of SNe with Galacto-centric radius are a plausible explanation for the lack of vertical filaments in the outer Galaxy. The number of SN remnants in the Milky Way decreases significantly for R_gal > 5 kpc (Case & Bhattacharya 1998; Green 2015). The distribution of stars in the Galaxy decreases exponentially with a characteristic radius R₀ that depends on the age of the population and is around 3.5 kpc for young stars (Drimmel & Spergel 2001; Frankel et al. 2018). Both of these observations indicate that for R_gal > 10 kpc, the number of SNe is significantly lower than in the inner Galaxy. At large R_gal not only the number of SNe is lower, but the area of the disk also is more extensive, so the SN rate per unit area falls off even faster with radius than the number of SNe. Consequently, due to the lack of events in the outer Galaxy, SNe’s energy and momentum input may not be the dominant process shaping the atomic gas.

The SNe do not only contribute to the formation of vertical filaments by lifting the H I gas, but also through the infall of the previously ejected gas back on the plane. Marasco & Fraternali (2011) studied the extraplanar H I gas of the Milky Way and found that it mainly consists of gas that is falling back to the MW after being ejected by SNe. The extraplanar H I is more conspicuous while falling because it cools while settling down. When it is lifted, it is hotter and so less prominent in the H I emission, as can be seen in the numerical simulations presented in Kim & Ostriker (2015) and Girichidis et al. (2016). It is difficult to determine the portion of the vertical filaments that correspond to outgoing or incoming H I gas. Still, the filament orientations provide a roadmap to study these effects in the Galactic plane and in Galactic-scale simulations that include the impact of stellar feedback.

Fig. 16 Mean intensity (〈I〉, top), projected Rayleigh statistic (V, middle), and mean orientation angle of the filaments (〈θ〉, bottom) for the tiles centered on l = 202.5 in the HI4PI and GALFA-H I. The left panel corresponds to the HI4PI observations analysis with a 18′ FWHM kernel for a 1.29-km s−1 velocity resolution. The center and right panels correspond to the GALFA-H I observations analysis with a 5′ FWHM kernel for 1.29-km s−1 and 0.184-km s−1 velocity resolutions, respectively.

Fig. 17 H I emission (top) and filamentary structures (bottom) toward a region dominated by vertical structures in the indicated radial velocity range.

Fig. 18 H I emission (top) and filamentary structures (bottom) for the region around the 〈z〉 dip toward ϕ ≈ 100° and 20 < Rgal < 25 kpc in Fig. 8.

5.1.2 Filaments and holes in other galaxies

Most the face-on nearby galaxies (2 ≲ d ≲ 15 Mpc) in The H I Nearby Galaxy Survey (THINGS, Walter et al. 2008) present holes in their inner regions, where the gas density and the star formation rate are higher, and elongated structures at larger radii. A detailed study of one of them, NGC 6946, a well studied spiral galaxy of about one-third of the Milky Way’s size located at 7.72 Mpc, revealed 121 H I holes in its inner regions and high-velocity gas associated with some of them, thus suggesting the presence of a fountain mechanism (Boomsma et al. 2008). It is not straightforward to demonstrate that the associations of vertical filaments in the inner Milky Way correspond to the holes seen in face-on nearby galaxies, particularly due to the difficulty in disentangling emission components and assigning reliable distances in our own Galaxy. However, there is significant evidence that toward some portions of the Galactic plane we are seeing the gas swept by the expansion of SN remnants, as it is the case of the H I shells identified in the range 10 kpc < R_gal 15 kpc toward QIV (McClure-Griffiths et al. 2002).

A recent study of the H I gas structure in a sample of 33 disk galaxies in THINGS revealed changes in the delta-variance spectrum in the inner and outer regions of the disks (Dib et al. 2021). The authors interpret these variations as the dominant effect of SN feedback in the inner parts of the galaxies and the prevalence of other large-scale dynamics in the outer parts. Certainly, SNe are not independent actors in the shaping the ISM. Recent numerical experiments indicate that feedback exclusively produced by SNe does not reproduce produce observed star cluster masses and warm gas volume-filling factors (Kim et al. 2018; Rathjen et al. 2021). Detailed modeling of the Galactic winds propelled by SNe indicate that cool and hot outflows can lift significant amounts of mass, energy, and metals (Kim et al. 2020). Our results indicate that the study of the H I filaments can shed light in the study of these Galactic winds and their distribution in the Milky Way.

5.2 Atomic filaments as a tracers of Galactic structure

5.2.1 Galactic warp

In addition to making the change in orientation with R_gal noticeable, the identification of filaments in the H I emission highlights the disk warp and other variations in the midplane height of the Milky Way. Figure 8 shows that the 〈z〉 features at R_gal > 20 kpc known as scalloping are highlighted by the filaments. This may be an indication that these features are not the product of global disk perturbations, as considered in Levine et al. (2006a), but can be the product of events such as gas inflow or the interaction with Milky Way satellites.

The most prominent variation in 〈z〉, aside from the warp, is the dip toward ϕ ≈ 100° and 20 < R_gal < 25 kpc. This feature, shown in Fig. 18, corresponds to a horizontal cloud at b ≈ −8° that contrasts with the location of most of the H I emission, which is concentrated at 0° < b < 10°. It is located in one of the few regions of the lυ-diagram that does not show a preferential filament orientation at 0° or 90°, as it is evidenced by the elongated structures oriented at 45°, among which is the giant atomic filament Cattail (Li et al. 2021). The filaments in this region seem to be pointing in the direction of the Magellanic stream and their location both in radial velocity and Galactic longitude seem to reinforce the potential association with the interwoven tail of filaments trailing the Magellanic Clouds in their orbit around the Milky Way (see D’Onghia & Fox 2016, and references therein). Although determining this association is beyond the scope of this work, it exemplifies the use of the H I filament properties as a tracer of the Galactic dynamics.

5.2.2 Filament scale height

Figure 11 indicates a difference of roughly a factor of two in the scale height traced by the bulk of the H I emission and just by the H I filaments at R_gal > 15 kpc. To explain this difference, we need to consider the physical conditions that may lead to the formation of these structures. Most likely the filaments are produced by large-scale compressive motions and thermal instability that produce overdensities in which the H I can cool rapidly (Wolfire et al. 2003). Thus, the filaments may correspond to a mixture of CNM and WNM surrounded by WNM. This scenario is supported by the relatively narrow emission line widths toward the filamentary structures, as shown in the example presented in Fig. 19.

Kalberla & Kerp (1998) estimated a scale height at the location of the Sun of about 150 pc for the for the CNM and 400 pc for the WNM, assuming hydrostatic equilibrium between the components. The authors inferred velocity dispersions of 4 and 24 km s⁻¹ for the CNM and WNM, which are roughly consistent with the combination of the turbulent and the thermal sound speeds of the gas in each phase. That result suggest that each component is separately in vertical equilibrium, so that the total turbulent and thermal midplane pressures balances the vertical weight of the ISM (see, for example, Ostriker & Shetty 2011). However, both the measured scale heights and the measured velocity dispersions are subject to relatively large uncertainties: for example, Heiles & Troland (2003) reports that the difference in scale heights of the two phases is closer to a factor of two, as result of their finding of large velocity dispersion for the CNM and smaller for the WNM.

Toward the inner Galaxy the scale height of the CNM has also been observed to be smaller than that of the WNM (Dickey & Lockman 1990; McClure-Griffiths & Dickey 2007). If we assume that the same conditions can be extended to the outer galaxy, the difference in scale heights reported in Fig. 11 is consistent with the expected difference for the two H I components. However, that conclusion is in tension with the results of the observations of H I absorption against strong continuum sources, which indicate that the scale height of the CNM is similar to that of the WNM, and both increase with Galactic radius in the outer disk (Dickey et al. 2009).

Liszt (1983) shows that one portion of the WNM must be associated with the CNM, while the other portion is more widely spread in the neutral medium between the CNM clouds. Dickey et al. (2009) argues that the scale height derived from the H I observed in absorption, which samples mostly the CNM clouds, and that from the H I emission, which samples both the CNM and WNM, are similar in the outer Galaxy, thus suggesting that the scale heights of the two H I components are not very different. One reason for the discrepancy with our results are the azimuthal variations that are evident in the radial profiles presented in Fig. 12.

Figure 7 of Dickey et al. (2009) shows that the best agreement between the scale heights sampled by absorption and emission is for the average over the azimuth range 0° < ϕ <180. In that range our analysis results in large dispersions, mostly due to the lower emission toward that portion of the Galaxy at R_gal > 20 kpc that results in fewer filaments. In the range 180° < ϕ <360, the agreement between the scale heights at R_gal > 17 kpc in Dickey et al. (2009) is less evident. It is precisely toward that azimuthal range that we find the clearest difference between the σ_z sampled by the filaments and by the bulk of the H I emission.

Our results are not conclusive toward 0° < ϕ < 180, where we do not find enough filaments to match the results of Dickey et al. (2009). However, we find a strong indication of different scale heights in the range 180° < ϕ < 360, where only 126 of 650 H I opacity observations available at the time of Dickey et al. (2009) are located. Whether or not the difference in the scale heights derived using the H I filaments are confirmed by the emission and opacity studies will depend on the analysis of future high-resolution surveys that increase the number of observations against continuum sources.

Fig. 19 H I emission (top) and selected spectra toward filamentary structures (bottom). The spectra correspond to the positions indicated by the crosses in the upper panel. The colored vertical bars in the bottom panel correspond to the three velocity channels presented in the corresponding colors in the top panel.

5.2.3 Column density in atomic filaments

We found that the ratio $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ is about 80% up to R_gal ≈ 15 kpc, as shown in Fig. 12. The reason for the decrease in $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ at R_gal > 15 kpc toward some lines of sight can be attributed to the lack of filamentary structure in the emission, particularly toward QI and QII. The explanation for the roughly constant values of $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ for R_gal < 15 kpc is potentially related to the almost constant CNM to WNM mass ratio throughout the Galaxy. The CNM fraction is almost constant with R_gal and corresponds to between 15 and 20% of the atomic gas in the outer disk (Dickey et al. 2009). The larger values of $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ can result from the filaments tracing not exclusively CNM but also some WNM associated with the cool gas. This, however, implies a very coarse average because on scales up to at least hundreds of parsecs the pressure and density of the ISM may be highly variable and far from equilibrium.

The decrease on the values of $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ happens at roughly the same R_gal at which we see the end of the regular trend in filament orientation illustrated in Fig. 7. This two observations can be interpreted as an indication of the change in the physical conditions of the atomic gas in the outer galaxy. These changes are not axisymmetric and can be caused by a variety of conditions that include the passage of the Milky satellites and variations in the ISM pressure and the stellar potential.

In principle, this observation implies that the formation of the structures that we observe as H I filaments is not affected by variations in star formation and stellar feedback with Galactocentric radius. Potentially converging flows compensate for the decrease in the SN rate with increasing distance from the Galactic center in maintaining the density enhancement that drives the warm gas into the cool structures that we observe as filaments. However, further interpretation calls for the detailed study of particular regions and the analysis of numerical simulations of Galactic disks.

5.3 Properties of atomic filaments

We have employed the word filaments to designate the elongated structures found in H I emission. We used their orientation to describe the structure of the emission and its anisotropy without assigning their existence to a particular kind of object in three-dimensional space. The fact that the preferential orientation is persistent across 1.29-km s⁻¹ velocity channels is a strong indication that this trend is produced by three-dimensional density structures in the ISM and not by fluctuations imprinted by the turbulent velocity field (see Clark et al. 2019, and references therein).

The deep examination of the nature of the filamentary structures in three dimensions, for example, to distinguish between filaments and sheets (Heiles 1997), is beyond the scope of this work. However, the characterization of the emission in terms of filaments highlights a large number of elongated structures that are reminiscent of the plumes in the H I emission toward the outskirts of nearby spiral galaxies (Walter et al. 2008). Recent studies have identified prominent elongated structures in the H I emission, for example, the 1.1-kpc-long Cattail filament with roughly 6.5 × 10⁴ M_⊙ (Li et al. 2021) and the 1.2-kpc-long Magdalena filament with approximately 7.2 × 10⁵ M_⊙ (Syed et al. 2022). Our results indicate that these elongated structures are the norm toward the outer Galaxy.

The Magdalena filament is singular given its orientation, parallel to the Galactic plane, and its distance from the H I emission midplane, roughly 500 pc. The kind of perturbations that can produce an object of these characteristics are currently unknown. However, they may be related to the horizontal orientation of the VHC complexes AC and H and as such the orientation of this structures can be used as an indicator of the dynamics of the Galactic disk and its interaction with the halo.

There are few detections of molecular gas toward the bands dominated by horizontal H I filaments in Fig. 1. Carbon monoxide (CO) observations indicate that there is barely any significant emission at these large radii (Dame et al. 2001; Miville-Deschênes et al. 2017). Toward the OSC, Dame & Thaddeus (2011) report ten detections in 220 pointings toward high-Hi density regions. Sun et al. (2015) reported the detection of 72 molecular clouds in the CO(1 → 0) emission toward the extension of the OSC into the second quadrant. Similar observations for the outer regions toward the outer third and fourth quadrants have not been made. Thus, it is hard at the moment to estimate if the H I filaments inherit their elongation and orientation into molecular clouds that they may help form within them, in contrast with the general trend for the clouds traced by CO in the inner Galaxy (Soler et al. 2021). Yet, the outer Galaxy reveals as an attractive laboratory to test our theories of molecular cloud formation and the coupling between the ISM phases.

6 Conclusions

We studied the orientation of the filamentary structure in the H I emission toward the Galactic plane. We have found that the orientation of the H I filaments with respect to the Galactic plane changes progressively with Galactocentric distance, from mostly perpendicular or having no preferred orientation in the inner Galaxy to mostly parallel up to R_gal ≈ 15 kpc and beyond. The radial dependence in the H I filament orientation strongly suggests that SNe’s energy and moment input is responsible for lifting the atomic gas in the inner Galaxy. Toward the outer Galaxy, the horizontal H I filaments are most likely the result of Galactic rotation and shear. Together these observations indicate that the atomic gas in the Milky Way follows a similar distribution as in nearby spiral galaxies, where holes potentially produced by SNe-propelled winds are found toward the inner portions, and arms and plumes dominate the outer ones.

The characterization of the H I structure in terms of filaments highlights the features in midplane height and dispersion around the midplane height that might be indicative of the local effects perturbing the Galactic disk, such as variations in the potential or the inflow of gas. We found that the scale height of the H I filaments is lower than that traced by the bulk of the H I emission by a factor of two. If the filaments are tracing cold gas clouds and the portion of warm gas associated with them, our results imply that the scale heights of the cold and warm gas are not close, as previously inferred from H I absorption observations.

We estimated that the amount of H I in filamentary structures is roughly constant with Galactocentric radius up to R_gal ≈ 15 kpc. This observation is potentially related to the mass ratio between the cold and warm H I phases. However, conclusively proving that hypothesis calls for further characterization of the H I phases, for example, through Gaussian decomposition and dedicated studies of the filament properties.

In this paper, we have presented a study of the characteristics of the H I emission. We did so without assuming that these characteristics define an idealized physical object in three-dimensional space, such as a spiral arm or a thin strand of gas. Through this generality, we identified a general anisotropy related to a physical quantity: the Galactocentric radius. Some of the emission features that we classified as filaments may correspond to thread-like gas clouds, but that was not the starting assumption of this study. Our conclusions were driven by a dominant trend found through the data analysis rather than adapting the data to a preconceived notion about the gas structure. Given the diversity of processes shaping the ISM, this data-driven approach offers a promising path to study the cycles of mass and energy that lead to the formation and destruction of star-forming molecular clouds. Even if the observed lack of correlation between star formation and neutral atomic hydrogen at low surface densities has traditionally excluded the H I observations from the global study of star-formation rates in galaxies (see, for example, Kennicutt 1989; Leroy et al. 2008; Schruba et al. 2011), the H I still carries crucial information on the dynamics that drive the diffuse gas into the condensations that become new stars.

Acknowledgements

J.D.S., R.K.S., and S.C.O.G. are funded by the European Research Council via the ERC Synergy Grant “ECOGAL – Understanding our Galactic ecosystem: From the disk of the Milky Way to the formation sites of stars and planets” (project ID 855130). R.S.K., H.B. and S.C.O.G. acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Center (SFB 881, Project-ID 138713538) “The Milky Way System” (subprojects A1, B1, B2 and B8) and from the Heidelberg cluster of excellence (EXC 2181 – 390900948) “STRUCTURES: A unifying approach to emergent phenomena in the physical world, mathematics, and complex data”, funded by the German Excellence Strategy. R.J.S. acknowledges funding from an STFC ERF (grant ST/N00485X/1) and HPC from the DiRAC facility (ST/P002293/1). The team at Z.A.H. acknowledges computing resources and data storage facilities provided by the State of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 35/1134-1 FUGG and INST 35/1503-1 FUGG, and they thank for computing time from the Leibniz Computing Center (LRZ) in project pr74nu. We thank the anonymous referee for the thorough review and appreciate the suggestions, which helped us improve the manuscript. J.D.S. thanks the following people who helped with their encouragement and conversation: Jonathan Henshaw, Jonas Syed, Eleonora Zari, Milena Benedettini, Robert Benjamin, and Thomas Henning. Part of the crucial discussions that lead to this work took part under the program Milky-Way-Gaia of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. The computations for this work were made at the Max-Planck Institute for Astronomy (MPIA) astro-node servers. This work has been written during a moment of strain for the world and its inhabitants. It would not have been possible without the effort of thousands of workers facing the COVID-19 emergency around the globe. Our deepest gratitude to all of them. Software: astropy (Astropy Collaboration 2018), SciPy (Virtanen et al. 2020), magnetar (Soler 2020).

Appendix A The Hessian analysis method

A.1 Studying angle distributions using circular statistics

Throughout this paper we have characterized the distribution of orientations using circular statistics, the branch of statistics that deals with observations that are directions (see, for example, Batschelet 1981; Mardia & Jupp 2009). The probability distributions against which we are comparing the data are angles. The most natural characteristic to describe the angle distributions is the mean, which is defined as $$\langle \theta \rangle =\text{arctan}\left(\frac{\sum w_i\text{sin}\left(2{\theta }_i\right)}{\sum w_i\text{cos}\left(2{\theta }_i\right)}\right),$$(A.1)

where the index i runs over the elements in a group of angle observations and w_i is the statistical weight corresponding to each of them.

The mean angle, 〈θ〉, is insufficient to fully describe the distribution of angles. For example, the angle distribution is ill-defined if the distribution is uniform or not unimodal. To test whether the angle distribution is uniform and unimodal one can use the Rayleigh statistic, which is defined as $$Z={\left[{\left(\sum w_i\text{cos}\left(2{\theta }_i\right)\right)}^2+{\left(\sum w_i\text{sin}\left(2{\theta }_i\right)\right)}^2\right]}^{1/2}.$$(A.2)

This value can be understood as the net displacement of a random walk in the complex plane with step sizes w_i. For a uniform distribution of angles, Z ≈ 0.

In our application there is one direction that produces an expected anisotropy in the filament orientation: the Galactic plane. That is why we used the projected Rayleigh statistic (V, Eq. 4) to test the orientation angle distribution against 0° and 90°, or parallel and perpendicular to the Galactic plane, respectively. The observations also reveal that these two orientations are predominant in the H I filaments.

Figure A.1 shows the results of our filament orientation analysis in terms of 〈θ〉 and Z. It is evident that the majority of the tiles have filamentary structures with mean orientations 〈θ〉 ≈ 0°. However, this observation alone is not enough to conclude the significance of this trend. The values of Z indicate that there is a significant tendency for the angle distributions to be unimodal. The highest significance is found in two bands located at the minimum and maximum υ_LSR toward the first two and the last two Galactic quadrants, respectively. These bands correspond to the high-V bands shown in Fig. 1.

A summary of the trends obtained using the circular statistic is shown in Fig. A.2. We found that the large majority of the tiles with significantly unimodal distributions, Z ≳ 3, correspond to large values of |V|. This confirms that 0° and 90° correspond to the most represented orientations in the H I filaments and justified the use of these two reference angles in the application of the calculation of V.

A.2 Selection of parameters in the Hessian analysis method

In the main body of this paper we have reported the results of a particular set of parameters in the Hessian analysis and the study of the filament orientations. In this section we present the parameter study that led to that selection. In what follows we detail the considerations related to the selection of the tile sizes, the derivative kernel size employed for the calculation of the Hessian matrix, and the thresholds in intensity and curvature used to select the filamentary structures. We considered the variations introduced by the changes in one of these parameters and leaving the other fixed to the values presented in Table. 2.

Fig. A.1 Same as Fig. 1, but for the filament orientation quantified by the mean orientation angle (〈θ〉, Eq. A.1, top) and the Rayleigh statistic (Z, Eq. A.2, bottom).

A.2.1 Tile sizes

We selected the |b| ≤ 10° range for the analysis aiming to resolve the warp and the flaring of the Galactic disk. However, we examined the dependence of the H I filament orientation trends on the tile size. Figure A.3 shows the V lυ-diagrams corresponding the square tiles centered on b = 0° and with side lengths of 3°, 5°, and 20°.

The 3° × 3° tiles cover the range included in the interferometric Galactic plane surveys THOR, |b| < 1.°25, and SGPS, |b| < 1.°4. The corresponding V lυ-diagram, shown in the top panel of Fig. A.3, indicates that the high-V band reported in the main body of the paper is not found in that range toward the first quadrant. This is due to the Galactic warp, which reaches midplane heights around 3.0 kpc for R_gal ≳ 15 kpc toward the first quadrant, as discussed in Sect. 5.2.

The 5° × 5° tiles cover the range included in the 1′-resolution interferometric Canadian Galactic Plane Survey (CGPS Taylor et al. 2003). The corresponding V lυ-diagram, shown in the middle panel of Fig. A.3, indicates that the high-V band in the first quadrant is present in this range, but its significance is relatively low. The high-V band toward the last two quadrants is clearly distinguishable in both the 3° × 3° and 5° × 5° tiles analysis, because the Galactic warp in that direction is lower. Finally, the bottom panel of Fig. A.3 shows the results corresponding to the 20° × 20° tiles.

The decrease in the tile size leads to a reduction of the maximum V values, as it is evident from the comparison of the panels in Figure A.3, which is the result of the lower number of pixels or, equivalently, independent filament orientations. However, the positive V in radial velocities corresponding to the outer galaxy and negative V toward the inner galaxy, which are highlighted as the main results of this paper, are persistent for the considered tile size selections. The selection of broader tiles washes away some of the azimuthal variations. Thus, we compromise between b-coverage and resolution in l by using rectangular tiles.

Figure A.4 illustrates that the selection of square tiles increases the significance in the signal produced by mostly vertical structures. However, this increase comes at the cost of reducing the number of tiles in Galactic longitude and significantly increasing the signal coming from the local emission. All things considered, we chose rectangular tiles to highlight the significance of the filamentary structure in the outer Galaxy and produce a more detailed face-on reconstruction of the filament orientation trends.

Fig. A.2 Two-dimensional histogram of the Rayleigh statistic (Z) and projected Rayleigh statistic (V) calculated for the filamentary structures in the HI4PI observations.

A.2.2 Tile aspect ratio

Square tiles are a natural selection to avoid any bias toward one of the preferential orientations studied in this analysis, that is, 0° or 90°, as discussed in Soler et al. (2020). Yet the use of rectangular tiles improves the resolution in l and the study of azimuthal effects. The 3°× 20° × 1.29 km s⁻¹ tiles selected for the analysis in the main body of the paper result from the compromise between azimuthal resolution and the V significance. To reach that selection, we empirically studied the effect of the tile aspect ratio on the values of V.

Figure A.5 presents the results obtained using 10° × 20° × 1.29 km s⁻¹ and 5°× 20° × 1.29 km s⁻¹ tiles. We found that the main features in the V lυ-diagram persist with the change in the tile aspect ratio. The increase in the l resolution resolves the negative-V patches in the bottom panel of Fig. A.3 into regions surrounded by V = 0 ranges. The rectangular tiles also constrain the l range to which the dominant features in the map are assigned. For example, the increase in the l resolution resolved the separation between the high-V band in the second quadrant and the HVC complex located around 155° < l < 190° and −340 < υ_lsr < −135km s⁻¹.

We found that for tiles with 3° coverage in l, the values of V are reduced below the level of significance in most of the l diagram. This effect is particularly acute for the V < 0, where the narrow tiles reduce the coverage of close-to-vertical structures and introduce a bias for V > 0. The high V-bands, which are mostly produced by a series of horizontal filaments distributed over a relatively broad range in b, are less affected by the tile aspect ratio.

Figure A.6 illustrates the effect of a broader tile selection in the trends reported in Fig. A.6. The 10° × 20° × 1.29 km s⁻¹ tiles maintain a high significance in the regions where the filaments are preferentially oriented perpendicular to the Galactic plane. This tile aspect ratio emphasizes the change in the orientation of the filaments, from mostly perpendicular to mostly parallel to the Galactic plane with increasing Galactocentric radius. However, the 3° × 20° × 1.29 km s⁻¹ tiles used to report the results in the mainbody of the paper increase the sampling in Galactic longitude, thus producing a more robust identification of the increasing values of V with R_gal.

A.2.3 Kernel sizes

We have chosen a 18′ FWHM derivative kernel to profit from the angular resolutions of the observations. The use of larger kernels implies necessarily a lowering of the significance of the filament orientations, given the reduction in the number of independent filament orientations. However, we tested the results of the analysis of 3° × 20° × 1.29 km s⁻¹ tiles using larger derivative kernels.

Figure A.7 show the V lυ-diagrams obtained using derivative kernels two and three times larger than the beam size, 32.′4 and 48.′6 FWHM, respectively. Besides the expected reduction in the V values, the larger derivative kernels displace the high-V bands to larger radial velocities This is due to the loss of significance produce by the reduction in the number of independent filament orientation measurements. Rather than producing more horizontal structures, which would increase the values of V, the smoothing with a larger derivative kernel reduces the preferential orientation in the filamentary structures.

Figure A.8 illustrates that the larger derivative kernels systematically reduce the values of V independently of the kinematic distance corresponding to each tile. Rather than artificially introducing high-V values at closer distances, the smoothing results in a loss of the significance in the relative orientation trends. This reinforces the argument that the high-V bands are not a product of the angular resolution of the HI4PI survey, which we proved using the 4′ resolution observations in GALFA-H I. For the sake of simplicity, we present only the results corresponding to the far distance, but the same is inferred for the nearby solution to the near kinematic distance.

Figure A.9 shows the ratio between the column density in filaments and in the bulk of the H I emission ($N_{\text{H}}^{\text{fil}}/N_{\text{H}}$) for the two derivative kernels. We found that the increase in the derivative kernel size does not substantially change the fraction of column density in filaments or the R_gal above which we found a drop in $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$. This result suggests that the drop in $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ is not an effect introduced by the distance. If distance was producing the drop, we would see a systematic shift in the R_gal values where $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ starts to decrease. The presumed shift in R_gal should be lower for the larger derivative kernel but this is not what is found in the data, where we found the drop at the same R_gal for the three considered kernel sizes.

Since we fixed the I S/N threshold to test the effect of the derivative kernel size, we also conclude that the drop in $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ is not simply caused by a drop in the emission. The derivative kernels effectively smooth the velocity channel maps, thus bringing a larger portion of the emission above the I S/N threshold. However, changing this smoothing by up to a factor of three above the beam size does not substantially change the trends in $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$. Therefore, we conclude that the drop in $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ is related to a transition toward less-structured H I emission rather than a change introduced by the decrease in H I emission with increasing R_gal.

A.2.4 Intensity threshold

We empirically evaluated the effect of the intensity threshold in the filament orientation analysis. Figure A.10 shows two examples of this selection above and below the threshold used for the results presented in the main body of this paper. The results corresponding to the 3σ_I intensity threshold show a vertical stripe around l ≈ 140°. These features are typical for low-intensity thresholds and are produced by features in low-intensity velocity channels. Most of these features disappear with a 5σ_I threshold, as shown in the right-hand side panel of Fig. 1.

We do not find significant differences in the results produced by a 9σ_I threshold. The increase in the I threshold tends to eliminate low-intensity velocity channels at the extremes of the velocity range rather than modifying the general trends in the V lυ diagram. This relatively low dependence on the intensity threshold is partly due to the stringent criterion in curvature, λ₋, which is dominant in determining the structures selected for the filament orientations analysis.

A.2.5 Curvature threshold

The lowest eigenvalue of the Hessian matrix, λ₋, plays a crucial role in the identification of the filamentary structure (see, for example, Planck Collaboration Int. XXXII 2016; Schisano et al. 2014). The filamentary structures correspond to the pixels where λ₊λ₋ < 0 and λ₋ < 0. In the observations, however, there is a potential contribution of the noise to the curvature that we account for by choosing a lower curvature ${\lambda }_{-}^{\text{thres}}$.

We present an example of the results obtained when using no curvature threshold, just λ₋ < 0, in Fig. A.11. The example tile shows a collection of low-I and relatively low curvature, λ₋ ≈ 0, filamentary structures. These faint filaments objects appear systematically over the whole map and do not show an evident preferred orientation.

The top panel of Fig. A.12 shows the net effect of the low-curvature structures. Globally, the V lυ-diagram presents the same high-V bands and the low V regions reported in Fig. 1. However, the |V| values are lower, which corresponds to a lower significance in the filament orientation trends. We also found a series of vertical stripes in the V lυ-diagram which are particularly prominent in the region υ_LSR < 0 toward the first two quadrants and υ_LSR > 0 toward the last two. These stripes are mitigated by applying a curvature threshold.

The bottom panel of Fig. A.12 shows the results obtained by using a curvature threshold based on the λ₋ values found in the lowest intensity-S/N channels toward the position of each tile. This selection increases the significance of the relative orientation trends, quantified by V, We selected the median of these per-tile λ₋ threshold as the value used for the results reported in the main body of this paper.

Fig. A.3 Longitude-velocity (lυ) diagrams of the H I filament orientation quantified by the projected Rayleigh statistic (V, right). Each of the pixel elements in the diagrams corresponds to a 3°× 3°× 1.29 km s−1 (top), 5° × 5°× 1.29 km s−1 (middle), and 20° × 20° × 1.29 km s−1 (bottom). These results correspond to the Hessian analysis performed using a 18′ FHWM derivative kernel. The overlaid curves correspond to the spiral arm features discussed in the main body of the paper.

Fig. A.4 Same as Fig. 7, but for tile sizes 20° × 20° × 1.29 km s−1.

Fig. A.5 Same as Fig. A.3, but for tile sizes 10° × 20° × 1.29 km s−1 and 5° × 20° × 1.29 km s−1, shown in the top and bottom panels, respectively.

Fig. A.6 Same as Fig. 7, but fortile sizes 10° × 20° × 1.29 km s−1.

Fig. A.7 Same as Fig. A.3, but for 3° × 20° × 1.29 km s−1 tiles analyzed using 32.′4 and 48.′6 FWHM derivative kernels, shown in the top and the bottom panels, respectively.

Fig. A.8 Same as Fig. 4, but for32.′4 and 48.′6 FWHM derivative kernels, shown in the top and the bottom panels, respectively.

Fig. A.9 Same as Fig. 12, but for 32.′4 and 48.′6 FWHM derivative kernels, shown in the top and the bottom panels, respectively.

Fig. A.10 Same as Fig. A.7, but for 3° × 20° × 1.29 km s−1 tiles and intensity thresholds of 3σ (top), 5σ (middle), and 9σ (bottom).

Fig. A.11 Example of the filamentary structures identified using the Hessian method with curvature thresholds λ− < 0 (top) and λ− < −1.8 Kdeg−2 (bottom). The three maps represent: Left. H I intensity after the smoothing by a Gaussian kernel with the size of the derivative kernel. Center Second eigenvalue of the Hessian matrix, λ−, calculated using Eq. 2. Left. Filament orientations, θ, calculated using Eq. 3. The radial velocity corresponds to a kinematic distance d = 12.7 kpc and Galactocentric radius Rgal = 16.1 kpc, as derived using the assumptions in App. B

Fig. A.12 Same as Fig. A.10, but for λ− < 0 (top) and a tile-dependent λ−-threshold (bottom).

Appendix B De-projection of the spectral observations

B.1 Construction of the face-on view

We de-projected the results our analysis results in the l, b, and υ_LSR grid into Galactocentric cylindrical coordinates R_gal, ϕ, and z. We calculated a Galactocentric radius, R_gal for each tile’s central position and radial velocity using the Kinematic Distance Utilities: KDUtils⁶ (Wenger et al. 2018). In practice, these routines solve $$R_{\text{gal}}=R_0\text{sin}\left(l\right)\frac{V\left(R_{\text{gal}}\right)}{{\upsilon }_p+V_0\text{sin}\left(l\right)},$$(B.1)

where R₀ and V₀ are the Galactocentric radius and tangent velocity of the solar orbit around the Galactic center, V(r) is the rotation curve, and υ_p is the radial velocity, which at the center of the tile is equal to υ_LSR. In this study we used the rotation curve derived in Reid et al. (2019).

We also evaluated the Galactocentric azimuth angle $$\varphi =\text{arccos}\left(\frac{R_0-d\text{cos}l}{R_{\text{gal}}}\right),$$(B.2)

where d is the distance given by $$d=R_0\text{cos}\left(l\right)\pm \sqrt{R_{\text{gal}}^2-R_0^2{\text{sin}}^2\left(l\right).}$$(B.3)

In the inner Galaxy (R_gal < R₀) this expression yields to two distances along the line of sight, the near and the far kinematic distances, located on either side of the tangent point. Figure B.1 illustrates the near- and far-distances obtained for the positions of the 3° × 20° × 1.29 km s⁻¹ tiles in the lυ-diagram.

The ambiguity introduced by the two solutions to Eq. B.3 makes the determination of kinematic distances very challenging in the inner Galaxy. Thus, we restricted the results of our de-projection to the range R_gal > 10 kpc to avoid the artifacts produced by the motions around the solar orbit. We also exclude the longitude ranges l < 15°, l > 345°, and 165° < l < 195°, where the radial velocities are not dominated by the circular motions and it is difficult to establish reliable distances (Levine et al. 2006b).

To test the potential effect of distance in the filament orientation trends, we presented the V distribution for the “far” heliocentric kinematic distances, d_far in Fig. 4. For the sake of completeness, we present the same diagram considering the “near” heliocentric distances in Fig. B.2, which differ from d_far only toward QI and QIV. We also present the unmasked face-on reconstruction of H I filament orientations obtained by using the far heliocentric distances in Fig. B.3.

B.2 Galactic warp and flaring

For each pixel in a tile, we calculate the height z with respect to b = 0° by computing $$z=d\text{arctan}\left(b\right).$$(B.4)

This approach assumes a single representative distance for each tile, that is, it does not account for the variation between b = 0° and b = ± 10°. The maximum differences in R_gal and d between the tile top and center are around 0.074 kpc and 0.076 kpc, respectively. These are found in tiles at large R_gal and d, where the variations in the two quantities due to the tile height are negligible compared to the uncertainties introduced by the resolution in radial velocity.

Our de-projection approach is different from the PPV grid interpolation and mapping into a fixed PPP grid followed in, for example, Levine et al. (2006b), Kalberla & Dedes (2008), and Koo et al. (2017). However, it does not yield significant differences in the results obtained for the midplane displacement (warp) and dispersion around the midplane (flaring), as shown in Sect. 4.2, and it is computationally more efficient. Our analysis does not require the continuous sampling of the PPP grid necessary for the Fourier mode analysis in Levine et al. (2006b) or the three-dimensional mass reconstruction presented in (Nakanishi & Sofue 2016), but rather aims to illustrate the variation in filament orientation with increasing R_gal. Given the one-to-one mapping from υ_LSR into R_gal in the outer galaxy under the assumption of circular motions, our approach is close to a discrete sampling of the PPP space.

B.3 Warp and flaring

For each one of the tiles de-projected into the Galactocentric cylindrical coordinates, we calculated the mid-plane displacement with respect to the b =0° position, or warp, by calculating the first-moment of the height z weight by the intensity, $$\langle z\rangle =\frac{\sum I_{ij}{z}_{ij}}{\sum I_{ij}},$$(B.5)

where the indices i and j run over the tile pixels.

We also calculated the dispersion around the mid-plane height, or flaring, by computing the second-moment of the height z weight by the intensity, $${\sigma }_z={\left[\frac{\sum I_{ij}{\left(z_{ij}-\langle z\rangle \right)}^2}{\sum I_{ij}}\right]}^{1/2}$$(B.6)

where 〈z〉 is calculated using Eq. B.5.

Fig. B.1 Near and far kinematic distances assigned to the 3° × 20° × 1.29 km s−1 tiles used to construct the lυ-diagrams in Fig. 1. These results correspond to the solutions of Eq. B.3 using the rotation curve presented in Reid et al. (2019).

Fig. B.2 Same as Fig. 4, but for the near heliocentric kinematic distances.

Fig. B.3 Same as Fig. 8, but without masking Rgal < 10kpc and assigning all kinematic heliocentric distances to the largest solution to Eq. B.3.

Appendix C Comparison between HI4PI and GALFA-H I

The only survey with enough coverage of the Galactic plane to expand the results of the HI4PI analysis to higher resolutions is GALFA-H I. With 4′ FWHM resolution, GALFA-H I provides a factor of four improvement with respect to HI4PI in the ranges 36.°5 < l < 70.°5 and 175.°3 < l < 209.°4 We applied the Hessian matrix analysis to the GALFA-H I data using the parameters presented in Table. 2. The analysis was applied to the data in its native velocity resolution, 0.184 km s⁻¹, and to the data smoothed in the spectral axis to match the HI4PI velocity resolution. As with the HI4PI data, the curvature threshold (${\lambda }_{-}^C$) was estimated empirically using noise-dominated velocity channels.

We present an example of the Hessian analysis applied to one 3° × 20° × 1.29 km s⁻¹ tile in the HI4PI and GALFA-H I observations in the top panels of Fig. C.1. This example tile illustrates the connection in the emission structures traced by the two surveys at different angular resolutions and projected into the same radial velocity resolution.

We also present an example of the Hessian analysis applied to the GALFA-H I observations in its native 0.184km s⁻¹ velocity resolution in the bottom panels of Fig. C.1. These tiles, which correspond to the same line of sight of those in the top panels, illustrate the similarity of the intensity structures in contiguous narrow velocity channels. This similarity, caused by the linewidth of the H I emission, is evident for the whole velocity range in Fig. 15.

We also tested the results of a larger derivative kernel and a corresponding ${\lambda }_{-}^C$, as detailed in Table C.1. We present the equivalent of Fig. 13 for the 8′ FWHM kernel in Fig. C.2. As expected from the lower number of gradients derived with a large kernel, the values of V are reduced. The equivalent of Fig. 15 for the 8′ FWHM kernel, shown in Fig. C.3, indicates a significant increase in the values of V with respect to those obtain for the HI4PI data.

Fig. C.1 Hessian matrix analysis applied to the HI4PI and the GALFA-H I observations. Top. Results obtained on a 3° × 20° × 1.29 km s−1 tile in the HI4PI and the GALFA-H I observations, shown in the left and right panels, respectively. Bottom. Results obtained on two 3° × 20° × 0.184 km s−1 tiles in the GALFA-H I observations. In each of the panels, the left-hand side map shows the H I emission convolved with a 2D Gaussian with the same diameter of the derivative kernel. The map in the center shows the second eigenvalues of the Hessian matrix, λ−, which highlights the filamentary structure. The right-hand side map shows the orientation angle of the filamentary structure.

Fig. C.2 Same as Fig. 13, but for an 8′ FWHM kernel.

Fig. C.3 Same as Fig. 15, but for an 8′ FWHM kernel.

Appendix D Filamentary structures in high velocity clouds

We present the result of the Hessian matrix analysis toward two high velocity cloud (HVC) complexes in Fig. D.1 and Fig. D.2.

Fig. D.1 Atomic hydrogen (H I) emission (top) and filamentary structures identified using the Hessian matrix method (bottom) toward the HVC complex H.

Fig. D.2 Same as Fig. D.1, but for the anticenter HVC complex.

References

All Tables

All Figures

Fig. 1 Longitude-velocity (lυ) diagrams of the mean intensity (〈I〉, left) and the filament orientation quantified by the projected Rayleigh statistic (V, right). Each of the pixel elements in the diagrams corresponds to a 3° × 20° × 1.29 km s−1 velocity channel map centered on b = 0°, which we call tile throughout out this paper. These results correspond to the Hessian analysis performed using an 18′ FHWM derivative kernel. Values of V > 0 (red) or V < 0 (blue) indicate a preferential orientation of the filaments parallel, θ = 0°, or perpendicular, θ = 90°, to the Galactic plane. The 3σ statistical significance for these two orientations corresponds to V > 2.87 or V < −2.87, respectively. The overlaid curves correspond to the main spiral arms features presented in Reid et al. (2016), the outer Scutum-Centaurus arm (OSC, Dame & Thaddeus 2011), and the extended outer arm (M-G2004, McClure-Griffiths et al. 2004).

In the text

Fig. 2 Mean H I emission intensity (〈I〉) and projected Rayleigh statistic (V, Eq. (4)) for the 3° × 20° × 1.29 km s−1 tiles. The colors indicate the tile percentage covered by the filamentary structures selected following the criteria described in Sect. 3.2. The horizontal dotted lines show V = ±2.87, which indicate the 3σ significance in the preferential orientation parallel, θ = 0°, or perpendicular, θ = 90°, to the Galactic plane, respectively. The vertical dashed line indicates the intensity threshold.

In the text

Fig. 3 Examples of the H I emission and filamentary structures toward the l and υLSR ranges with predominantly horizontal filamentary structures. The top image in each panel corresponds to the H I emission intensity in three velocity-channel maps presented in red, green, and blue colors. The bottom image corresponds to the second eigenvalue of the Hessian matrix, λ−, which highlights the filamentary structures in each velocity-channel map. The white disk in the lower right corner corresponds to the size of the derivative kernel used to calculate the Hessian matrix.

In the text

Fig. 4 Kinematic heliocentric distance (d) and projected Rayleigh statistic (V, Eq. (4)) for the 3° × 20° × 1.29 km s−1 tiles. We report the largest (far) solution to the heliocentric distance equation, as detailed in Appendix B. From top to bottom, the panels correspond to each Galactic quadrant and to the whole Galactic longitude range. The colors indicate the Galactic longitude of the tiles. The horizontal dashed lines indicate V = ±2.87, which correspond to the 3σ significance in the preferential orientation around θ = 0° or θ = 90°, respectively. The solid line in the bottom panel corresponds to the mean V values in 1-kpc dfar bins.

In the text

Fig. 5 Mean intensity (〈I〉, top), projected Rayleigh statistic (V, middle), and derivative kernel size (Δ, bottom) for 3° × 20° × 1.29 km s−1 toward the indicated positions. The two lines in the middle panel indicate the results obtained with a derivative kernel with a fixed angular size and a kernel with a fixed physical size, as indicated on the legend. The kernel with a fixed physical size is exclusively used in the radial velocity range corresponding to kinematic distances where the angular size of the kernel is below 1°.

In the text

Fig. 6 Face-on view of the mean H I emission intensity (〈I〉; left) and the filament orientation quantified by the projected Rayleigh statistic (V; right). Values of V > 2.87 or V < −2.87 indicate a statistical significance above 3σ in the preferential orientation around θ = 0° or θ = 90°, respectively. The yellow circular marker indicates the position of the Sun. The overlaid curves correspond to the spiral arm models identified in Reid et al. (2019), shown in black color, and their extrapolation to larger Galactocentric radii, shown in green.

In the text

Fig. 7 Galactocentric radius (Rgal) and projected Rayleigh statistic (V, Eq. (4)) for the 3°× 20° × 1.29 km s−1 tiles. From top to bottom, the panels correspond to each Galactic quadrant and to the whole Galactic longitude range. The colors indicate the Galactic longitude of the tiles. The horizontal dashed lines indicate V = ±2.87, which indicate the 3cc significance in the preferential orientation around θ = 0° or θ = 90°, respectively. The vertical dotted line indicates the radius of the solar orbit. The solid line in the bottom panel corresponds to the mean V values in 1-kpc Rgal bins. The vertical solid lines mark the central locations of the H I shells identified in the LAB survey (Ehlerová & Palouš 2013) and the SGPS (McClure-Griffiths et al. 2002) in the range |b| < 10°, shown in gray and green colors, respectively.

In the text

	Fig. 8 Face-on view of the H I midplane height (〈z〉) reconstructed from all the emission in the range |b| < 10° (left) and the emission in filamentary structures (right).
In the text

Fig. 9 Midplane height (〈z〉) and Galactocentric radius (Rgal) for the 3° × 20° × 1.29 km s−1 tiles. The left-hand side panel corresponds to all the H I emission in the |b| < 10° range. The right-hand side panel corresponds to the emission in the structures classified as filaments. The colors indicate the Galactic longitude, l. The vertical dotted line indicates the radius of the solar orbit around the Galactic center.

In the text

	Fig. 10 Face-on view of the dispersion around the midplane height (σz) reconstructed from all the emission in the range |b| < 10° (left) and the emission in filamentary structures (right).
In the text

	Fig. 11 Same as Fig. 9, but for the dispersion around the midplane height (σz).
In the text

Fig. 12 Radial profiles of the ratio between the column density in the filamentary structures $N_{\text{H}}^{\text{fil}}$ and the total column density (NH). From top to bottom, the panels correspond to each Galactic quadrant and to the whole Galactic longitude range. The colors indicate the Galactic longitude of the tiles. The vertical dotted line indicates the radius of the solar orbit. The solid line in the bottom panel indicates the mean values of $N_{\text{H}}^{\text{fil}}/N_{\text{H}}$ in 1-kpc Rgal bins.

In the text

Fig. 13 Longitude-velocity (lυ) diagrams of the filament orientation quantified by the projected Rayleigh statistic (V, right) in the two portions of the Galactic plane covered in the range |b| ≤ 10° by both HI4PI and GALFA-H I. TheGALFA-H I observations were analyzed at 1.29 and 0.184 km s−1 resolution using a 5′ FHWM derivative kernel. The curves correspond to the same spiral arm features introduced in Fig. 1.

In the text

	Fig. 14 Examples of the H I emission and filamentary structures identified in the HI4PI and the GALFA-H I data, shown in the top and bottom panels respectively. The disk in the lower right corner of each panel indicates the size of the derivative kernels used for the calculation of the Hessian matrix, 18′ and 5′ for the HI4PI and the GALFA-H I observations, respectively.
In the text

	Fig. 15 Projected Rayleigh statistic (V) calculated for the filamentary structure in the HI4PI and GALFA-H I observations. Each point corresponds to a 3° × 20° × 1.29 km s−1 tile. The gray dashed lines indicate V = 2.87 and V = −2.87, which correspond to a 3σ significance in the preferential orientation around θ = 0° or θ = 90°, respectively.
In the text

Fig. 16 Mean intensity (〈I〉, top), projected Rayleigh statistic (V, middle), and mean orientation angle of the filaments (〈θ〉, bottom) for the tiles centered on l = 202.5 in the HI4PI and GALFA-H I. The left panel corresponds to the HI4PI observations analysis with a 18′ FWHM kernel for a 1.29-km s−1 velocity resolution. The center and right panels correspond to the GALFA-H I observations analysis with a 5′ FWHM kernel for 1.29-km s−1 and 0.184-km s−1 velocity resolutions, respectively.

In the text

	Fig. 17 H I emission (top) and filamentary structures (bottom) toward a region dominated by vertical structures in the indicated radial velocity range.
In the text

	Fig. 18 H I emission (top) and filamentary structures (bottom) for the region around the 〈z〉 dip toward ϕ ≈ 100° and 20 < Rgal < 25 kpc in Fig. 8.
In the text

	Fig. 19 H I emission (top) and selected spectra toward filamentary structures (bottom). The spectra correspond to the positions indicated by the crosses in the upper panel. The colored vertical bars in the bottom panel correspond to the three velocity channels presented in the corresponding colors in the top panel.
In the text

	Fig. A.1 Same as Fig. 1, but for the filament orientation quantified by the mean orientation angle (〈θ〉, Eq. A.1, top) and the Rayleigh statistic (Z, Eq. A.2, bottom).
In the text

	Fig. A.2 Two-dimensional histogram of the Rayleigh statistic (Z) and projected Rayleigh statistic (V) calculated for the filamentary structures in the HI4PI observations.
In the text

Fig. A.3 Longitude-velocity (lυ) diagrams of the H I filament orientation quantified by the projected Rayleigh statistic (V, right). Each of the pixel elements in the diagrams corresponds to a 3°× 3°× 1.29 km s−1 (top), 5° × 5°× 1.29 km s−1 (middle), and 20° × 20° × 1.29 km s−1 (bottom). These results correspond to the Hessian analysis performed using a 18′ FHWM derivative kernel. The overlaid curves correspond to the spiral arm features discussed in the main body of the paper.

In the text

	Fig. A.4 Same as Fig. 7, but for tile sizes 20° × 20° × 1.29 km s−1.
In the text

	Fig. A.5 Same as Fig. A.3, but for tile sizes 10° × 20° × 1.29 km s−1 and 5° × 20° × 1.29 km s−1, shown in the top and bottom panels, respectively.
In the text

	Fig. A.6 Same as Fig. 7, but fortile sizes 10° × 20° × 1.29 km s−1.
In the text

	Fig. A.7 Same as Fig. A.3, but for 3° × 20° × 1.29 km s−1 tiles analyzed using 32.′4 and 48.′6 FWHM derivative kernels, shown in the top and the bottom panels, respectively.
In the text

	Fig. A.8 Same as Fig. 4, but for32.′4 and 48.′6 FWHM derivative kernels, shown in the top and the bottom panels, respectively.
In the text

	Fig. A.9 Same as Fig. 12, but for 32.′4 and 48.′6 FWHM derivative kernels, shown in the top and the bottom panels, respectively.
In the text

	Fig. A.10 Same as Fig. A.7, but for 3° × 20° × 1.29 km s−1 tiles and intensity thresholds of 3σ (top), 5σ (middle), and 9σ (bottom).
In the text

Fig. A.11 Example of the filamentary structures identified using the Hessian method with curvature thresholds λ− < 0 (top) and λ− < −1.8 Kdeg−2 (bottom). The three maps represent: Left. H I intensity after the smoothing by a Gaussian kernel with the size of the derivative kernel. Center Second eigenvalue of the Hessian matrix, λ−, calculated using Eq. 2. Left. Filament orientations, θ, calculated using Eq. 3. The radial velocity corresponds to a kinematic distance d = 12.7 kpc and Galactocentric radius Rgal = 16.1 kpc, as derived using the assumptions in App. B

In the text

	Fig. A.12 Same as Fig. A.10, but for λ− < 0 (top) and a tile-dependent λ−-threshold (bottom).
In the text

	Fig. B.1 Near and far kinematic distances assigned to the 3° × 20° × 1.29 km s−1 tiles used to construct the lυ-diagrams in Fig. 1. These results correspond to the solutions of Eq. B.3 using the rotation curve presented in Reid et al. (2019).
In the text

	Fig. B.2 Same as Fig. 4, but for the near heliocentric kinematic distances.
In the text

	Fig. B.3 Same as Fig. 8, but without masking Rgal < 10kpc and assigning all kinematic heliocentric distances to the largest solution to Eq. B.3.
In the text

Fig. C.1 Hessian matrix analysis applied to the HI4PI and the GALFA-H I observations. Top. Results obtained on a 3° × 20° × 1.29 km s−1 tile in the HI4PI and the GALFA-H I observations, shown in the left and right panels, respectively. Bottom. Results obtained on two 3° × 20° × 0.184 km s−1 tiles in the GALFA-H I observations. In each of the panels, the left-hand side map shows the H I emission convolved with a 2D Gaussian with the same diameter of the derivative kernel. The map in the center shows the second eigenvalues of the Hessian matrix, λ−, which highlights the filamentary structure. The right-hand side map shows the orientation angle of the filamentary structure.

In the text

	Fig. C.2 Same as Fig. 13, but for an 8′ FWHM kernel.
In the text

	Fig. C.3 Same as Fig. 15, but for an 8′ FWHM kernel.
In the text

	Fig. D.1 Atomic hydrogen (H I) emission (top) and filamentary structures identified using the Hessian matrix method (bottom) toward the HVC complex H.
In the text

	Fig. D.2 Same as Fig. D.1, but for the anticenter HVC complex.
In the text