Bayesian inference of three-dimensional gas maps: II. Galactic HI

The 21-cm emission from atomic hydrogen (HI) is one of the most important tracers of the structure and dynamics of the interstellar medium. Thanks to Galactic rotation, the line is Doppler shifted and, assuming a model for the velocity field, data from gas line surveys can be deprojected along the line of sight. However, given our vantage point in the Galaxy, such a reconstruction suffers from a number of ambiguities. Here, we argue that those can be cured by exploiting the spatial coherence of the gas density that is implied by the physical processes shaping it. We have adopted a Bayesian inference framework that allows reconstructing the three-dimensional map of HI and quantifying its uncertainty. We employ data from the HI4PI compilation to produce three-dimensional maps of Galactic HI. The reconstructed density shows structure on a variety of scales. In particular, some spurs and spiral arms can be identified with ease. We discuss the morphology of the surface mass density and the radial and vertical profiles. The reconstructed three-dimensional HI densities are available at https://doi.org/10.5281/zenodo.5956696.


Introduction
Hydrogen is the most abundant element in the universe by far. All heavier elements are formed from it, either by primordial or stellar nucleosynthesis. In the Galaxy, most of the hydrogen is locked up in stars, with its gaseous fraction only amounting to about 10 % of its stellar mass Sparke & Gallagher (2006). This interstellar gas, however, forms the reservoir from which new stars are formed. Its distributions, properties, and dynamics are therefore of paramount importance for understanding and modelling star formation. On larger scales, observations of the gas density in the Galaxy have important consequences for questions of galaxy structure and galaxy evolution (Kewley et al. 2019). In highenergy astrophysics, the gas distribution is an important input for the modelling of diffuse emission in gamma-rays which is primarily due to the production of neutral pions by interactions of non-thermal cosmic rays with the interstellar gas Tibaldo et al. (2021).
More than half of the gaseous hydrogen in the Galaxy is in the form of atomic hydrogen (Hi), which dominates the cold phases of the interstellar medium (ISM). The mass of molecular hydrogen (H 2 ) is about half of that and is located mostly in molecular clouds. Ionised hydrogen (Hii) is found in the warm and hot phases of the ISM, that is, mostly in the halo. Because of the low temperature of the phases carrying the Hi, most of it will be in the atomic ground state, such that there is no emission by transition from excited states. Instead, the hyperfine transition of the hydrogen ground state has become a standard diagnostic tool for a broad range of astrophysical disciplines. Since the electric dipole transition is forbidden, the magnetic dipole transition has a very long lifetime of ∼ 3.5 × 10 14 s, but the low collision rates in the ISM guarantee that the higher state can de-excite radiatively. In addition, the high column densities available on Galactic and extragalactic scales guarantee that observable emission is accumulated. The emission occurs at a frequency of 1.4204058 GHz or a wavelength of ∼ 21 cm, which gives the emission line its name. Because of the long lifetime, the intrinsic line width is very narrow such that the 21 cm line lends itself to velocity spectroscopy.
Ever since the discovery of the 21 cm line in space (Muller & Oort 1951;Ewen & Purcell 1951), surveys of the Galactic Hi distribution have been carried out with increasingly higher sensitivity and larger sky coverage. One concrete example of such surveys is the Leiden/Argentine/Bonn Survey (LAB Survey, Kalberla et al. 2005), which over the last 20 years has been the de facto source of information on Hi in the Galaxy. The LAB Survey has revealed the flared, warped, and multi-arm spiral structure of the Hi disk (Levine et al. 2006a,b). Recently, a new all-sky survey of the 21 cm line emission that combines different data sets for both low and high Galactic latitude data was published by the HI4PI collaboration (HI4PI Collaboration et al. 2016). This survey provides the brightness temperature obtained from the 21 cm line emission in a three-dimensional grid of Galactic longitude ℓ, Galactic latitude b, and the radial velocity v LSR with respect to the local standard of rest. For a given gas flow model, that is, a model of v LSR = v LSR (r) for the Milky Way, the threedimensional gas density can in principle be reconstructed from the gas brightness temperature (Kulkarni et al. 1982; Article number, page 1 of 17 arXiv:2202.02341v2 [astro-ph.GA] 25 May 2023 A&A proofs: manuscript no. higift Burton & te Lintel Hekkert 1986;Nakanishi & Sofue 2003;Kalberla & Dedes 2008;Marasco et al. 2017).
There are a few complications that degrade the quality of the deprojection from an ℓbv-data cube of brightness temperature to the xyz-cube of gas density. First, there are two distance solutions inside the solar circle. It is not clear how we can determine the fractions of the emission that are at near and far distances. This effect is known as the kinematic distance ambiguity. In addition, the line of sight close to the Galactic centre (longitude ℓ ≃ 0 • ) and anti-centre (longitude ℓ ≃ ±180 • ) directions exhibit little to no radial velocity and therefore lack kinematical resolution. This means that all emission piles up around v LSR = 0 km s −1 and cannot be deprojected along the line of sight. Finally, peculiar velocities, that is, random motions of gas in addition to the large-scale gas flow, for instance due to stellar winds, supernova explosions, or spiral structures, perturb the smooth mapping of distance to radial velocity. These perturbations become visible as artefacts in the deprojected gas maps. In this case, the distribution of gas is smeared out along the line of sight, leading to the well-known finger-of-god effect.
Historically, the pioneering work of Westerhout (1957) and Oort et al. (1958) provided parts of the face-on Galactic Hi density map using the 21 cm line data. Ever since, a number of studies have tried to improve on these early studies, oftentimes circumventing some of these abovementioned difficulties by limiting the scope of the deprojection. Using data from the LAB survey, Levine et al. (2006a,b) attempted to uncover both the spiral structure and the radial and vertical distributions. For the velocity model, they assumed circular rotation with a constant velocity of 220 km s −1 . They limited themselves to regions outside the solar circle where there is no kinematic distance ambiguity and lines of sight far enough away from ℓ ≃ 0 • and ℓ ≃ ±180 • . However, the effect of peculiar velocities is still very much visible in terms of finger-of-god effects. This limited the authors in their attempt to identify large-scale structures; in order to identify spiral structures, some more intricate analysis was necessary. One way to limit the influence of low-intensity features that are rather spread out in velocity is to focus only on the densest regions. Such a study has recently been performed by Koo et al. (2017), also employing LAB data. They also constrained themselves to regions outside the solar circle and to sight-lines with finite radial velocities in a circular rotation model. For each sight-line, local maxima of the velocity spectrum above a certain threshold were identified and resolved into individual cloudlets. The resulting density of these emission regions was much better resolved than the overall density, but this comes at the price of only reconstructing part of the Hi density. The authors estimate that their surface density map contains ≲ 10% of the total Hi gas in the Galaxy.
Attempts to reconstruct most of the Galactic Hi density have been presented by Nakanishi & Sofue (2003) and Nakanishi & Sofue (2016). Again based on LAB data, the authors attempted to break the kinematic distance ambiguity inside the solar circle by employing a variant of the double-Gaussian method (e.g. Clemens et al. 1988). If we assume a parametrised vertical profile and the same height of the gas distribution at the near and far distance, the latitudinal profile at a certain distance should be the sum of two contributions, one narrower and one wider in latitude. When this is fit with the parametrised profile, the emission observed at a certain v LSR can be broken down into contributions at the near and far distance. Emission close to the Galactic centre or anticentre directions still needed to be excluded. While there are some differences between the earlier and later analysis, both reconstructions show a deficit in gas density in the solar circle. In addition, the density is markedly smeared out along the sight lines.
Two other studies are noteworthy here, even though they do not attempt a direct reconstruction of the threedimensional gas density from line survey data. Focussing on regions inside the solar circle, Marasco et al. (2017) assumed spherical symmetry for the gas distribution, at least separately in the half-planes corresponding to positive and negative longitudes. Identifying the tangent points in velocity diagrams, that is, the maximum velocities with significant emission for a given line of sight, they could identify the emissivity for all radii r ≤ R ⊙ where R ⊙ is the distance of the Sun from the Galactic centre. While only giving some (azimuthal) average, the radial profiles are still useful for studying some radial dependences. A different approach altogether has been presented by Jóhannesson et al. (2018), who presented a parametric model of Hi (and H 2 ) fitted to gas line data. While a large number of model parameters can be determined in this way, the model is not flexible enough to reproduce all of the observed emission. The gas density therefore should be considered a lower bound on the actual gas density in the Galaxy.
All these previous attempts ignore an important physical constraint that the reconstructed gas density needs to satisfy and that can be exploited in reconstructing the gas density. All processes determining the gas distribution lead to spatial correlations, be it on large scales (spiral arms from density waves; e.g. Shu 2016, disk profile from gravitational collapse) or small scales (turbulence of the ISM; e.g. Kolmogorov 1941). Some theoretical arguments and numerical simulations suggest that the statistics of the threedimensional gas density is log-normal (Nordlund & Padoan 1999;Ostriker et al. 2001). Such a log-normal field can be characterised by the two-point correlation in coordinate space. Under the additional assumption of (spatial) stationarity and isotropy, the two-point function in harmonic space becomes a function of the wave number only, the power spectrum. Reconstructing the gas density under the constraint of a known or to be determined correlation structure is a Bayesian inference problem. The Bayesian nature allows us to reconstruct not only the gas density, but also to obtain an estimate of its uncertainty. This method has been successfully applied to a deprojection of the threedimensional distribution of H 2 in the Galaxy (Mertsch & Vittino 2021).
For the problem at hand, a number of the abovementioned short comings can be overcome when the correlations are taken into account. For instance, the kinematic distance ambiguity can be resolved by demanding that the reconstructed gas density globally satisfy the correlation structure. These correlations can also help where no distance information is available, that is in the directions of ℓ = 0 • and ℓ = ±180 • . Given that the gas density model is probabilistic, the uncertainty will naturally be larger in spatial regions for which the data are less constraining. Finally, the data model employed takes the presence of noise into account, thus providing some denoising of the observational data.
In the following, we discuss our method and the results of our application to HI4PI data. In Sect. 2, we briefly re-view the data used, introduce the gas flow models employed, and explain the Bayesian reconstruction framework. Our results are presented in Sect. 3 and are compared with the results of previous analyses. We provide a short summary and outlook in Sect. 4.

Survey data
The reconstruction of the Hi density was carried out using HI4PI data (HI4PI Collaboration et al. 2016) that were merged from data of the 21 cm line spectra of two individual surveys, that is, the Effelsberg-Bonn Hi Survey (EBHIS; Kerp et al. 2011;Winkel et al. 2016) and the Galactic All-Sky Survey (GASS; McClure-Griffiths et al. 2010;Kalberla et al. 2010;Kalberla & Haud 2015). Both surveys have rather similar angular and spectral resolutions (see the survey parameters presented in Tbl. 1 of HI4PI Collaboration et al. 2016). Since EHBIS covers the declination range of δ ≤ 1 • and GASS has data for δ ≥ −5 • , the combined HI4PI survey provides the full-sky data. We downloaded and assembled the prepared data in the form of HEALPix maps with N side = 1024 (Górski et al. 2005) from the Strasbourg astronomical Data Center 1 .
For the spatial resolution we aim for, a coarser representation of the survey data is sufficient. We, therefore, degraded the HEALPix maps from N side = 1024 to N side = 128 for the reconstruction. More importantly, even though the survey data consists of 933 maps, that correspond to the velocity range from −600 km s −1 to 600 km s −1 , we focused only on the velocity range −320 ≲ v LSR ≲ 320 km s −1 as in the reconstruction of H 2 presented in Mertsch & Vittino (2021). We also masked out the combination of the latitude and longitude ranges −60 • ≤ b < −20 • and 270 • ≤ ℓ < 330 • in order to remove contamination from the Magellanic Clouds.
We show as an example the degraded data cube of 21 cm line emission integrated over latitude, also known as the ℓ-v-diagram, in Fig. 1. The all-sky map of column density is shown in Fig. 2.

Gas flow model
The distribution of Hi brightness temperature in ℓ, b, and v LSR depends on both the three-dimensional density of the gas and its velocity along the line of sight. This means that the reconstruction of the gas density requires knowledge about the gas motion as briefly mentioned above. Even though it is commonly assumed that the gas follows purely circular motions around the Galactic centre, this assumption is problematic for the reconstruction in that it does not provide any kinematic resolution towards the Galactic centre and anti-centre directions. In the inner Galaxy, we know, however, that the gas also has radial motions due to the presence of the Galactic bar (Blitz & Spergel 1991). The modelling of non-circular motions in this region is quite uncertain and we follow Mertsch & Vittino (2021) in considering two different gas flow models in order to somewhat bracket the systematic uncertainty. Consequently, a brief discussion for the two gas flow models is in order.
1 http://cdsarc.u-strasbg.fr/ftp/J/A+A/594/A116/HPX/ The first model is based on the smoothed particle hydrodynamics simulation by Bissantz et al. (2003) (hereafter BEG03). The BEG03 velocity profile is shown in the top panel of Fig. 3. We also note that this gas flow model was extended beyond 8 kpc using a flat rotation curve similar to the approach by Pohl et al. (2008) (see also Mertsch & Vittino 2021, for more details).
We also considered a modified version of the semianalytic model for gas-carrying orbits in the potential dominated by the Galactic bar from Sormani et al. (2015) (referred to as SBM15; see Appendix A of Mertsch & Vittino 2021). The SBM15 velocity map seen by an observer moving with the velocity of the local standard of rest at a distance 8.15 kpc from the Galactic centre Reid et al. (2019) and at an angle of 20 • with respect to the major axis of the bar is shown in the middle panel of Fig. 3.
The difference between the BEG03 and SBM15 models is shown in the bottom panel of Fig. 3. We note that the two models do not agree very well within 2 kpc from the Galactic centre. Furthermore, there is no perturbation due to spiral structures in the SBM15 model resulting in the difference from around 10 km/s to 30 km/s in certain regions in comparison to the BEG03 model. It should be noted that purely circular motion is assumed in both models beyond about 4 kpc from the Galactic centre and the differences as large as 30 km/s outside the solar circle are due to the different adopted rotation curves. Differences between a velocity model and the true gas flow will manifest as gas placed at an incorrect distance along the line of sight, which is a concern common to any deprojection from the ℓbv cubes alone. While either of our models likely exhibits deviations from the true gas flow, the differences between the gas densities that are reconstructed with either model can illustrate the dependence of features along the line of sight on the gas flow model to a certain degree.

Bayesian inference
The framework adopted here is similar to the one presented in Mertsch & Vittino (2021) for the reconstruction of the H 2 density from the Galactic CO line survey. We treated the reconstruction of the Hi density map as a high-dimensional Bayesian inference problem and sought the posterior distribution of the gas density for the given survey data using a method called Metric Gaussian Variational Inference (MGVI) (Knollmüller & Enßlin 2019). The MGVI is a method from information field theory (IFT; Enßlin 2019) that relies on approximating the posterior distribution by a parametrised Gaussian distribution with the covariance taken to be the inverse Fisher information metric. The parameters of the approximate distribution are determined by repeatedly minimising the "distance" measured by the Kullback-Leibler divergence (Kullback 1968) between this distribution and the exact posterior. In practice, this is performed in an iterative fashion by alternating between computing the covariance for the current mean and updating the mean obtained from the minimisation of the KL divergence for the current covariance. The result is no explicit representation of the posterior distribution, but an ensemble of samples drawn from it that can be used to study the Hi density and its uncertainty.    All-sky column density of Hi obtained from HI4PI data using Eq. 8 for the velocity range −320 ≲ vLSR ≲ 320 km s −1 . The map is presented in Galactic coordinates using a Mollweide projection. The graticule lines are spaced by 10 • .

Data model
The radio data of the 21 cm emission line provide the brightness temperature T (ℓ, b, v), and we would like to reconstruct the volume emissivity ε(x, y, z). We assumed that we are always in the optically thin regime, rendering the re-lation linear, where we also defined the response matrix R = R xyz lbv . The discretised version of the response matrix to transform from emissivity ε αβγ ≡ ε(x α , y β , z γ ) to brightness temperature Peculiar motions can lead to deviations between the velocities described by the large-scale velocity field and the actual velocities. While this cannot be predicted, we can allow for some deviations by smearing the response matrix with a Gaussian kernel with variance σ 2 v , The velocity dispersion σ v has been estimated to be 9.2 km/s for the first quadrant (0 • < l < 90 • ) and 9.0km/s for the fourth quadrant (270 • < l < 360 • ) (Malhotra 1995). We adopted σ v = 10 km/s in line with Nakanishi & Sofue (2003). It should also be noted that the measured brightness temperature would also be altered by noise term which here we denote by n. We modelled n as an additive Gaussian white noise p(n) = G(n, N ) with the diagonal covariance in harmonic spaceÑ = 2πσ 2 n δ(k−k ′ ) as in Mertsch & Vittino (2021). We approximated the noise level of the HI4PI data (see Tab. 1 in HI4PI Collaboration et al. 2016) with σ n = 50 mK. The final data model for the brightness temperature is then The likelihood function can then be written as (see also Mertsch & Vittino 2021 for more details):

Prior
In order to infer the posterior distribution, we have to specify also the prior. We modelled the Hi emissitivity as a lognormal distributed random field meaning s(r) = ln [ε(r)/ε 0 ] follows a normal distribution. We further assumed s(r) to be statistically homogeneous and isotropic, such that the two-point correlation function in the harmonic domain could be characterised simply as ⟨s (k) is the power spectrum. To allow for enough flexibility, we modelled P (k) as the product of a (broken) power law and a random field (in log k), Here, ϕ y and ϕ m represent the slope and the variance which are random variables following the Gaussian distributions G(µ y , σ y ) and G(µ m , σ m ) respectively, F −1 denotes the inverse Fourier transform from the variable t to its conjugate log(k), and τ (t) is the above-mentioned random field (in log(k)). The metaparameters were fixed as follows: µ y = −13, σ y = 1, µ m = −4, σ m = 0.2, a = 1, and t 0 = 0.1. We reconstructed the random field τ (t) at the same time as the signal itself and, thus, determined P (k) together with the gas density. Where the data were constraining enough, the data were able to deviate from the prior structure; where they were not the gas density was reconstructed with a realisation of the log-normal field prescribed by the power spectrum.
We note that the assumption of a homogeneous random field cannot be fully realised as the Hi is localised in a disk of finite extent in radius. Ignoring this can lead to biases, most severely affected by the quick fall-off of gas density with distance z from the disk. We therefore scaled the signal field ε xyz with an exponential profile in z, exp (−|z|/z h ), where the scale-height z h varies with galactocentric radius r as The scale of radial variation is the result of an earlier analysis (Kalberla & Dedes 2008), and we fixed z 0 h = 150 pc. While this z-profile does not preclude the possibility that other profiles are found by the reconstruction, we find that z 0 h = 150 pc in particular leads to consistent results.
The inference method was implemented using the NIFTy package (Selig et al. 2013;Steininger et al. 2019) 2 , specifically NIFTy5 (Arras et al. 2019a), which was previously applied to a variety of problems ranging from the reconstruction of the three-dimensional dust density in the Galaxy from reddening data (Leike & Enßlin 2019) to radio interferometry (Arras et al. 2019b).

Reconstructed density
The conversion from the velocity-integrated brightness temperature T B to the Hi column density is normally defined (Draine 2011 where N Hi is the total column density for the line of sight (ℓ, b) and X Hi ≃ 1.823 × 10 18 cm −2 (K km s −1 ) −1 . We note that Eq. 8 is applicable only in the optically thin limit which means that the derived column density strictly only provides a lower limit for the true column density in the disk (Nakanishi & Sofue 2003;Jóhannesson et al. 2018). The density can then be obtained from the reconstructed volume emissivity as n Hi (r) = X Hi ε(r).

Results and discussion
The main results of our analysis are the three-dimensional density maps of Galactic Hi for the BEG03 and SBM15 gas flow models. These maps are presented on a Cartesian grid with (n x , n y , n z ) = (768, 768, 64) stretching from −24 kpc to 24 kpc in the x-and y-directions and from −2 kpc to 2 kpc in the z-direction. Therefore, the spatial resolution for these reconstructions is 62.5 pc, that is the same as for the H 2 reconstruction of Mertsch & Vittino (2021) and the highest resolution among all three-dimensional reconstructions of Hi line surveys to date. We have made our maps accessible on zenodo 3 . We present in Fig. 4 the mean and standard deviation of the Hi surface mass density for both the BEG03 (top panels) and SBM15 (bottom panels) gas flow models. It is clear that the survey data are successfully reconstructed into localised clusters of gas following the spiral structures (fitted independently with data of molecular masers associated with very young high-mass stars) from Reid et al. (2019). Generally, the maps agree on the presence of largescale features such as spiral arms, but their exact position and other properties differ between the gas flow models. This is expected given the disagreement between the radial velocities (see Fig. 3) which upon deprojection results in different distances to emission features. Some of these differences are also present in the reconstructions of H 2 , and they have been linked to local extrema in the gas flow models (Mertsch & Vittino 2021).
For example, for the BEG03 model there seems to be a spur-like feature stretching from (x, y) ≃ (−8, 2) kpc to (x, y) = (−4, 6) kpc, and it is rather tempting to identify this feature with the Scutum-Centaurus (Sct-Cen) spiral arm. However, this is probably local emission with velocities close to v LSR ≈ 0 km/s that is misreconstructed at large distances from the Solar System: Brightness temperature spectra peak at v ≳ 5 km s −1 for lines of sight with −70 • ≲ l ≲ −20 • . The gas flow models, however, give negative v LSR inside the solar circle along these directions. Part of the emission is therefore placed outside the solar circle, while the large latitudinal spread indicates that the emission is local. Another noticeable difference that can be spotted is the small arc stretching from (x, y) = (−5, 0) kpc to (x, y) = (−1, −4) kpc along the Sagittarius-Carina (Sgr-Car) arm in the BEG03 reconstructed map. This arc does not appear in the SBM15 reconstructed map which again might be linked to the arc in the same position in the map showing the difference in velocity profiles (see the bottom panel of Fig.  3).
Another version of the reconstructed mean density overlaid with velocity contours from the gas flow models and a longitude grid is shown in Fig. 5. Since the dynamical range of the gas density is rather large, as expected for a log-normal field, we have presented this map on a logarithmic colour scale to better highlight features that are not visible on a linear scale. For the reconstruction based on the BEG03 gas flow model, it is apparent that some of the bright features are aligned with regions of strong gradients in the velocity field, for example the arc stretching from (x, y) = (−5, 0) kpc to (x, y) = (−1, −4) kpc that was already mentioned above.
So far, we have discussed localised structures in the two models based on the reconstructed mean. However, some of the localised structures with relatively high surface mass density might also have rather large uncertainties. Unlike previous methods, the Bayesian inference method readily provides an estimate of the uncertainty. It might be more informative to judge the validity of specific features by comparing the mean µ with its uncertainty σ. We therefore followed Mertsch & Vittino (2021) and defined the signal- to-noise ratio (S/N) as µ/σ. The surface mass density S/N overlaid with the spiral arms from Reid et al. (2019) is presented in Fig. 6 for the BEG03 (top panel) and SBM15 (bottom panel). We shall now proceed with a more detailed discussion on some notable features of the reconstructed maps of the BEG03 and SBM15 models: -Most of the structures with a S/N roughly above 3 could be associated with spiral arms, especially for the Sct-Cen, Sgr-Car, Local and Perseus arms. We note, however, that little Hi has been reconstructed along the Norma arm in the inner Galaxy for either model (see Fig. 4). -Interestingly, even though significant structures with relatively high S/N are present around the Local arm in both maps in Fig. 6, the Local arm seems to be better reconstructed in the BEG03 model (see the top left panel of Fig. 4). The corresponding structure in the SBM15 reconstructed map seems to be more scattered around the Local arm. However, we caution again that this might be an artefact of the spiral structure already present in the velocity model for the case of the BEG03 model. -The BEG03 reconstructed map seems to extend to larger galactocentric radii than the SBM15 one. Some clusters of Hi gas outside the solar circle are reconstructed at different distances in the two models due to the different rotation curves adopted beyond 4 kpc from the Galactic centre. Because of the rather small velocity gradient, this easily translates into differences of the order of a kiloparsec and thus affects the association with spiral arms. We could see, for example, the arc of Hi gas in the BEG03 reconstructed map stretching between (x, y) ≃ (10, 15) kpc and (x, y) = (12, 2) kpc, which is a little off the Norma arm (see the top left   Fig. 4). In the SBM15 reconstructed map, this arc instead ends up on the Norma arm.
We note that the comparison between the surface mass densities reconstructed with both gas models might not reflect the full systematic uncertainty due to our ignorance of the true velocity field. One possible way of estimating deviations between the gas flow models and the true velocity field makes use of the velocities measured for molecular masers in certain high-mass star-forming regions (Reid et al. 2019). While deviations from a circular rotation model are clustered around zero, the spread can be several tens of km s −1 (Tchernyshyov & Peek 2017;Tchernyshyov et al. 2018) (see also Fig. 6 of Reid et al. 2019).
It has also been suggested that when gas spectra were combined with differential dust extinction data, the Galactic velocity field can be reconstructed (Tchernyshyov & Peek 2017). When this is applied, localised devitations of the reconstructed velocities from a circular rotation curve as large as 30 km s −1 have been inferred. As stressed in Sect. 2.2, these devitations can lead to gas placed at incorrect distances along the line of sight, in particular, along directions of small velocity gradients. While these studies are very promising, they require a rather tight correlation between gas and dust, that is, essentially a constant gas-todust ratio. We have turned this problem around and compared our reconstructed gas maps with gas maps converted from dust maps assuming this constant gas-to-dust ratio. In Appendix A, we present the results and show that some features in our gas maps can be clearly identified with features in extinction maps while others cannot.
In Fig. 7, we again show the surface mass density in the BEG03 and SBM15 models together with previous reconstructions from Nakanishi & Sofue (2003)  All surface density maps are presented with the same dynamical range, except for the one by Koo et al. (2017). This is because Koo et al. (2017) focused essentially on regions of dense Hi concentrations and, thus, provided the surface mass density map only for the densest structures. A few comments on the similarities and differences among our reconstructions and previous ones are in order: -Both the BEG03 and SBM15 reconstructed maps have spiral structures beyond the solar circle that are similar to all earlier analyses. For example, segments of the Sgr-Car arm in the region with y > 0 or segments of the Sct-Cen arm appear in all of these reconstructed maps. We note, however, that not many localised structure with high surface mass density being reconstructed inside the solar circle for the maps of Nakanishi & Sofue (2003) and Nakanishi & Sofue (2016). In our case, we clearly see more high surface mass-density clusters of gas inside the solar circle. -As discussed above, peculiar motion of the gas can result in the finger-of-god effect, which manifests itself as elongated structures along the lines of sight in most of the previous reconstructed maps (see e.g. Nakanishi & Sofue 2003;Levine et al. 2006b;Nakanishi & Sofue 2016).
In the Bayesian framework, these structures seem less apparent (see the top panels of Fig. 7) thanks to both the consideration of correlations and the treatment of velocity dispersion in the response matrix. -Finally, our reconstructed maps show clusters of Hi gas in regions with galactocentric radii between 15 and 20 kpc along the Galactic centre and anti-centre directions that were excluded in all previous studies. However, we caution that most of them have much lower S/N than other features, and their distances are, given the lack of association with known spiral structure, most likely unreliable.
We now examine the Hi profile in the z-direction. We present in Fig. 8 the weighted mean of the Hi disk vertical offset defined as (Jóhannesson et al. 2018) z 0 (x, y) = dz z n Hi (x, y, z) dz n Hi (x, y, z) , where we approximated the integrals with sums over the bins in z-direction. It is clear that the midplane of the Hi disk is displaced from z = 0 for galactocentric radii r ≳ 12 kpc for both the BEG03 (upper panel) and SBM15 (lower panel) models. This feature is sometimes referred to as the Galactic warp and has also been found in many of the previous studies (see e.g. Burton & te Lintel Hekkert 1986;Diplas & Savage 1991;Nakanishi & Sofue 2003;Levine et al. 2006a;Nakanishi & Sofue 2016;Jóhannesson et al. 2018). In our case, the Hi disk beyond 12 kpc from the Galactic centre bends towards negative z in the region between φ ≃ 30 • to φ ≃ 150 • and positive z in the region between φ ≃ 210 • to φ ≃ 330 • for both velocity models. We note, however, that the exact values of z 0 are slightly different between the two models in certain parts of the Galaxy. For example, the outer disk (r > 12 kpc) seems to be more warped in the BEG03 map than in the SBM15 map in the regions between φ = 60 • and φ = 90 • and between φ = 300 • and φ = 330 • . Some summary statistics of the Hi gas distribution in galactocentric radius can also be extracted from our reconstructed maps. We present in the upper panel of Fig. 9 the surface mass density averaged in galactocentric rings r i ≤ r < r i+1 where r i = i∆r with ∆r = 0.5 kpc. In the SBM15 model, the reconstructed Hi surface mass density increases slightly from Σ Hi ≃ 1 M ⊙ pc −2 around the Galactic centre to the maximum value of Σ Hi ≃ 6 M ⊙ pc −2 at r ≃ 12 kpc. Similarly, the Hi surface mass density from the BEG03 model also has Σ Hi ≃ 1 M ⊙ pc −2 around the Galactic centre and the maximum value is Σ Hi ≃ 7 M ⊙ pc −2 . The BEG03 profile, however, seems to have a few more local maxima and the Hi surface mass density is slightly higher than that of SBM15 for rings between r = 4 kpc and r = 8 kpc. We note that the surface mass density with these galactocentric radii appears to be compatible with data from Wouterloot et al. (1990) (filled circles in Fig. 9) and also comparable to the commonly quoted value of 5 M ⊙ pc −2 for the inner Galaxy (Dickey & Lockman 1990;Kalberla & Kerp 2009). Beyond r ≃ 12 kpc, both models exhibit a continuously decreasing Hi surface mass density up to r ≃ 20 kpc, with values slightly lower than data from Wouterloot et al. (1990)  (2008). Finally, it is worth mentioning that, in the region r ≳ 20 kpc, the Hi surface mass density again increases slightly, but we caution that this is likely due to the mis- In the bottom panel of Fig. 9, we show the spread of the reconstructed Hi disk around the midplane as a function of galactocentric radius r. Specifically, we fitted a Gaussian profile, as is oftentimes assumed, and display its full width at half maximum (FWHM), which is related to the standard deviation σ as FWHM = 2 √ 2 ln 2σ ≃ 2.4σ. The warping means that the maximum position of the vertical profile in z varies with x and y, therefore, we fit the Gaussian profile for limited ranges of galactocentric radius r and azimuth φ and then averaged over the bins in φ for each r bin. The width increases steadily from the inner to the outer Galaxy. Unlike earlier work (Wouterloot et al. 1990;Kalberla & Dedes 2008), however, we do not find an exponential behaviour and instead the FWHM almost saturates beyond r ∼ 15 kpc. We stress that while some assumption on the gas being confined to the plane was necessary in order to model the gas density as a homogeneous random field, see Sect. 2.3.2, the reconstruction algorithm has the freedom to deviate from the specific profile assumed if this is preferred by the data. We note that the reconstructed widths beyond r ∼ 15 kpc differ significantly between the different samples from the posterior distribution, which is reflected in the relatively large error bars.

Summary and perspectives
We have presented new three-dimensional Hi maps with a spatial resolution of 62.5 pc deprojected from the HI4PI 21 cm line survey (HI4PI Collaboration et al. 2016). These maps were reconstructed using two gas flow models, namely the BEG03 and SBM15 models. Both of them consider the effects of the Galactic bar. The former was derived from a simulation of the gas flow in a predetermined potential (Bissantz et al. 2003) and the latter is based on a model for gas-carrying orbits in the bar potential (Sormani et al. 2015).
We assumed the optically thin limit such that the velocity-integrated brightness temperature and the Hi gas column density could be related by a simple relation. The deprojection was then carried out using MGVI (Knollmüller & Enßlin 2019). This Bayesian framework has previously been adopted for the CO line survey in Mertsch & Vittino (2021) and allows evaluating not only the mean gas density but also its uncertainty. We can therefore distinguish between statistically significant structures and noise artefacts. While we have assumed the existence of correlations in configuration space, the specific form of the power spectrum was not imposed but instead determined together with the gas density during the reconstruction.
We have made our mean Hi maps and their uncertainty available 3 . We examined the surface mass density of deprojected maps and found that there are structures that are compatible with spiral arms fitted from parallax measurements of masers in Reid et al. (2019). Some radial profiles out of the three-dimensional gas distribution were also investigated and apparently agree well with earlier analyses for regions within roughly 12 kpc from the Galactic centre. The warping of the Hi disk, which has been identified in some of the previous studies, is observed in these reconstructions as well.
In the future, we plan to improve this analysis in two main ways: First, we will relax the optically thin assumption, which will provide more reliable reconstructions of the Hi density, in particular, in the inner Galaxy. Second, we will use more flexible velocity models than those used here. In particular, we will constrain these models with the velocity information from molecular masers Reid et al. (2019). If this is successful, even a joint reconstruction of gas density and velocity field could be attempted, as suggested by Mertsch & Vittino (2021).
In addition, we note that dust reddening data might provide additional constraints for the three-dimensional distribution of gas. In fact, the IFT framework has been successfully applied to derive a three-dimensional view of nearby dust (Leike & Enßlin 2019;Leike et al. 2022), and results agree rather well with other analyses within ∼ 300 pc around the Solar System (Lallement et al. 2018;Green et al. 2019). If the correlation between dust and gas was known with sufficient accuracy, dust reddening data could be used as an additional input for the gas reconstruction. Moreover, by combining dust data with emission lines of HI and CO, a three-dimensional model for the line-of-sight velocity could be derived. Tchernyshyov & Peek (2017), however, pointed out that a combined analysis requires matching current reddening data, which are local (within a few kiloparsecs from the Solar System), with emission line data, which presumably are global (the intensity is contributed by all the gas along the line of sight). In addition, the conversion between ∞ −∞ dz dE(B − V )/ds should provide the rough estimate for the surface density of dust, we could see that the surface density of dust in the case of Green et al. (2019) decreases at large heliocentric distance. More investigations might be required to better understand the differences between these analyses. We caution that given the current disagreement in the dust reddening maps, it might not be straightforward to adopt these data sets as constraints for gas models.
Even though we would at this point not adopt the dust reddening data as a part of the input for our Bayesian analysis because of the issues mentioned above, it might be interesting to provide some qualitative comparisons between the gas map converted from dust reddening with our gas reconstructions. We took the three-dimensional maps of the dust reddening E(B − V ) from both Green et al. (2019) and Leike et al. (2022) and converted them into three-dimensional maps of hydrogen atoms (both atomic and molecular) using the simple relation provided in Peek The relation between the differential dust reddening and the gas density becomes We note that this relation has implicitly assumed a relation between the reddening and the amount of dust and, more importantly, a dust-to-gas ratio that might vary at different positions in our Galaxy (Tchernyshyov & Peek 2017). The surface mass density of hydrogen atoms converted from dust reddening data is presented in the top panels of Fig. A.2. We shall refer to the maps derived from Green et al. (2019) and Leike et al. (2022) respectively as the GSZ19 and LAE22 maps. The reconstruction by Nakanishi & Sofue (2016) is also shown in Fig. A.2 (bottom panel) together with our reconstructions from the BEG03 and SBM15 gas flow models (middle left and right panels). Note that we combined the Hi gas maps from this work with the H 2 gas maps from Mertsch & Vittino (2021) to produce the BEG03 and SBM15 maps of all hydrogen atoms. We also overlay the spiral structures fitted with data of molecular masers associated with very young high-mass stars, (Reid et al. 2019).
One can actually spot many structures which appear in both in the GSZ19 and in LAE22 maps and our reconstructions. We have marked some of these structures with red ellipses, for example, the segment of gas following the Local arm (roughly stretching between (8, 0) kpc and (7, −3) kpc in the gas map, region R1) or the clump of gas at around (10, 0.3) kpc (region R2). Some other similarities exist as well, but they might be less obvious and differ slightly among these maps. For instance, the clump of gas at (6, −1.5) kpc on the Sgr-Car arm (region R3) in the GSZ19 map seems to stretch along the line of sight similar to the corresponding structures in the BEG03 and SBM15 maps but it has lower surface density than the corresponding structures in our reconstructions. The LAE22 map shows, in contrast, that the structure in the R3 region is clustered around the Sgr-Car arm and does not stretch along the line-of-sight, but it has a higher surface density that is comparable to the surface density of the corresponding structure in our reconstructions. These differences and similarities seem to hold true also for the clump of gas around (9, −3) kpc on the Perseus arm (region R4). We also see that some structures are specific to one of these maps, especially on very small scales (≲ 1 kpc). Such differences might be related to the differences in the threedimensional reconstructed dust maps and the uncertainties of the gas flow model. We believe, however, that the majority of structures in the Green et al. (2019) and Leike et al. (2022) maps could be identified in the gas maps for the region within roughly 3 kpc from the Solar System.
At larger heliocentric distance, the conversion might be not straightforward between dust and gas. Yet, there are a few clusters of gas density that are likely misplaced in the BEG03 map, possibly also in the SBM15 map. The clearest such example is possibly in the region B2 at a distance of ∼ 6 kpc, see Fig. A.2. In the BEG03 map, there is significant gas density along longitude 134 • . The SBM15 map is based on a different rotation curve and places this emission at a much closer distance of 3 -4 kpc. The maser distance along this direction is ∼ 2 kpc (Reid et al. 2019) instead. . Surface mass density of hydrogen atoms reconstructed from HI4PI data using the BEG03 and SBM15 gas flow models (middle left and right panels respectively). Surface mass density of hydrogen atoms reconstructed by Nakanishi & Sofue (2016) (bottom panel). The Galactic spiral structure fitted from data of molecular masers as in Reid et al. (2019) is also overlaid on all the plots. The red (blue) ellipses mark some representative regions in which structures from our gas maps are seen (are not seen) in the dust maps.