The history of dynamics and stellar feedback revealed by the HI filamentary structure in the disk of the Milky Way

We present a study of the filamentary structure in the emission from the neutral atomic hydrogen (HI) at 21 cm across velocity channels in the 40"-resolution observations in The HI/OH/Recombination (THOR) line survey of the inner Milky Way. Using the Hessian matrix method in combination with tools from circular statistics, we find that the majority of the filamentary structures in the HI emission are aligned with the Galactic plane. Part of this trend can be assigned to long filamentary structures that are coherent across several velocity channels. However, we also find ranges of Galactic longitude and radial velocity where the HI filamentary structures are preferentially oriented perpendicular to the Galactic plane. These are located (i) around the tangent point of the Scutum spiral arm and the terminal velocities of the Molecular Ring, around $l\approx 28${\deg} and $v_{\rm LSR}\approx 100$ km/s, (ii) toward $l\approx 45${\deg} and $v_{\rm LSR}\approx 50$ km/s, (iii) around the Riegel-Crutcher cloud, and (iv) toward the positive and negative terminal velocities. Comparison with numerical simulations indicates that the prevalence of horizontal filamentary structures is most likely the result of the large-scale Galactic dynamics and that vertical structures identified in (i) and (ii) may arise from the combined effect of supernova (SN) feedback and strong magnetic fields. The vertical filamentary structures in (iv) can be related to the presence of clouds from extra-planar HI gas falling back into the Galactic plane after being expelled by SNe. Our results indicate that a systematic characterization of the emission morphology toward the Galactic plane provides an unexplored link between the observations and the dynamical behaviour of the interstellar medium, from the effect of large-scale Galactic dynamics to the Galactic fountains driven by SNe.


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 (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 Corresponding author, e-mail: soler@mpia.de at 21 cm using the observations with broadest dynamic range in spatial scales toward the Galactic plane available to this date.
The structure of the Hi emission in small velocity intervals has revealed a multitude of filamentary (slender or threadlike) structures, 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(Heiles , 1984McClure-Griffiths et al. 2002). These observations constitute clear evidence of the injection of energy into the ISM by supernova explosions (see for example, Cox & Smith 1974;McKee & Ostriker 1977;Mac Low & Klessen 2004).
The study of the Hi structure has been 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 3,000 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 applying the unsharp mask (USM), another technique from machine vision to enhance the contrast 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 −1 . Both of these studies find a significant correlation between the elongation of these filamentary structures and the orientation of the local interstellar magnetic field, which may be the product of magneticallyinduced velocity anisotropies, 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;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. 2006b), 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 4 K or cold neutral medium (CNM) with T ≈ 10 2 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 hairlike CNM filaments aligned with the ambient magnetic field toward the Riegel-Crutcher cloud. However, there has been no 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 survey of the inner Milky Way (THOR, 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 in-terferometric observations made in VGPS and THOR, respectively. We focus on the statistics of a 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 Sec. 2. We present the method used for the characterization of the topology of the Hi emission in Sec. 3 and comment on the results of our analysis in Sec. 4. In Sec. 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 Sec. 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  projects. Finally, we present our conclusions in Sec. 7. We reserve additional analysis 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 the analysis of the MHD simulations from FRIGG and CloudFactory, respectively.

Atomic hydrogen emission
For the main analysis in this paper we use the Hi positionposition-velocity 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 position-position-velocity (PPV) cubes, I(l, b, v), are set in Galactic coordinates and a Cartesian projection in a spatial grid with 10 × 10 pixels and 1.5 km s −1 velocity channels. For details on the calibration and imaging of these data we refer to Bihr et al. (2015); Beuther et al. (2016); 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. We re-projected the GALFA-Hi into the THOR-Hi spatial and spectral grid in two steps. First, we smoothed and re-gridded the data in the spectral dimension by using the tools in the spectral-cube package in astropy (Astropy Collaboration et al. 2018). Second, we projected the observations into the same spatial grid of the THOR-Hi data by using the reproject tools, also from astropy.

Carbon monoxide (CO) emission
We compared the Hi emission observations with the 13 CO (J = 1 → 0) observations from The Boston University-Five College Radio Astronomy Observatory Galactic Ring Survey (GRS, Jackson et al. 2006). The GRS survey has 46 angular resolution with an angular sampling of 22 . In this particular region, it covers the velocity range −5 ≤ v LSR ≤ 135 km s −1 at a resolution of 0.21 km s −1 . It has a typical root mean square (RMS) sensitivity of 0.13 K. We re-projected and re-gridded the GRS data using the same procedure followed with the GALFA-Hi. We also used the 12 CO (J = 1 → 0) presented in (Dame et al. 2001) for the comparison with the results of our analysis, in particular for the illustration of the molecular ring structure across Galactic longitude and velocity.

Catalogs of Hii 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 catalogue 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 . We also referred to the catalogs of Galactic supernova remnants presented in Anderson et al. (2017) and Green (2019).

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 threedimensional space, but rather aims to study the charateristics of elongated structures in the position-position-velocity 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 and FilFinder (Koch & Rosolowsky 2015), as we show in App. C. Our implementation of this method is as follows.

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 estimate the first and second derivatives with respect to the local coordinates (x, y) and build the Hessian matrix, where H xx ≡ ∂ 2 I/∂x 2 , H xy ≡ ∂ 2 I/∂x∂y, H yx ≡ ∂ 2 I/∂y∂x, and H yy ≡ ∂ 2 I/∂y 2 are the second-order partial derivatives and x and y refer to the Galactic coordinates (l, b) as x ≡ l cos b and y ≡ b, so that the x-axis is parallel to the Galactic plane. 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 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 −1 to match the THOR-Hi spectral resolution. The results obtained with different derivative kernel sizes and ∆v selections are presented in App. A. The two eigenvalues (λ ± ) of the Hessian matrix are found by solving the characteristic equation, 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. The eigenvector corresponding to λ − defines the orientation of intensity ridges with respect to the Galactic plane, which is characterized by the angle This angle 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.

Selection of the filamentary structures
We apply the Hessian matrix analysis to each velocity channel map, as illustrated in Fig. 1, but then select 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 select regions where I(l, b) > 5σ I , where σ I is approximately 4 K and is estimated from the standard deviation of I in the velocity channels with the lowest mean I. 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 App. 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 velocity channels with the lowest mean intensity and use the median of their minimum value of λ − as the threshold value, λ C − . This selection minimizes the effect of artifacts in the noise-dominated channels without any loss of generality. We exclusively consider regions of each velocity channel map where λ − < λ C − , given that the most filamentary structure has more negative values of λ − .

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 θ .
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, θ 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 (r ≈ 0) or multimodal (V ≈ 0).

Is there a preferred orientation of the Hi filaments?
The mean resultant vector is defined as 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 i j is the statistical weight of each angle θ i j . We account for the spatial correlations introduced by the telescope beam by choosing w i j = (δx/∆) 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.

3.3.2.
What is the significance of the 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 non-uniform and peaked at 0 • or 90 • . It is defined as which follows the same conventions introduced in Eq. (4). The value of V represents the likelihood test against a von Mises distribution, which is the circular normal distribution centered on 0 • (see for example, Jow et al. 2018, and references therein).

Measurements of V
√ 2 indicate a significant detection of structures parallel to the Galactic plane, while measurements of V − √ 2 indicate a significant detection of structures perpendicular to it.
3.3.3. Is there another mean orientation different to parallel or perpendicular to the Galactic plane?
For the sake of completeness, we also computed the mean orientation angle, defined as θ ≡ 0.5 arctan i j w i j sin(2θ i j ) i j w i j cos(2θ i j ) .
This quantity highlights orientations that are not considered in the statistical test V, that is, 0 • or 90 • . However, given the broad distribution of angles theta θ, shown in Fig. 2, θ has large standard deviations Figure 3 presents an example of the values of the three quantities that we used to characterize the orientation angle distributions, across velocity channels toward the test tile presented in Fig. 1. For reference, we also include the values estimated from the structure of the noise. We note that the values of r, V, and θ change across velocity channels and, except for very few channels, their values differ from those corresponding to the structure of the sidelobe noise.
We also note that 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 will further discuss in Sec. 4.
The relatively large positive values of V in Fig. 3 imply that the Hi filamentary structures are mostly parallel to the Galactic plane ( θ ≈ 0 • ). This tendency changes across velocity channels and it is significantly different from that of the noise, confirming that the side lobes are not a major contributor to the reported trends. These trends in V and θ are independent of the average intensity in each velocity channel and they are an exclusive result of the distribution of orientation angles of the filamentary structures, a fact that we will further discuss in Sec. 4.

Longitude-velocity diagrams
We applied the Hessian matrix analysis 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 App. 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 where the selection criteria introduced in Sec. 3.2 discard that region from the analysis or, in the case of V, the results of the orientation 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| √ 2 in the middle panel of Fig. 4, indicate that the Hi filaments are preferentially orientated either parallel or perpendicular to the Galactic plane. This is further confirmed by the values of θ , 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 App B. Figure 4 shows that the most significant trend is for the Hi filaments to be parallel to the Galactic plane and if the 90 • orientation is significant for this analysis, it is because it is clearly grouped in specific ranges of l and v LSR .

Comparison to previous observations
We compare 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 the THOR data 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 is very 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 by the presence of negative values of V in the THOR-Hi, which can be linked to the increase in the angular resolution of the observations by a factor of 6.
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 −1 , 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 App. B

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 that there are no significant differences in the trends corresponding to v LSR < 0 km s −1 and v LSR > 0 km s −1 , although the latter velocity range is more prone to velocity crowding due to the mapping of at least two positions into the same radial velocity, under the assumption of circular motions (Reid et al. 2014). 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 −1 velocity range.
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 Sec. 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 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 −1 and l ≈ 36 • and v LSR ≈ 40 km s −1 . 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 Sec. 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 θ values found toward the maximum and minimum velocities in Fig. 4, which deviates from the general trend V 0 and θ ≈ 0 • . However, velocity crowding does not provide a conclusive explanation for the prevalence of vertical Hi filaments, V 0 and θ ≈ 90 • , found around the maximum and minimum velocities toward 55. • 0 l 65. • 0.

Discussion
We articulate the discussion of the physical phenomena related to the results presented in Sec. 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.

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 paral- lel 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 trends appear to be concentrated 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 width of the derivative kernel (see App. A for a discussion on the selection of the kernel size).
The predominantly positive values of V observed at v LSR < 0 km s −1 can be associated with a distance effect. Assuming circular motions, v LSR < 0 km s −1 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 −1 , 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 −1 , 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 line-widths of the Hi emission across this filament 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). 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 −1 , marked as ROI A in Fig. 7. This position has been previously singled out in the study of the Hi emission toward the Galactic plane, 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, which appear 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 preferencial 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 Hi filamentary structures in the range l 22 • and v LSR 70 km s −1 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 effect 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 centre (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 distinguished 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 • , has also 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 filaments in ROI B, where the effect of the Galactic dynamics is less evident. 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 −1 . 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 − √ 2 and θ ≈ 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 −1 ) and W51 (l = 49. • 4, b = −0. • 3, v LSR ≈ 60 km s −1 ), two of the most promi-  nent 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 towards 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 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 Sec. 6.

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 −1 . 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). 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 selfabsorption (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 sharing 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); Kavars et al. (2005); and Wang et al. (2020c); 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.

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 to 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 2006Fraternali 2017).   Figure 14 shows that the structure of the emission at −83 and 46 km s −1 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 clearly 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 −1 identified in the VGPS (Stil et al. 2006a). 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. The vertical lines indicate the v LSR that corresponds to the emission shown in the other two panels.

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 13 CO (Goodman et al. 2014;Ragan et al. 2014;Wang et al. 2015;Zucker et al. 2015;Abreu-Vicente et al. 2017;Wang et al. 2020b). To establish a link between the Hi filaments and the GMFs, without detailing individual objects, we also applied the Hessian analysis to the 13 CO (J = 1 → 0) emission observations in the GRS survey. Following the selection criteria presented in Sec. 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 compila- tion presented in Zucker et al. (2018), we find that most of the 13 CO filamentary structures are parallel to the Galactic plane. However, we found no evident correlation between the orientation of the filamentary structures in both tracers.
There can be several reasons for this general lack of correlation between the Hi and the 13 CO filamentary structures. First, the linewidths of the 13 CO emission are narrower than those of the Hi and it is possible that we are washing away part of the orientation of the filaments by projecting 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 −1 and the interquartile range is 0.68 < σ v < 1.89 km s −1 , thus re-gridding the data to 1.5 km s −1 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 √ 2 indicates that there is a preferential orientation of the filamentary structures in 13 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 13 CO, as quantified in Soler et al. (2019), the filamentary structure in the Hi emission is not necessarily related to that of the 13 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 13 CO can appear completely decoupled ). This does not discard the local correlation between the morphology of both tracers toward filamentary structures, as reported in Wang et al. (2020c) and Syed et al. (2020b), but it shows that the orientation of the filamentary structures are generally different.

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 explore the effect of SNe, magnetic fields, and galactic rotation in the multiphase medium in two families of numerical simulations. First, we consider 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 consider 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.

Initial conditions
The FRIGG simulations use 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, 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 use 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, of 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 figure 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 0 of about 3 µG, and one with "very high" magnetization, B 0 ≈ 12 µG. Both simulations have an initial particle density of n 0 = 1.5 cm −1 . 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 (http://www.galactica-simulations.eu). This selection guarantees that the simulations have reached a quasi-steady state and does not affect the reported results. Details on the construction of the synthetic observations from these simulations are presented in App. D.

Hessian analysis results
The prevalence of longer filamentary structure with higher magnetization has been reported in previous studies (Hennebelle 2013;Seifried & Walch 2015;. 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.
The results of the Hessian analysis of the Hi emission from the FRIGG simulations are summarized 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 would fragment if the magnetic field were weaker. Subsequently, the gas layer consisting of ejected clouds falls back on the plane and is stretched 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 towards 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 observations.

Initial conditions
We consider the effect of the Galactic dynamics by using the CloudFactory simulations, presented in Smith et al. (2020). These simulations, run 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 10 M , 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 years, 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 years, 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. Details of the construction of the synthetic observations from these simulations are presented in App E.

Hessian analysis results
The results of the Hessian analysis of the Hi emission from the CloudFactory simulations 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| √ 2 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 towards quantifying the effect of SN feedback in the Galactic plane.

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 −1 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. Comparison with the numerical simulation of the Galactic dynamics and chemistry in the CloudFactory project indicate that elongated and non-selfgravitating 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 −1 and toward l ≈ 36 • and v LSR ≈ 40 km s −1 . 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 are signatures of the effect of the history of star formation and stellar feedback in the currect structure of the atomic gas in the Galactic plane.
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 the gas together, creating favourable 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 robust 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.    We calculated the curvature threshold values λ C − introduced in Sec. 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 λ 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.
To characterize the spatial structure of the noise, we use the noise map of the continuum emission at 1.42 GHz, σ 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 σ 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 σ I 1.42 rarely corresponds to that found in the Hi emission, as we show in the example presented in 3. But that is not necessarily the case in all the 2 • × 2 • regions, which motivates our masking scheme.

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 introduced with a very small kernel and the loss of information resulting from using a very coarse one. 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 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-neighbour 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.
The lv diagrams of the V and θ obtained with the 80 and 160 FWHM derivative kernels are shown in Fig. A.4 and Fig. 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. Filament selection
One of the main differences introduced by the selection of the derivative kernel size is the percentage of the map that is selected in the 2 • × 2 • and 1.5-km s −1 tiles across l and v LSR , as illustrated in Fig. A.6. The significant changes in the percentage of the maps covered by filamentary structures indicates that most of them are not resolved at the spatial scales corresponding to the 120 and 160 derivative kernels.
The selected percentage across l and v LS R roughly follows the same distribution of the mean intensity, shown in Fig. 4, which potentially indicate 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-Hi
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 −1 tile in both surveys in Fig. B.1.     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 Rolling Hough Transform, 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. 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 ρ = x cos θ + y sin θ, (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  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 in Fig. 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 C.2. However, this fundamental difference does not produce a significant difference in the distribution of orientations of the filamentary structures.     The output PPV cubes contain 512 2 pixels. The front domain of the simulations is located 3.2 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 pre-