© The Authors 2025

1 Introduction

Tracing the flows of matter and energy in the interstellar medium (ISM) is instrumental for understanding the structure and evolution of the Milky Way and other galaxies (see, for example, Ballesteros-Paredes et al. 2020; Tacconi et al. 2020; Chevance et al. 2023). Reconstructing the flow of gas in the ISM enables the identification of the mechanisms that make gas available for star formation and allows us to quantify the impact of galactic dynamics, magnetic fields, stellar winds, outflows, and supernovae (SNe) in the ISM (e.g., Field & Saslaw 1965; Tomisaka 1984; Koyama & Inutsuka 2000; Blitz et al. 2007; Hennebelle et al. 2008; Klessen & Hennebelle 2010; Dobbs et al. 2014; Krumholz et al. 2014; Schmidt et al. 2016; Hennebelle & Inutsuka 2019). In this work, we present a reconstruction of the ISM line-of-sight (LOS) motions in the Solar neighborhood using state-of-the-art models of the three-dimensional dust density distribution, which we associate with the velocity information from the neutral atomic hydrogen (H I) and carbon monoxide (CO) line emission using a machine vision tool to quantify the morphological similarities between distance and velocity channels (Soler et al. 2019).

Classic studies of the ISM distribution and kinematics rely on the assumption of circular motions around the Galactic center to compute what is usually called “kinematic distances” (see, for example, Oort et al. 1958; Roman-Duval et al. 2009; Sofue 2011; Wenger et al. 2018; Hunter et al. 2024). Kinematic distances depend on the Galactic rotation model and are affected by the kinematic distance ambiguity (see, for example, Eilers et al. 2020; Reid 2022). Moreover, the motion of the Galactic ISM is not purely rotational; non-circular streaming motions are induced by the Galactic bar, the spiral arms, and stellar feedback (Burton 1971; Gómez 2006; Moisés et al. 2011). The latter effect is expected to be prominent in the Solar neighborhood, where the motions introduced by supernovae (SNe)-blown bubbles and other local effects are potentially dominant over Galactic rotation as seen from the local standard of rest (LSR; see, for example, Zucker et al. 2023, and references therein)

In recent years, the simultaneous rise of Gaia and deep wide-field photometric surveys, such as the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS Chambers et al. 2016) and the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Allende Prieto et al. 2008), has sparked a revolution in our understanding the Galactic ISM in 3D. The remarkable stellar parallax and reddening observations provided by these surveys have triggered unprecedented reconstructions of the interstellar dust distribution in the three spatial dimensions (see, for example, Green et al. 2019; Leike et al. 2020; Lallement et al. 2022). These “3D dust” maps have enabled a variety of ISM studies, which include the characterization of the Local Bubble (Pelgrims et al. 2020; Zucker et al 2022; O’Neill et al. 2024), SNe-blown cavities (Bracco et al. 2023; Liu et al. 2024), and molecular cloud (MC) envelopes (Zucker et al. 2021; Mullens et al. 2024), among others.

In this work, we associated the state-of-the-art 3D dust models presented by Edenhofer et al. (2024) with the LOS velocity information provided by H I and CO surveys toward the Galactic plane. Following the pioneering work in Tchernyshyov & Peek (2017); Tchernyshyov et al. (2018), we adopted the term “kinetic tomography” to describe our reconstruction of the density and LOS velocity distribution. Our approach is based on the morphological similarity quantified by the histogram of oriented gradients (HOG) method (Soler et al. 2019). Using HOG we linked the position-position-distance (PPD) 3D dust density cubes and the position-position-velocity (PPV) H I and CO line emission cubes. The result is a position-position-distance-velocity (PPDV) hypercube, which we employ to characterize the local ISM motions along the line of sight. Line emission from ISM tracers and 3D dust have previously been combined to investigate the ISM in four dimensions (see, for example, Ivanova et al. 2021; Duchêne et al. 2023), but not with the level of reconstruction in the Edenhofer et al. (2024) maps, and never before with the HOG method.

We used the H I line emission at 21-cm wave-length, which traces both the cold, pre-molecular state before star formation and the warm, diffuse ISM before and after star formation (see, McClure-Griffiths et al. 2023, for a recent review). Since its discovery, H I has been instrumental in studying the diffuse ISM in the Milky Way (Ewen & Purcell 1951; Muller & Oort 1951; Pawsey 1951). Early observations of the H I absorption against radio continuum sources revealed the presence of narrow, few-km s⁻¹-wide, spectral features (e.g., Hagen et al. 1955; Clark 1965). In emission, these narrow features appear on top of broader 10 to 20-km s⁻¹-wide features (e.g., Matthews 1957). These observations inspired the formulation of a “two phase” H I model (Field et al. 1966; McKee & Ostriker 1977), in which at the pressure of the ISM, the heating and cooling processes naturally lead to two thermally stable states: a dense cold neutral medium (CNM; T ≈ 50 K and n ≈ 50 cm⁻³) immersed in a diffuse warm neutral medium (WNM; T ≈ 8000 K and n ≈ 0.3 cm⁻³).

Part of our understanding of the H I multiphase structure comes from studying absorption toward continuum sources. (see, for example, Strasser et al. 2007; Stanimirović et al. 2014; Murray et al. 2018). Comparisons between 21-cm H I emission and absorption measurements indicate that, in the vicinity of the Sun, the WNM has roughly the same column density as the CNM (Falgarone & Lequeux 1973; Liszt 1983). Crucial additional information about the distribution of the CNM in and around MCs comes from the portion of the CNM sampled by the extended absorption of background H I emission by cold foreground H I , which is generically known as H I self-absorption (HISA; Heeschen 1955; Gibson et al. 2000; Seifried et al. 2020). Observations indicate that the HISA distribution is related to that of the CO emission, although the CNM it reveals corresponds to less than 15% than the total gas mass (see, for example, Gibson et al. 2005; Krčo et al. 2008; Wang et al. 2020; Syed et al. 2020).

We also used the CO (J = 1 → 0) emission at 2.6 mm wave-length, which is the archetypal tracer of the cooler and denser ISM (Wilson et al. 1970). With its low excitation energy and critical density, CO provides an irreplaceable proxy for H₂, which is the most abundant molecule in the Galaxy but is much harder to observe (Combes 1991; Bolatto et al. 2013; Heyer & Dame 2015). The focus on the molecular material traced by CO in star formation (SF) studies was cemented by the observed correlation between the SF and CO surface densities in nearby galaxies, which contrasts with the lack of correlation between SF and H I at low surface densities in those objects (see, for example, Kennicutt 1998; Wong & Blitz 2002; Leroy et al. 2008; Schinnerer & Leroy 2024).

Soler et al. (2023) presented a pilot study of the HOG-based kinetic tomography using the 3D dust density reconstruction from Leike et al. (2020) and archival H I and CO line emission toward the Taurus molecular cloud (MC). The authors found anti-correlation between the dust density and the H I emission, which uncovers the CNM associated with the MC. They also found a pattern in LOS velocities and distances consistent with converging gas motions in the Taurus MC, with the cloud’s near side moving at higher velocities than the far side. This result is consistent with the kinematic imprint of the MC location at the intersection of two bubble surfaces, the Local Bubble (Pelgrims et al. 2020; Zucker et al. 2022) and the Per-Tau shell (Bialy et al. 2021).

In this paper, we applied the HOG method to study the atomic and molecular gas motions toward the nearby Milky Way’s disk, defined as the Galactic latitude range |b| ≤ 5° within the 1.25-kpc extent of the Edenhofer et al. (2024) 3D dust reconstruction. The presentation of the data, analysis methods, and results is organized as follows. In Sect. 2, we present the 3D dust models and the H I and CO observations used in the analysis. Section 3 summarizes the main aspects of the HOG method implementation for analyzing the Galactic plane. We describe the morphological correlation between the 3D dust models and the line emission in Sect. 4. In Sect. 5, we detail the streaming motions obtained with the HOG method and the energy and momentum densities derived from them. Section 6 discusses our results and their implications for understanding the ISM dynamics in the Solar neighborhood. We present our conclusions of this work in Sect. 7. We complement the main results of this work with the analysis shown in a set of appendices. Appendix A presents details on the HOG method’s error propagation and selection of parameters. We consider the effect of HISA in our results and the potential of the HOG method to identify these features in Appendix B. In Appendix C, we evaluate the effects of the fixed angular resolution and distance in the HOG results. Appendix D presents the HOG analysis of synthetic line emission and 3D density from a multiphase magnetohydrodynamic (MHD) MC simulation. Appendix E compares our results with the distances and line-of-sight velocities for the five maser sources within the studied volume. Finally, Appendix F presents the distance-velocity mapping for a few regions of interest, including areas toward the Galactic center and anticenter, where uncertainties in the kinematic distance estimates are considerable.

Fig. 1 Examples of the 3D dust, 12CO, and H I data combined in this paper. Top: Nucleon density derived from the 3D extinction model presented by Edenhofer et al. (2024) for three distance bins indicated in the figure. Middle: 12CO line emission from the Dame et al. (2001) survey for the three LOS velocity intervals indicated in the figure. Bottom: H I 21-cm line emission observations presented by HI4PI Collaboration (2016).

2 Data

We limit our analysis to the band of sky within |b| < 5°, which is the area covered by the most extensive spectroscopic survey of CO emission for the same energy transition and the same isotopologue (Dame et al. 2001). This region is covered by the most wide-ranging MC studies (Rice et al. 2016; Miville-Deschênes et al. 2017). As our results later indicate, the proximity to the Galactic midplane also facilitates the interpretation of the reconstructed LOS velocity pattern. We complement the CO observations with the whole-sky 3D dust models from Edenhofer et al. (2024) and H I emission observations introduced in HI4PI Collaboration (2016). We show an example of these three data sets in Fig. 1.

In this paper, we use v_LOS to refer to the radial velocity inferred from the emission line Doppler shift relative to the LSR, usually identified as v_LSR in the H I and CO observations. We introduce this convention to distinguish it from the general motion relative to the LSR, for which the three components of the velocity vector can be inferred in the Solar vicinity from the measurement of stellar motions (see, for example, Miret-Roig et al. 2022; Ratzenböck et al. 2023).

The selected b range implies that we are studying the volume covered by the revolution of a trapezoid with base sizes 12.16 and 220.4 pc located at 69 and 1250 pc from the LSR, as illustrated in Fig. 2. The limits along the LOS are set by the boundaries of the Edenhofer et al. (2024) 3D dust model. At the furthest point in this analysis, we cover a region roughly within the characteristic CO scale height around the Solar circle (~100 pc, Heyer & Dame 2015).

We consider the line emission in the range −25 < v_LOS < 25 km s⁻¹, which corresponds to the expected amplitude of LOS velocities for a Galactic rotation model for heliocentric distances within 1250 pc plus or minus a few kilometers per second. This restricted input range mitigates the spurious chance correlations produced by the limited angular resolution in the 3D dust reconstruction, as detailed in Appendix A. However, it also introduces a generic limitation in our analysis, as some local ISM structures can have significant morphological correlations beyond the input v_LOS range. Thus, in practice, the restricted v_LOS input range limits the amplitude of reconstructed LOS motions for the local density structures. Until future higher-resolution 3D dust reconstructions break the ambiguities that produce chance correlations and expand the reliable input v_LOS range, our HOG-based reconstruction should be considered lower limits for the streaming motions and other derived quantities.

Fig. 2 Rendering of the 3D dust density distribution from the extinction models presented in Edenhofer et al. (2024) for the region |b| < 5° considered in this paper. The yellow sphere represents the position of the Sun. The associated movie is available online.

2.1 3D dust density distribution

The primary dataset enabling our study is the reconstruction of the 3D dust density between 69 and 1250 pc from the Sun presented in Edenhofer et al. (2024). This model uses the extinction estimates from Zhang et al. (2023b), which are primarily based on the Gaia satellite’s BP/RP spectra, with a spectral resolution (R) of roughly between 30 and 100 (Gaia Collaboration 2023). Zhang et al. (2023b) forward-modeled the extinction, distance, and intrinsic parameters of each star given the combination of the Gaia spectra (Carrasco et al. 2021; De Angeli et al. 2023; Montegriffo et al. 2023) and infrared photometry from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), and unWISE, a processed catalog based on the observations from NASA’s Wide-field Infrared Survey Explorer (WISE) mission (Wright et al. 2010; Schlafly et al. 2019). The model was trained using a subset of stars observed with higher spectral resolution, R ~ 1800, available from the Large Sky Area Multi-object Fibre Spectroscopic Telescope (LAMOST; Wang et al. 2022; Xiang et al. 2022).

The Edenhofer et al. (2024) 3D dust maps achieve a compromise between angular resolution and volume coverage, alleviating the limitations of previous 3D extinction models, which can be roughly divided into two groups. First, there are Cartesian reconstructions, which commonly feature fewer artifacts produced by the smearing of extinction structures along the line of sight (“fingers of god”) but which are either limited to covering a small volume at high resolution (see, for example, Leike & Enßlin 2019; Leike et al. 2020) or a large volume at a low resolution (see, for example, Capitanio et al. 2017; Lallement et al. 2019; Vergely et al. 2022). Second, there are spherical reconstructions, which have a much higher angular resolution and probe large volumes of the Galaxy but feature more pronounced finger-of-god artifacts (see, for example, Chen et al. 2019; Green et al. 2019). Other approaches, such as using many small reconstructions (e.g., Leike et al. 2022), an analytical approach (e.g., Rezaei Kh. et al. 2017), or inducing-point methods (e.g., Dharmawardena et al. 2022), have so far been unsuccessful in modeling dust distributions at high resolution over large volumes without artifacts.

Edenhofer et al. (2024) obtains a spherical-coordinate reconstruction beyond 1 kpc while still resolving nearby dust clouds at parsec-scale resolution by implementing a new Gaussian process (GP) methodology to incorporate smoothness in a spherical coordinate system, mitigating fingers-of-god artifacts. The authors modeled the 3D distribution of differential extinction for stars in the Zhang et al. (2023b) catalog, assuming that the dust extinction distribution is spatially smooth. The posterior of their extinction model was reconstructed using variational inference and Gaussian variational inference (MGVI, Knollmüller & Enßlin 2019; Leike & Enßlin 2019). We refer to Edenhofer et al. (2024) for further data modeling and reconstruction details.

The result of the Edenhofer et al. (2024) model is a set of 12 samples drawn from the variational posterior of the 3D dust extinction distribution; that is, a set of 12 possible distributions consistent with the observations within the uncertainties. Each sample gives the value of the dust density in voxels arranged as a series of HEALPix¹ spheres. There are 516 of these spheres, each corresponding to a distance logarithmically spaced grid between 69 to 1250 pc. Each sphere has a resolution parameter N_side = 256, corresponding to an angular size of around 13′7 for each voxel. We analyzed each of the 3D dust extinction samples, considering them independent realizations of 3D dust extinction, and reported the mean trends, as described in detail in Appendix A.

Using the tools in the Python healpy package, we produced a Cartesian projection of the 3D dust extinction covering the range −180° < l < 180° and −5° < b < 5° with a pixel size δl = δb = $7_{\cdot }^{\mathrm{\prime }}5$. The reconstruction was made in terms of the unitless extinction defined by Zhang et al. (2023b), which we transformed into V-band extinction (A_V) by multiplying with the 2.8 factor derived from their extinction curve. We then transformed A_V into hydrogen nucleon column density (N_H) by multiplying by the mean extinction per H nucleon factor 5.8 × 10²¹ cm⁻²/mag from Bohlin et al. (1978) and dividing by the total-to-selective extinction ratio, R_V ≈ 3.1 (Savage & Mathis 1979). Finally, we computed the hydrogen nucleon density (n_H) by dividing by the width of each distance channel.

Figure 3 shows a face-on view of the 3D dust distribution across the studied regions. The pixelization in this representation corresponds to the distance cells in which we divide the studied volume for the analysis presented in Sect. 3. For reference, we indicate in the figure the position of the large-scale features identified in the Solar neighborhood (Zucker et al. 2023). They are the Sagittarius spur (Kuhn et al. 2021), the Cepheus Spur (Pantaleoni González et al. 2021), the Split (Lallement et al. 2019), the Radcliffe Wave (Alves et al. 2020), and the Local Spiral Arm, as defined in Reid et al. (2019). Although our analysis does not assume any prior information on the presence of these features, they are relevant for interpreting our results.

Fig. 3 Heliocentric distribution of the mean nucleon density derived from the 3D dust extinction models presented in Edenhofer et al. (2024) for the sky region |b| ≤ 5°. The superimposed curves indicate large-scale features in the Solar neighborhood, as reported in Zucker et al. (2023). For the Radcliffe Wave, we indicate the full extent of the structure with the dashed white line and the segments within |b| ≤ 5° in green. The white disks indicate the locations of MC complexes referenced in the text.

2.2 Carbon monoxide (CO) emission

We employed the ¹²CO(J = 1 → 0) emission maps presented by Dame et al. (2001), which cover the whole Galactic plane and have an angular resolution comparable to the 3D dust data. This dataset is a combination of the observations obtained over two decades with two 1.2-meter-aperture telescopes: one at Columbia University in New York City, and later in Cambridge, Massachusetts, and one at the Cerro Tololo Inter-American Observatory in Chile. These observations have an angular resolution of $8_{\cdot }^{\mathrm{\prime }}5$ at 115 GHz , the frequency of the ¹²CO(J = 1 → 0) line. This study used the dataset covering the whole Galactic plane within the Galactic latitude range |b| ≤ 5° with 1.3-km s⁻¹-wide spectral channels. We used the raw dataset, which is not interpolated in the spatial coordinates or the spectral axis.

The noise level throughout the Dame et al. (2001) data is not uniform, as it comprises the combination of surveys with different instruments acquired at various times. This is discussed in Appendix A of Miville-Deschênes et al. (2017), where the authors identified three peaks in the noise distribution at 0.06 , 0.10 , and 0.19 K. We conservatively adopt the latter as the global noise level for this dataset.

We used the astropy reproject package to project this data into the same spatial grid of the 3D dust (Astropy Collaboration 2018). We also applied the astropy spectral-cube² package to project the spectral axis of these observations into that of the H I observations introduced below. This spectral reprojection does not affect the results of our analysis; it was performed for convenience in treating this large dataset. An example of the resulting CO emission data is shown in the middle panel of Fig. 1.

2.3 Neutral atomic hydrogen (H I) emission

We employed the publicly available H I 21-cm-wave-length line observations in the H I 4π (H I4PI) survey (HI4PI Collaboration 2016). This survey is based on data from the Effelsberg-Bonn H I Survey (EBHIS, Kerp et al. 2011) and the Galactic All-Sky Survey (GASS, McClure-Griffiths et al. 2009; Kalberla et al. 2010). It comprises observations over the whole sky in the range −600 < v_LOS < 600 km s⁻¹ for declination δ > 0° and −470 < v_LOS < 470 km s⁻¹ for δ < 0°, as observed with the Effelsberg 100-m radio telescope in Bad Münstereifel, Germany and the 64-m radio telescope at Parkes, New South Wales, Australia.

The HI4PI observations have complete spatial sampling, thus overcoming the central issue of pioneering whole-sky H I observations in the Leiden/Argentine/Bonn (LAB) survey (Kalberla et al. 2005). The final HI4PI data product is a set of whole-sky H I maps with an angular resolution of ${16}_{\cdot }^{\mathrm{\prime }}2$ and sensitivity of 43 mK per 1.29-km s⁻¹ velocity channel.

We used the data distributed in FITS-format binary tables containing lists of spectra sampled on a HEALPix grid with N_side = 1024, which corresponds to a pixel angular size of $3_{\cdot }^{\mathrm{\prime }}44$ (Górski et al. 2005). We arranged these spectra in all-sky HEALPix maps corresponding to the 1.29-km s⁻¹-wide velocity channels. Using the tools in the Python healpy package, we produced a Cartesian projection for each one of these velocity channels into the same spatial grid used for the 3D dust data. An example of the resulting H I emission data is shown on the bottom panel of Fig. 1.

3 Methods

3.1 Gradient comparison

We matched the H I and CO emission across velocity channels and the 3D dust cube using the histogram of oriented gradients (HOG) method introduced in Soler et al. (2019). This method is based on characterizing the similarities in the emission distribution using the orientation of its gradients. In the HOG method, a LOS-velocity channel map from the line emission and a distance channel from the 3D dust are morphologically similar if their gradients are mainly parallel and dissimilar if they are randomly oriented. The distribution of angles between gradient vectors is evaluated using the projected Rayleigh statistic (V), a statistical test of non-uniformity in an angle distribution around a particular direction (see Jow et al. 2018, and references therein).

We implemented HOG as follows. We split the sky band within |b| < 5° into 36 10° × 10° regions, identified by the index k. For each region, we used the PPV cube derived from the line emission by tracer $X,I_{ijkp}^X$, and the nucleon density PPD cube derived from the 3D dust reconstruction n_ijkq. The i and j indexes run over the pixels in the sky coordinates, which in our case are Galactic longitude (l) and latitude (b), and the indexes p and q run over LOS velocity and distance channels, respectively. We calculated the relative orientation angles between the emission and density gradients, $${\theta }_{ijkpq}^X=\arctan \left(\frac{\parallel \mathrm{\nabla }{I}_{ijkp}^{\mathrm{X}}\times \mathrm{\nabla }{n}_{ijkq}\parallel }{\mathrm{\nabla }{I}_{ijkp}^{\mathrm{X}}\cdot \mathrm{\nabla }{n}_{ijkq}}\right),$$(1)

where ∇ is the differential operator corresponding to the gradient in the spatial sky coordinates l and b and ||. . .|| indicates the vector norm.

We computed the gradients using Gaussian derivatives resulting from the image’s convolution with the spatial derivative of a two-dimensional Gaussian function. The width of the two-dimensional Gaussian determines the area of the vicinity over which the gradient is calculated. Varying the width of the Gaussian derivative kernel enables the sampling of different scales and reduces the effect of noise in the pixels (see, Soler et al. 2013, and references therein).

3.2 Quantifying the morphological correlation

We used Eq. (1) to calculate the relative orientation angles in the range (−π, π], thus accounting for the direction of the gradients. The values of ${\theta }_{ijkpq}^X$ are only meaningful in regions where both $\left|\mathrm{\nabla }{I}_{ijkp}^{\mathrm{X}}\right|$ and |∇n_ijkq| are greater than zero or above thresholds that are estimated according to the noise properties of each emission and 3D dust cube, as further discussed in Appendix A. We synthesized the information contained in ${\theta }_{ijkpq}^X$ by summing over the spatial coordinates, indexes i and j, using the projected Rayleigh statistic (Jow et al. 2018), which we defined as $${\left(V_{\mathrm{d}}\right)}_{kpq}=\frac{\sum _{ij}{w}_{ijkpq}\cos \left({\theta }_{ijkpq}\right)}{{\left(\sum _{ij}{w}_{ijkpq}/2\right)}^{1/2}}.$$(2)

This definition differs from that used by Soler et al. (2019), where the gradients’ orientation (and not their direction) was considered. To distinguish our definition from that provided by Soler et al. (2019), we introduced the subscript “d”, for direction. This modification improves the significance when comparing 3D dust and CO, where only parallel gradients are meaningful in the column density range considered in the 3D dust reconstruction. A comparison between the two metrics is presented in Appendix A.1.

The projected Rayleigh statistic (either V or V_d) tests nonuniformity in a distribution of angles around a particular direction. In Eq. (2), the angles of interest are θ₀ = 0° and 180°, such that V_d > 0 or V_d < 0 correspond to clustering around those angles (Durand & Greenwood 1958). Values of V_d > 0, which imply that the gradients are primarily parallel, quantify the significance of the morphological similarity between the line emission, I^X, and the density, n, for a pair distance-v_LOS channels. Values of V_d < 0, which imply that the gradients are mostly antiparallel, are potentially relevant toward HISA features, which are characterized by an anti-correlation between the dust density and the H I emission intensity.

The projected Rayleigh statistic’s null hypothesis is that the angle distribution is uniform. In the particular case of independent and uniformly distributed angles and for a large number of samples, values of V_d ≈ 1.64 and 2.57 correspond to the rejection of the null hypothesis with a probability of 5% and 0.5%, respectively (Batschelet 1972). Thus, a value of V_d ≈ 2.87 is roughly equivalent to a 3σ confidence interval. Similarly to the χ²-test probabilities, V_d and its corresponding null hypothesis rejection probability are reported in the classical circular statistics literature as tables of “critical values”, as, for example, in Batschelet (1972), or computed in the circular statistics packages, such as circstats in astropy (Astropy Collaboration 2018).

3.2.1 Selection of HOG parameters

There are two main parameters in HOG calculation. The first is the pixel vicinity for the finite-differences gradient calculation, that is, the angular scale at which the gradients are calculated. The second is the statistical weight in Eq. (2), which accounts for the statistical dependence of the gradients.

The derivative kernel size (Δ) sets the area over which the finite differences derivative is calculated. Varying the size of the derivative kernel enables the sampling of different angular scales and reduces the effect of noise in the pixels. For this particular application, we chose Δ = 30′, which roughly Nyqvist-samples the observations and provides 40 × 40 independent relative orientation measurements per tile. The results of different kernel size selections are discussed in Appendix A.1. The size of the derivative kernel sets a common angular scale for comparing the line emission and 3D dust channels. Thus, we did not smooth the input data to the same angular resolution for the HOG computation.

We accounted for the spatial correlations introduced by the scale of the derivative kernel by choosing the statistical weights w_ijkpq = (δx/Δ)², where δx is the angular size of the pixel and Δ is the derivative kernel’s FWHM. For pixels where the gradient norm is negligible or can be confused with the signal produced by noise, we set w_ijkpq = 0. If all the gradients in an image pair are negligible, Eq. (2) has an indeterminate form that is treated numerically as Not a Number (NaN).

3.2.2 Handling of errors in the HOG morphological correlation and LOS velocity estimation

We estimated the errors in V_d by applying the HOG method to the line emission cubes and the twelve 3D dust distribution realizations produced by Edenhofer et al. (2024) from the variational posterior. The result is twelve V_d values for each line emission and 3D dust channel map. Additionally, we analyzed 100 realizations of the line emission maps produced with Monte Carlo sampling, assuming a Gaussian distribution centered on the observed value in each channel and a standard deviation equal to the noise level of each data set. We reported the mean value of V_d for the twelve 3D dust cubes and the 100 realizations.

The V_d standard deviation estimated from the 12 × 100 samples, σ_Vd, is a measurement of variance in the morphological correlation. In Appendix A.2, we show that σ_Vd is led by the variance among the twelve 3D dust realizations rather than by the noise in the line emission observations. Given that each of the twelve realizations is equally valid, large σ_Vd value reflects an ambiguity in the 3D dust morphology that limits the significance of the morphological matching with a line emission channel. Thus, we exclude distance-v_LOS pairs with |V_d/σ_Vd| < 3.0 from the calculations of the representative v_LOS in each distance channel, as described in Sect. 3.3. We also exclude distance-n_LOS pairs with V_d below the chance correlation thresholds obtained when flipping one of the input maps in the vertical, horizontal, and diagonal directions, as described in Appendix A.2.3.

Fig. 4 Example of the HOG method morphological correlation between the 3D dust density reconstruction and the line emission observations for the region 30° < l < 40° and |b| < 5°. Each panel corresponds to the correlation between distance channels, from the 3D dust density, and velocity channels, from the H I and CO line emission, shown in the left and right panels, respectively. The correlation metric is the direction-sensitive projected Rayleigh statistic (Vd, Eq. (2)). Values of Vd ≈ 0 correspond to a random orientation between the gradients of the two tracers, thus indicating low morphological correlation. Values of Vd > 2.87 correspond to mostly parallel gradients, thus indicating a significant morphological correlation. Values of Vd < 2.87 correspond to mostly antiparallel gradients. The white markers indicate the vLOS with the highest Vd for each distance channel. The green dashed lines correspond to the distance and vLOS with the highest Vd for the whole velocity and distance ranges. The three cyan lines represent the kinematic distances for the Galactic longitude at the center and the region’s edges, according to the Reid et al. (2019) Galactic rotation model.

3.3 Dissecting the nearby Galactic plane

We carried out the Galactic plane analysis by splitting the observations within the |b| < 5° range into 36 contiguous 10° × 10° tiles centered on b = 0°. For each tile, we computed Eq. (2) using the HOGcorr_ima routine in the publicly available astroHOG package³. The result is a V_d matrix representing the correlation between each distance and velocity-channel map toward a tile, as illustrated in Fig. 4.

The V_d matrix is a map of the relationship between d and v_LOS for each studied region, something which is usually estimated by assuming circular motions around the Galactic center to compute kinematic distances (see, for example, Oort et al. 1958; Sofue 2011; Wenger et al. 2018; Hunter et al. 2024). In the HOG method, this mapping is based on the morphological similarity between the line emission and the 3D dust density distribution for each d and v_LOS pair. Thus, d and v_LOS are not linked by a one-to-one relation but rather by a distribution of V_d, which accounts for the correlation across LOS velocities introduced by the emission linewidth and the coherence expected in 3D dust structure for contiguous distance channels.

We employed the V_d matrices to assign a v_LOS to each distance channel in a tile. We used the critical value 2.87 and the V_d standard deviation, (σ_Vd)_vk, as significance thresholds for the LOS velocity assignment to each distance channel. If the maximum value of V_d is below either of these values, no v_LOS is assigned to a distance channel. If more than one velocity channel has V_d above these thresholds, we assign the v_LOS corresponding to the maximum V_d. The selection by maximum V_d excludes distance channels potentially dominated by HISA from the v_LOS reconstruction. This choice is motivated by analysis in Appendix B, which shows that a selection based on |V_d| instead of V_d only favors negative values in a small number of the distance channels and that the interpretation of V_d < 0 exclusively as the product of HISA is not always correct. We reserve the specific study of HISA using HOG for a subsequent publication.

The simplification of the velocity field in our HOG application aims to identify a representative LOS velocity for the bulk of the dust in each distance slice along a 10° × 10° area. Thus, we estimated the average LOS motions of volumes between 12 × 12 × 0.4 pc³ and 220 × 220 × 7.0 pc³. We acknowledge that this selection extracts just a segment of the turbulent energy cascade, which connects galactic motions to the smallest scales in the ISM. However, our goal is reconstructing the bulk motions of density structures revealed by the 3D dust model rather than describing the full complexity of the ISM velocity field across scales, which has been considered in other works (see, for example, Elmegreen & Scalo 2004; Hennebelle & Falgarone 2012, and references therein).

Figure 4 illustrates an example of the LOS velocity selection from a V_d matrix, where the markers represent the assigned v_LOS for each distance channel. In this example, the representative LOS motions for CO roughly follow the expected v_LOS from Galactic rotation, indicated by the solid lines. In contrast, the H I v_LOS shows deviations up to tens of km s⁻¹ from the CO v_LOS and the LOS velocities expected from the Galactic rotation model. In the following sections, we discuss this and other dynamical effects on the reconstructed velocity field after presenting the global HOG results.

Fig. 5 Maximum morphological correlation between the 10° × 10° distance channels and the H I (left) and CO (right) line emission in the range −25 < vLOS < 25 km s−1, as quantified by the direction-sensitive projected Rayleigh statistic (Vd; Eq. (2)).

4 Results

4.1 Correlation between 3D dust and line emission morphology

The fundamental hallmark of the HOG analysis is quantifying the similarity in the plane-of-the-sky distribution of the 3D dust extinction and the H I and CO line emission across distance and velocity channels. Figure 5 shows the maximum values of the correlation metric across v_LOS, max(V_d)_v, for the distance channels in the 36 tiles in which we divided the |b| < 5° sky region. That is, each pixel in Fig. 5 corresponds to the value of Eq. (2) for the particular case of p = p*, where p* is the velocity channel with the maximum value of V_d.

The values of max(V_d)_v reported in Fig. 5 indicate whether or not a dust parcel has a match in the H I or CO emission in the range −25 < v_LOS < 25 km s⁻¹. By construction, max(V_d)_v highlights configurations where H I and 3D dust density gradients are preferentially parallel. Although this choice is biased against regions potentially dominated by HISA features, where V_d < 0, it minimizes the spurious signal produced by chance correlation between voids in the 3D dust and H I emission, as further discussed in Appendix B.

Figure 5 shows that most of the distance channels display a significant correlation (V_d > 2.87) between the 3D dust morphology and the line emission. This result is remarkable given that it comes from two pairs of independent datasets, H I and 3D dust and CO and 3D dust. Although the association of 3D dust and line emission components has been previously studied toward particular regions in the Solar vicinity (see, for example, Piecka et al. 2024; Rybarczyk et al. 2024), this is the first time this relation is reported using a quantitative measure of the morphological similarity and covering such an extensive region.

An alternative representation of the correlation between the line emission and the 3D dust is presented in Fig. 6, where we show the maximum values of V_d across distances, max(V_d)_d, for the v_LOS channels in the 36 tiles in which we divided the |b| < 5° sky region. That is, each pixel in Fig. 6 corresponds to the value of Eq. (2) for the particular case of q = q*, where q* is the distance channel with the maximum value of V_d. Figure 6 can be understood as a different projection of the PPDV space mapped with the HOG method, although some of the line emission signal in the |v_LOS| ≤ 25 km s⁻¹ may come from density structures beyond the 3D dust reconstruction distance range. Therefore, Fig. 6 should be understood as a visual representation of whether a line emission channel has a morphological match in the 3D dust distribution between 69 and 1250 pc , rather than as a full reconstruction of the PPV space from PPD information.

It is evident on the right panel of Fig. 6 that the CO emission has the most considerable morphological correlation with the 3D dust within the velocity ranges expected from the circular motion around the Galactic center, which we estimated using the Reid et al. (2019) model. In contrast, the H I emission is highly correlated with the 3D dust beyond the LOS velocity limits expected from pure circular motions. We further discuss the prevalence of these streaming motions in Sect. 4.2.

Figure 7 shows an example of the H I and CO line emission and the density distribution for the tile with the highest morphological correlation between CO and the 3D dust, as identified by the highest V_d in the right-hand-side panel of Figure 5. The noticeable similarity between the CO emission and the density distribution visually confirms the HOG results. Note that the fact we see CO emission associated with gas with a density n_H ~ 10–20 cm⁻³, which we would typically expect to be atomic, is likely an indication that the gas is clumped on scales smaller than the resolution of the 3D dust map at this distance, as confirmed by higher-angular-resolution CO surveys (see, for example, Jackson et al. 2006; Benedettini et al. 2020; Ma et al. 2021).

The H I emission in Fig. 7 shows the typical shadows produced by the cold H I observed on top of the warm background emission, which is a characteristic signature of a HISA feature (see, for example, Gibson et al. 2000). We presume that the lower H I emission within the CO contours in the central panel of Fig. 7 is a HISA based on the morphological similarity between the CO emission and HISA identified in other observations of the Galactic plane (see, for example, Gibson et al. 2005; Soler et al. 2019; Wang et al. 2020). However, a confirmation of the HISA for this example and throughout the Galactic plane requires further analysis of the H I spectra, which is beyond the scope of this paper (see, for example Syed et al. 2023). The presumed HISA in this particular example does not produce a prominent V_d < 0 in the HOG comparison between H I and 3D dust, as shown in Fig. 4. This is most likely due to the combination of antiparallel gradients in the HISA feature and the parallel gradients outside it, resulting in V > 0.

Figure 8 shows the H I and CO line emission and density distribution for the tiles with the highest V_d in the H I and 3D dust comparison, as identified in the left-hand-side panel of Fig. 5. In contrast with the region in Fig. 7, there is no CO counterpart to the matching structures in H I and 3D dust. This example illustrates that the morphological correlation traced by H I and CO is not necessarily identical because of the different structures each tracer shows.

We further compared the morphological correlations of the 3D dust and the two gas tracers in Fig. 9. The scatter plot contrasts the values of max(V_d)_v for H I and CO reported in Fig. 5. There are some tiles where max(V_d) is high for the two gas tracers. However, the correlation between the 3D dust and one gas tracer does not indicate its correlation with the other.

We found that the highest max(V_d)_v are not concentrated around near or far distances. We interpret this observation as an indication that no significant bias is introduced by distance in the morphological correlation, as further considered in Appendix C. The highest max(V_d)_v is found for H I , as expected from the largest angular extend of H I emission.

Figure 9 shows a few tiles where max(V_d)_v is negative, most of them in H I. This suggests that despite our selection of the maximum V_d in each tile, there are regions where the H I emission gradients are mostly antiparallel to the 3D dust gradients. However, their significance is low, V_d < −2.87, and they are related to chance correlation rather than HISA, as discussed in Appendix B.

Fig. 6 Maximum morphological correlation between the 10° × 10° LOS velocity channels and the dust density in the range 69 < d < 1250 pc, as quantified by the direction-sensitive projected Rayleigh statistic (Vd; Eq. (2)). The cyan curves represent the vLOS expected for d = 69 and 1250 pc according to the Reid et al. (2019) Galactic rotation model, shown by dotted and dashed lines, respectively.

Fig. 7 Example of the distance and CO velocity channel with high morphological correlation. It corresponds to the highest Vd in the comparison between the CO line emission and the 3D dust presented in the right panel of Fig. 5. From left to right, the panels show the CO and H I line emission at the indicated vLOS and the nucleon density inferred from the 3 D dust model at the corresponding distance.

Fig. 8 Same as Fig. 7, but for high morphological correlation between H I and 3D dust. This example corresponds to the highest Vd in the comparison between the H I line emission and the 3D dust density presented in the left panel of Fig. 5.

Fig. 9 Comparison of the maximum Vd values for H I and CO reported in Fig. 5.

Fig. 10 Expected line-of-sight velocity (right) from the Reid et al. (2019) Galactic rotation model across the cells in our analysis.

4.2 Local ISM motions

We used the regions with high morphological correlation (V_d > 2.87) to assign a prevalent velocity to the dust parcels within the 3D dust reconstruction, following the procedure described in Sect. 3.3. Assuming circular rotation around the Galactic center implies a likely LOS velocity pattern in our studied region. For example, Fig. 10 shows the expected v_LOS pattern derived from the state-of-the-art Galactic rotation model presented in Reid et al. (2019).

Figure 11 shows the velocity field reconstructed from 3D dust and line emission correlation identified using the HOG method. The regions without assigned v_LOS correspond to distance channels with V_d < 2.87 for all velocity channels in either tracer. The first compelling result from the velocity reconstruction is the large-scale similarity between the velocity fields in Fig. 11 and the theoretical expectation shown in Fig. 10, that is, the correspondence between the motions away from the Sun in the first and third Galactic quadrants and toward the Sun in the second and fourth quadrants. The correspondence of this quadrupolar pattern in the model and our reconstruction is evident in both gas tracers, despite the large portion of the first quadrant for which the morphological correlation with H I does not result in an unambiguous v_LOS assignment.

Appendix D presents the HOG analysis of synthetic line emission from a simulated MC, demonstrating that our approach correctly reproduces the mean v_LOS of the numerical cloud. Dispersions around the central v_LOS are dominated by the emission linewidth. Larger velocity excursions registered with the HOG are produced by the MC’s complex dynamic structure rather than by a systematic effect introduced by the method.

We quantified the departures from the LOS motions produced by pure Galactic rotation by subtracting the expected LOS velocity from the Reid et al. (2019) model, $v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{\mathrm{R}19}$, from the reconstructions presented in Fig. 11, $v_{\text{LOS }}^{\mathrm{H}\mathrm{O}\mathrm{G}}$. In what follows, we refer to the quantity $v_{\text{LOS }}^{\text{R19 }}-v_{\text{LOS }}^{\text{HOG }}$ as streaming motions.

The streaming motions distribution, shown in Fig. 12, is not strictly symmetric but is centered around zero, indicating that, on average, the departures from pure Galactic rotation are relatively small. The mean values streaming motions are roughly 0.2 and 1.1 km s⁻¹ for H I and CO , respectively. These values are below the 1.3-km s⁻¹ spectral resolution of the line emission data, so they are within the LOS velocity determination uncertainties. The characteristic amplitude of the streaming motions, quantified by the standard deviation of $v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{\mathrm{R}19}-v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{\mathrm{H}\mathrm{O}\mathrm{G}}$, is roughly 10.8 and 6.6 km s⁻¹, for H I and CO , respectively. Given the |v_LOS| < 25 km s⁻¹ input range limitation, large excursions from Galactic rotation are currently excluded from the reconstruction, so these dispersions may represent lower limits of the actual streaming motions.

The streaming motions map, presented in Fig. 13, does not show a prevalent large-scale pattern of streaming motions that may indicate the effect of compression toward the Local Arm. Nor does it reveal global agreement in the radial velocity pattern sampled by the H I and CO line emission, indicating that the residual motions due to Galactic rotation are relatively small. In some positions, for example, at d ~ 800 pc for 80 < l < 90°, the H I shows diverging motions that the CO does not match. However, in general, the radial velocities do not reveal the concentration of converging or diverging motions but rather a succession of compressions and rarefactions along the line of sight, as expected from the behavior of a turbulent medium. Given the scales considered in this study, averaging of cells that range between 12 pc × 12 pc × 0.4 pc and 218 pc × 218 pc × 7 pc in size, it is plausible that the reported motions result from a turbulent energy cascade introduced at kiloparsec scales by global galactic motions (see, for example, Colman et al. 2022, and references therein).

Fig. 11 Line-of-sight velocity derived from the H I (left) and CO (right) line emission associated to the distance channels in the 3D dust reconstruction using the HOG method $(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{\mathrm{H}\mathrm{O}\mathrm{G}})$. The colored circles correspond to the position and vLOS for the five high-mass star-forming regions with VLBI parallax observations within the studied volume.

Fig. 12 Histograms for the differences between the LOS velocity estimated from the H I and CO line emission, $v_{\text{LOS }}^{\mathrm{H}\mathrm{O}\mathrm{G}}$, as reported in Fig. 11, and the expected LOS velocity according to the Galactic rotation model presented by Reid et al. (2019), $v_{\text{LOS }}^{\mathrm{R}19}$, presented in Fig. 10.

Fig. 13 Differences between the line-of-sight velocity derived from the H I (left) and CO (right) line emission associated with the distance channels in the 3D dust reconstruction using the HOG method, $v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{\mathrm{H}\mathrm{O}\mathrm{G}}$, and the expected LOS velocities according to the Galactic rotation model presented by Reid et al. (2019), $v_{\text{LOS }}^{\mathrm{R}19}$.

5 Physical quantities derived from the local ISM motions

The results presented in Sec 4 indicate a significant correlation between the density traced by the 3D extinction models and the line emission. We used this fact to estimate the distribution of departures from circular motions, also known as streaming motions. Using these streaming motions, we calculated the distributions of kinetic energy density, momentum density, and mass flow rates associated with each gas tracer, as reported in Fig. 14. We calculated these physical quantities as follows.

5.1 Effective densities

In the standard, line-emission-based PPV approach, the density and LOS velocity of an object are obtained by isolating a line component, identifying its centroid, and integrating the emission (see, for example, Miville-Deschênes et al. 2017). Thus, the density and the LOS motion are defined for the same PPV object. Our analysis defines the object in PPD space as a slice of 3D dust density for a particular sky area. We assigned a v_LOS to that object using the HOG method through morphological matching. However, not all the 3D dust slice is represented by the v_LOS from one tracer. For example, the HI and CO emission can match different portions of a 3D dust slice, and each portion can have different velocities, as illustrated in Fig. 4. Thus, we introduce an “effective density” corresponding to the portions of the 3D dust slice with a significant morphological match with either gas tracer. This approach has the advantage of separating the H I and the CO motions and assigning them a portion of the 3D dust instead of using the whole slice, which would lead to an overestimation of the energy and momentum densities.

We employed the map of relative orientation angles between gradients for each 3D density slice and v_LOS-channel pair. We split this map into 9 × 9 blocks, computed V_d for each block, and calculated the mass contained in the blocks with V_d > 1.0, which roughly corresponds to a 1-σ significance for V_d. Given the limited number of independent gradient vectors in each block, the threshold on V_d is necessarily less rigid than that we used for the whole tile. Currently, the HOG analysis is limited by the angular resolution in the gas observations and the 3D dust reconstruction, which ultimately define the number of independent gradient vectors used in the morphological comparison. More stringent values lead to null values of the effective density for most of the distance channels.

The effective density in each density slice is determined by dividing the total mass in the blocks with V_d > 1.0 by the volume of the slice. Explicitly, for each distance channel q in the portion of the Galactic plane identified with the subindex k, we identified a velocity channel p* that has the highest V_d for the gas tracer X. For that distance channel, we compute the effective nucleon density, $$n_{k{p}^{\ast }q}^{\mathrm{e}\mathrm{f}\mathrm{f},X}=\frac{\sum _{ab}^{N_A,N_B}\left({\omega }_{ab}\sum _{i^{\mathrm{\prime }}{j}^{\mathrm{\prime }}}^{N_{i^{\mathrm{\prime }}},N_{j^{\mathrm{\prime }}}}{n}_{i^{\mathrm{\prime }}{j}^{\mathrm{\prime }}}{\mathcal{V}}_{i^{\mathrm{\prime }}{j}^{\mathrm{\prime }}kq}\right)}{{\mathcal{V}}_{kq}},$$(3)

where n is the nucleon density derived from the 3D dust distribution model, the indices a and b run over the blocks, i′ and j′ run over the pixels within the block, and _i′j′kq is the volume occupied by the i′-th and j′-th pixel at the distance $d_{kq}^2$. These definitions imply that the total volume in the distance channel q is $${\mathcal{V}}_{kq}={\displaystyle \sum _{i,j}^{N_i,N_j}}{\mathcal{V}}_{ijkq}={\displaystyle \sum _{i,j}^{N_i,N_j}}{d}_{kq}^2\tan (\mathrm{\Delta }l)\tan (\mathrm{\Delta }b)\mathrm{\Delta }{d}_q,$$(4)

where Δl and Δb are the angular sizes of the pixel and Δd_q is the distance channel width along the LOS. The statistical weight ω_ab is equal to one if V_d > 1.0 in the block identified by the indices a and b and zero otherwise.

The top-left panel of Fig. 14 shows the effective density distribution across distance cells. By construction, the effective densities for the gas tracers are below the total dust mean. While the mean density from all the dust is around 0.51 cm⁻³, the effective densities are around 0.27 and 0.14 cm⁻³ for H I and CO, respectively. The lower value for CO reflects that tracer’s lower volume filling factor, which implies a smaller contribution than that from H I to the kinetic energy and momentum for each volume element in our reconstruction.

We note that the effective densities calculated using Eq. (3) are not the specific density of a gas tracer but a combination of its volume filling factor and a selection based on its correlation with the 3D dust density. For example, CO is associated with densities roughly above 100 cm⁻³, but the Edenhofer et al. (2024) 3D dust models do not have the spatial resolution to detail the regions where these densities are prevalent. The value of n^eff does not aim to estimate the specific volumetric density of the gas but rather identify the portion of the 3D dust slice associated with the motion of the gas tracer.

The fundamental uncertainties in estimating the densities come from the limitation of the 3D extinction reconstruction, described in detail in section 6 of Edenhofer et al. (2024). We converted the extinction into nucleon density using the value from Bohlin et al. (1978), which recent observations of H I emission and interstellar reddening indicate is too low by a factor of 1.5 (Lenz et al. 2017). Variations from the diffuse environment at high Galactic latitude to the denser environment in the plane are also relevant, adding a factor of a few to the uncertainty. Recent observations also reveal spatial variability of R_V characterized by a standard a standard deviation of around 0.25 (Zhang et al. 2023a).

The HOG analysis of synthetic line emission from a simulated MC presented in Appendix D indicates that the reconstructed H I densities are roughly within a factor of two from the actual values. The spatial averaging in Eq. (3) results in smearing the density peaks across the area of a block, effectively limiting the maximum reconstructed density values. For CO , the comparison with the numerical model is limited by the extent of the synthetic diffuse CO emission, which is a complex chemical modeling problem beyond this work’s scope. However, the HOG method reproduces the CO density profile, although it over-estimates the CO density by a factor of a few. In either case, the limitations in the HOG reconstruction are below the variance across the studied region.

Fig. 14 Histograms of the physical quantities derived from the HOG analysis. Top left: Effective energy density, neff, that is assigned to each gas tracer, as defined in Eq. (3). Top right: Kinetic energy density, Ek, calculated with Eq. (5). Bottom left: Radial momentum, p, calculated with Eq. (7). Bottom right: Mass flow rate, $\dot{M}$, calculated with Eq. (8).

5.2 Kinetic energy

We estimated the average kinetic energy from streaming motions for each distance channel q as $${⟨E_{\mathrm{K}}^X⟩}_{kq}=\frac{1}{2}{\rho }_{k{p}^{\ast }q}^{\mathrm{e}\mathrm{f}\mathrm{f},X}{[{\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^X\right)}_{k{p}^{\ast }q}-{\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{R19}\right)}_{kq}]}^2,$$(5)

where X corresponds to the velocity tracer, either H I or CO, ${\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{R19}\right)}_{kq}$ is the expected LOS velocity from the Reid et al. (2019) Galactic rotation model for the Galactic plane position k and the distance channel q, and ${\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^X\right)}_{k{p}^{\ast }q}$ is the LOS velocity from the line emission channel p* with the highest morphological correlation (V_d) with the distance channel q. We calculated the effective volumetric density, $${\rho }_{k{p}^{\ast }q}^{\mathrm{e}\mathrm{f}\mathrm{f},X}=\mu m_p{n}_{k{p}^{\ast }q}^{\mathrm{e}\mathrm{f}\mathrm{f},X},$$(6)

using Eq. (3), the proton mass, m_p, and μ is the mean nucleon mass factor, which for solar metallicity is 1.402 (Asplund et al. 2009; Draine 2010).

The E_K distribution is shown in the top-right panel of Fig. 14. We found that the average E_K is around 0.11 and 0.04 eV/cm³ for H I and CO, respectively. These estimates are comparable to the 0.22 eV/cm³ obtained for the turbulent kinetic energy in the local ISM, obtained with a mean nucleon density n = 30 cm⁻³ and velocity v = 1 km s⁻¹, or mean values ⟨n_H⟩ = 1 cm⁻³ and σ_v = 5.5 km s⁻¹ (Draine 2010). Given that we are only sampling one component of the velocity field, the actual E_K value can be up to a factor of three higher.

The E_K standard deviations are around 1.95 and 0.53 eV/cm³ for H I and CO, respectively. These values, roughly a factor of four above the mean, suggest significant fluctuations of this quantity within the 2.5-kpc-diameter ISM parcel considered in this paper. The maximum E_K values are around 100 and 30 eV/cm³ for H I In general, the highest E_K are associated with the highest n^eff, suggesting that density is the dominant factor in our reconstruction of the kinetic energy distribution.

Figure 15 presents the E_K maps obtained from the H I and CO emission. We found no evident similarity in the E_K distribution in the two gas tracers. However, for both tracers the highest E_K are found in the regions at d < 300 pc.

The E_K radial profile, presented in the top panel of Fig. 16, shows values up to ten times over the mean for regions outside d ≈ 300 pc and up to 20 times for regions inside that range. The E_K peaks at d < 300 pc do not seem to result from a concentration of high V_d due to a bias in the morphological correlation toward near distances, as discussed in Appendix C. However, the angular resolution in the 3D dust reconstruction may limit the mapping of E_K enhancements to further distances. One potential explanation for the E_K peak emerging at d ≈ 100 pc is the LSR location within the Local Bubble. We are observing the gas motions from a location near the center of that cavity so that the peak is at a similar radius in all directions. This is less likely to be the case for material that is further away from us.

The E_K azimuthal profile, obtained by adding the energy along the LOS, is presented in the bottom panel of Fig. 16, It shows peaks toward l ≈ 40 to 90° for H I and CO. This range covers the general directions along the axis of two large-scale structures identified in 3D dust reconstructions: the Split (Lallement et al. 2019) and the Radcliffe wave (Alves et al. 2020), shown in Fig. 3. Additional overdensities are found toward the fourth quadrant, roughly between l ≈ 300 to 340°.

Fig. 15 Kinetic energy density (Ek) derived from the velocity field reconstruction in Fig. 13 and n derived from the 3D dust extinction modeling, as detailed in Eq. (5).

Fig. 16 Kinetic energy density, EK, radial and azimuthal profiles normalized to the mean values over the studied region.

Fig. 17 Momentum density (p) obtained with the LOS velocity reconstructions in Fig. 13 and Eq. (7).

5.3 Momentum density and momentum

We estimated the LOS momentum density from the streaming motions for each distance channel q as $${\left(p_{\mathrm{K}}^X\right)}_{kq}={\rho }_{k{p}^{\ast }q}^{\mathrm{e}\mathrm{f}\mathrm{f},X}[{\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^X\right)}_{k{p}^{\ast }q}-{\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{R19}\right)}_{kq}],$$(7)

where the LOS velocities ${\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^X\right)}_{k{p}^{\ast }q}$ and ${\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{R19}\right)}_{kq}$ are defined in the same way as in Eq. (5). The effective density for tracer X, ${\rho }_{k{p}^{\ast }q}^{\text{eff,}X}$, is estimated using Eq. (6).

The distribution of ${\left(p_{\mathrm{K}}^X\right)}_{kq}$ is shown in the bottom-left panel of Fig. 14. The volume-weighted mean is around 5.9 and −0.6 × 10⁻³ M_⊙ km s⁻¹/pc³ for H I and CO, respectively. However, there are large excursions from these mean trends, characterized by standard deviations of roughly 0.8 and 0.2 M_⊙ km s⁻¹/pc³. This implies that the local fluctuations are much more important than the mean values.

Figure 17 shows the p maps derived from the radial motions sampled by the two gas tracers. The peak momentum density in H I , found toward l ≈ 345° at d ≈ 190 pc, corresponds to a momentum around 6.1 × 10³ M_⊙ km s⁻¹. This value can be interpreted by comparison with the radial momentum input from a single supernova remnant (SNR), identified as roughly between 1 and 15 × 10⁴ M_⊙ km s⁻¹ in numerical simulations of one SNR expansion in an inhomogeneous density field for different ambient medium’s densities, metallicities (Martizzi et al. 2015). These values should be taken just as a reference; determining whether SNRs produce the reconstructed LOS momentum requires considering stochastic parameters such as initial conditions, efficiency in the conversion between thermal and kinetic energy, and the projection of momenta along the LOS. Our reconstruction aims to characterize the distribution of LOS momentum density and identify momentum sources in the local galactic plane. Reconstruction of formation scenarios for individual regions is beyond the scope of this work.

It is apparent from Fig. 17 that there is a concentration of high p values around d ≈ 200 pc, both in H I and CO. For regions beyond that heliocentric radius, the highest |p| values appear toward 40 < l < 90°, roughly corresponding to the directions of the Split and the Radcliffe wave. There are also high p values in H I and CO toward l ≈ 270° for 800 < d < 1100 pc, roughly around the Vela C location indicated in Fig. 3, and toward 300 < l < 360° for d > 800 pc.

5.4 Mass flow rates

Finally, we considered the mass flow rates, $\dot{M}$, corresponding to the motion of the material in each distance channel across its radial width. For that purpose, we calculated $${\dot{M}}_{kq}={\rho }_{k{p}^{\ast }q}^{\mathrm{e}\mathrm{f}\mathrm{f},X}{\mathcal{V}}_{kq}[{\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^X\right)}_{k{p}^{\ast }q}-{\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{R19}\right)}_{kq}]/\mathrm{\Delta }{d}_q,$$(8)

where the radial velocities ${\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^X\right)}_{k{p}^{\ast }q}$ and ${\left(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{R19}\right)}_{kq}$ are defined in the same way as in Eq. (5), the effective density for tracer X, ${\rho }_{k{p}^{\ast }q}^{\text{eff,}X}$, is estimated using Eq. (6), _kq is the distance channel volume, and Δd_q is the distance channel width.

The $\dot{M}$ distribution, shown on the bottom-right panel of Fig. 14, is approximately centered around zero, implying that the LOS net mass flow for the studied region is close to zero. The $\dot{M}$ standard deviations are around 1.1 × 10⁻³ and 0.4 × 10⁻³ M_⊙/yr for H I and CO, respectively. The azimuthal $\dot{M}$ profiles are presented in Fig 19. They correspond to the mean and the standard deviations of $\dot{M}$ along the LOS for H I and CO. The mean mass flows, displayed in Fig 19, show net matter displacements of roughly couple 10⁻³ M_⊙/yr for lines of sight toward the first Galactic quadrant and 250 < l < 360°. This could be interpreted as residual velocity from the Galactic rotation if not because it has the opposite sign of the expected LOS motions indicated on the right-hand-side panel of Fig. 11. Further interpretation of these mass flows should consider that they may be a product of the LSR position and the boundaries of the 3D dust reconstruction.

The standard deviation of $\dot{M}$, shown in Fig 19, accounts for the amplitude of the mass flow fluctuations along the line of sight. Two prominent peaks in this quantity traced by H I appear toward l ≈ 40 and 80°, roughly the directions along the Split and the Radcliffe wave. Additional peaks are found toward l ≈ 270 and 340°, coinciding with the locations of momentum overdensities in Fig. 17. These peaks do not correspond to the $\dot{M}$ profiles derived from CO.

Fig. 18 Mass flow rates ($\dot{M}$) obtained with the LOS velocity reconstructions in Fig. 13 and Eq. (8).

Fig. 19 Azimuthal profiles of mean mass flow rates and its standard deviation.

6 Discussion

6.1 3D dust and line emission correlation

Figure 5 shows the significant similarity in the morphology of the structures traced by the 3D extinction modeling and the H I and CO line emission. Previous works demonstrated that dust and gas are well-mixed in the diffuse ISM, as inferred by the tight empirical linear relation between the H I column density and dust extinction and emission (see, for example, Burstein & Heiles 1978; Boulanger et al. 1996). However, our results present the first quantitative estimation of the global likeness between the extinction structures and the line emission observations.

We found that roughly 67% and 15% of the 10° × 10° distance channels in Fig. 5 show a significant morphological correlation, max(V_d)_v > 2.87, between the 3D dust and the H I and the CO emission, respectively. Around 9% and 1% of the distance channels present max(V_d)_v < 0. Additional analysis of the max(V_d)_v distributions are presented in Appendix C.

The lower correlation with CO can be explained by the lower volume filling factor of that tracer (Kwan & Sanders 1986; Kalberla & Kerp 2009). Clouds traced by CO are expected to be smaller and more concentrated than H I structures. Hence, CO has fewer independent CO gradients for the HOG morphological comparison on the same sky area and at the common angular scale set by the derivative kernel size.

The lack of correlation in the rest of the channels can be attributed to several factors, which we separated into those related to the HOG method and those intrinsic to the dust and gas distribution in the local ISM. Among the former factors, we can distinguish two cases. First, there can be structures beyond 1.25 kpc with a morphological match with the emission within the input range. But even if they become available with a future 3D dust reconstruction, including them would not affect the results in Fig. 5 and only potentially fill a few gaps in the lv-diagrams presented in Fig. 6. Second, there can be line emission structures in channels beyond the v_LOS input range with a morphological match with the dust within the 3D reconstruction domain. We considered this scenario in Appendix A.3.3, where we reported the results for the input range −120 < v_LOS < 120 km s⁻¹. The expanded v_LOS range emphasizes higher streaming motion amplitudes at the inevitable price of an increased chance correlation. However, it does not patch the low V_d channels in Fig. 5.

The low morphological correlation for certain portions of the studied volume can also be due to the intrinsic dust and gas distribution in the local ISM, which we consider in three main effects. First, there is low extinction in local ISM volumes, such as found in the regions around 600 < d < 800 pc toward −45° < l < 30° and 1000 < d < 1200 pc toward 250° < l < 290°, to mention a couple of examples. Second, line emission is absent across some LOS velocity ranges toward particular lines of sight. For example, toward the GSH 238+00+09 superbubble (Heiles 1998), where there is very little CO, as illustrated in the right panel of Fig. 6. Identifying the exact physical conditions that produce the local ISM under-densities and a lack of emission toward particular regions is beyond the scope of this work. However, the maps in Fig. 5 and Fig. 6 provide a global context for identifying the conditions toward specific regions and their implications for the ISM dynamics at smaller scales.

There is a third global factor that reduces the correlation between the 3D dust and the line emission. Observations of other ISM tracers, such as gamma-rays and the [C II] and OH line emission, reveal that not all of the ISM material is traced by H I and CO emission (see, for example, Grenier et al. 2005; Langer et al. 2014; Pineda et al. 2013; Busch et al. 2021). Including this material in the distance-velocity correlation is not straightforward. In the case of [C II], there is no map with enough angular resolution and coverage to obtain significant spatial correlations using HOG. In the case of OH , the emission distribution involves the complex physical conditions that lead to the emission from that radical, which calls for a separate dedicated study (see, for example, Rugel et al. 2018; Dawson et al. 2022). We acknowledge that those and other line tracers should be accounted for in future reconstructions of the local ISM, adding richness to the results we present in this paper. Yet, our analysis has shown we have enough distance channels where the correlation between the dust distribution and the line emission is significant enough to identify prevalent LOS velocities for different portions of the local ISM, which we discuss in the following section.

6.2 Local ISM motions

We used the morphological correlation across distances and LOS velocities to assign a v_LOS to the distance channels in each of the 10° × 10° regions in which we divided the studied volume. This assignment simplifies the velocity field by identifying a single value to ISM parcels with volumes between 10² to 10⁵ pc³, corresponding to distance channels at the inner and outer edges of the 3D dust model. Thus, our velocity reconstruction represents the local ISM dynamics on spatial scales of tens of parsecs in the plane of the sky and around a parsec in the LOS direction.

The velocity assignment is unambiguous in most of the 10° × 10° regions, as illustrated in Fig. 4 and Fig. 23. There are distance ranges for which the assigned v_LOS shows an oscillatory pattern. We do not interpret these as actual oscillations since they can arise from the prevalence of the correlation measured by V_d alternating between different structures in the same sky area. We do not discard the possibility of oscillations, such as the ones recently identified in works such as Henshaw et al. (2020). Still, its identification requires a dedicated analysis beyond this work’s scope.

The maps of reconstructed LOS velocities shown in Fig. 11 indicate that the local LOS velocity field carries the imprint of the Galactic rotation pattern in both tracers, as inferred from the comparison with the expected v_LOS in the Reid et al. (2019) Galactic rotation model, shown in Fig. 10. The large-scale pattern of H I and CO emission is deeply influenced by Galactic rotation, as shown in figure B.2 of HI4PI Collaboration (2016) and figure 3 of Dame et al. (2001). However, the peculiar velocities of MCs within 3 kpc cannot be explained purely by Galactic rotation (see, for example, Stark & Brand 1989).

The radial velocity of stars from the Gaia DR2 catalog within 1 kpc of the Galactic mid-plane and roughly 4 kpc from the Sun, reported by Katz et al. (2019), show a global pattern consistent with the LOS-projected differential rotation of the stars around the Galactic center, as observed from the Sun. The HOG method enabled mapping this rotational pattern in the gas in the Solar neighborhood by anchoring the LOS velocity to a distance using the morphological similarity between the dust and gas tracers. Figure 11 strongly suggests that despite its collisional nature and being exposed to the influence of local conditions, such as the energy and momentum input from stellar feedback, the gas within 1.25 kpc follows the pattern of Galactic rotation.

We subtracted the Reid et al. (2019) model of Galactic rotation from the LOS velocity pattern in Figure 11 to characterize the deviations from pure Galactic rotation in the local ISM (Burton 1971). We refer to the resulting LOS motions as streaming motions to avoid attributing them to fluctuations in the large-scale turbulent velocity field or local fluctuations introduced, for example, by SN explosions. Separating the mean and the turbulent velocity fields is not straightforward, even when all velocity field components are available. Thus, our characterization of the streaming motions reported in Figure 13 aims to relate them to the known structure of the local ISM rather than a generic separation that may be hindered by the limitations of our reconstruction.

Figure 12 shows that the distribution of LOS streaming motions roughly centered around zero, within the uncertainties imposed by the velocity channel widths in the input data. The standard deviations of the rotation-subtracted LOS velocities are 10.8 and 6.6 km s⁻¹ for H I and CO , respectively. Previous studies have identified similar deviations from Galactic rotation. Clemens (1985) used the Massachusetts-Stony Brook Galactic plane CO survey to identify streaming motions with an amplitude of around 5 km s⁻¹ on length scales larger than 0.22R₀, where R₀ is the radius of the solar orbit around the Galactic center, so that 0.22R₀ corresponds to a scale of between 1.87 and 2.2 kpc for the two values considered in that work. Brand & Blitz (1993) reported streaming motions with mean values of 12 km s⁻¹ in their modeling of the spectro-photometric distances for a sample H II regions and the LOS velocities for their associated MCs. The authors found a quadrupolar LOS velocity pattern, qualitatively comparable to those in Fig. 11, and report that the molecular gas is streaming past the LSR at about 3.8 km s⁻¹, which are within a factor two we obtained for CO in the Solar vicinity. Our reported velocity dispersions in H I are also comparable to those found in The H I Nearby Galaxy Survey (THINGS; Walter et al. 2008; Schmidt et al. 2016). However, the physical conditions in their sample are not in general equivalent to those in the solar vicinity.

The difference in the amplitude of the streaming motions estimated with the HOG method for H I and CO can be related to the properties of the global turbulent motion in the media dominated by each tracer. Dense, CO-dominated clouds form in regions where the velocity field is convergent, corresponding to regions where shocks dissipate turbulent kinetic energy in a supersonically turbulent flow (see, for example, Koyama & Inutsuka 2000; Audit & Hennebelle 2005; Heitsch et al. 2006). Therefore, it is natural that the velocity dispersions in those regions should be smaller. Moreover, the observed linear relation between size and velocity dispersion implies that smaller regions sampled by CO have lower dispersion than the more extended H I clouds (Larson 1981; Heyer et al. 2009).

6.3 Comparison to previous works

The most direct precedent of our global reconstruction of v_LOS as a function of distance is that presented in Tchernyshyov & Peek (2017), from here on TP17. Other previous reconstructions have focused on individual spiral arms, using ISM flow models to invert emission line observations, (see, for example, Shane 1972; Foster & MacWilliams 2006) or individual regions (see, for example, Traficante et al. 2014). Reid et al. (2016) used the Very Long Baseline Interferometry (VLBI) observations of maser sources associated with young massive stars (see Reid et al. 2009; Honma et al. 2012, and references therein) in combination with what the authors considered spiral arm signatures in CO and H I surveys to produce distance probability density functions as a function of v_LOS. None of these approaches uses the information available in extinction-based 3D dust density reconstructions.

TP17 used as input maps of reddening as a function of l, b, and d (PPP) and maps of H I and CO line emission as a function of l, b, and v (PPV) to produce the ISM distribution as a function of l, b, d, and v (PPDV). Their method is based on minimizing the difference between the observed PPV cube, from H I or CO line emission, and a model PPV cube, obtained by modeling the density distribution in PPDV as the product of a density function and a Gaussian distribution of v centered on the LOS velocity from a flat Galactic rotation curve, and projecting into PPV by summing across distance channels. The product of this minimization is a map of the LOS velocity residual, which is regularized by adding a term penalizing significant differences between the LOS velocities of contiguous pixels.

The HOG-based kinetic tomography differs from that in TP17 in three main aspects. First, the TP17 input is the radial difference of the 3D interstellar dust reddening presented in Green et al. (2015). These maps cover the distances between 63 pc and 63 kpc in steps of half a distance modulus (roughly 12.6 pc), but are limited to the Galactic longitude range 5° < l < 243°. We employed the 3D dust reconstruction from Edenhofer et al. (2024), which is limited to distances between 69 pc and 1.25 kpc but uses an improved source selection based on observations not available at the time of Green et al. (2015).

Second, the method in TP17 can be considered as an inversion of the usual kinematic distance method, assuming rotational motions around the Galactic center as a prior. The HOG does not assume a prior on the velocity distribution. This means that the two methods would converge where v_LOS is close to the expectation from rotation around the Galactic center, within the 5 km s⁻¹ LOS velocity dispersion assumed for the regularization in TP17, but diverge where it is not. Thus, the HOG-based kinetic tomography favors the identification of unusual velocity offsets, assuming that the plane-of-the-sky morphological matching represents the relation between the dust and gas distributions.

Third, TP17 is a pixel-by-pixel analysis where additional information is introduced in the model by assuming similar v_LOS for adjacent pixels. The HOG method uses the gradient orientation, a non-local quantity calculated over the angular scale set by the Gaussian derivative kernel size. Then, it accumulates the statistics over the 10° × 10° region to obtain the relation between d and v_LOS. Thus, TP17 obtains a higher-resolution reconstruction of PPDV than the HOG method at the expense of having to assume priors for the LOS velocity and spatial and velocity coherence. The higher angular resolution in TP17 favors the comparison with the VLBI observations. Still, it can be biased toward identifying the most massive object along a LOS and associating it with the most massive object for the same position in the line emission data. This would correctly assign a distance and velocity to a dense object like a giant molecular cloud (GMC), but produce a less accurate solution for less-concentrated neutral gas, which is more exposed to the acceleration by SN feedback and deviations from pure Galactic rotation (Iffrig & Hennebelle 2015; Walch et al. 2015).

Figure 20 presents the LOS velocity reconstruction obtained in TP17 for the region within 1250 pc from the Sun. From visual comparison with Fig. 11, it is evident that TP17 also recovers the quadrupolar v_LOS pattern resulting in circular motions around the Galactic center. However, the amplitude of the reconstructed v_LOS in TP17 is much higher than the values we report in Fig. 11, reaching maxima and minima of up to ±40 km s⁻¹. These large residual velocities are beyond the range of our HOG-based reconstruction, given our input range −25 < v_LOS < 25 km s⁻¹. They are recovered using the broader input v_LOS range considered in Appendix A.3.3, but at the expense of a higher chance correlation that biases the HOG results toward higher v_LOS.

There are only five VLBI maser observations within the ISM volume considered in this paper. Figure 21 shows the difference between the v_LOS reconstructions presented in Fig. 11 and TP17 and the observations for those five sources. For our reconstruction, we considered the difference between the nearest distance channel with an assigned v_LOS and the position of the maser. The v_LOS differences between the masers and the TP17 reconstruction are positive and negative and characterized by a standard deviation of roughly 6.7 km s⁻¹. The most significant discrepancy is for the source G176.51+00.20, where the difference is around −14 km s⁻¹. The v_LOS differences between the masers and our reconstructions are also in either direction. They are characterized by standard deviations around 13.5 and 11.3 km s⁻¹ for H I and CO , respectively. The most considerable difference is for source G14.33-00.64, where the lack of morphological correlation between CO and 3D dust limits the v_LOS reconstruction.

The HOG-based v_LOS reconstruction requires the averaging over large sky areas. Thus, its results are expected to show less agreement with the masers results than the TP17 reconstruction, which has an angular resolution much closer to the VLBI observations. Yet, Fig. 21 shows that the discrepancy between the maser and HOG v_LOS for CO is less than a few km s⁻¹ for most of the sources, and in the case of G090.21+02.11 it is closer than TP17. The fact that the HOG v_LOS for H I is systematically further from the maser values reflects that it samples more diffuse gas, which is not necessarily anchored to the denser gas producing the maser emission. Given the low number of maser sources to anchor both reconstructions, the comparison with the VLBI results is limited. Thus, we compare our results to those in TP17. We reserve further comparison between the HOG and the masers for Appendix E.

The differences in the HOG and TP17 v_LOS reconstruction methods provide several potential causes for the discrepancies in the results reported in Fig. 11 and Fig. 20. It is plausible that the high v_LOS values in TP17 for the range d < 1250 kpc are due to features not resolved in the Edenhofer et al. (2024) 3D extinction models or washed away by the spatial averaging required to compute the HOG. Peek et al. (2022) presented a detailed study of a collection of CO clouds in the range 135° < l < 160° using the 3D dust reconstructions from Green et al. (2019) and identified a few objects within d < 1250 kpc with v_LOS < −30 km s⁻¹; namely the clouds identified with numbers 3397, 3430, 4338, 5935, 6078, and 6142 in the Miville-Deschênes et al. (2017) catalog. The angular size of these clouds ranges between $0_{\cdot }^{\circ }23$ and $0_{\cdot }^{\circ }67$, so it is unlikely that they are the dominant feature of the average across each 10° × 10° distance channels in which the HOG is computed. In that case, the discrepancy is introduced by the different scales sampled by each method.

The HOG-based kinetic tomography employs the information contained in the plane-of-the-sky distribution of the tracers to distinguish the motions in H I and CO. It provides additional information to break the degeneracies in the inversion problem tackled with the TP17. Yet, it is limited by the resolution of the 3D dust reconstruction and the spatial averaging. Thus, an optimal approach is the combination of both methods, using the pixel-by-pixel inversion in TP17 and the HOG-based spatial correlation to regularize the solutions and solve the LOS reconstruction ambiguities. This merging of techniques is not straightforward and falls beyond the scope of this paper, but should bridge between our results and those in TP17.

One of the advantages of HOG-based kinetic tomography is the separate treatment of H I and CO as tracers of different ISM dynamics. Although the global LOS velocity reconstructions reported in Fig. 11 are qualitatively similar, the streaming motions obtained are different, as illustrated in Fig. 13. The differences between the streaming motions from H I and CO, reported in Fig. 22, are roughly around zero, indicating that, on average, the two tracers are co-moving in the studied volume. The dispersion of the LOS velocity difference distribution is around 15 km s⁻¹, suggesting that there are significant local velocity offsets in the material sampled by each tracer.

Several physical conditions can lead to a velocity offset between H I and CO. Large-scale shocks, from spiral density waves or SNe , would couple more efficiently to the diffuse material sampled by H I, accelerating it with respect to the denser material, traced by CO (see, for example, Gatto et al. 2015; Iffrig & Hennebelle 2015; Izquierdo et al. 2021).

There are some tiles where the H I shows diverging velocities that are not matched by the CO motions. A significant example of this behavior is found at d ≈ 800 pc toward the region 80 < l < 90°. In this particular example, which is located right at the axis of the Radcliffe wave, the distance-LOS velocity mapping produced with the HOG, shown in Fig. 23, indicates that the structure traced by CO closely follows Galactic rotation. In contrast, the structures traced by H I appear offset from Galactic rotation by roughly 10 km s⁻¹ in either direction. The center of these diverging motions appears to be close to the location of the North America nebula (NGC 7000; ${81}_{\cdot }^{\circ }7$ < l < ${86}_{\cdot }^{\circ }6$, $-2_{\cdot }^{\circ }2$ < b < $0_{\cdot }^{\circ }1$, d ≈ 731 to 878 pc; Zucker et al. 2020). Previous observations have reported expanding motions in the CO MCs around H II region W80 (Bally & Scoville 1980). Clustering of the young stellar objects (YSOs) in this region reveals six spatio-kinematic groups, three expanding with velocity gradient between 0.3 and 0.5 km s⁻¹ pc⁻¹. However, on a global scale, the relative motions of the groups do not appear either divergent or convergent (Kuhn et al. 2020). Determining whether or not this and other star-forming regions and their stellar content are responsible for the velocity offset in H I is further considered in the context of the energy and momentum estimates discussed in the following section.

Fig. 20 Line-of-sight velocity reconstruction from Tchernyshyov & Peek (2017) for the distance range considered in this paper.

Fig. 21 Line-of-sight velocity reconstruction from Tchernyshyov & Peek (2017) for the distance range considered in this paper.

Fig. 22 Histogram of the difference between LOS velocities across distance channels reconstructed using H I and CO emission.

Fig. 23 Same as Fig. 4, but for the tile 80 < l < 90°.

6.4 Energy and momentum in the Solar neighborhood

6.4.1 Kinetic energy density

We used our streaming motion estimates to calculate the kinetic energy densities, E_k, across tiles in the studied volume. These correspond to just one component of the velocity field. However, until new revolutionary techniques or observational methods are invented, this is the only component of the gas velocity we can measure.

The mean E_k values of 0.11 and 0.04 eV/cm³ obtained for H I and CO are comparable to the 0.22 eV/cm³ obtained by assuming standard conditions for the Solar neighborhood (Draine 2010). Assuming equipartition between the three components of the velocity vector, these estimates should be roughly within a factor of three from the total values.

It is a known fact that previous generic estimates of the kinetic energy density are close to the thermal energy density, estimated from the distribution of thermal pressure throughout the ISM, 0.49 eV/cm³ (Jenkins & Tripp 2011). They are also within a factor of a few from the cosmic ray energy density, 1.39 eV/cm³ (Webber & Yushak 1983), and the magnetic energy density derived from the H I Zeeman effect observations, 0.89 eV/cm³ (Heiles & Crutcher 2005). Finally, they are also close to the energy densities in starlight and far-infrared radiation from dust, 0.54 and 0.31 (see, chapter 12 in Draine 2010). Our analysis provides a new and independent estimate of the kinetic energy density, demonstrating that this near equipartition is also found in the local ISM.

Still, a novelty introduced by our analysis is the mapping of the distribution of this quantity, which shows maxima a factor of ten and more above the mean value for particular regions, as illustrated in Fig. 15. The distribution of E_k across heliocentric distances, presented in the top panel of Fig. 16, displays overdensities concentrated at lower distances. This is most likely produced by the sampling of different angular scales across distances. That is, nearby distance channels have higher spatial resolutions, which resolve dense structures on spatial scales that are not sampled in the most distance channels. This results in a concentration of higher densities at a lower heliocentric radius that is reflected in our E_k estimates.

The potential distance bias is marginalized in the E_k azimuthal profiles, presented in the bottom panel of Fig. 16. There it is evident that E_k overdensities of a factor of a few are found toward 40° < l < 90° in H I and CO, around l ≈ 170° and 270° in CO, and 320° < l < 340° in H I. The first l range coincides with the position of the Radcliffe wave and the Split, as illustrated in Fig. 3.

The Radcliffe wave comprises the majority of nearby starforming regions and contains about three million solar masses of gas (Alves et al. 2020). Its existence was not among the premises of our study, but it is easily distinguishable in the H I E_k distribution presented in Fig. 15. Konietzka et al. (2024) used the motion of young stellar clusters to identify oscillatory motions with an amplitude of 14 km s⁻¹ along the Radcliffe wave. The energy overdensities we identified toward this structure are likely related to the oscillatory motions. Konietzka et al. (2024) argued that “a superposition of feedback-driven structures could reproduce the observed wave-length and amplitude of the wave”. The expanding motions in H I identified in Fig. 23 potentially be remnants of such an input from stellar feedback, but further comparison calls for a study of the Radcliffe wave in it full extent and not just in the region within |b| ≤ 5°.

The l ≈ 270° energy-overdensity corresponds to the location of the Vela C MC, also shown in Fig. 3, that is representative of the Vela GMC complexes, which host intermediate-mass star formation, up to early B-type, late O-type stars (Murphy & May 1991; Netterfield et al. 2009; Massi et al. 2019; Hottier et al. 2021). The region around 320° < l < 340° covers the position of the Ara OB1 association (l ≈ ${336}_{\cdot }^{\circ }7$; b ≈ $-1_{\cdot }^{\circ }6$), a stellar group extended across 1° on the plane of the sky, where the velocity distribution of H I has been suggested to correspond to an expanding structure (Arnal et al. 1987). It is unlikely that these individual star-forming regions would produce the energy over-densities we find in H I. However, they may be guideposts to larger-scale energy input into the ISM.

Mechanical and radiative energy input from high-mass stars is one of the primary sources of kinetic energy and ISM turbulence (see, for example, Hennebelle & Falgarone 2012; Klessen & Glover 2016). The largest contribution from high-mass stars to interstellar turbulence most likely comes from SN explosions, although recent studies indicate that the mechanical energy injected by the supernovae alone is not sufficient to explain the kinetic energy of the ionized gas in nearby galaxies (Egorov et al. 2023). Numerical simulations of a single supernova remnant (SNR) expanding in an inhomogeneous density field for different ambient medium’s densities, metallicities, and structures indicate thermal energy inputs of up to 10⁵¹ erg (Martizzi et al. 2015; Iffrig & Hennebelle 2015). The efficiency of energy transfer from supernova blast waves to the interstellar gas depends on many factors, but it has been estimated to be around 10% (Norman & Ferrara 1996; Thornton et al. 1998).

The total amount of kinetic energy in the studied volume is around 1.7 × 10⁵² erg. The total amount of kinetic energy excess relative to a background level, calculated by subtracting the mean from each reconstruction element, is 1.2 × 10⁵² erg. Given the energy input and energy transfer efficiency described above, it would take roughly 120 SNe to introduce that energy into that volume. Recently, Swiggum et al. (2024) used Gaia astrometry and other large-scale spectroscopic surveys to show that around 57% of the young stellar clusters within one kiloparsec from the Sun arise from three distinct spatial volumes and these families produced more than 200 SNe in the past 30 Myr , including the sources of the Local Bubble and GSH 238+00+09. This number of SNe is roughly consistent with the 15 Myr⁻¹ SNe for gas at surface densities compared to the local ISM, around 10 M_⊙/pc² (Walch et al. 2015). Given the 10.8 km s⁻¹ velocity dispersion, turbulent mixing would have spread this input energy covering scales of roughly 400 pc in 36 Myr. Thus, SN feedback is a plausible candidate for the source of the energy overdensities reported in Fig. 15, but it does not entirely explain its distribution.

In the absence of a direct tracer to locate the SN energy input in the local ISM, our best proxy for the potential location of SNe is the distribution of the current population of high-mass stars. One of the state-of-the-art catalogs of O-, B-, and A-type (OBA) stars in the Solar vicinity is presented in Zari et al. (2021). Although this catalog is not complete, due, for example, to the saturation of Gaia for bright sources, it provides a global view of the potential stellar energy input into the ISM. Figure 24 shows the distribution of OBA stars projected into the same grid used for the streaming motion reconstructions. It is evident that the energy overdensities found toward 320° < l < 340° in Fig. 15 are in the vicinity of an OBA-type stellar overdensity, making plausible the argument of SNe or other types of stellar feedback as the source of the kinetic energy accumulation. However, this argument does not hold in the region toward 40° < l < 90°.

If the kinetic energy overdensities reported in Fig. 15 do not correspond to overdensities in the population of high-mass stars, we can consider three alternative scenarios. The first is that the stellar population producing the SN feedback is not currently located at the position of the energy overdensities. This would mean we see the imprint of a past SNe not traced by the current distribution of high-mass stars. However, it is improbable that the cluster families identified in Swiggum et al. (2024) originate from the Radcliffe Wave. Most likely, they primarily formed in compact regions 30 to 60 Myr ago, before what we call the Radcliffe Wave existed. There are signs that the cluster families were responsible for forming the Radcliffe Wave and the Split by driving superbubbles from the Alpha Per and Cr135 families, as illustrated in Figure 2 of Swiggum et al. (2024).

A second scenario is that the extinction within the Split and the Radcliffe Wave hides the stellar population from which the SN feedback originated. Redding by dust extinction complicates the identification of OBA stars, making catalogs incomplete toward high-extinction regions. Zari et al. (2021) identified in their figure 14 that the dust density in the range 40° < l < 90° is higher than the stellar density, a fact that they interpreted as the result of the decoupling of gas and stars after the passage of the spiral arm. Testing whether this is the case is problematic because it would involve measuring the plane-of-the-sky gas motions with the same accuracy as we do for stellar motions, which is beyond the reach of current observational techniques. Yet, it is unlikely that the dust in the Split and the Radcliffe Wave are hiding a prominent overdensity of high-mass stars such as that shown in Fig. 24 for distances between 900 and 1250 kpc toward the fourth Galactic quadrant.

A third scenario is one in which the energy overdensities are not exclusively produced by SNe but by large-scale driving (see, for example, Elmegreen 2024). Large-scale energy sources such as spiral arms, magneto-rotational instability (MRI), gravity-driven turbulence, tidal interactions, and inflows can drive turbulence on scales larger than one kiloparsec (see, for example, Wada 2008; Bournaud et al. 2010; Colman et al. 2022). Although the properties of turbulence in the inertial range do not depend on the injection mechanism, large-scale turbulence implies that an additional factor drives the distribution of energy beyond the domain of our reconstruction.

Tchernyshyov et al. (2018) presents a kinematic-tomography-based study of Galactic-scale numerical simulations with a fixed background potential, corresponding to the spiral density wave model (Lin & Shu 1964), and a dynamically evolving background potential, corresponding to a transient and recurrent spiral structure model (Sellwood & Carlberg 1984), without feedback. The results of the comparison with their Milky Way kinetic tomography favor the dynamic spiral structure simulation but show that both spiral arm models produce streaming motion structures that are coherent over a few kiloparsecs (see Figure 5 in Tchernyshyov et al. 2018). Whether or not these structures are compatible with the Split or the Radcliffe way is a matter of discussion. Still, they show that large-scale motions can produce spatially coherent energy overdensities without the necessity of supernova feedback. Yet, as our momentum estimates further indicate, separating a purely SN-driven scenario from a large-scale forcing scenario is challenging.

Fig. 24 Distribution of O-, B-, and A-type stars in the Zari et al. (2021) catalog limited to the range b < 5° presented in Cartesian coordinates (left) and projected into the polar grid introduced in Fig. 3 (right).

6.4.2 Momentum density

We reported our estimates of radial momentum density, p, in Fig. 17. As in the streaming motions maps, p^HI reveals multiple reversals along the line of sight, the most prominent of which is found around the location of the North America MC, along the Radcliffe wave. However, the are significant concentration of relatively high |p^HI| for d < 700 pc toward 290° < l < 360° and 700° < d < 1100 pc toward l ≈ 270°. The accumulations of relatively high |p^CO| values are found toward 40° < l < 80° for d < 900 pc and for d < 700 pc toward 290° < l < 360°.

Matter accumulation along the Radcliffe Wave may produce star formation that injects energy and momentum into the ISM. However, the initial accumulation of matter that formed the Radcliffe Wave is probably the product of the large-scale Galactic dynamics. Thus, we are again faced with the dichotomy between an SN-dominated and large-scale-forcing-dominated scenario, as we were in the case of the energy distribution over the whole studied volume. Without further evidence, we must conclude that our reconstruction of the momentum density distribution is most likely a combination of both effects.

6.4.3 Mass flow rates

Finally, we considered the mass flows related to the reconstructed streaming motions. The distributions of $\dot{M}$, presented in Fig. 18, indicate a series of converging and diverging mass flows along the line of sight with standard deviations around 3.0 × 10⁻³ and 1.2 × 10⁻³ M_⊙/yr. It is beyond the scope of this work to determine whether there is a known MC complex at the interphase of all the H I converging flows, as suggested by some of the MC formation scenarios (see, for example, Koyama & Inutsuka 2000). Reconstructing the specific flows in and around MC calls for dedicated studies of each region, as done with the HOG method for the Taurus MC in Soler et al. (2023). The map in Fig. 18 and the azimuthal profiles in Fig. 19 provide a framework to distinguish specific flows from the global fluctuations and establish a statistical comparison with the results of numerical simulation of the Galactic environment.

A rough interpretation of the $\dot{M}$ standard deviation for H I , 3.0 × 10⁻³ M_⊙/yr, can be made considering that these mass displacements correspond to stochastic fluctuations in the density and velocity fields. If we assume that the mass displacements continue stochastically with an amplitude equal to the $\dot{M}$ standard deviation, they lead to an accumulation or dissipation of a mass comparable to that of a low-mass star formation region such as the Taurus MC (M ≈ 2.4 × 10⁴ M_⊙; Goldsmith et al. 2008) in about 8 Myr. This oversimplification ignores the effect of gravity, stellar feedback, and the bias introduced by exclusively measuring radial motions with respect to the LSR. Still, it provides a rough estimate of the timescales implied by the derived $\dot{M}$.

Following the same line of argumentation, the $\dot{M}$ standard deviation for H I implies that a GMC like Orion A (M ≈ 7.57 × 10⁵ M_⊙; Lombardi et al. 2011) would assemble or disperse in roughly 250 Myr. This is considerably larger than the cloud dispersal times of a few Myr estimated from various observations and statistical considerations (see, for example, Elmegreen 2000; Kruijssen et al. 2019). If we assume that the $\dot{M}$ standard deviation for H I exclusively represents converging flows, the times are equally large compared to the estimated cloud formation times (Fukui et al. 2009; Mac Low et al. 2017, see, for example). Thus, we conclude that the typical $\dot{M}$ fluctuations we identified in our analysis do not directly correspond to GMC-forming or dispersing flows but to fluctuations in the ISM density and velocity fields. However, identifying the flows that lead to the formation and dissipation of star-forming regions like Taurus or Orion requires a focused study beyond the scope of the global averages considered in this paper.

7 Conclusions

In this paper, we studied the morphological relation between a model of three-dimensional dust extinction and the H I and CO line emission. We identified a significant correlation between dust and gas that we employed to study the ISM velocities across the region within |b| < 5° and d < 1250 pc. We produced maps of the LOS motions, streaming motions, kinetic energy and moment densities, and mass flow rate distributions across the studied volume.

We found that the material within 1250 pc from the Sun carries the global kinematic imprint of circular motions around the Galactic center. Yet, we also identified average departures from Galactic rotation, streaming motions, of around 10.8 and 6.6 km s⁻¹ for the material sampled by H I and CO , respectively.

We estimated the kinetic energy density corresponding to streaming motions. Even if we only have access to one component of the gas velocity field, we obtained mean values compatible with previous estimates, which confirms the near-equipartition with other forms of energy in the ISM. However, our result is novel because it maps that quantity across the studied volume, producing a reference point for reconstructions and deeper comparison with the results of numerical experiments.

We found energy overdensities of up to a factor of ten above the mean concentrated in the Radcliffe Wave, the Split, and regions around the Vela C and Ara star-forming regions. We provided evidence that these overdensities are most probably produced by the combined effect of large-scale forcing and SN explosions. Marginalizing the impact of one effect in favor of the other is unlikely to fully describe the physical conditions in the local 1.25-kpc radius. Yet, we did not find evidence that the local spiral arm affects the studied region’s energy or momentum density distributions.

We found that the energy overdensities correspond to regions with momentum densities and mass flow rates above the mean. It is difficult to assign these concentrations to individual feedback events or accumulations of material, as establishing those processes would require modeling the gas’s past configuration with the only velocity component available to observations. However, the synergy of this kind of analysis with the reconstruction of the stellar dynamics enabled by Gaia and other observations and the novel techniques of data modeling promises to produce a complete picture of the ISM evolution in the Solar neighborhood and beyond.

Data availability

A movie associated to Fig. 2 is available at: https://www.aanda.org

Acknowledgements

We thank the anonymous referee for carefully reading our manuscript and providing insightful comments and suggestions. JDS thanks the following people who helped with their encouragement and conversation: João Alves, Bruce Elmegreen, Torsten Eßlin, Álvaro Hacar, Sebastian Hutschenreuter, Claire Murray, Martin Piecka, Eugenio Schisano, Mark Heyer, and Oleg Egorov. JDS also thanks Daniel Seifried and Stefanie Walch-Gassner for providing the numerical simulations and synthetic observations to test the HOG method. RSK also thanks the 2024/25 Class of Radcliffe Fellows for the stimulating discussions. Some of the crucial discussions that led to this work took place under the program Milky-Way-Gaia of the PSI2 project, which is funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. This project was funded by the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130). RSK also acknowledges financial support from the German Excellence Strategy via the Heidelberg Cluster of Excellence (EXC 2181 – 390900948) “STRUCTURES”, and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 500O2206). The computations for this work were performed at the Max-Planck Institute for Astronomy (MPIA) astro-node servers. RSK is also grateful for computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Württemberg through bwHPC and the German Science Foundation (DFG) through grants INST 35/1134-1 FUGG and 35/1597-1 FUGG, and also for data storage at SDS @hd funded through grants INST 35/1314-1 FUGG and INST 35/1503-1 FUGG. Software: astropy (Astropy Collaboration 2018), SciPy (Virtanen et al. 2020), astroHOG (Soler 2020).

Appendix A Insights into the HOG method

A.1 Direction- and orientation-sensitive Rayleigh statistic

The original implementation of the HOG method presented by Soler et al. (2019) used the projected Rayleigh statistic, defined as $$V_{lm}=\frac{\sum _{ij}{w}_{ij,lm}\cos \left(2{\theta }_{ij,lm}\right)}{{\left(\sum _{ij}{w}_{ij,lm}/2\right)}^{1/2}}.$$(A.1)

In this definition, the angles between the gradient vectors, θ_ij,lm, are doubled, effectively comparing their orientation and not their direction. This means that the result of Eq. (A.1) is the same for parallel or antiparallel gradient vectors. This selection was based on the applications of this metric to the studies of magnetic field orientations, where the directions of interest are 0° and 90° (see Jow et al. 2018, and references therein). However, in our application there is additional information in the direction of the gradient vectors.

The orientation between the dust density gradients and CO emission is expected to be parallel in matching structures. The same expectation applies to the H I emission, with the exception of H ISA, where the dust density gradients and HI emission are antiparallel. Thus, the directions of interest are 0° and 180°. For that reason, we applied a direction-sensitive projected Rayleigh statistic (V_d; Eq. 2), which corresponds to Eq. (A.1) without the doubling of the angle.

Fig. A.1 compares V and V_d for a test region. The distinction between parallel and antiparallel gradients increases the significance of the distance-v_LOS pairs with morphological correlation, remarkably increasing the contrast in comparing the 3D dust density and the CO emission. It also improves the comparison between the 3D dust density and the H I emission by excluding outlying distance-v_LOS pairs that are antiparallel, as shown by their V_d < 0 values. These outliers may be produced by H ISA structures or chance correlation. The detailed study of these features requires further H ISA identification using spectral methods (see, for example, Gibson et al. 2005; Syed et al. 2023), which are beyond the scope of this work.

Using V_d also reduces the fluctuations in the distance-v_LOS mapping for the different regions in which we split the Galactic plane, as illustrated in the example shown in Fig. A.1. This regularity enabled the assignment of v_LOS to the dust density channels that are at the core of our study, which we discuss in App. C In what follows, we present studies of the effects of noise and segmentation assuming V_d as the core metric of our analysis.

A.2 Error estimation

There are two primary sources of error in the input data used in the HOG analysis. First, the uncertainties in the observations of the H and CO line emission. Second, the uncertainties in the reconstruction of the 3D dust. In what follows, we describe how these uncertainties were propagated through the HOG analysis.

A.2.1 Line emission uncertainties

We used Monte Carlo sampling to propagate the uncertainties in the line emission observations following the procedure described in appendix B of Soler et al. (2019). For each gas tracer, we produced N_MC velocity channel maps $I_{ij,l}^{\ast }$ generated from random draws of Gaussian probability distribution centered on the observed values, I_ij,l, and standard deviation equal to the noise level of the observation, σ_I. This data model assumes that the noise is constant throughout the line emission observations, which is adequate for the combination of surveys in H I4PI but less so for the CO data. However, given that we assume conservative values for the noise level in the Dame et al. (2001) observations, we are most likely overestimating the observational uncertainties in the CO in exchange for a homogeneous treatment of the data.

Fig. A.1 Comparison of the results obtained with V and Vd for the mean 3D dust and the line emission toward the tile covering the Galactic longitude range 30° < l < 40°.

Figure A.2 shows an example of the line emission and 3D dust density correlation for N_MC = 10 and 100. For comparison, we used the mean 3D dust density cube. The negligible differences between the bottom panel of Fig. A.1 and Fig. A.2 suggest that the noise in the line emission observations is not a crucial uncertainty factor in the HOG method results. This is expected if we consider that the noise in the line emission observations produces a spurious gradient orientation that is unlikely to be correlated with the structure in the 3D dust maps.

Figure A.3 presents the V_d standard deviation for the example presented in Fig. A.2. As expected from the negligible variations in the results obtained with N_MC = 0, 10, and 100 , the values of σ_Vd are less than 1% below the mean values obtained with the Monte Carlo sampling. In essence, the maps in Fig. A.3 reflect the significance of the emission detection across v_LOS. Thus, the lower σ_Vd for the H I observations are the result of the lower noise level in those observations, in comparison with the CO composite survey.

Figure A.4 presents the global results of the HOG method applied to the mean 3D dust maps and the line emission observations with and without Monte Carlo sampling. The negligible effect of the noise level in the line emission observations reported for a particular example in Fig. A.3 can be extended to the whole Galactic longitude range. This indicates that the noise in the line emission surveys does not produce any features that introduce a spurious morphological correlation in the HOG comparison with the 3D dust. The mapping technique and instrumental variations undoubtedly introduce spatial features related to the noise level, as can be seen in figure B.1 of Riener et al. (2020) for the CO observations of the Galactic plane. However, given that these noise features are unlikely to be found in the 3D dust reconstruction, they do not acutely affect the HOG analysis results.

Fig. A.2 Same as the bottom Fig. A.1, but for NMC = 10 and 100 , shown in the top and bottom panels respectively.

Fig. A.3 Standard deviation of Vd obtained for the comparison between the mean 3D dust cube and the Monte Carlo sampling of the line emission observations with for NMC = 10 and 100 , shown in the top and bottom panels respectively.

A.2.2 3D dust reconstruction uncertainties

The dust extinction modeling result presented in Edenhofer et al. (2024) is a set of 12 samples from the variational posterior of the 3D dust distribution. We computed the correlation with the line emission data for each sample to propagate the reconstruction uncertainties into the HOG analysis. We reported the mean value of V_d. This data model assumes that each sample is a plausible reconstruction of the 3D dust distribution and that the variance across the set represents the reconstruction’s uncertainties.

Fig. A.5 presents the V_d standard deviation calculated from the 12 samples of the 3D dust reconstruction toward the test region used in Fig. A.1. The values of σ_Vd are significantly larger than those found with the Monte Carlo sampling of the uncertainties in the line emission maps, thus indicating that the variance in the 3D dust reconstructions is the dominant source of error in the HOG results. To deal with this effect, we report in the main body of this paper the mean values of V_d obtained from the set of 3D dust samples. We establish the significance of the HOG results for each distance-v_LOS pair by comparing the mean V_d with the standard deviation σ_Vd. We consider different levels of contrast to σ_Vd when assigning a prevalent velocity to a 3D dust channel and estimating other physical quantities from the HOG results, as we detail in App. C.

A.2.3 Chance correlation

A significant effect to consider in the morphological matching using the HOG method is the chance correlation, which is the increase in V_d due to the accidental similarity between images. For example, most landscape photographs taken by humans on the surface of the Earth are correlated by chance due to the presence of the horizon. In the same way, the vertical symmetry of the Galactic plane imposes a general level of morphological correlation in all tracers. This effect can be minimized by employing the highest possible number of independent gradient vectors to resolve individual characteristics in the two images and produce contrast above the chance correlation level. In our application, the number of independent gradient vectors is fixed by the 3D dust reconstruction’s angular resolution. Thus, we employ additional tests to marginalize the effect of change correlation.

Figures A.6 and A.7 present the V_d values obtained for the test tile in Fig. 4 when flipping the emission PPV cubes with respect to the 3D dust PPD cube in the vertical (JK01), horizontal (JK10), and diagonal (JK11) directions. The JK label comes from Soler et al. (2019), where this operation was introduced as “jackknife” tests. The flipped cubes have the same global statistical properties, but their morphology should differ from the originals. That difference level is quantified by the V_d values obtained with the flipped sets. For example, if the line emission is distributed like a large symmetric blob, it would produce a similar V_d in the three flip tests, thus suggesting that it has insufficient unique features to produce an unambiguous morphological match.

The root-mean-square of the jackknife tests, ${\sigma }_{V_{\mathrm{d}}}^{\mathrm{J}\mathrm{K}}$, reveals a significant number of d-v_LOS pairs where the chance correlation can be acute, as illustrated in Fig. A.6 and Fig. A.7. This is expected, given the multiscale correlations in the ISM that can mimic similarity when observed or reconstructed at low angular resolution. We account for this effect in v_LOS reconstruction by evaluating the analysis results against the jackknife tests’ amplitude, as described in Sect. 3.

Fig. A.4 Same as Fig. 5, but for the mean 3D dust density and the line emission data without and with Monte Carlo sampling, shown in the top and bottom panels, respectively.

A.3 Selection of analysis parameters

The sky portion chosen for our study of the morphological correlation between the 3D dust and the H I and CO line emission is mainly driven by the coverage of the Dame et al. (2001) composite survey, |b| < 5°. However, we were flexible in choosing the spatial scale at which we made the comparison, which was set by the derivative kernel size and the segmentation of the Galactic plane. We discuss these two analysis choices in the following sections.

A.3.1 Kernel size

The Δ = ${30}_{\cdot }^{\mathrm{\prime }}0$ FWHM derivative kernel used for the HOG results in the main body of this paper was mainly driven by the angular resolution of the data. That kernel size is roughly twice the resolution of the H I4PI data and covers four telescope beams in the Dame et al. (2001) composite survey, which was only beam-width sampled in many regions. Introducing a larger derivative kernel implies a loss in V_d significance, given the lower number of independent gradient vectors for the same sky area. However, we tested Δ = ${60}_{\cdot }^{\mathrm{\prime }}0$, and ${90}_{\cdot }^{\mathrm{\prime }}0$ kernels to determine the effect of that choice in our results.

Figure A.8 shows the effect of the two kernel sizes for the same example region considered in Fig. A.1. The distance-v_LOS pairs for which we obtain the highest correlation seem unaffected by the kernel size. This can be explained by the persistence of the 3D dust density and line emission structures producing the highest V_d across these angular scales; for example, a blob producing high V_d at the ${30}_{\cdot }^{\mathrm{\prime }}0$ scale is unlikely to disappear at ${60}_{\cdot }^{\mathrm{\prime }}0$ and ${90}_{\cdot }^{\mathrm{\prime }}0$. However, the increase in the kernel size leads to a decrease in the maximum values of V_d. Thus, we set for Δ = ${30}_{\cdot }^{\mathrm{\prime }}0$ FWHM to maintain a reasonable high V_d level in the 10° × 10° tiles throughout the |b| < 5° band.

Fig. A.5 Standard deviation of Vd obtained for comparing the 12 posterior samples of the 3D dust cube and line emission observations without Monte Carlo sampling.

Fig. A.6 Standard deviation of Vd obtained for comparing the 12 posterior samples of the 3D dust cube and line emission observations without Monte Carlo sampling.

A.3.2 Galactic plane segmentation

We chose to report the result of a segmentation of the |b| < 5° band into 10° × 10° tiles starting at l₀ = 0°. This segmentation is convenient in the Galactic coordinates reference frame but arbitrary. Thus, we considered two additional segmentations starting at the reference Galactic longitudes l₀ = $-2_{\cdot }^{\circ }5$ and $-5_{\cdot }^{\circ }0$, which are also arbitrary but we use to illustrate and estimate the effect of our division of the Galactic plane in the quantities derived with the HOG method.

Figure A.9 shows the V_d values for the two additional segmentations. The values obtained for individual regions differ, which is expected given that a high-correlation feature may be split between tiles in one segmentation or another. However, there is a considerable similarity in the global correlation pattern. We used the variations among these segmentations to quantify the variance in the physical quantities estimated in App. C

We evaluated the effect of the Galactic plane segmentation in the v_LOS assignment by considering the results obtained for three l divisions. Figure A.10 presents the distribution of the v_LOS assigned to the dust tiles using the H I and CO emission. We found that the reconstructed v_LOS have a global distribution primarily unaffected by the sky segmentation. The discrepancies between the three v_LOS distributions indicate the suppression of some velocity features toward particular lines of sight, which are worth exploring in a focused study of specific regions but do not significantly change the conclusions of our work.

Fig. A.7 Same as Fig. A.6, but for CO.

Fig. A.8 Same as bottom panel of Fig. A.1, but for derivative kernel sizes Δ=${60}_{\cdot }^{\circ }0$ and ${90}_{\cdot }^{\circ }0$, shown in the top and bottom panels, respectively

Fig. A.9 Same as the top panel of Fig. 5, but Galactic longitude segmentations starting on l = $-2_{\cdot }^{\circ }5$ and $-5_{\cdot }^{\circ }0$, as shown in the top and bottom panels respectively.

A.3.3 Line-of-sight velocity range

In the main body of this paper, we chose the range −25 < v_LOS < 25 km s⁻¹ for the line emission data input. This selection was initially motivated by the maximum LOS velocity expected from the Reid et al. (2019) Galactic rotation model within the heliocentric distances d < 1.25 kpc, which is 21.6 km s⁻¹. This restricted input v_LOS range aims to mitigate the effect of chance correlation in the morphological comparison since, in principle, line emission at higher v_LOS is not necessarily related to the local material. This section tests that assumption by considering the results obtained with broader v_LOS input ranges.

Figure A.11 presents the distribution of the maximum V_d for the comparison between the 3D dust and the H I and CO line emission in the range −80 < v_LOS < 80 km s⁻¹. Visual comparison between Fig. 5 and Fig. A.11 demonstrates that the extension of the v_LOS input range does not increase the morphological correlation, quantified by max(V_d)_VLOS. This result suggests that the low-correlation regions, max(V_d)_VLOS ≈ 0, are not the consequences of the limited v_LOS-input but a result of the 3D dust distribution across distance channels.

Figure A.12 shows the results of expanded LOS velocity input ranges in the distance-v_LOS diagrams introduced in Fig. 4. For the region in that example, a significant morphological correlation exists for H I beyond the |v_LOS| < 25 km s⁻¹ limit. Except for a few outliers, the most significant morphological correlation for CO is within the range |v_LOS| < 25 km s⁻¹. This result implies that the v_LOS input limit excludes some significant departures from circular motions, effectively limiting the amplitude of the reconstructed streaming motions and the energy and momentum estimates. However, the HOG method also indicates contamination by chance correlation when employing a broader v_LOS range.

Fig. A.10 Histogram of the vLOS assigned to the 3D dust tiles in three Galactic longitude segmentations, characterized by the position of the right edge of the first tile, l0.

Figure A.13 presents the distribution of the streaming motions derived with the extended input range −120 < v_LOS < 120 km s⁻¹. The histogram shows a background signal in H I that extends beyond the range |v_LOS| < 25 km s⁻¹, particularly for negative v_LOS. This signal is unlikely to have a natural origin, as it would imply a significant velocity gradient in the solar vicinity, in tension with the observations of other tracers and the conditions for the Galactic disk stability. It is more plausible that this signal results from accidental correlation remanent after the selection described in App. A.2.3. The CO emission appears less affected by this additional chance correlation, so we use it as a reference to set the input v_LOS limits of the reconstruction. For the main body of the paper, we chose an input v_LOS range centered on zero and spanning across three times the standard deviation of the CO streaming motions estimated with the extended for the range −120 < v_LOS < 120 km s⁻¹, which is 52.1 km s⁻¹. Hence the resulting |v_LOS| < 25 km s⁻¹ range. We reserve the study of broader departures from circular motion to a subsequent, less general study.

Fig. A.11 Same as Fig. 5, but for the input LOS velocity ranges −120 < vLOS < 120 km s−1.

Fig. A.12 Same as Fig. 4, but for the input LOS velocity range −120 < vLOS < 120 km s−1.

Fig. A.13 Same as Fig. 12, but for the input LOS velocity range −120 < vLOS < 120 km s−1.

Appendix B H I self absorption and the HOG results

The plots in Fig. 5 and Fig. 6 present the maximum values of V_d identified across v_LOS channels for the volume cells in which we divided the studied region. This selection highlights structures where the density and line emission gradients are parallel, thus highlighting a positive correlation between emission intensity and density. However, there may be information in the minimum values of V_d, which in the case of significant negative values may highlight a HISA feature, where the gradients for the dust density and H I emission are antiparallel. To evaluate this possibility, we considered the values of V_d corresponding to the maximum |V_d| across distance channels, as we report in Fig. B.1. That selection accentuates distance channels dominated by antiparallel density and emission gradients.

Figure B.1 presents the V_d values corresponding to the maximum |V_d| across distance channels in the studied volume. Visual comparison with Fig 5 shows that, in general, the negative V_d are only dominant in regions with relatively low max(V_d). That result implies the positive correlation between 3D dust density and H I emission dominates throughout the studied volume. However, we performed additional tests to determine whether dust structures are associated with HISA in the distance channels with significant V_d < 0.

Figure B.2 show the morphological correlation between 3D dust distance channels and line-emission v_LOS channels, indicating the v_LOS with the maximum |V_d| at each distance after the selection criteria described in Sect. 3. As expected from the results in Fig. B.1, in most distance channels, the maximum |V_d| coincides with the maximum V_d. However, in this example, the distance channel with the minimum V_d corresponds to a v_LOS that would have been excluded in selecting positive correlation by maximum V_d.

The top panel of Fig. B.3 presents the line emission and 3D dust density for the channels with the minimum V_d identified in Fig. B.2. It is apparent from the images that there is a decrease in the H I intensity toward the region with high CO emission on the central-left part of the maps. Inspection of the H I spectrum confirms the presence of a few-km s⁻¹-wide dip that usually characterizes HISA (Gibson et al. 2000). A large fraction of the V_d signal comes from the morphological coincidence between the depression in the H I emission and the 3D structures toward the lower portion of the map. The spectra toward that portion of the map are not straightforward to interpret but show some promising indications of HISA. However, a generalization of this result is complicated.

The bottom panel of Fig. B.3 displays the line emission and 3D dust density for the distance-v_LOS pair with the lowest V_d in the studied volume. The gradient in 3D dust density coincides with the H I emission decrease. In this case, however, the depression in the H I emission is extended tens of km s⁻¹, thus it is unlikely to come from HISA. Our results can be used to indicate the potential HISA locations, but their confirmation requires additional spectral analysis beyond the scope of this work. Consequently, we report only the results for the positive correlation between H I emission and 3D dust in the main body of the paper.

Fig. B.1 Same as Fig. 5, but reporting the Vd values corresponding to the maximum |Vd|.

Fig. B.2 Same as Fig. 4, but also indicating the maximum |Vd| for H I in each distance channel (magenta markers). The dashed vertical and horizontal lines indicate the distance and vLOS with the minimum Vd toward this region.

Fig. B.3 Two examples of regions where the H I emission gradients are mainly antiparallel to the 3D dust density gradients, as identified by the minimum values of Vd.

Appendix C Effects of fixed angular resolution and distance

In this section, we consider the effect of the fixed angular resolution introduced by the derivative kernel size in the V_d distribution and its potential effect in the results of the HOG method. Our selection of a 0.5° derivative kernel implies that we are sampling scales around 0.6 pc at the front of the 3D dust reconstruction, d = 69 pc, and 10.9 pc at the back, d = 1250 pc. We evaluate the potential effect of this feature in the reported V_d values by considering its distribution across distance channels.

Figure C.1 presents the distribution of maximum V_d values reported in Fig. 5, but divided in distance bins. The trends obtained for the 50th and 99th percentile of this quantity indicate that the fixed angular resolution does not produce a systematic fall in the max(V_d) values, but that these rather increase with distance. This increasing trend is most likely the result of the inclusion of additional ISM structures in the most distant channels compensating for the lack of spatial resolution. These results are not unexpected, considering the self-similar nature of the ISM, which implies that the correlation between H I, CO, and 3D dust is extended across the scales considered in this study. Figure C.1 indicates that the highest max(V_d) values are not clustered around particular distances but are distributed along the line of sight, thus indicating that HOG is not biased by an effect related to the distance in the 3D dust model.

We also considered the potential bias introduced by the mean density, ⟨n_H⟩, in a tile in the HOG results. It is expected that max(V_d) should have higher values toward tiles with higher densities due to the larger number of significant density gradients for the same area. Yet, a morphological correlation should not be exclusively found toward high-⟨n_H⟩ tiles but distributed across ⟨n_H⟩. This is what we found in the distribution of V_d as a function of ⟨n_H⟩ reported in Fig. C.2. We found that the highest V_d are not exclusively found in high-⟨n_H⟩ tiles, thus suggesting the HOG is finding a correlation that is not only related to the amount of dust in a tile but to the coincidence between its distribution and that of the line emission tracer in the plane of the sky.

The definition in Eq. (2) implies that V_d depends on the number of gradient pairs. This is useful in our application because V_d highlights regions with a high number of parallel gradients but does not provide information about the portion of the maps that is producing the morphological correlation. To estimate the percentage of gradients producing the morphological correlation identified with V_d, we used the normalized projected Rayleigh statistic introduced in Mininni et al. (2025), $${\chi }_{kpq}=\frac{{\left(V_{\mathrm{d}}\right)}_{kpq}}{{\left(V_{\mathrm{d},max}\right)}_{kpq}}=\frac{{\left(V_{\mathrm{d}}\right)}_{kpq}}{{\left(2\sum _{ij}{w}_{ijkpq}\right)}^{1/2}}.$$(C.1)

The quantity of V_d,max represents the maximum value of Eq. (2), which corresponds to all gradient pairs being parallel, θ_ijkpq = 0°, or the comparison of identical images. In the particular case w_ijkpq = 1, (V_d,max)_kpq = (2N_kpq)^1/2, where N_kpq is the number of relative orientation angles between gradients in distance channel p and velocity channel q. The values of χ are bound to the range −1 and 1 and roughly correspond to the percentage of the gradient pairs representing the parallel or antiparallel direction, for V_d > 0 and V_d < 0, respectively.

Figure C.3 shows the maximum χ across v_LOS channels for each distance channel in our reconstruction. This quantity can be considered an analog of max(V_d) in Fig. 5, but marginalizing over the difference in the number of independent gradients across tiles. Visual comparison between Fig. 5 and Fig. C.3 indicates that the max(V_d) and max(χ) distributions are very similar in H I , as expected from the homogeneous number of independent gradients across distance and v_LOS channels. However, in CO, some portions of the b ≤ 5° range are not uniform, producing a dissimilar number of independent gradients across Galactic longitudes and v_LOS. This effect is evident in the differences between max(V_d) and max(χ) for that tracer toward the third Galactic quadrant and other low CO emission regions. In those regions, the normalization implied in Eq. C.1 may emphasize a morphological correlation represented by just a few independent gradients. Thus, we settle for V_d as a morphological metric rather than χ.

Fig. C.1 Two-dimensional histograms showing the distribution of Vd and distances for the data presented in Fig. 5. The solid white lines correspond to the 50th and 99th percentile of the max(Vd) values in each distance bin.

Fig. C.2 Two-dimensional histograms showing the distribution of Vd and ⟨nH⟩ for the data presented in Fig. 5. The solid white lines correspond to the 50 th percentile of the max(Vd) values in each density bin.

Fig. C.3 Same as Fig. 5, but for the normalized projected Rayleigh statistic, χ, Eq. (C.1).

Appendix D The HOG method in synthetic observations

We conducted a HOG study of the synthetic H I and CO observations of a simulated MC in the Simulating the life-Cycle of molecular clouds (SILCC) Zoom project (Seifried et al. 2017, 2020). It corresponds to a 125-pc-side cube extracted for the kiloparsec-sized stratified boxes in the SILCC project (Walch et al. 2015; Girichidis et al. 2016), modeled at 0.06-pc resolution using adaptive mesh refinement in combination with a chemical network to follow heating, cooling and the chemistry of H I and the formation of H₂ and CO including (self-) shielding. We used their MC1-MHD model, one of the four clouds studied through H I synthetic observations in Seifried et al. (2022). We refer to the aforementioned publications for further details on the simulations and synthetic observations.

We used the synthetic H I and ¹²CO(1-0) from the SILCC-Zoom MC1-MHD simulation to mimic observations comparable to those employed in this paper. We placed the simulated cube at a distance of 709 pc parsec to cover the range l < ±5°. We smoothed each observation using a 2D-Gaussian point spread function matching the corresponding beam and used an identical grid to that employed in the analysis of observational data, as described in Sect. 2. As a proxy for the 3D dust PPD cube, we used the 3D density cube from the simulation but smoothed employing a 2D-Gaussian kernel to match the angular resolution of the Edenhofer et al. (2024) reconstruction and projected it into a regular distance grid with a channel width of 0.24 pc . The resulting maps are presented in Fig. D.1.

We applied the HOG analysis to the MC1-MHD synthetic observations using the same parameters described in Sect. 3. The results, shown in Fig. D.2, confirm the capabilities of HOG to identify the line emission from an object in the PPD density cube. We employed this result to evaluate the density and v_LOS reconstructions presented in Sect. 5.

Figure D.3 shows a comparison of the densities across LOS distance channels obtained directly from the simulation and reconstructed with the HOG using Eq. (3). The reconstructed densities for H I are roughly within a factor of two from their values in the simulation and recover the global density profile along the LOS. The values reconstructed with the HOG are smaller than the actual values, which is an expected consequence of estimating the densities using block averaging. Reconstructions and line emission observations with higher angular resolutions enable finer segmentation and better estimates. For CO, the HOG method tends to overestimate the density, possibly by smearing the concentrated emission along the block area. However, this comparison is limited by the amount of extended, diffuse CO emission in this model, which is still a subject of active discussion (e.g., Glover & Clark 2012; Levrier et al. 2012; Godard et al. 2023; Beitia-Antero et al. 2024)

Regarding v_LOS reconstruction, the HOG method correctly reproduces the global LOS motion in the MC1-MHD simulation, which is 0 km s⁻¹. The width of the V_d signal in the v_LOS axis of Fig. D.2 corresponds to the emission linewidth coming from this isolated object (Seifried et al. 2022). Excursions from the central v_LOS are most likely related to the internal motions of the cloud and its interaction with the environment. Further identification of the internal MC dynamics using the HOG is beyond the scope of this work. However, this numerical experiment indicates that the recovered v_LOS dispersion is related to the H I and CO emission linewidths.

Fig. D.1 Projected synthetic 12CO(1-0) and H I line emission and 3D dust reconstruction from the SILCC Zoom MC1-MHD simulation placed at 709 pc from the LSR.

Fig. D.2 Same as Fig. 4, but for the synthetic observations of the SILCC-Zoom MC1-MHD simulation.

Fig. D.3 Comparison between the total, H I , and CO densities across distance slices obtained directly from the simulation cube and reconstructed with HOG using Eq. (3).

Appendix E Comparison with parallax and proper motions of star-forming regions

The best estimates of the connection between the line of sight distances and velocities throughout the Galaxy come from VLBI observations of maser sources associated with young massive stars (see, for example, Reid et al. 2009; Honma et al. 2012, and references therein). We considered the compilation of parallaxes and proper motions for about 200 masers presented by Reid et al. (2019), which is the result of observations from the National Radio Astronomy Observatory’s Very Long Baseline Array (VLBA), the Japanese VLBI Exploration of Radio Astrometry (VERA) project, the European VLBI Network, and the Australian Long Baseline Array. Only five observations in this catalog are within the Galactic volume considered in this paper. We summarize their properties in Table E.1.

Figure E.1 shows the result analysis for H I and CO emission and the 3D dust in 5° × 5° regions centered on the position of each source in Table table:masers. The region’s size is selected to achieve minimum significance levels in V_d. However, it implies that the volume sampled with the HOG can include ISM volumes not associated with the maser. We found a coincidence between the HOG analysis results and the VLBI observations for sources G090.21+02.32 and G109.87+02.11, both for CO. The HOG-derived distance and LOS velocity for G121.29+00.65 in CO are within 200 pc and a few kilometers per second from the VLBI observations. For G176.51 +00.20 there are d − v_LOS pairs with high V_d in the vicinity of the maser, but they do not correspond to the highest V_d toward that region. Finally, G014.33-00.64 does not show a significant correlation close to the d and v_LOS for the maser source.

More than revealing a profound truth about the ISM, the general lack of agreement between the HOG and VLBI results indicates that the two methods measure different ISM features. The VLBI observes high-density regions where the physical conditions favor the maser emission while the HOG matches features in the emission across larger scales. The two methods coincide in some cases, for example, in the CO emission toward G090.21+02.32 and G109.87+02.11. The cases in which the two methods disagree reveal regions where the dynamics sampled by the diffuse, extended gas emission differ from those sampled by the denser gas. Thus, the HOG reveals a velocity structure across distance channels inaccessible to the maser observations throughout the Galaxy.

Fig. E.1 Same as Fig. 4, but for a 5° × 5° area centered on the position of maser sources within |b| < 5°. They are, from top to bottom, G014.33–00.64, G090.21+02.32, G109.87+02.11, G121.29+00.65, and G176.51+00.20. The magenta disks indicate the parallax-derived distance and vLOS for the maser sources.

Appendix F Distance-LOS velocity correlation toward the Galactic center and anticenter

The assumption of circular motions around the Galactic center to derive distances from the observed emission across v_LOS is particularly catastrophic in directions where the LOS component of the velocity is small, for example, toward the Galactic center and anticenter (see, for example, Hunter et al. 2024, and references therein). Therefore, we report the HOG-derived distance-v_LOS mapping for four reference regions in Fig. F.1.

Toward the Galactic center, the analysis shows that the highest V_d are distributed in dust parcels located at d ≈ 950 and 1125 pc It is apparent that the signal from the dust clouds in that distance range is smeared across a broad v_LOS range. This is unlikely the result of the cloud dynamics but rather a manifestation of the mapping of one density structure into multiple velocities in line emission due to the large velocity gradients along the LOS (see, for example, Beaumont et al. 2013). The presence of substantial v_LOS gradients and the high densities sampled across this LOS also implies that a significant portion of the line emission may originate from distances beyond the range of the reconstruction. However, testing this hypothesis requires deeper LOS 3D dust reconstructions, which are not yet available.

Toward the Galactic anticenter, the highest V_d values appear distributed across multiple distances in a limited v_LOS range. This is an archetypal example of the effect known as “velocity crowding”, which is the superposition of multiple PPD structures into a few channels in PPV (see, for example, Ballesteros-Paredes & Mac Low 2002). For CO, the highest V_d appears around v_LOS ≈ 10 and −10 km s⁻¹ across distances between 69 and 600 pc , although the radial velocity expectation is close to zero. There are also considerable offsets from the expected v_LOS in H I. For both gas tracers, the amplitude of the departure from the expected v_LOS is comparable to the global velocity offsets reported in Fig. 22.

Fig. F.1 Same as Fig. 4, but for four 10° × 10° regions of interest centered on b = 0°, including the Galactic center (|l| < 5°) and anticenter (175° < |l| < 185°).

References

All Tables

All Figures

Fig. 1 Examples of the 3D dust, 12CO, and H I data combined in this paper. Top: Nucleon density derived from the 3D extinction model presented by Edenhofer et al. (2024) for three distance bins indicated in the figure. Middle: 12CO line emission from the Dame et al. (2001) survey for the three LOS velocity intervals indicated in the figure. Bottom: H I 21-cm line emission observations presented by HI4PI Collaboration (2016).

In the text

	Fig. 2 Rendering of the 3D dust density distribution from the extinction models presented in Edenhofer et al. (2024) for the region |b| < 5° considered in this paper. The yellow sphere represents the position of the Sun. The associated movie is available online.
In the text

Fig. 3 Heliocentric distribution of the mean nucleon density derived from the 3D dust extinction models presented in Edenhofer et al. (2024) for the sky region |b| ≤ 5°. The superimposed curves indicate large-scale features in the Solar neighborhood, as reported in Zucker et al. (2023). For the Radcliffe Wave, we indicate the full extent of the structure with the dashed white line and the segments within |b| ≤ 5° in green. The white disks indicate the locations of MC complexes referenced in the text.

In the text

Fig. 4 Example of the HOG method morphological correlation between the 3D dust density reconstruction and the line emission observations for the region 30° < l < 40° and |b| < 5°. Each panel corresponds to the correlation between distance channels, from the 3D dust density, and velocity channels, from the H I and CO line emission, shown in the left and right panels, respectively. The correlation metric is the direction-sensitive projected Rayleigh statistic (Vd, Eq. (2)). Values of Vd ≈ 0 correspond to a random orientation between the gradients of the two tracers, thus indicating low morphological correlation. Values of Vd > 2.87 correspond to mostly parallel gradients, thus indicating a significant morphological correlation. Values of Vd < 2.87 correspond to mostly antiparallel gradients. The white markers indicate the vLOS with the highest Vd for each distance channel. The green dashed lines correspond to the distance and vLOS with the highest Vd for the whole velocity and distance ranges. The three cyan lines represent the kinematic distances for the Galactic longitude at the center and the region’s edges, according to the Reid et al. (2019) Galactic rotation model.

In the text

	Fig. 5 Maximum morphological correlation between the 10° × 10° distance channels and the H I (left) and CO (right) line emission in the range −25 < vLOS < 25 km s−1, as quantified by the direction-sensitive projected Rayleigh statistic (Vd; Eq. (2)).
In the text

	Fig. 6 Maximum morphological correlation between the 10° × 10° LOS velocity channels and the dust density in the range 69 < d < 1250 pc, as quantified by the direction-sensitive projected Rayleigh statistic (Vd; Eq. (2)). The cyan curves represent the vLOS expected for d = 69 and 1250 pc according to the Reid et al. (2019) Galactic rotation model, shown by dotted and dashed lines, respectively.
In the text

Fig. 7 Example of the distance and CO velocity channel with high morphological correlation. It corresponds to the highest Vd in the comparison between the CO line emission and the 3D dust presented in the right panel of Fig. 5. From left to right, the panels show the CO and H I line emission at the indicated vLOS and the nucleon density inferred from the 3 D dust model at the corresponding distance.

In the text

	Fig. 8 Same as Fig. 7, but for high morphological correlation between H I and 3D dust. This example corresponds to the highest Vd in the comparison between the H I line emission and the 3D dust density presented in the left panel of Fig. 5.
In the text

	Fig. 9 Comparison of the maximum Vd values for H I and CO reported in Fig. 5.
In the text

	Fig. 10 Expected line-of-sight velocity (right) from the Reid et al. (2019) Galactic rotation model across the cells in our analysis.
In the text

Fig. 11 Line-of-sight velocity derived from the H I (left) and CO (right) line emission associated to the distance channels in the 3D dust reconstruction using the HOG method $(v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{\mathrm{H}\mathrm{O}\mathrm{G}})$. The colored circles correspond to the position and vLOS for the five high-mass star-forming regions with VLBI parallax observations within the studied volume.

In the text

	Fig. 12 Histograms for the differences between the LOS velocity estimated from the H I and CO line emission, $v_{\text{LOS }}^{\mathrm{H}\mathrm{O}\mathrm{G}}$, as reported in Fig. 11, and the expected LOS velocity according to the Galactic rotation model presented by Reid et al. (2019), $v_{\text{LOS }}^{\mathrm{R}19}$, presented in Fig. 10.
In the text

Fig. 13 Differences between the line-of-sight velocity derived from the H I (left) and CO (right) line emission associated with the distance channels in the 3D dust reconstruction using the HOG method, $v_{\mathrm{L}\mathrm{O}\mathrm{S}}^{\mathrm{H}\mathrm{O}\mathrm{G}}$, and the expected LOS velocities according to the Galactic rotation model presented by Reid et al. (2019), $v_{\text{LOS }}^{\mathrm{R}19}$.

In the text

	Fig. 14 Histograms of the physical quantities derived from the HOG analysis. Top left: Effective energy density, neff, that is assigned to each gas tracer, as defined in Eq. (3). Top right: Kinetic energy density, Ek, calculated with Eq. (5). Bottom left: Radial momentum, p, calculated with Eq. (7). Bottom right: Mass flow rate, $\dot{M}$, calculated with Eq. (8).
In the text

	Fig. 15 Kinetic energy density (Ek) derived from the velocity field reconstruction in Fig. 13 and n derived from the 3D dust extinction modeling, as detailed in Eq. (5).
In the text

	Fig. 16 Kinetic energy density, EK, radial and azimuthal profiles normalized to the mean values over the studied region.
In the text

	Fig. 17 Momentum density (p) obtained with the LOS velocity reconstructions in Fig. 13 and Eq. (7).
In the text

	Fig. 18 Mass flow rates ($\dot{M}$) obtained with the LOS velocity reconstructions in Fig. 13 and Eq. (8).
In the text

	Fig. 19 Azimuthal profiles of mean mass flow rates and its standard deviation.
In the text

	Fig. 20 Line-of-sight velocity reconstruction from Tchernyshyov & Peek (2017) for the distance range considered in this paper.
In the text

	Fig. 21 Line-of-sight velocity reconstruction from Tchernyshyov & Peek (2017) for the distance range considered in this paper.
In the text

	Fig. 22 Histogram of the difference between LOS velocities across distance channels reconstructed using H I and CO emission.
In the text

	Fig. 23 Same as Fig. 4, but for the tile 80 < l < 90°.
In the text

	Fig. 24 Distribution of O-, B-, and A-type stars in the Zari et al. (2021) catalog limited to the range b < 5° presented in Cartesian coordinates (left) and projected into the polar grid introduced in Fig. 3 (right).
In the text

	Fig. A.1 Comparison of the results obtained with V and Vd for the mean 3D dust and the line emission toward the tile covering the Galactic longitude range 30° < l < 40°.
In the text

	Fig. A.2 Same as the bottom Fig. A.1, but for NMC = 10 and 100 , shown in the top and bottom panels respectively.
In the text

	Fig. A.3 Standard deviation of Vd obtained for the comparison between the mean 3D dust cube and the Monte Carlo sampling of the line emission observations with for NMC = 10 and 100 , shown in the top and bottom panels respectively.
In the text

	Fig. A.4 Same as Fig. 5, but for the mean 3D dust density and the line emission data without and with Monte Carlo sampling, shown in the top and bottom panels, respectively.
In the text

	Fig. A.5 Standard deviation of Vd obtained for comparing the 12 posterior samples of the 3D dust cube and line emission observations without Monte Carlo sampling.
In the text

	Fig. A.6 Standard deviation of Vd obtained for comparing the 12 posterior samples of the 3D dust cube and line emission observations without Monte Carlo sampling.
In the text

	Fig. A.7 Same as Fig. A.6, but for CO.
In the text

	Fig. A.8 Same as bottom panel of Fig. A.1, but for derivative kernel sizes Δ=${60}_{\cdot }^{\circ }0$ and ${90}_{\cdot }^{\circ }0$, shown in the top and bottom panels, respectively
In the text

	Fig. A.9 Same as the top panel of Fig. 5, but Galactic longitude segmentations starting on l = $-2_{\cdot }^{\circ }5$ and $-5_{\cdot }^{\circ }0$, as shown in the top and bottom panels respectively.
In the text

	Fig. A.10 Histogram of the vLOS assigned to the 3D dust tiles in three Galactic longitude segmentations, characterized by the position of the right edge of the first tile, l0.
In the text

	Fig. A.11 Same as Fig. 5, but for the input LOS velocity ranges −120 < vLOS < 120 km s−1.
In the text

	Fig. A.12 Same as Fig. 4, but for the input LOS velocity range −120 < vLOS < 120 km s−1.
In the text

	Fig. A.13 Same as Fig. 12, but for the input LOS velocity range −120 < vLOS < 120 km s−1.
In the text

	Fig. B.1 Same as Fig. 5, but reporting the Vd values corresponding to the maximum |Vd|.
In the text

	Fig. B.2 Same as Fig. 4, but also indicating the maximum |Vd| for H I in each distance channel (magenta markers). The dashed vertical and horizontal lines indicate the distance and vLOS with the minimum Vd toward this region.
In the text

	Fig. B.3 Two examples of regions where the H I emission gradients are mainly antiparallel to the 3D dust density gradients, as identified by the minimum values of Vd.
In the text

	Fig. C.1 Two-dimensional histograms showing the distribution of Vd and distances for the data presented in Fig. 5. The solid white lines correspond to the 50th and 99th percentile of the max(Vd) values in each distance bin.
In the text

	Fig. C.2 Two-dimensional histograms showing the distribution of Vd and ⟨nH⟩ for the data presented in Fig. 5. The solid white lines correspond to the 50 th percentile of the max(Vd) values in each density bin.
In the text

	Fig. C.3 Same as Fig. 5, but for the normalized projected Rayleigh statistic, χ, Eq. (C.1).
In the text

	Fig. D.1 Projected synthetic 12CO(1-0) and H I line emission and 3D dust reconstruction from the SILCC Zoom MC1-MHD simulation placed at 709 pc from the LSR.
In the text

	Fig. D.2 Same as Fig. 4, but for the synthetic observations of the SILCC-Zoom MC1-MHD simulation.
In the text

	Fig. D.3 Comparison between the total, H I , and CO densities across distance slices obtained directly from the simulation cube and reconstructed with HOG using Eq. (3).
In the text

	Fig. E.1 Same as Fig. 4, but for a 5° × 5° area centered on the position of maser sources within |b| < 5°. They are, from top to bottom, G014.33–00.64, G090.21+02.32, G109.87+02.11, G121.29+00.65, and G176.51+00.20. The magenta disks indicate the parallax-derived distance and vLOS for the maser sources.
In the text

	Fig. F.1 Same as Fig. 4, but for four 10° × 10° regions of interest centered on b = 0°, including the Galactic center (|l| < 5°) and anticenter (175° < |l| < 185°).
In the text