© J. D. Soler et al. 2020

1 Introduction

The diffuse neutral atomic hydrogen (HI) is the matrix within which star-forming clouds reside and the medium that takes in the energy injected by stellar winds, ionizing radiation, and supernovae (SNe, see, for example, Kulkarni & Heiles 1987; Dickey & Lockman 1990; Kalberla & Kerp 2009; Molinari et al. 2014). The observation of its distribution and dynamics provides a crucial piece of evidence to understand the cycle of energy and matter in the interstellar medium (ISM; for a review see Ferrière 2001; Draine 2011; Klessen & Glover 2016). In this paper, we present a study of the spatial distribution of the emission by HI at 21 cm, using the observations with the broadest dynamic range in spatial scales toward the Galactic plane available to date.

The structure of the HI emission in small velocity intervals has revealed a multitude of filamentary (slender or thread-like) structures, which were first identified in the earliest extended observations (see, for example, Weaver & Williams 1974; Heiles & Habing 1974; Colomb et al. 1980). Many of these filaments are curved arcs that appear to be portions of small circles on the sky. In some cases, the diameters of these arc structures change with velocity in the manner expected for expanding shells (Heiles 1979, 1984; McClure-Griffiths et al. 2002). These observations constitute clear evidence of the injection of energy into the ISM by SN explosions (see, for example, Cox & Smith 1974; McKee & Ostriker 1977; Mac Low & Klessen 2004).

The study of the HI structure has been made possible with the advent of single-dish surveys, such as the Galactic All-Sky Survey (GASS, McClure-Griffiths et al. 2009), the Effelsberg-Bonn HI Survey (EBHIS, Kerp et al. 2011), and the Galactic Arecibo L-Band Feed Array HI survey (GALFA-HI, Peek et al. 2018). Using the GALFA-HI observations of 3000 square degrees of sky at 4′ resolution in combination with the rolling Hough transform (RHT), a technique from machine vision for detecting and parameterizing linear structures, Clark et al. (2014) presented a pioneering work on the systematic analysis of HI filamentary structures. Using the EBHIS and GASS observations to produce a whole-sky HI survey with a common resolution of 30′ and in applying the unsharp mask (USM), which is another technique from machine vision that can enhance thecontrast of small-scale features while suppressing large-scale ones, Kalberla et al. (2016) presented a study of the filamentary structure of the local Galactic HI in the radial velocity range |v_LSR| < 25 km s⁻¹. Both of these studies find a significant correlation between the elongation of the HI filaments and the orientation of the local interstellar magnetic field, which 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).

Higher resolution HI observations can only be achieved by using interferometric arrays. The Galactic plane has been observed at a resolution of up to 1′ in the Canadian Galactic Plane Survey (CGPS, Taylor et al. 2003), the Southern Galactic Plane Survey (SGPS, McClure-Griffiths et al. 2005), and the Karl G. Jansky Very Large Array (VLA) Galactic Plane Survey (VGPS, Stil et al. 2006a), as well as intermediate Galactic latitudes in the Dominion Radio Astrophysical Observatory (DRAO) HI Intermediate Galactic Latitude Survey (DHIGLS, Blagrave et al. 2017). Although these surveys are limited in sensitivity compared to the single-dish observations, they have been instrumental in the study of the multiphase structure of HI, through the absorption toward continuum sources (see, for example, Strasser et al. 2007; Dickey et al. 2009) and the absorption of background HI emission by cold foreground HI (generically known as HI self-absorption, HISA; Heeschen 1955; Gibson et al. 2000).

Much of the HI is observed to be either warm neutral medium (WNM) with T ≈10⁴ K or cold neutral medium (CNM) with T ≈10² K (Heiles & Troland 2003). Detailed HISA studies of the CGPS observations reveal a population of CNM structures organized into discrete complexes that have been made visible by the velocity reversal of the Perseus arm’s spiral density wave (Gibson et al. 2005a). Using a combination of observations obtained with the Australia Telescope Compact Array and the Parkes Radio Telescope, McClure-Griffiths et al. (2006) reported a prominent network of dozens of hair-like CNM filaments aligned with the ambient magnetic field toward the Riegel-Crutcher cloud. However, there has not been a dedicated systematic study of the characteristics of these or other kinds of elongated structures in the HI emission toward the Galactic plane.

In this paper, we present a study of the elongated structures in the high-resolution observations of HI emission in the area covered by The HI/OH/recombination-line (THOR) survey of the inner Milky Way (Beuther et al. 2016; Wang et al. 2020a). We use the 40′′ -resolution HI maps obtained through the combination of the single-dish observations from the Robert C. Byrd Green Bank Telescope (GBT) with the VLA D- and C-array interferometric observations made in VGPS and THOR, respectively. We focus on the statistics of the following particular property of the elongated HI structures: its relative orientation with respect to the Galactic plane across radial velocities and Galactic longitude.

This paper is organized as follows. We introduce the observations in Sect. 2. We present the method used for the characterization of the topology of the HI emission in Sect. 3 and comment on the results of our analysis in Sect. 4. In Sect. 5, we discuss the observational effects, such as the mapping of spatial structure into the velocity space or “velocity crowding” and the HISA in the interpretation of our results. In Sect. 6, we explore the relationship between our observational results and the dynamical processes included in a set of numerical simulations of magnetohydrodynamic (MHD) turbulence taken from the From intermediate galactic scales to self-gravitating cores (FRIGG, Hennebelle 2018) and the CloudFactory (Smith et al. 2020) projects. Finally, we present our conclusions in Sect. 7. We reserve additional analyses for a set of appendices as follows. Appendix A presents details on our implementation of the Hessian analysis, such as the selection of derivative kernels, noise masking, and spatial and velocity gridding. Appendix B provides further details on the study of the filamentary structures in the GALFA-HI observations. Appendix C presents a comparison between our analysis method and the RHT and FilFinder methods. Appendices D and E expand on the analysis of the MHD simulations from FRIGG and CloudFactory, respectively.

2 Observations

2.1 Atomic hydrogen emission

For the main analysis in this paper we used the HI position-position-velocity (PPV) cube introduced in Wang et al. (2020a), which we call THOR-HI throughout this paper. It corresponds to a combination of the single-dish observations from the GBT and the VLA C- and D-array configurations in THOR and VGPS. The resulting data product covers the region of the sky defined by 14. °0 ≤ l ≤ 67. °0 and |b| ≤ 1. °25 and has an angular resolution of 40′′. The THOR-HI PPV cube, I(l, b, v), was set in Galactic coordinates and a Cartesian projection in a spatial grid with 10′′ × 10′′ pixels and 1.5 km s⁻¹ velocity channels. Ithas a typical root mean square (RMS) noise of roughly 4 K per 1.5 km s⁻¹ velocity channel. Fordetails on the calibration and imaging of these data we refer to Bihr et al. (2015), Beuther et al. (2016), and Wang et al. (2020a).

We also used the GALFA-HI observations described in Peek et al. (2018) to establish a comparison with the highest resolution single-dish HI observations. GALFA-HI has a 4′ FHWM angular resolution. At this resolution, the RMS noise of GALFA-HI and THOR-HI per 1.5 km s⁻¹ velocity channel is 126 and 358 mK, respectively. We reprojected the GALFA-HI into the THOR-HI spatial and spectral grid in two steps. First, we smoothed and regridded the data in the spectral dimension by using the tools in the spectral-cube package¹ in astropy (Astropy Collaboration 2018). Second, we projected the observations into the same spatial grid of the THOR-HI data by using the reproject package², also from astropy.

2.2 Carbon monoxide (CO) emission

We compared the HI emission observations with the ¹³CO (J = 1 → 0) observations from The Boston University-Five College Radio Astronomy Observatory Galactic Ring Survey (GRS, Jackson et al. 2006). The GRS extends over most of the area covered by THOR-HI, but it is limited to the velocity range −5 ≤ v_LSR ≤ 135 km s⁻¹. It has 46′′ angular resolution with an angular sampling of 22′′ and a velocity resolution of 0.21 km s⁻¹. Its typical root mean square (RMS) sensitivity is 0.13 K. We reprojected and regridded the GRS data using the same procedure followed with the GALFA-HI. We also used the ¹²CO (J = 1 → 0) presented in Dame et al. (2001).

2.3 Catalogs of H II regions, masers, and supernova remnants

For the study of the relation between the orientation of the HI filamentary structure and star formation, we used the catalog of HII regions from the WISE observations presented in Anderson et al. (2014). Additionally, we referred to the OH masers that are also part of the THOR observations (Beuther et al. 2019). We also used the catalogs of Galactic supernova remnants presented in Anderson et al. (2017) and Green (2019).

3 Method

Filaments, fibers, and other denominations of slender threadlike objects are common in the description and study of the ISM (see, for example, Hacar et al. 2013; André et al. 2014; Clark et al. 2019). In this work, we refer to the elongated structures in the HI emission maps across velocity channels. We characterize these structures using the Hessian matrix, a method broadly used in the study of HI (Kalberla et al. 2016) and other ISM tracers (see, for example, Polychroni et al. 2013; Planck Collaboration Int. XXXII 2016; Schisano et al. 2014).

The Hessian method uses the eigenvalues of the Hessian matrix at each pixel to classify them as filament-like or not. The Hessian matrix for a given pixel is constructed by convolving the local image patch with a set of second order Gaussian derivative filters. Different variances of the Gaussians can be used to find filaments of various widths. This approach does not imply that the identified structures are coherent objects in three-dimensional space, but rather aims to study the charateristics of elongated structures in the PPV space sampled by the HI emission.

The Hessian matrix method requires a relatively low computing time, which allows for a fast evaluation of the large set of HI observations. It also allows for the repeated testing that is required to assess the impact of the side-lobe noise, a process that would be prohibitively time consuming using more sophisticated methods (for example DisPerSE, Sousbie 2011). It also yields similar results to the RHT (Clark et al. 2014) and FilFinder (Koch & Rosolowsky 2015), as we show in Appendix C. Our implementation of this method is as follows.

3.1 The Hessian matrix method

For each position of an intensity map corresponding to v_LSR =v and a velocity interval Δv, I(l, b, v) ≡I(l, b, v±Δv), we estimatethe first and second derivatives with respect to the local coordinates (x, y) and build the Hessian matrix, $$H(x,y)\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, and H_yy ≡∂²I∕∂y² are the second-order partial derivatives and x and y refer to the Galactic coordinates (l, b) as x ≡lcosb and y ≡b, so that the x-axis is parallel to the Galactic plane. The partial spatial derivatives are obtained by convolving I(l, b) with the second derivatives of a two-dimensional Gaussian function with standard deviation w. Explicitly, we use the gaussian_filter function in the open-source software package SciPy. In the main text of this paper we present the results obtained by using a derivative kernel with 120′′ full width at half maximum (FWHM), which corresponds to three times the values of the beam FWHM in the THOR-HI observations. We also select Δv = 1.5 km s⁻¹ to match the THOR-HI spectral resolution. The results obtained with different derivative kernel sizes and Δv selections are presented in Appendix A.

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

Both eigenvalues define the local curvature of the intensity map. In particular, the minimum eigenvalue (λ₋) highlights filamentary structures or ridges, as illustrated in Fig. 1. Throughout this paper we use the word curvature to describe the topology of the scalarfield I and not the arc shape of some of the filaments.

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}\arctan \left[\frac{H_{xy}+H_{yx}}{H_{xx}-H_{yy}}\right].$$(3)

In practice, we estimated an angle θ_ij for each one of the m×n pixels in a velocity channel map, where the indexes i and j run over the pixels along the x- and y-axis, respectively. 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 Selection of the filamentary structures

We applied the Hessian matrix analysis to each velocity channel map, as illustrated in Fig. 1, but then selected the regions that are significant for the analysis by following three criteria. The first selection of filamentary structures is based on the noise properties of the HI observations. We selected regions where I(l, b) > 5σ_I, where σ_I is approximately 4 K and is estimated from the standard deviation of I in noise-dominated velocity channels.

The second selection criterion addresses the fact that the noise in the interferometric data is affected by the artifacts resulting from residual side lobes with amplitudes that vary depending on the sky position. The side lobes can introduce spurious filamentary structures in the HI emission around continuum sources, which are seen in absorption. To mitigate this effect, we mask regions of the map by using a threshold on the continuum emission noise maps introduced in Wang et al. (2018), as detailed in Appendix A. For the sake of illustration, we include the orientation of the noise features characterized with the Hessian method in the examples presented in Figs. 2 and 3, which correspond to a 2° × 2° tile toward the center of the observed Galactic longitude range.

The final selection criterion is based on the values of the eigenvalue λ₋. Following the method introduced in Planck Collaboration Int. XXXII (2016), we estimate this quantity in velocity channels dominated by noise. For that purpose, we select the five noise-dominated velocity channels and estimate the minimum value of λ₋ in each one of them. We use the median of these five λ₋ values as thethreshold value, ${\lambda }_{-}^C$. The median aims to reduce the effect of outliers, but in general the values of λ₋ in the noise-dominated channels are similar within a 10% level and this selection does not imply any loss of generality. We exclusively consider regions of each velocity channel map where λ₋ <${\lambda }_{-}^C$, given that the most filamentary structure has more negative values of λ₋.

Fig. 1 Examples of the Hessian analysis applied to three different velocity-channel maps toward a section of the Galactic plane observed by THOR. Left: intensity map. Center: map of the eigenvalue of the Hessian matrix identified as λ− in Eq. (2), which is used to characterize the filamentary structures in the intensity map. The overlaid gray map shows the filamentary structure obtained from the Hessian matrix analysis of the continuum noise maps, see Appendix A for details. Right: map of the orientation angles evaluated using Eq. (3).

3.3 Characterization of the filament orientation

Once the filamentary structures are selected, we use the angles estimated using Eq. (3) to study their orientation with respect to the Galactic plane, as illustrated in Fig. 1. The histograms presented in Fig. 2 show the variation of the preferential orientation across velocity channels. For a systematic evaluation of this variation, we use three quantities commonly used in the field of circular statistics (see, for example, Batschelet 1981): the mean resultant vector (r), the projected Rayleigh statistic (V), and the mean orientation angle $<\theta >$, which are calculated using the functions in the magnetar package³.

These statistical tools provide different and complementary information on the distribution of orientation angles, for example, r indicates if the distribution is flat or there is a preferred orientation. The values of V indicate the significance of that preferred orientation with respect to two directions of interest, namely, parallel or perpendicular to the Galactic plane, that is, 0 and 90°. Finally, $<\theta >$ indicates if the preferred orientation is different to the reference directions indicated in the computation of V, although this can be an ill-defined quantity if the angle distribution is random or multimodal (r ≈ 0). Figure 3 presents an example of the values of these three quantities across velocity channels toward the tile presented in Fig. 1.

Fig. 2 Histograms of orientation angles (θ, as defined in Eq. (3)) of the filamentary structures obtained by applying the Hessian matrix method to the velocity channels maps presented in Fig. 1. The angles θ = 0° and θ = 90° correspond to the directions parallel and perpendicular to the Galactic plane, respectively.

3.3.1 Quantifying a preferred orientation of the H I filaments

The mean resultant vector is defined as $$r=\frac{{\left({\left[{\sum }_{ij}^{m,n}{w}_{ij}\cos (2{\theta }_{ij})\right]}^2+{\left[{\sum }_{ij}{w}_{ij}\sin (2{\theta }_{ij})\right]}^2\right)}^{1/2}}{{\sum }_{ij}{w}_{ij}},$$(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. We accountfor the spatial correlations introduced by the telescope beam by choosing w_ij =${(\delta x/\Delta )}^2$, where δx is the pixel size and Δ is the diameter of the derivative kernel that we use to calculate the gradients, if it is larger than the beam size. If it is smaller than the beam size, the main scale of spatial correlation is that of the beam and consequently Δ should correspond to the beam size. We note that the statistical weights cancel out in Eq. (4), but we include them for the sake of completeness and because they are crucial to evaluate the significance of V (see, for example, Fissel et al. 2019; Soler 2019; Heyer et al. 2020).

The mean resultant vector, r, is a descriptive quantity of the angle distribution that can be interpreted as the percentage of filaments pointing in a preferential direction. If r ≈ 0, the distribution of angles is either uniform or on average does not have a well-defined mean orientation. Consequently, r ≈ 1 only if almost all of the angles are very similar.

The uncertainty on r is related to the uncertainty on the orientation angles estimated with Eq. (3). Given the selection of the filamentary structures based on the I signal-to-noise ratio (S/N), presented in Sect. 3.2, and the additional spatial smoothing introduced by the derivative kernel, the uncertainty on r is expected to be very small. This assumption is confirmed by the results of the Monte Carlo sampling introduced in Appendix A.1, which leads to the imperceptible error bars in Fig. 3.

Fig. 3 Mean spectrum and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels toward the region presented in Fig. 1. Top: mean HI intensity (I). Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: circular mean ($<\theta >$) and standard deviation (σV) of the orientation angles, represented by the line and the shaded area, respectively. The error bars correspond to the observation uncertainties and the error propagation based on the Monte Carlo sampling introduced in Appendix A.1. The dashed horizontal lines in each panel indicate the values of those quantities in the 1.4 GHz noise maps used to characterize the effect of the side lobes and continuum sources in the HI maps. The vertical gray lines correspond to the velocity channels shown in Fig. 1.

3.3.2 The significance of the preferred orientations parallel or perpendicular to the Galactic plane

We also considered the projected Rayleigh statistic (V), a test to determine whether the distribution of angles is nonuniform 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). It is defined as $$V=\frac{{\sum }_{ij}^{n,m}{w}_{ij}\cos (2{\theta }_{ij})}{\sqrt{{\sum }_{ij}^{n,m}{w}_{ij}/2}},$$(5)

which follows the same conventions introduced in Eq. (4).

In the particular case of independent and uniformly distributed angles, and for large number of samples, values of V =1.64485 and 2.5752 correspond to statistical significance levels of 5 and 0.5%, respectively (Durand & Greenwood 1958; Batschelet 1972). Values of V ≈ 1.71 and 2.87 are roughly equivalent to 2- and 3-σ confidence intervals. We note, however, that the correlation across scales in the ISM makes it very difficult to estimate the absolute statistical significance of V. In our application, we accounted for the statistical correlation brought in by the beam by introducing the statistical weights w_ij. Throughout this paper we considered values of V >3 or V < − 3 as a compelling indication of a preferred orientation of 0° or 90°, respectively.

The uncertainty on V can be estimated by assuming that each orientation angle θ_ij derived from Eq. (3) is independent and uniformly distributed, which leads to the bounded function $${\sigma }_V^2=\frac{2{\sum }_{ij}^n{w}_{ij}{\left[\cos (2{\theta }_{ij})\right]}^2-V^2}{{\sum }_{ij}{w}_{ij}},$$(6)

as described in Jow et al. (2018). In the particular case of identical statistical weights w_ij =w, ${\sigma }_V^2$ has a maximum value of w. We also estimated ${\sigma }_V^2$ by directly propagating the Monte Carlo sampling introduced in Appendix A.1. This method produces slightly higher values than those found using Eq. (6), most likely because it accounts for the correlation between the orientation angles in the map. Thus, we report the results of the latter method in the error bars in Fig. 3.

3.3.3 Mean orientations different to parallel or perpendicular to the Galactic plane

For the sake of completeness, we also computed the mean orientation angle, defined as $$<\theta >\equiv 0.5\arctan \left[\frac{{\sum }_{ij}{w}_{ij}\sin (2{\theta }_{ij})}{{\sum }_{ij}{w}_{ij}\cos (2{\theta }_{ij})}\right].$$(7)

This quantity highlights orientations that are not considered in the statistical test V, that is, different to 0° or 90°.

As in the case of r, the selection of the structures introduced in Sect. 3.2 guarantees very small uncertainties on the values $<\theta >$. There is, however, important information in the standard deviation of the angle distribution, which corresponds to the circular variance $${\sigma }_V\equiv \sqrt{-2\ln (r)}$$(8)

(see, for example Batschelet 1981), which we also report in Fig. 3.

We found that the values of r, V, and $<\theta >$ change across velocity channels, as we show in Fig. 3 for the example region presented in Fig. 1. The values of r are on average around 0.1, which means that the signal of a preferential orientation in the HI filamentary structure corresponds to roughly 10% of the selected pixels in each velocity channel map. A clear exception are the channels close to the maximum and minimum v_LSR, where the values of r are large, a result that we considered in more detail in Sect. 4.

The values of V >0 in Fig. 3 imply that the HI filamentary structures are mostly parallel to the Galactic plane. This result is further confirmed by the values of $<\theta >$, which are close to 0°. We found that roughly half of the velocity channels have values of V >3, which corresponds to a 3-σ confidence interval. The trends in V and $<\theta >$ from the HI emission are significantly different from that of the continuum noise, which we used as a proxy of the artifacts introduced by strong continuum sources. They are also independent of the average intensity in each velocity channel, which implies that they are a result of the spatial distribution of the emission and not the amount of it.

4 Results

4.1 Longitude-velocity diagrams

We applied the Hessian matrix analysis with a derivative kernel of 120′′ FWHM to the whole area covered by the THOR-HI observations by evaluating 2° × 2° nonoverlapping tiles for which we estimated r, V, and θ across velocity channels. The 2° × 2° area provides enough spatial range for the evaluation of the filamentary structures. This selection does not significantly affect our results, as further described in Appendix A.

The results of the Hessian analysis of the tiles are summarized in the longitude-velocity (lv) diagrams presented in Fig. 4. The empty (white) portions of the lv diagrams correspond to tiles and velocity channels discarded by the selection criteria introduced in Sect. 3.2. The color map for V was chosen such that the white color also corresponds to V ≈ 0, where the results of this statistical test are inconclusive.

We note that the values of r around the maximum and minimum velocities are particularly high with respect to most of the other tiles. The fact that these high values of r appear in velocity channels where $<I>$ is low indicates that the circular statistics are dominated by just a few filamentary structures that are above the I threshold. For example, if there is just one straight filament in a tile, r = 1. If there are ten straight filaments in a tile, all of them would have to be parallel to produce r = 1. Thus, the high values of r are a product of the normalization of this quantity rather than due to a lack of significance in other areas. Despite this feature, r is a powerful statistic that indicates that, on average across l and v_LSR, the preferential orientations correspond to around 17% of the HI filamentary structures.

The values of |V | ≫ 0 in the middle panel of Fig. 4, indicate that the HI filaments are preferentially oriented either parallel or perpendicular to the Galactic plane. This is further confirmed by the values of $<\theta >$, shown in the bottom panel of Fig. 4. However, it does not imply that other orientations different from 0 and 90° are not present in some tiles and velocity channels or that the distribution of the orientation of HI filaments is bimodal, as detailed in Appendix B. Figure 4 shows that the most significant trend is for the HI filaments to be parallel to the Galactic plane. The 90° orientation is significant because it is clearly grouped in specific ranges of l and v_LSR.

4.2 Comparison to previous observations

We compared the results of our study with the analysis of the GALFA-HI observations, which do not cover the entire area of the THOR survey but provide a benchmark dataset to test the possible artifacts introduced by the side lobe noise in the interferometric observations. Figure 5 presents an example of the Hessian analysis applied to the same region and velocity channels in GALFA-HI and THOR-HI. Visual inspection of these and other regions confirms that the structures selected in THOR-HI have a counterpart in GALFA-HI. In general, structures that appear monolithic in GALFA-HI are decomposed into smaller components, and filamentary structures appear to be narrower in THOR-HI. Moreover, the general orientation of the structures isvery similar, further confirming that the filamentary structures in THOR-HI are not produced by the side lobe noise.

The direct comparison of the V values in both datasets, shown in Fig. 6, provides a quantitative evaluation of the differences in the orientation of the filaments in GALFA-HI and THOR-HI. It is evident that there is a linear correlation between the values of V in both datasets, as expected from the common scales in the observations. This trend is offset from the direct one-to-one relation, for example, the tiles that present V ≈ 2 in GALFA-HI show V ≈ 4 in THOR-HI. This increase in the values of V is related tothe higher number of independent gradients that is obtained with the higher resolution THOR-HI observations, as detailed in Appendix B. The factor of two is consistent with the square root of the ratio between the statistical weights used for each dataset, which following Eq. (5) is equivalent to the ratio between the spatial scales at which we evaluated the filamentary structures, namely, the kernel size of 120′′, and the GALFA-HI 240′′ beam.

The positive values of V in Fig. 6 indicate that the general trend across l and v_LSR is roughly the same in both datasets, that is, filamentary structures in the HI emission run mostly parallel to the Galactic plane. Unfortunately, the region of the lv-diagram with the most prominent grouping of V < 0 in THOR-HI, l ≈ 27° and v_LSR ≈ 100 km s⁻¹, is not covered by the GALFA-HI observations. However, the general agreement obtained in the region covered by both surveys indicates that the global results of the Hessian analysis of THOR-HI are not the product of potential artifacts in the interferometric image reconstruction. We provide further details on the comparison between the GALFA-HI and THOR-HI in Appendix B.

Fig. 4 Longitude-velocity (lv) diagrams corresponding to the mean intensity and the results of the Hessian analysis of the velocity-channel maps in the region covered by the THOR observations. Each pixel in the diagrams represents the quantities calculated in 2° × 2° tile in a 1.5 km s−1-velocity channel as follows. Top: mean HI intensity (I). Middle top: mean resulting vector (r). Middle bottom: projected Rayleigh statistic (V). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures, which is only well-defined if r > 0. The overlaid curves correspond to selected spiral arms from the model presented in Reid et al. (2016). The vertical line indicates thecentral l of the 2° × 2° tile presented inFig. 1 and whose circular statistics are shown in Fig. 3.

4.3 General trends

The general trend revealed by the Hessian analysis of the HI emission across velocity channels is the presence of filamentary structures that tend to be preferentially parallel to the Galactic plane, as illustrated in Fig. 4. Below, we detail some observational considerations that are worth highlighting.

First, there is a significant preferential orientation of the filamentary structures in the HI. This holds despite the fact that many of the density structures in position-position-position (PPP) space are crammed into the same velocity channel in PPV, an effect called velocity crowding (see, for example, Beaumont et al. 2013). Figure 6 shows similar orientation trends in the ranges v_LSR < 0 km s⁻¹ and v_LSR > 0 km s⁻¹, although there is a larger spread in V in the range v_LSR < 0 km s⁻¹. This is presumably related to the velocity crowding due to the mapping of at least two positions into the same radial velocity, which is more pronounced at v_LSR > 0 km s⁻¹, under the assumption of circular motions (Reid et al. 2014). The velocity crowding is most likely averaging the orientation of the HI filaments and eliminating the outliers seen at v_LSR < 0 km s⁻¹. This averaging is not producing a random distribution of orientation but rather a tighter distribution around 0°. If velocity crowding had completely randomized the orientation of the HI filaments, we would see a clearer tendency toward V ≈ 0 in the v_LSR > 0 km s⁻¹ velocity range, which is clearly not the case in Fig. 6.

Given that we are using square tiles in l and b for our analysis, it is unlikely that the prevalence of horizontal structures is biased by the aspect ratio of the analysis area, although several filaments can extend from one tile to another. Additionally, given that we are using the same curvature criterion for all velocity channels in a particular l-range, as discussed in Sect. 3.2, we can guarantee that the filamentary structures are selected independently of the intensity background in a particular velocity channel.

Second, the preferential orientation of the HI filamentary structures is persistent across contiguous velocity channels, as illustrated in the example shown in Fig. 3. This implies that there is a coherence in the structure of these filamentary structures in PPV space and they are not dominated by fluctuations of the velocity field on the scale of a few kilometers per second. This is a significant fact given the statistical nature of our study, which characterizes the curvature of the emission rather than looking for coherent structures in PPV.

Third, there are exceptions to the general trend and they correspond to HI filaments that are perpendicular to the Galactic plane. These are grouped in Galactic longitude and velocity around l ≈ 27° to 30° and v_LSR ≈ 100 km s⁻¹ and l ≈ 36° and v_LSR ≈ 40 km s⁻¹. The fact that these velocity and Galactic longitude ranges correspond to those of dynamic features in the Galactic plane and active sites of star formation is further discussed in Sect. 5.

Fourth, there is a deviation of the relative orientation of the HI filaments from the general trend, V > 0, in the channels corresponding to the terminal velocities, v_t. Inside the solar circle, each line of sight has a location where it is tangent to a circular orbit about the Galactic center. At this location, known as the tangent point, the projection of Galactic rotational velocity, onto the local standard of rest (LSR) is the greatest, and the measured v_LSR is called the terminal velocity. Assuming that azimuthal streaming motions are not too large, the maximum velocity of HI emission can be equated to the terminal velocity, v_t (see, for example, McClure-Griffiths & Dickey 2016).

The velocity range around v_t is the most affected by the effect of velocity crowding, that is, a velocity channel in the PPV space corresponds to multiple positions in the PPP space. Thus, velocity crowding is a plausible explanation for the V and $<\theta >$ values found toward the maximum and minimum velocities in Fig. 4, which deviates from the general trend V ≫ 0 and $<\theta >$ ≈ 0°. However, velocity crowding does not provide a conclusive explanation for the prevalence of vertical HI filaments, V ≪ 0 and $<\theta >$ ≈ 90°, found around the maximum and minimum velocities toward 55. °0 ≲l ≲ 65. °0.

Fig. 5 GALFA-HI (top) and THOR-HI (bottom) observations and their respective filaments from the Hessian analysis toward the same region of the Galactic plane in the indicated velocity channels.

Fig. 6 Projected Rayleigh statistic V calculated from Hessian analysis of the THOR-HI and GALFA-HI observations using a 120′′ FWHM derivative kernel. Each point corresponds to a 2° × 2° × 1.5 km s−1 velocity-channel tile fully covered by both surveys. The vertical and horizontal dotted lines indicate |V | = 3, which roughly corresponds to a 3-σ confidence level in the estimation of a preferential orientation of 0° or 90°. The difference in the range of V is the result of the larger significance obtained with the higher-resolution observations.

5 Discussion

We articulate the discussion of the physical phenomena related to the results presented in Sect. 4 as follows. First, we discuss the general trend in the orientation of HI filamentary structures. which is being parallel to the Galactic plane. Then, we focus on the areas of the lv-diagram dominated by HI structures perpendicular to the Galactic plane, identified as regions of interest (ROIs) in the lv-diagram presented in Fig. 7.

Fig. 7 Projected Rayleigh statistic (V) corresponding to the orientation of the filamentary structures identified with the Hessian matrix method in the HI emission in 2° × 2° tiles and across velocity channels. The overlaid curves correspond to selected spiral arms from the model presented in Reid et al. (2016). The black contours corresponds to the 97 and 99 percentile of the 12CO emission in the observations presented in Dame et al. (2001) integrated over the Galactic longitude range |b| ≤ 1°. The green stars correspond to the HII regions in the WISE catalog of Galactic HII regions (Anderson et al. 2014). The two most prominent regions of high-mass star formation in this portion of the Galaxy, W49 and W51, are indicated by the corresponding labels. The magenta crosses correspond to the 1720-MHz OH masers in the catalog (Beuther et al. 2019), which are typically excited by shocks. The black crosses correspond to the compact HI clouds at high forbidden velocities identified in the VGPS (Stil et al. 2006b). The ellipses labeled with letters indicate the regions of interest (ROIs) discussed in Sect. 5.2.

5.1 Filamentary structures parallel to the Galactic plane

Most of the filamentary structures identified with the Hessian method in the HI emission across velocity channels are parallel to the Galactic plane. This is hardly a surprise considering that this is the most prominent axis for an edge-on disk, but it confirms that in the observed area other effects that break this anisotropy, such as the expansion of SN remnants, do not produce significant deviations from this preferential direction. Most of the deviations from this general trend appear to beconcentrated in the l and v_LSR ranges that we discuss in the next sections.

However, not all of the HI filamentary structures that are parallel to the Galactic plane appear to have the same morphology. Figure 8 shows an example of three velocity channels with high values of V, which correspond to being significantly dominated by HI filaments parallel to the Galactic plane. We note a few narrow structures that appear to be randomly oriented, which can be related to the fibers identified in Clark et al. (2014), and even some prominent vertical filaments, but the prevalent structures are parallel to the Galactic plane. These extend over several degrees in l and appear to have a width of at least 0. °5, although they are decomposed in smaller segments that correspond to the 120′′ width of the derivative kernel (see Appendix A for a discussion on the selection of the kernel size).

The predominantly positive values of V observed at v_LSR < 0 km s⁻¹ can be associated with a distance effect. Assuming circular motions, v_LSR < 0 km s⁻¹ roughly corresponds to distances larger than 10 kpc, where the apparent size of the HI disk fits in the THOR-HI b range. This would appear as a concentration of brightness that could bias r to higher values and V to more positive values, as illustrated in Fig. 3. Yet we observe vertical HI filaments at v_LSR < 0 km s⁻¹, as shown in Fig. 5.

At shorter distances, the horizontal structures in HI emission are reminiscent of the filamentary structures parallel to the Galactic plane found using unsupervised machine learning on the Gaia DR2 observations of the distance and kinematics of stars within 1 kpc (Kounkel & Covey 2019). The fact that the same work identifies that the youngest filaments (<100 Myr) are orthogonal to the Local Arm may also be a relevant hint to establish a link between the HI filaments and the process of star formation. However, most of the structures identified in Kounkel & Covey (2019) are located at intermediate Galactic latitude and are outside of the region of THOR-HI.

We show another example of a very long filamentary structure that is coherent across velocity channels in Fig. 9. Given its uniqueness, we have named it Magdalena, after the longest river in Colombia. The Magdalena (“Maggie”) filament extends across approximately 4° in Galactic longitude. With a central velocity v_LSR ≈−54 km s⁻¹, assuming circular motions, Maggie would be located at approximately 17 kpc from the Sun and 500 pc below the Galactic plane with a length that exceeds 1 kpc. The linewidths obtained with a Gaussian decomposition of the HI emission across this filament are around 13 km s⁻¹ FWHM (Syed et al. 2020a). This indicates that it is mostly likely a density structure, such as those identified in Clark et al. (2019), rather than fluctuations imprinted by the turbulent velocity field (velocity caustic, Lazarian & Pogosyan 2000; Lazarian & Yuen 2018).The physical processes that would produce such a large structure, if the kinematic distances provide an accurate estimate of its location, is still not well understood but can provide crucial insight into the dynamics of the atomic gas in the Galactic disk and halo. We present a detailed study of Maggie in an accompanying paper (Syed et al. 2020a).

Fig. 8 THOR-HI observations (top) and their respective filaments from the Hessian analysis (bottom) toward l and vLSR ranges dominated by filamentary structures parallel to the Galactic plane.

Fig. 9 Magdalena HI filament. Top: atomic hydrogen (HI) emission in three velocity channels. Middle: filaments obtained using the Hessian technique toward the corresponding channels. Bottom: spectra toward the positions indicated by the crosses in the top panel. The vertical lines indicate the vLSR that corresponds to the emission shown in the other two panels.

5.2 Regions of interest (ROIs) in the lv-diagram

5.2.1 ROI A: l ≈ 27° and v_LSR ≈ 100 km s⁻¹

We found that the most prominent association of tiles with HI filaments perpendicular to the Galactic plane is located around l ≈ 27° and v_LSR ≈ 100 km s⁻¹, marked as ROI A in Fig. 7. This position has been previously singled out in the study of the HI emission toward the Galacticplane, particularly in observations that indicate the presence of voids several kiloparsecs in size centered approximately on the Galactic center, both above and below the Galactic plane (Lockman & McClure-Griffiths 2016). These voids, whichappear to map the boundaries of the Galactic nuclear wind, are evident in the sharp transition at galactocentric radius of around 2.4 kpc from the extended neutral gas layer characteristic of much of the Galactic disk to a thin Gaussian layer with 125 pc FWHM.

It is plausible that this reported thinning of the HI disk is related to the change in the preferential orientation of the HI filaments at l ≈ 27°. Visual inspection of a few velocity channels within the l and v_LSR ranges of the ROI A, shown in Fig. 10, indicate that there is indeed a thinning of the Galactic plane for l ≲ 22°. But the HIfilamentary structures in the range l ≲ 22° and v_LSR ≳ 70 km s⁻¹ are mostly parallel to the Galactic plane, as shown by the positive V values in that range in Fig. 7.

It is also possible that the filaments perpendicular to the Galactic plane are prominent just because of a geometrical effect. The tangent point of the Scutum arm is close to l ≈ 30° (Reid et al. 2016), thus, the filaments that are parallel to this spiral arm will appear shortened in the plane of the sky and their orientation will be either random or dominated by the passage of the spiral arm. However, if that was the case, we should also see a significant variation around l ≈ 50°, toward the tangent of the Sagittarius arm. Such an effect was seen in the strong excess of Faraday rotation toward that position (Shanahan et al. 2019), but this is not seen in the orientation of the HI filament, as illustrated in Fig. 7.

A different effect that singles out the position of the ROI A in l and v_LSR is the Galactic feature identified as the Molecular Ring around 4 kpc from the Galactic center (Cohen & Thaddeus 1977; Roman-Duval et al. 2010). This structure forms two branches with tangent points at l ≈ 25° and 30° close to the positive terminal velocity curve, which can be reproduced in numerical simulations that include the gravitational potential of the Galactic bar, thus indicating that the presence of this feature does not depend on local ISM heating and cooling processes (Fux 1999; Rodriguez-Fernandez & Combes 2008; Sormani et al. 2015). It is currently unclear whether this structure is really a ring, or emission from many spiral arms crowded together, as suggested by the numerical simulations. The coincidence between the tangent points of the Galactic Ring and the ROI A suggests that there is an imprint of the Galactic dynamics in the vertical structure of the atomic hydrogen. The fact that a similar effect is not observed in the tangent of the Sagittarius arm implies that this is not a geometrical effect or the result of the passage of a spiral arm.

One observational fact that distinguishes ROI A is the large density of HII regions and sources of RRLs, also shown in Fig. 7. There is observational evidence of a significant burst of star formation near the tangent of the Scutum arm, in the form of a large density of protoclusters around W43 (Motte et al. 2003) and multiple red supergiant clusters (Figer et al. 2006; Davies et al. 2007; Alexander et al. 2009). This star formation burst can be associated to an enhancement of SN feedback forming a “Galactic fountain” (Shapiro & Field 1976; Bregman 1980; Fraternali 2017; Kim & Ostriker 2018), whose remants shape the vertical HI filaments.

The relation between the star-formation and the vertical HI filaments is further suggested by the observed asymmetry in the density of HI clouds in the lower Galactic halo toward the tangent points in the first and the fourth quadrants of the Milky Way reported in Ford et al. (2010). There, the authors show that there are three times more disk-halo clouds in the region 16. °9 <l < 35. °3 than in the region 324. °7 <l < 343. °1 and their scale height is twice as large. Our results indicate that the potential origin of this population of disk-halo clouds, which are found in the range |b| < 20°, also has a significant effect in the structure of the HI gas at lower Galactic latitude. Given the symmetry in the Galactic dynamics of the tangent point in the first and the fourth quadrant, the difference between these two regions appear to be linked to the amount of star formation and SN feedback. This hypothesis is reinforced by the prevalence of vertical HI filamentsin ROI B, where the effect of the Galactic dynamics is less evident.

Fig. 10 THOR-HI (top) observations and their respective filaments from the Hessian analysis (bottom) toward the Galactic longitude range of the region of interest (ROI) A, as defined in Fig. 7. The yellow circles correspond to the positions and sizes of the supernova remnants in the Green (2019) catalog. The cyan circles correspond to the positions and sizes of the HII regions in the Anderson et al. (2014) catalog with vLSR within ±10 km s−1 of the central velocity indicated in the figure. The magenta crosses mark the position of the THOR OH 1720-MHz masers in that velocity range.

Fig. 11 Same as Fig. 10, but for the ROI B.

5.2.2 ROI B: l ≈ 36° and v_LSR ≈ 40 km s⁻¹

The second region where we found a significant association of HI filaments perpendicular to the Galactic plane is located around l ≈ 36° and v_LSR ≈ 40 km s⁻¹. Figure 11 shows that the filamentary structure in the HI emission toward ROI B forms an intricate network where many orientations are represented. However, the values of V ≪0 and $<\theta >$ ≈ 90° indicate that statistically, these filaments are preferentially perpendicular to the Galactic plane. We also can identify some prominent vertical filaments around l ≈ 38°.

The region of interest B coincides with a large density of HII regions, as illustrated in Fig. 7. However, this is not enough to establish a causality between the presence of HII regions and a preferential orientation in the HI filamentary structure. As a matter of fact, we did not find a prevalence of vertical HI filaments toward W49 (l = 43. °2, b = 0. °0, v_LSR ≈ 11 km s⁻¹) and W51 (l = 49. °4, b =−0. °3, v_LSR ≈ 60 km s⁻¹), two of the most prominent regions of high-mass star formation in the Galaxy (Urquhart et al. 2014). Both W49 and W51 are relatively young and have a small ratio of SN remnants relative to the number of HII regions (Anderson et al. 2017). The absence of a preferred vertical orientation in the HI filament toward these two regions suggests that the effects observed toward ROI A and B are not associated with very young star-forming regions, where most of the massive stars are still on the main sequence.

A plausible explanation for the prevalence of vertical filaments in both ROI A and B is the combined effect of multiple SNe. In this case, however, it does not correspond to the walls of expanding bubbles, such as those identified in Heiles (1984) or McClure-Griffiths et al. (2002), but rather the preferential orientation produced by multiple generations of bubbles, which depressurize when they reach the scale height and produce structures that are perpendicular to the Galactic plane. There are at least six 0. °5-scale supernova remnants (SNRs) toward ROI B, including Westerhout 44 (W44, Westerhout 1958), as shown in Fig. 11. There is also a strong concentration of OH 1720-MHz masers toward ROI B, which are typically excited by shocks and found toward either star-forming regions or SNRs (see, for example, Frail et al. 1994; Green et al. 1997; Wardle & Yusef-Zadeh 2002). Most of these OH 1720-MHz masers appear to be associated with W44 and would not trace the effect of older SNe that may be responsible for the vertical HI filamentary structures that are seen, for example, around l ≈ 38. °0 and l ≈ 38. °8 in Fig.7. Thus, we resort to numerical experiments to evaluate this effect statistically in Sect. 6.

5.2.3 ROI C: Riegel-Crutcher cloud

Among the prominent features in the orientation of the HI filamentary structure presented in Fig. 7, one that is of particular interest is that found around 18° <l < 26° and 0 <v_LSR < 10 km s⁻¹. This location corresponds to the Riegel-Crutcher (RC) cloud, a local (d = 125 ± 25 pc, Crutcher & Lien 1984) CNM structure that extends approximately 40° in Galactic longitude and 10° in latitude (Riegel & Jennings 1969; Crutcher & Riegel 1974). Figure 12 reveals features in the emission around l ≈ 22° that are reminiscent of the strands found in the most studied portion of the RC cloud (|l| < 5° and |b| < 5°, McClure-Griffiths et al. 2006; Dénes et al. 2018), but with a different orientation. The analysis of the filamentary structure presented in this paper servesas an alternative way of defining the RC cloud, through the contrast in the orientation of the filamentary structures with respect to that in the surrounding velocity channels.

Many of the structures in the RC cloud are seen as shadows against an emission background. It is common at low Galactic latitudes that cold foreground HI clouds absorb the emission from the HI gas behind. This effect is often called HI self-absorption (HISA), although it is not self-absorption in the standard radiative transfer sense, because the absorbing cloud may be spatially distant from the background HI emission, but share a common radial velocity (Knapp 1974; Gibson et al. 2000).

HISAs are often identified by a narrow depression in the HI spectrum and as a coherent shadow in the emission maps. The systematic identification of HISA features and the evaluation of the completeness of the census of cold HI that it provides is a complex process, as described in Gibson et al. (2005b) and Kavars et al. (2005) or Wang et al. (2020b) and Syed et al. (2020b) in the particular case of THOR-HI. The HI filament orientations provide a complementary method to identify HISA structures, by quantifying the difference in orientation of the HISA with respect to that of the intensity background. Figure 13 shows an example of the HISA identification in the spectrum and in the orientation of the filamentary structure toward the molecular cloud GRSMC 45.6+0.3 studied in Jackson et al. (2002). The spectra and the values of V across velocity channels toward this region indicate that the HISA feature is evident in both the intensity and the morphology of the HI emission.

Fig. 12 Same as Fig. 10 for the ROI C.

5.2.4 ROI D: terminal velocities

The final exception to the general trend of HI filaments parallel to the Galactic plane is found at the maximum and minimum v_LSR with a large prominence in multiple velocity channels at l > 56°. Part of the emission at the extremes of the radial velocities has been identified as a separate HI component from that in the Galactic disk and belongs to the Galactic Halo (Shane 1971; Lockman 2002). This component is 1–2 kpc above the Galactic plane, and it is usually called “extraplanar” HI gas, which avoids potential confusion with the Galactic Halo as intended in, for example, cosmological simulations (see, for example, Fraternali & Binney 2006, 2008; Fraternali 2017).

Figure 14 shows that the structure of the emission at − 83 and 46 km s⁻¹ is clearly vertical and different from that at intermediate velocities. It is also different from the structure in ROIs A and B, shown in Figs. 10 and 11, where the vertical HI structures are part of an intricate network of filamentary structures. In the ROI D, most of the vertical filamentary structures in HI emission are clearly separated from other structures and notably extend to high b, further suggesting that they are HI clouds from the halo. At least one of the tiles with preferentially vertical filaments corresponds to a cloud above the maximum velocity allowed by Galactic rotation, that at l = 60. °8 and v_LSR = 58.0 km s⁻¹ identified in the VGPS (Stil et al. 2006b).

Marasco & Fraternali (2011) studied the extraplanar HI gas of the Milky Way, and found that it is associated with SN feedback and mainly consists of gas that is falling back to the MW after being ejected by SNe. The reason why the extraplanar HI is more conspicuous falling down than going up is because it cools while settling down, whereas it is hotter and so less visible when going up, as can be seen in the numerical simulations presented in Kim & Ostriker (2015) and Girichidis et al. (2016). So potentially, the vertical filaments around ROI D are falling back after being ejected by SNe instead of material going up, as is possibly the case for ROI A and ROI B. The time and spatial delay between SN explosion and gas falling back might also explain why we do not find that every region associated with SNe has vertical filaments.

Fig. 13 Selected example of a prominent HI self-absorption (HISA) structure. Top: atomic hydrogen (HI) emission in three velocity channels. Middle: projected Rayleigh statistic (V, Eq. (5)) corresponding to the orientation of the filamentary structures in the HI emission in the corresponding channels. Bottom: spectra toward the positions indicated by the crosses in the top panel. The vertical lines indicate the vLSR that corresponds to the emission shown in the other two panels.

5.3 Relation to the structure of the molecular gas

Filamentary structures designated giant molecular filaments (GMF) have been previously identified toward the Galactic plane in the emission from molecular species, such as ¹³CO (Goodman et al. 2014; Ragan et al. 2014; Wang et al. 2015, 2020c; Zucker et al. 2015; Abreu-Vicente et al. 2017). To establish a link between the HI filaments and the GMFs, without detailing individual objects, we also applied the Hessian analysis to the ¹³CO (J = 1 → 0) emission observations in the GRS survey. Following the selection criteria presented in Sect. 3.2, we estimated the orientation of the filamentary structures identified using the Hessian method in the GRS observations projected into the same spectral axis of THOR-HI.

Figure 15 shows a comparison of the values of V in THOR-HI and GRS observations. In agreement with the GMF compilation presented in Zucker et al. (2018), we find that most of the ¹³CO filamentary structures are parallel to the Galactic plane. However, we found no evident correlation between the orientation of the filamentary structures in the two tracers.

There can be several reasons for this general lack of correlation between the HI and the ¹³CO filamentary structures. First, the linewidths of the ¹³CO emission arenarrower than those of the HI and it is possible that we are washing away part of the orientation of the filaments byprojecting both data sets into the same spectral grid. This effect, however, may not be dominant. The Gaussian decomposition of the GRS presented in Riener et al. (2020) indicates that the mean velocity dispersion (σ_v) is approximately 0.6 km s⁻¹ and the interquartile range is 0.68 <σ_v < 1.89 km s⁻¹, thus regridding the data to 1.5 km s⁻¹ resolution would not completely alter the morphology of most of the emission. The fact that we found a large number of emission tiles with V ≫ 0 indicates that there is a preferential orientation of the filamentary structures in ¹³CO, parallel to the Galactic plane, and this preferential orientation is not washed away by the integration over a broad spectral range.

Second, although there is a morphological correlation of the HI and the ¹³CO, as quantified in Soler et al. (2019), the filamentary structure in the HI emission is not necessarily related to that of the ¹³CO. In general, the much larger filling factor of the HI makes it unlikely that most of its structure is related to that of the less-filling molecular gas. Moreover, when evaluating comparable scales, the HI and the ¹³CO can appear completely decoupled (Beuther et al. 2020). This does not discard the local correlation between the morphology of the two tracers toward filamentary structures, as reported in Wang et al. (2020b) and Syed et al. (2020b), but it shows that the orientation of the filamentary structures are generally different.

Fig. 14 THOR-HI emission for two velocity channels in the region of interest (ROI) D (left and right) and a channel toward the same l range but dominated by filaments parallel to the Galactic plane (middle).

Fig. 15 Comparison between the values of the projected Rayleigh statistic V calculated using the THOR and 13CO GRS observations in the indicated Galactic-longitude ranges. Both are calculated using a derivative kernel with 120′′ FWHM. Each point in the scatter plot corresponds to the values estimated using the Hessian analysis in the same 2° × 2° tile and 1.5-km s−1 velocity channel in both datasets.

6 Comparison to MHD simulations

Among the plethora of physical processes that can be responsible for the preferential orientation of HI filamentary structures reported in this paper, we explored the effect of SNe, magnetic fields, and galactic rotation in the multiphase medium in two families of state-of-the-art and readily-available numerical simulations. First, we considered the simulations in the “From intermediate galactic scales to self-gravitating cores” (FRIGG) project, introduced in Hennebelle (2018), which we used to explore the effects of SN feedback and magnetization in the orientation of the HI structures. Second, we considered the Cloud Factory simulation suite, which is introduced in Smith et al. (2020) and is designed to study SN feedback effects while also including the larger-scale galactic context in the form of the galactic potential, the differential rotation of the disk, and the influence of spiral arms.

Fig. 16 Synthetic observations of the HI emission from the FRIGG simulations with a standard (3 μG) and very high (12 μG) initial magnetizations, shown in the left and right panels, respectively. The colors represent the emission in the three indicated radial velocity bins with the same average emission.

6.1 Supernovae and magnetic fields

6.1.1 Initial conditions

The FRIGG simulations used the RAMSES code and take place in a stratified 1-kpc side box with SN explosions and MHD turbulence. This is a standard configuration that can be found in other works, such as, de Avillez & Breitschwerdt (2007), Walch et al. (2015), and Kim & Ostriker (2017). It includes the cooling and heating processes relevant to the ISM, which produce a multiphase medium. They reproduce the vertical structure of the Galactic disk, which results from the dynamical equilibrium between the energy injected by the SNe and the gravitational potential of the disk. In particular, we used the set of simulations described in Iffrig & Hennebelle (2017), which have different levels of magnetization that we analyze to assess the role of the magnetic field in the orientation and characteristics of the HI filaments. These simulations have a resolution that is limited to 1.95 pc, however, this is enough for a first glimpse at the orientation of the structures formed under the general 1-kpc scale initial conditions.

Iffrig & Hennebelle (2017) report that the efficiency of the SNe in driving the turbulence in the disk is rather low, on the order of 1.5%, and strong magnetic fields increase it by a factor of between two and three. It also reports a significant difference introduced by magnetization in the filamentary structures perpendicular to the Galactic plane, illustrated in their Fig. 1. To quantify the differences introduced by the magnetization in the morphology of the emission from atomic hydrogen, we compared one snapshot in the simulation with “standard” magnetization, initial magnetic field B₀ of about 3 μG, and one with “very high" magnetization, B₀ ≈ 12 μG. Both simulations have an initial particle density of n₀ = 1.5 cm⁻¹. The initial magnetic field strengths are chosen around the median magnetic field strength 6.0 ± 1.8 μG observed in the CNM (Heiles & Troland 2005). We selected snapshots at 75 and 81 Myr for the standard and very high magnetization cases, respectively, both of which are available in the simulation database Galactica⁴. This selection guarantees that the simulations have reached a quasi-steady state and does not affect the reported results. We present the synthetic observations in Fig. 16 and describe details of their construction in Appendix D.

6.1.2 Hessian analysis results

The prevalence of longer filamentary structure with higher magnetization has been reported in previous studies (Hennebelle 2013; Seifried & Walch 2015; Soler & Hennebelle 2017). It is related to the effect of strain, which means that these structures simply result from the stretch induced by turbulence, and the confinement by the Lorentz force, which therefore leads them to survive longer in magnetized flows. However, their orientation in this kind of numerical setup has not been systematically characterized.

We applied the Hessian analysis with a 120′′ FWHM derivative kernel to the synthetic HI emission from the FRIGG simulations. We summarized the results in Fig. 17. Our first significant finding is that the standard magnetization case reproduces some of the filamentary structures parallel to the Galactic plane that are broadly found in the observations, but these do not show the same significance in terms of the values of V. This means that the initial conditions in the standard case reproduce some of the horizontal filaments just with the anisotropy introduced by the vertical gravitational potential. However, this setup is missing the Galactic dynamics that are the most likely source of the stretching of structures in the direction of the plane.

Our second finding is that the magnetization does not constrain the filamentary structures to the plane, but rather maintains the coherence of the structures that are blown in the vertical direction by the SNe, as show in the high magnetization case. When the bubble blown by a SN reaches the scale height and depressurizes, the magnetic field maintains the coherence in its walls, which wouldfragment if the magnetic field were weaker. Subsequently, the gas layer consisting of ejected clouds falls back on the plane and isstretched along the field lines.

The aforementioned results suggest that the magnetic field may play a significant role in the prevalence of vertical structures in the regions indicated in Fig. 7. Filamentary structures have been observed in radio continuum toward the Galactic center and their radio polarization angles indicate that these structures follow their local magnetic field (see, for example, Morris & Serabyn 1996; Yusef-Zadeh et al. 2004). Studies of radio polarization at higher Galactic latitude indicate a correspondence between the depolarization canals and the HI filamentary structures, which suggests that the filamentary structures share the orientation of the magnetic field (Kalberla et al. 2017). Thus, it is tempting to think that a similar effect can be responsible for the orientation of the vertical filaments in the THOR-HI observations.

Fig. 17 Synthetic HI emission intensity and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels in the synthetic observations presented in Fig. 16. Top: HI emission averaged over the whole map and toward the indicated position. Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

6.2 Supernovae and Galactic dynamics

6.2.1 Initial conditions

We considered the effect of the Galactic dynamics by using the CloudFactory simulations, presented in Smith et al. (2020). These simulations, ran using the AREPO code, consist of a gas disk inspired by the Milky Way gas disk model of McMillan (2017) and focus on the region between 4 and 12 kpc in galactocentric radius. The simulations start with a density distribution of atomic hydrogen that declines exponentially at large radii. Molecular hydrogen forms self-consistently as the gas disk evolves. We used a 1-kpc-side box region within the large-scale setup, with a mass resolution of 10M_⊙, and gas self-gravity.

We compared two simulation setups. In one, SN were placed randomly in the galactic disk at a fixed rate of 1 per 50 yr, chosen to match the value appropriate for the Milky Way. The other setup combined a random SN component with a much smaller rate of 1 per 300 yr, designed to represent the effect of type Ia SNe, with a clustered SN component whose rate and location were directly tied to the rate and location of star formation in the simulation. Following the terminology of Smith et al. (2020), we refer to these simulations as potential-dominated and feedback-dominated, respectively.

Smith et al. (2020) reports the alignment of filamentary structures in the disk by spiral arms and the effect of differential rotation. The authors also note that clustered SN feedback randomize the orientation of filaments and produce molecular cloud complexes with fewer star-forming cores. To quantify these effects in the HI emission, we studied the synthetic observation of one snapshot in the potential- and feedback-dominated simulations. We present the synthetic observations in Fig. 18 and describe details on their construction in Appendix E.

The simulations of the whole disk were performed at a lower resolution with cell masses of 1000 M_⊙ for 150 Myr until the simulated ISM reached a steady state. This initial setup step was followed by a progressive resolution increase in a corotating 1-kpc sized box. The snapshot used for analysis is chosen when the gas within the high resolution box has had time to perform two spiral arm passages, enough for the simulated ISM to adapt to the new resolution and physics. Unfortunately, the simulations did not progress much further in time such that an analysis of the simulation at different ISM conditions is not possible at the moment.

Fig. 18 Synthetic observations of the HI emission from the potential and feedback dominated CloudFactory simulations, presented in the left and right panels, respectively. The colors represent the emission in the three indicated radial velocity bins, each with the same average emission.

6.2.2 Hessian analysis results

We applied the Hessian analysis with a 120′′ FWHM derivative kernel to the synthetic HI emission from the CloudFactory simulations. The results are summarized in Fig. 19. The most significant outcome of this study is that the Galactic dynamics in these simulations naturally produce filamentary structures parallel to the Galactic plane across velocity channels, which are comparable to those found in the GALFA-HI and THOR-HI observations. These filamentary structures are coherent across several velocity channels and correspond to overdensities that are clearly identifiable in the density cubes from the simulation, thus, they are not exclusively the product of fluctuations in the velocity field.

The clustered SNe in the feedback-dominated simulation produce structures that resemble clumpy filaments in the synthetic HI PPV cube, as shown in Fig. 18. These structures do not show a significant preferential orientation, as illustrated by the values of |V | ≲ 0 in the corresponding panel of Fig. 19. This confirms and quantifies the randomization of the structures described in Smith et al. (2020).

Both the potential-dominated and feedback-dominated cases considered in this numerical experiment correspond to extreme cases. The fact that the potential-dominated simulation does not produce a significant number of vertical HI filaments suggests that these are most likely related to the effect of clustered SNe. The fact that the SN feedback erases all the anisotropy introduced by the dynamics in the direction of the Galactic plane indicates that the prevalence of vertical filaments is an indication that, at least in a few specific locations, SN feedback has a dominant effect in the structure of the ISM. Therefore, the observation of this vertical HI filaments is a promising path toward quantifying the effect of SN feedback in the Galactic plane.

7 Conclusions

We presented a study of the filamentary structure in the maps of the HI emission toward inner Galaxy using the 40′′ -resolution observations in the THOR survey. We identified filamentary structures in individual velocity channels using the Hessian matrix method and characterized their orientation using tools from circular statistics. We analyzed the emission maps in 2° × 2° tiles in 1.5-km s⁻¹ velocity channels to report the general trends in orientation across Galactic longitude and radial velocity.

We found that the majority of the filamentary structures are aligned with the Galactic plane. This trend is in general persistent across velocity channels. A comparison with the numerical simulation of the Galactic dynamics and chemistry in the CloudFactory project indicate that elongated and nonself-gravitating structures naturally arise from the galactic dynamics and are identified in the emission from atomic hydrogen.

Two significant exceptions to this general trend of HI filaments being parallel to the Galactic plane are grouped around l ≈ 37° and v_LSR ≈ 50 km s⁻¹ and toward l ≈ 36° and v_LSR ≈ 40 km s⁻¹. They correspond to HI filaments that are mostly perpendicular to the Galactic plane. The first location corresponds to the tangent point of the Scutum arm and the terminal velocities of the Molecular Ring, where there is a significant accumulation of HII regions. The second position also shows a significant accumulation of HII regions and supernova remnants.

Comparison with numerical simulations in the CloudFactory and FRIGG projects indicate that the prevalence of filamentary structures perpendicular to the Galactic plane can be the result of the combined effect of SN feedback and magnetic fields. These structures do not correspond to the relatively young (<10 Myrs) structures that can be identified as shells in the HI emission, but rather to the cumulative effect of older SNe that lift material and magnetic fields in the vertical direction. Thus, their prevalence in the indicated regions is the potential signature of an even earlier history of star formation and stellar feedback in the current structure of the atomic gas in the Galactic plane. Our observational results motivate the detailed study of more snapshots and different numerical simulations to firmly establish the relation between these relics in the atomic gas and the record of star formation.

Another exception to the general trend of HI filaments being parallel to the Galactic plane is found around the positive and negative terminal velocities. Comparison with previous observations suggests that these structures may correspond to extraplanar HI clouds between the disk and the halo of the Milky Way. A global explanation for the vertical HI filaments is that the combined effect of multiple SNe creates a layer of gas consisting of ejected clouds, some of which are falling back on the plane. Such clouds would naturally tend to be vertically elongated and coherent due to the effect of the magnetic fields. Galactic dynamics may be responsible for creating the observed vertical filaments only in an indirect way: it helps bringing thegas together, creating favorable conditions for SNe to cluster together, explode, and create the vertical structure.

The statistical nature of our study unveils general trends in the structure of the atomic gas in the Galaxy and motivates additional high-resolution observations of the HI emission in other regions of the Galaxy. Further studies of the nature and the origin of the HI filamentary structures call for the identification of other relevant characteristics, such as their width and length, as well as the physical properties that can be derived using other complementary ISM tracers. Our results demonstrate that measuring the orientation of filamentary structures in the Galactic plane is a promising tool to reveal the imprint of the Galactic dynamics, stellar feedback, and magnetic fields in the observed structure of the Milky Way and other galaxies.

Fig. 19 Synthetic HI emission intensity and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels in the synthetic observations presented in Fig. 18. Top: HI emission averaged over the whole map and toward the indicated position. Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

Acknowledgements

J.D.S., H.B., and Y.W. acknowledge funding from the European Research Council under the Horizon 2020 Framework Program via the ERC Consolidator Grant CSF-648505. H.B., J.S., S.C.O.G., and R.K. acknowledge support from the Deutsche Forschungsgemeinschaft via SFB 881, “The Milky Way System” (subprojects A1, B1, B2, and B8). S.C.O.G. and R.K. also acknowledge support from the Deutsche Forschungsgemeinschaft via Germany’s Excellence Strategy EXC-2181/1-390900948 (the Heidelberg STRUCTURES Cluster of Excellence). R.K. also acknowledges funding from the European Research Council via the ERC Synergy Grant ECOGAL (grant 855130) and the access to the data storage service SDS@hd and computing services bwHPC supported by the Ministry of Science, Research and the Arts Baden-Württemberg (MWK) and the German Research Foundation (DFG) through grants INST 35/1314-1 FUGG as well as INST 35/1134-1 FUGG. N.S. acknowledges support by the Agence Nationale de la Recherche (ANR) and the DFG through the project “GENESIS” (ANR-16-CE92-0035-01/DFG1591/2-1). We thank the anonymous referee for the thorough review. We highly appreciate the comments, which significantly contributed to improving the quality of this paper. J.D.S. thanks the following people who helped with their encouragement and conversation: Peter G. Martin, Marc-Antoine Miville-Deschênes, Josh Peek, Jonathan Henshaw, Steve Longmore, Paul Goldsmith, and Hans-Walter Rix. 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. 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: matplotlib (Hunter 2007), astropy (Greenfield et al. 2013), magnetar (Soler 2020), RHT (Clark et al. 2020).

Appendix A: Details of the Hessian method implementation

Fig. A.1 Hessian analysis applied to a velocity-channel map dominated by noise toward the same region presented in Fig. 1. Left: HI intensity map. Center: map of the eigenvalue of the Hessian matrix identified as λ− in Eq. (2), which is used to characterize the filamentary structures in the intensity map. The overlaid gray map shows the filamentary structure obtained from the Hessian matrix analysis of the continuum noise maps, see Appendix A for details. Right: map of the orientation angles evaluated using Eq. (3).

Fig. A.2 Same as Fig. A.1, but for the noise map of the continuum emission at 1.42 GHz (${\sigma }_{I_{1.42}}$), which we use as a proxy of the structure introduced by the side-lobes in the imaging of the HI emission.

A.1 Accounting for the effect of noise

A.1.1 Error propagation through Monte Carlo sampling

We estimated the effect of the noise in the circular statistics introduced in Sect. 3.3 by using Monte Carlo sampling. We generated 10⁴ draws I_i(l, b, v) from a Gaussian probability distribution function described by the mean value I(l, b, v) and the variance (σ_I, where I(l, b, v) are the observed intensities in the PPV cubes and σ_I is estimated from noise-dominated velocity channels. As a zeroth-order approximation, we assumed that the noise is constant in each velocity-channel map and across the velocity channels. This assumption is well justified by the uniform uv coverage and total integration times throughout the THOR observations (Beuther et al. 2016; Wang et al. 2020a).

From each of the generated I_i(l, b, v) we calculated r, V, and $<\theta >$. Assuming that the errors ineach of these quantities have a Gaussian distribution, we reported the mean value and the standard deviation for each on of them. Given the selection by I signal-to-noise introduced in Sect. 3.2, the estimated uncertainties obtained through this method are very small.

A.1.2 Curvature threshold

We calculated the curvature threshold values ${\lambda }_{-}^C$ introduced in Sect. 3.2 by considering a velocity channel with very low HI in each 2° × 2° region. For that channel, we estimated the mean intensity and the maximum curvature λ₋, as defined in Eq. (2), which we assign to be ${\lambda }_{-}^C$ for that particular position. Figure A.1 shows an example of this procedure for the region presented in Fig. 1. Given that the selected velocity channel is dominated by noise, the filamentary structures cover the whole maps and present values of λ₋ close to 0, that is, very low curvature. The orientation of these filaments is rather homogeneous and it is not indicative of the spatial distribution of the noise.

A.1.3 Spatial structure of the noise

To characterize the spatial structure of the noise, we use the noise map of the continuum emission at 1.42 GHz, ${\sigma }_{I_{1.42}}$, presented in Fig. A.2. The noise map of the continuum emission serves as a proxy for the linear structures that can be potentially introduced in the HI maps by continuum sources in absorption. By masking the HI emission based on the filamentary structures found in ${\sigma }_{I_{1.42}}$, we exclude the strong continuum sources and the side-lobe features around them from the Hessian analysis. We note that in general the orientation of the filamentary structure in ${\sigma }_{I_{1.42}}$ rarely corresponds to that found in the HI emission, as we show in the example presented in Fig. 3. But that is not necessarily the case in all the 2° × 2° regions, which motivates our masking scheme.

Fig. A.3 Same as Fig. 1, but for derivative kernels of 80′′ (top) and 160′′ FWHM (bottom), which correspond to two and four times the native resolution of the THOR-HI observations, respectively.

Fig. A.4 Longitude-velocity diagram of the projected Rayleigh statistic V estimated using the indicated derivative kernel sizes.

A.2 Derivative kernel size

In the main body of this paper, we have shown the results for a particular selection of the derivative kernel size with 120′′ FWHM. This selection, which sets the scale at which the filamentary structures are evaluated, was selected empirically by reaching a compromise between the spurious filamentary structures obtained with a kernel size close to the resolution of the observations (40′′) and the loss of information resulting from using a very coarse one.

Fig. A.5 Longitude-velocity diagram of the mean orientation angle $<\theta >$ of the HI filamentary structures estimated using the Hessian method and the indicated derivative kernel sizes.

Figure A.3 shows an example of two different kernels sizes applied to the same velocity channel map toward the region presented in Fig. 1. The kernel that corresponds to twice the native resolution of the observations (80′′ FWHM) highlights a large number of narrow filamentary structures, but it is very sensitive to the features of the HI imaging. Some of these features are the result of the artifacts from the interferometric data and are common when considering next-neighbor derivatives. The coarser 160′′ FWHM kernel, shows a much clearer contrast in terms of λ₋, but washes away some of the structures in the intensity map. This selection may need further investigation for the study of other filament properties, such as the width, but it does not critically affect the results of the orientation study.

Fig. A.6 Longitude-velocity diagram showing the percentage of each velocity-channel map that is selected in the Hessian analysis. Each pixel corresponds to a 2° × 2° tile.

The lv diagrams of the V and $<\theta >$ obtained with the 80′′ and 160′′ FWHM derivative kernels are shown in Figs. A.4 and A.5. The results are in general similar to those presented in Fig. 4. However, it is evident that the 80′′ FWHM kernel appears noisier in both V and θ, most likely related to the effect of the spatial features shown in Fig. A.3. The 160′′ FWHM kernel shows lower maximum levels of V, but the main regions of interest in Fig. 7 are still clearly identifiable in the lv diagrams. A.3:

A.3 Filament selection

In addition to the I signal-to-noise considerations, the selection of the filaments is based on the curvature threshold ${\lambda }_{-}^C$, which is calculated independently for each derivative kernel size. The area of each tile that is covered by the selected filaments changes with the derivative kernel selection. For example, if we picked a very large derivative kernel, we would obtain a very large filament, which is the Galactic plane itself, that covers most of the tile. As the kernel size is reduced one expects that the area of the tile covered by filaments is also reduced, as we are selecting smaller scales and sampling narrower filaments. This expectation is confirmed for the THOR-HI observations in Fig. A.6.

The fact that the percentage of the tile covered by filaments is in general changing with the size of the derivative kernel is an indication that most of them are not resolved at the spatial scales considered in Fig. A.6. If there was a particular scale at which the filaments are resolved, the map coverage would remain constant for derivative kernel sizes below that scale. This is due to the fact that the curvature of a resolved structure, quantified by the value of λ₋, should be independent of the sampled scale. This motivates further observations at higher resolution to study the width of the structures.

The percentage of a tile that is identified as filaments changes across l and v_LSR and roughly follows the same distribution of the mean intensity, shown in Fig. 4. This suggests that more filamentary structures are found in the highest HI intensity tiles. For the smallest derivative kernel, 80′′ FWHM, the selected filamentary structures correspond to up to 25% of the area of the tiles but it can be up to 80% in the case of the 160′′ FWHM kernel. The selected percentage does not show any evident correlation with the orientation of the filamentary structures.

Appendix B: Comparison with GALFA-H I

With an angular resolution of 4′, GALFA-HI is the highest resolution single-dish observation that we can use to evaluate potential artifacts introduced by the inteferometer and the validity of our masking scheme in the analysis of the THOR-HI data. We present an example of the Hessian analysis of a 2° × 2° and 1.5-km s⁻¹ tile in both surveys in Fig. B.1. Figure B.3 confirms that the Hessian analysis of the GALFA HI and THOR-HI observations gives very similar results when using the same derivative kernel size (240′) and the same criteria for filament selection.

Fig. B.1 Examples of the results of the Hessian analysis applied to the GALFA HI and THOR-HI observations, shown in the top and bottom panels respectively.

Fig. B.2 Longitude-velocity diagram of the projected Rayleigh statistic V corresponding to the results of the Hessian analysis of the GALFA-HI and THOR-HI observations, shown in the top and bottom. Both results were obtained with the same derivative kernel size, indicated in the figure.

Fig. B.3 Projected Rayleigh statistic V calculated from the Hessian analysis of the GALFA-HI and THOR-HI observations using a 240′′ FWHM derivative kernel. Each point corresponds to a 2° × 2° × 1.5 km s−1 velocity-channel tile fully covered by both surveys. The vertical and horizontal dotted lines indicate |V | = 3, which roughly corresponds to a 3-σ confidence level in the estimation of a preferential orientation of 0° or 90°.

Fig. B.4 Histogram of the mean orientations of the HI filaments identified using the Hessian analysis with a 240′′ FWHM derivative kernel applied to 2° × 2° × 1.5 km s−1 velocity-channel tiles in the GALFA-HI and THOR-HI observations.

Figure B.2 indicates that the vertical filaments identified in THOR-HI around l ≈ 36° and v_LSR ≈ 40 km s⁻¹ are not prominent in GALFA-HI. It also indicates that the deviations from the general trend of V > 0 around the maximum and minimum v_LSR are commonto both datasets, but the higher sensitivity GALFA-HI allows for their observation in more velocity channels.

Appendix C: Other filament finding methods

Throughout this paper, we have used an algorithm for the identification of filamentary structures that is based on the Hessian matrix. This is not the only filament finder algorithm, although it is considerably faster than other alternatives. Additionally, the definition of its parameters (that is, the derivative kernel size and the curvature threshold) can be readily made in a self-consistent fashion.

In this section, we consider two alternative algorithms, FilFinder and the RHT, and show that their results are consistent with those found using the Hessian matrix. These two algorithms are computationally more demanding, but they do not offer a significant advantage in the study of the orientation of the filamentary structures. However, they are very powerful tools for studying properties such as the filament length and width that may be of interest in a follow-up study. C.1:

C.1 FilFinder

FilFinder is a Python package for extraction and analysis of filamentary structure in molecular clouds introduced in Koch & Rosolowsky (2015). It segments filamentary structure by using adaptive thresholding. This thresholding is made over local neighborhoods, allowing for the extraction of structure over a large dynamic range. Using the filament mask, the length, width, orientation and curvature are calculated. Further features include extracting radial profiles along the longest skeleton path, creating filament-only model images, and extracting values along the skeleton.

However, one limiting restriction in the implementation of FilFinder is the size of the map. When applying FilFinder to the 2° × 2° tiles in THOR-HI, the memory requirements made it impractical to use it without careful masking of the maps. We used a mask based on the curvature λ₋ obtained from the Hessian analysis to produce the example presented in Fig. C.1. Different masking schemes, such as those based on the intensity, will introduce bias toward filaments in channels with low intensity background.

Figure C.1 shows that FilFinder highlights most of the elongated structures found with the Hessian algorithm. It also finds significant connectivity among them, resulting in very long curved filaments. We conclude from this experiment with FilFinder that the masking scheme would constitute the main source of discrepancy between this and other methods. Given that we obtained the best results using a mask that is based on the curvature, FilFinder does not produce a different result than that in the Hessian method. Needless to say, FilFinder constitutes a powerful tool to study other properties of the HI filaments, such as their width or length, but those studies are beyond the scope of this work.

C.2 The rolling Hough transform method

The rolling Hough transform (RHT, Clark et al. 2014) is a tool for quantifying the linearity and spatial coherence of HI structures. In contrast to the Hessian matrix analysis, which is based on second order spatial derivatives, the RHT is a mapping between the image space (x, y) and the space defined by the transformation $$\rho =x\cos \theta +y\sin \theta ,$$(C.1)

whose coordinates are (ρ, θ) Although the procedure described in Clark et al. (2014) introduces additional steps aimed to select and evaluate image space features at a particular scale, at the core of the method is the transformation of each of the image points into a straight line in a parameter space. The output of the RHT is the function R(θ, x, y), which contains information on the directionality (θ) of the image features in the position (x, y).

In this section, we evaluate the difference in the results of the Hessian matrix and the RHT method by applying both techniques to the same area of the THOR HI observations. Figure C.2 shows the results of both analysis in a particular velocity channel. The leftmost panel shows the intensity maps and the middle panel show the values of the curvature eigenvalues of the Hessian matrix (λ₋), introduced in Eq. (2). The rightmost panel shows the RHT backprojection $$R(x,y)={\displaystyle \int R}(x,y,\theta ),$$(C.2)

which is a visualization of the linear structures identified by the RHT. At first glimpse, and without any selection based on the signal-to-noise ratio of the intensity map, both methods seem to trace the same structures.

Figure C.3 shows the HI filament orientation obtained using the Hessian method and the RHT in the tiles presented inFig. C.2. The great similarity in the histograms of the orientation angles obtained with the two methods confirms that there is no significant difference between the global results obtained with either of them. This is further confirmed in Fig. C.4, which shows the results for both methods across velocity channels toward the same region.

The description of the filamentary structure in the Hessian and RHT methods is different. While the Hessian matrix offers a characterization of the topology of the 2D scalar field, intensity map, the RHT describes the filamentary structures in that scalar field as straight line segments, as shown in Fig. C.2. However, this fundamental difference does not produce a significant difference in the distribution of orientations of the filamentary structures.

Fig. C.1 Example of the FilFinder application to a pair of velocity channel maps in THOR-HI. Left: intensity map at the indicated line-of-sight velocity. Middle: map of the eigenvalue λ− obtained with a 120′ FWHM derivative kernel. This map is used as a mask for FilFinder. Right: FilFinder filaments overlaid on the respective HI intensity map smoothed to 120′ FWHM.

Fig. C.2 Comparison of the results of the Hessian and the RHT methods for two different velocity-channel maps toward a region within the area covered by THOR-HI. Left: intensity map at the indicated line-of-sight velocity, vLSR. Middle: mapof the eigenvalue λ− of the Hessian matrix. Right: map of the RHT backprojection.

Fig. C.3 Histograms of the orientation of the filamentary structures identified using the Hessian matrix and the RHT methods in the velocity-channel maps presented in Fig. C.2.

Fig. C.4 Circular statistics used for the study of the orientation of the HI filamentary structures found using the Hessian and the RHT methods across velocity channels toward the region presented in Fig. C.2. Top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

Appendix D: FRIGG simulations

We produced synthetic HI emission observations following the procedure described in Appendix C of Soler et al. (2019) for a line of sight (LOS) that is perpendicular to the initial direction of the magnetic field. We followed the general radiative transfer equations presented in Draine (2011). We calculated the emissivity $$j(v_{\text{LSR}})=\frac{3}{16\pi }h{\nu }_0{A}_{10}{n}_{\text{H}}\varphi (v_{\text{LSR}})$$(D.1)

and absorption coefficient $$\kappa (v_{\text{LSR}})=\frac{3}{32\pi }\frac{h{c}^2}{k_{\text{B}}{\nu }_0}{A}_{10}\frac{n_{\text{H}}}{T_{\text{s}}}\varphi (v_{\text{LSR}}),$$(D.2)

where ν₀ = 1.4204 GHz is the frequency at the line center, A₁₀ = 2.876 × 10⁻¹⁵ s⁻¹ is the Einstein coefficient of the line, n_H is the HI number density, and ϕ(v_LSR) is the normalized line profile as a function of velocity. We assumed a Gaussian velocity distribution $$\varphi (v_{\text{LSR}})=\frac{1}{\sqrt{\pi }\Delta v}\frac{c}{{\nu }_0}{\text{e}}^{-{[(v_{\text{LSR}}-v)/\Delta v]}^2},$$(D.3)

where Δv ≡${(2k_{\text{B}}{T}_{\text{k}}/m_{\text{H}})}^{1/2}$ is the Doppler width of the line. Given that the FRIGG simulations do not include a chemical network, we assumed n_H to be 60% of the gas density in each cell. This rough assumption, which is based on the average values in the Milky Way at galactocentric distances R < 20 kpc (Draine 2011), does not have a critical effect in the results of our study of this numerical experiment.

We assumed that the excitation temperature (spin temperature, T_s) of the 21-cm line is the same as the kinetic temperature of the gas, T_k. This assumption leads to an over estimation of T_s in warm gas(see, for example Kim et al. 2014, and the references therein), but it is sufficient for the current exploration of the simulations. We used parallel ray-tracing assuming a far-field geometry. We calculated the absorption coefficient using Eq. (D.2) and the contribution to the local optical depth from the ith simulation element along the LOS, which is τⁱ (v_LSR) =κ(v_LSR;s_i)Δs, where s_i =iΔs and Δs = 1.95 pc. The total optical depth along an LOS is $$\tau (v_{\text{LSR}})={\sum }_{i=1}^N{\tau }^i(v_{\text{LSR}}),$$(D.4)

where N is the total number of LOS element in a given l, b direction. The synthetic HI emission is the result of the sum of the brightness temperature of each element and the foreground attenuation,such that $$T_{\text{B}}(v_{\text{LSR}})={\sum }_{i=1}^N\left[T_{\text{s}}^i\left(1-{\text{e}}^{-{\tau }^i(v_{\text{LSR}})}\right)\exp \left(-{\sum }_{j=1}^{i-1}{\tau }^i(v_{\text{LSR}})\right)\right].$$(D.5)

We calculated the emission between − 30 and 30 km s⁻¹ in velocity channels of 0.5 km s⁻¹ width, which we smoothed and regridded to match the velocity resolution of the THOR-HI observations. The output PPV cubes contain 512² pixels in the spatial dimensions, which we convolve with a two-dimensional Gaussian function with 40′′ FWHM to match the angular resolution of the THOR-HI observations.

The front domain of the simulations is located 7.5 kpc away to match the kinematic distance of the features found at l ≈ 20° and v_LSR = 120 km s⁻¹. Given that these simulations do not include a chemical network, we use thresholds in density and temperature to define the regions that produce HI emission. Explicitly, we consider the emission from the simulation voxels with particle densities n < 50 cm⁻³ and temperatures in the range 20 <T <3× 10⁴ K, which correspond to the general range of physical conditions where atomic hydrogen is expected (see, for example, Walch et al. 2015). Although these thresholds provide a very rough approximation to the real conditions in the ISM, they are sufficient for the numerical experiment that we are proposing here.

Fig. D.1 Synthetic observations of the HI emission from the FRIGG simulations with a standard (3 μG) and very high (12 μG) initial magnetizations, shown in the left and right panels, respectively. These snapshots are taken 14 Myr before those shown in Fig. 16. The colors represent the emission in the three indicated radial velocity bins with the same average emission.

Fig. D.2 Synthetic HI emission intensity and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels in the synthetic observations presented in Fig. D.1. Top: HI emission averaged over the whole map and toward the indicated position. Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

The synthetic HI observations, shown in Fig. 16, reveal a plethora of filamentary structures in the emission from both simulations. In the simulation with standard magnetization, these filaments form an intricate network that extends up to roughly 125 pc above and below the midplane. In the simulation with very high magnetization, these filaments are much longer, they appear perpendicular to the plane both below and above it, and extend further that in the standard magnetization case. In both cases, these filaments mostly correspond to density structures, as shown in the integrated density plots in Iffrig & Hennebelle (2017), rather that to fluctuations in the velocity field. For the sake of comparison with the results presented in Sect. 6.1, we present in Figs. D.1 and D.2 the results of the synthetic observations and Hessiananalysis of the standard and high-magnetization simulations in snapshots taken at 50 and 66 Myr, respectively.

Appendix E: Cloud Factory simulations

We used the cubes of atomic hydrogen resulting from the time-dependent chemical model in the CloudFactory simulations to perform synthetic observations using the lime radiative transfer code (Brinch & Hogerheijde 2010). The modifications involve the optimization of the code and its adaptation to deal with simulated unstructured meshes in the second paper of the CloudFactory series (Izquierdo et al. 2020). Prior to running lime, we use the grid and rt tools of the sf3dmodels package to read, process and provide format to the simulation snapshots (Izquierdo et al. 2018).

The HI 21 cm hyperfine line has a low emission coefficient A₁₀ = 2.876 × 10⁻¹⁵ s⁻¹ and hence it is easily populated by collisions even in the diffuse interstellar medium of the Galaxy (Purcell & Field 1956).For this reason, it is sufficient to assume an LTE approximation to solve the radiation transport in this scenario, though lime provides as well an iterative non-LTE implementation suitable for different gas tracers and physical conditions. For the level populations, the statistical weights of the upper and lower spin states are g₁ = 3 (triplet) and g₀ = 1 (singlet state), respectively. We are assuming that the spin temperature (T_s) of the HI tracks the gas kinetic temperature (T_k), which is an acceptable approximation for T_k < 5 ×10³ K (see, for example, Kim et al. 2014). For T_k > 5 ×10³ K, this assumption leads to an over estimation of the spin temperature and the linewidths of the HI emission that does not significantly affect the results of our study.

We calculated the emission between − 50 and 50 km s⁻¹ in velocity channels of 0.5 km s⁻¹ width, which we smoothed and regridded to match the velocity resolution of the THOR-HI observations. The output PPV cubes contain 512² pixels in the spatial dimensions, which we convolve with a two-dimensional Gaussian function with 40′′ FWHM to match the angular resolution of the THOR-HI observations. As in the FRIGG simulations, the front domain of the simulations is located 7.5 kpc away.

We produced two edge-on views of the simulation snapshots to explore variations in the relative orientation of HI filamentary structures. The line of sight orthogonal to the Galactic rotation is presented in Fig. 18. The line of sight tangent to the Galactic rotation is presented in Fig. E.1. The line of sight selection does not produce a significant difference in the orientation of the HI filaments, as illustrated in Fig. E.2.

Fig. E.1 Synthetic observations of the HI emission from the potential and feedback dominated CloudFactory simulations, presented in the left and right panels, respectively. These results correspond to a line of sight at 90° from that considered in the synthetic observations presented in Fig. 18. The colors represent the emission in the three indicated radial velocity bins, each of them with the same average emission.

Fig. E.2 Synthetic HI emission intensity and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels in the synthetic observations presented in Fig. E.1. Top: HI emission averaged over the whole map and toward the indicated position. Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

References

All Figures

Fig. 1 Examples of the Hessian analysis applied to three different velocity-channel maps toward a section of the Galactic plane observed by THOR. Left: intensity map. Center: map of the eigenvalue of the Hessian matrix identified as λ− in Eq. (2), which is used to characterize the filamentary structures in the intensity map. The overlaid gray map shows the filamentary structure obtained from the Hessian matrix analysis of the continuum noise maps, see Appendix A for details. Right: map of the orientation angles evaluated using Eq. (3).

In the text

	Fig. 2 Histograms of orientation angles (θ, as defined in Eq. (3)) of the filamentary structures obtained by applying the Hessian matrix method to the velocity channels maps presented in Fig. 1. The angles θ = 0° and θ = 90° correspond to the directions parallel and perpendicular to the Galactic plane, respectively.
In the text

Fig. 3 Mean spectrum and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels toward the region presented in Fig. 1. Top: mean HI intensity (I). Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: circular mean ($<\theta >$) and standard deviation (σV) of the orientation angles, represented by the line and the shaded area, respectively. The error bars correspond to the observation uncertainties and the error propagation based on the Monte Carlo sampling introduced in Appendix A.1. The dashed horizontal lines in each panel indicate the values of those quantities in the 1.4 GHz noise maps used to characterize the effect of the side lobes and continuum sources in the HI maps. The vertical gray lines correspond to the velocity channels shown in Fig. 1.

In the text

Fig. 4 Longitude-velocity (lv) diagrams corresponding to the mean intensity and the results of the Hessian analysis of the velocity-channel maps in the region covered by the THOR observations. Each pixel in the diagrams represents the quantities calculated in 2° × 2° tile in a 1.5 km s−1-velocity channel as follows. Top: mean HI intensity (I). Middle top: mean resulting vector (r). Middle bottom: projected Rayleigh statistic (V). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures, which is only well-defined if r > 0. The overlaid curves correspond to selected spiral arms from the model presented in Reid et al. (2016). The vertical line indicates thecentral l of the 2° × 2° tile presented inFig. 1 and whose circular statistics are shown in Fig. 3.

In the text

	Fig. 5 GALFA-HI (top) and THOR-HI (bottom) observations and their respective filaments from the Hessian analysis toward the same region of the Galactic plane in the indicated velocity channels.
In the text

Fig. 6 Projected Rayleigh statistic V calculated from Hessian analysis of the THOR-HI and GALFA-HI observations using a 120′′ FWHM derivative kernel. Each point corresponds to a 2° × 2° × 1.5 km s−1 velocity-channel tile fully covered by both surveys. The vertical and horizontal dotted lines indicate |V | = 3, which roughly corresponds to a 3-σ confidence level in the estimation of a preferential orientation of 0° or 90°. The difference in the range of V is the result of the larger significance obtained with the higher-resolution observations.

In the text

Fig. 7 Projected Rayleigh statistic (V) corresponding to the orientation of the filamentary structures identified with the Hessian matrix method in the HI emission in 2° × 2° tiles and across velocity channels. The overlaid curves correspond to selected spiral arms from the model presented in Reid et al. (2016). The black contours corresponds to the 97 and 99 percentile of the 12CO emission in the observations presented in Dame et al. (2001) integrated over the Galactic longitude range |b| ≤ 1°. The green stars correspond to the HII regions in the WISE catalog of Galactic HII regions (Anderson et al. 2014). The two most prominent regions of high-mass star formation in this portion of the Galaxy, W49 and W51, are indicated by the corresponding labels. The magenta crosses correspond to the 1720-MHz OH masers in the catalog (Beuther et al. 2019), which are typically excited by shocks. The black crosses correspond to the compact HI clouds at high forbidden velocities identified in the VGPS (Stil et al. 2006b). The ellipses labeled with letters indicate the regions of interest (ROIs) discussed in Sect. 5.2.

In the text

	Fig. 8 THOR-HI observations (top) and their respective filaments from the Hessian analysis (bottom) toward l and vLSR ranges dominated by filamentary structures parallel to the Galactic plane.
In the text

	Fig. 9 Magdalena HI filament. Top: atomic hydrogen (HI) emission in three velocity channels. Middle: filaments obtained using the Hessian technique toward the corresponding channels. Bottom: spectra toward the positions indicated by the crosses in the top panel. The vertical lines indicate the vLSR that corresponds to the emission shown in the other two panels.
In the text

Fig. 10 THOR-HI (top) observations and their respective filaments from the Hessian analysis (bottom) toward the Galactic longitude range of the region of interest (ROI) A, as defined in Fig. 7. The yellow circles correspond to the positions and sizes of the supernova remnants in the Green (2019) catalog. The cyan circles correspond to the positions and sizes of the HII regions in the Anderson et al. (2014) catalog with vLSR within ±10 km s−1 of the central velocity indicated in the figure. The magenta crosses mark the position of the THOR OH 1720-MHz masers in that velocity range.

In the text

	Fig. 11 Same as Fig. 10, but for the ROI B.
In the text

	Fig. 12 Same as Fig. 10 for the ROI C.
In the text

Fig. 13 Selected example of a prominent HI self-absorption (HISA) structure. Top: atomic hydrogen (HI) emission in three velocity channels. Middle: projected Rayleigh statistic (V, Eq. (5)) corresponding to the orientation of the filamentary structures in the HI emission in the corresponding channels. Bottom: spectra toward the positions indicated by the crosses in the top panel. The vertical lines indicate the vLSR that corresponds to the emission shown in the other two panels.

In the text

	Fig. 14 THOR-HI emission for two velocity channels in the region of interest (ROI) D (left and right) and a channel toward the same l range but dominated by filaments parallel to the Galactic plane (middle).
In the text

Fig. 15 Comparison between the values of the projected Rayleigh statistic V calculated using the THOR and 13CO GRS observations in the indicated Galactic-longitude ranges. Both are calculated using a derivative kernel with 120′′ FWHM. Each point in the scatter plot corresponds to the values estimated using the Hessian analysis in the same 2° × 2° tile and 1.5-km s−1 velocity channel in both datasets.

In the text

	Fig. 16 Synthetic observations of the HI emission from the FRIGG simulations with a standard (3 μG) and very high (12 μG) initial magnetizations, shown in the left and right panels, respectively. The colors represent the emission in the three indicated radial velocity bins with the same average emission.
In the text

Fig. 17 Synthetic HI emission intensity and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels in the synthetic observations presented in Fig. 16. Top: HI emission averaged over the whole map and toward the indicated position. Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

In the text

	Fig. 18 Synthetic observations of the HI emission from the potential and feedback dominated CloudFactory simulations, presented in the left and right panels, respectively. The colors represent the emission in the three indicated radial velocity bins, each with the same average emission.
In the text

Fig. 19 Synthetic HI emission intensity and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels in the synthetic observations presented in Fig. 18. Top: HI emission averaged over the whole map and toward the indicated position. Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

In the text

Fig. A.1 Hessian analysis applied to a velocity-channel map dominated by noise toward the same region presented in Fig. 1. Left: HI intensity map. Center: map of the eigenvalue of the Hessian matrix identified as λ− in Eq. (2), which is used to characterize the filamentary structures in the intensity map. The overlaid gray map shows the filamentary structure obtained from the Hessian matrix analysis of the continuum noise maps, see Appendix A for details. Right: map of the orientation angles evaluated using Eq. (3).

In the text

	Fig. A.2 Same as Fig. A.1, but for the noise map of the continuum emission at 1.42 GHz (${\sigma }_{I_{1.42}}$), which we use as a proxy of the structure introduced by the side-lobes in the imaging of the HI emission.
In the text

	Fig. A.3 Same as Fig. 1, but for derivative kernels of 80′′ (top) and 160′′ FWHM (bottom), which correspond to two and four times the native resolution of the THOR-HI observations, respectively.
In the text

	Fig. A.4 Longitude-velocity diagram of the projected Rayleigh statistic V estimated using the indicated derivative kernel sizes.
In the text

	Fig. A.5 Longitude-velocity diagram of the mean orientation angle $<\theta >$ of the HI filamentary structures estimated using the Hessian method and the indicated derivative kernel sizes.
In the text

	Fig. A.6 Longitude-velocity diagram showing the percentage of each velocity-channel map that is selected in the Hessian analysis. Each pixel corresponds to a 2° × 2° tile.
In the text

	Fig. B.1 Examples of the results of the Hessian analysis applied to the GALFA HI and THOR-HI observations, shown in the top and bottom panels respectively.
In the text

	Fig. B.2 Longitude-velocity diagram of the projected Rayleigh statistic V corresponding to the results of the Hessian analysis of the GALFA-HI and THOR-HI observations, shown in the top and bottom. Both results were obtained with the same derivative kernel size, indicated in the figure.
In the text

Fig. B.3 Projected Rayleigh statistic V calculated from the Hessian analysis of the GALFA-HI and THOR-HI observations using a 240′′ FWHM derivative kernel. Each point corresponds to a 2° × 2° × 1.5 km s−1 velocity-channel tile fully covered by both surveys. The vertical and horizontal dotted lines indicate |V | = 3, which roughly corresponds to a 3-σ confidence level in the estimation of a preferential orientation of 0° or 90°.

In the text

	Fig. B.4 Histogram of the mean orientations of the HI filaments identified using the Hessian analysis with a 240′′ FWHM derivative kernel applied to 2° × 2° × 1.5 km s−1 velocity-channel tiles in the GALFA-HI and THOR-HI observations.
In the text

	Fig. C.1 Example of the FilFinder application to a pair of velocity channel maps in THOR-HI. Left: intensity map at the indicated line-of-sight velocity. Middle: map of the eigenvalue λ− obtained with a 120′ FWHM derivative kernel. This map is used as a mask for FilFinder. Right: FilFinder filaments overlaid on the respective HI intensity map smoothed to 120′ FWHM.
In the text

	Fig. C.2 Comparison of the results of the Hessian and the RHT methods for two different velocity-channel maps toward a region within the area covered by THOR-HI. Left: intensity map at the indicated line-of-sight velocity, vLSR. Middle: mapof the eigenvalue λ− of the Hessian matrix. Right: map of the RHT backprojection.
In the text

	Fig. C.3 Histograms of the orientation of the filamentary structures identified using the Hessian matrix and the RHT methods in the velocity-channel maps presented in Fig. C.2.
In the text

Fig. C.4 Circular statistics used for the study of the orientation of the HI filamentary structures found using the Hessian and the RHT methods across velocity channels toward the region presented in Fig. C.2. Top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

In the text

	Fig. D.1 Synthetic observations of the HI emission from the FRIGG simulations with a standard (3 μG) and very high (12 μG) initial magnetizations, shown in the left and right panels, respectively. These snapshots are taken 14 Myr before those shown in Fig. 16. The colors represent the emission in the three indicated radial velocity bins with the same average emission.
In the text

Fig. D.2 Synthetic HI emission intensity and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels in the synthetic observations presented in Fig. D.1. Top: HI emission averaged over the whole map and toward the indicated position. Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

In the text

Fig. E.1 Synthetic observations of the HI emission from the potential and feedback dominated CloudFactory simulations, presented in the left and right panels, respectively. These results correspond to a line of sight at 90° from that considered in the synthetic observations presented in Fig. 18. The colors represent the emission in the three indicated radial velocity bins, each of them with the same average emission.

In the text

Fig. E.2 Synthetic HI emission intensity and circular statistics used for the study of the orientation of the HI filamentary structures across velocity channels in the synthetic observations presented in Fig. E.1. Top: HI emission averaged over the whole map and toward the indicated position. Middle top: mean resulting vector (r), which indicates if the distribution of orientation angles is flat (r ≈ 0) or sharply unimodal (r ≈ 1). Middle bottom: projected Rayleigh statistic (V), which indicates if the distribution of orientation angles is clearly peaked around 0° (V ≫ 0) or 90° (V ≪ 0). Bottom: mean orientation angle ($<\theta >$) of the HI filamentary structures.

In the text