© The Authors 2025

1 Introduction

The gas in the Milky Way exists in a number of different phases, characterised by their typical number densities and temperatures (e.g. Ferrière 2001). The element dominating throughout by number, mass, and volume filling fraction is hydrogen. Apart from the hot ionised medium, hydrogen occurs mostly electrically neutral, either in atomic (HI) or molecular form (H₂). Whereas H₂ exhibits very pronounced structures on small scales, HI is much more diffuse. Both show a structure on large, Galactic scales. For instance, HI and H₂ accumulate in Galactic discs with scale heights of 𝒪(100) pc and otherwise have been argued to form similar structures as the stellar population, for instance, spiral arms. The gaseous disc is not flat though, but warping, meaning that the position of the midplane differs from z = 0 kpc in Galactic coordinates. The vertical scale height is also not constant but flaring, that is, it varies from 100 pc at the solar position to values a factor of a few larger at further distances from the Galactic centre.

Information on the density and distribution of gas can be chiefly obtained from surveys of emission lines. For atomic hydrogen, this emission is due to the hyperfine transition at a frequency of ~1.4 GHz or a wavelength of ~21 cm. Due to the symmetry of the molecule, molecular hydrogen is a radiatively inefficient emitter; instead, molecular carbon monoxide (CO), which is largely believed to be well-mixed with H₂, can be used as a tracer (Bolatto et al. 2013). CO emits at a number of frequencies with the most important one being the 1 → 0 transition of ¹²CO at ~115 GHz or 0.26 cm. It should be noted that this emission line is often optically thick, hampering its usefulness in learning about the properties of dense cloud cores.

While radio surveys (Kalberla et al. 2005; HI4PI Collaboration 2016; Dame et al. 2001; Dame & Thaddeus 2022) of these emission lines nowadays offer sub-degree resolution, the distance resolution is comparatively bad. In fact, on Galactic scales, distance information can only be obtained from the Doppler shift of the emission lines if the emitting gas cloud has a relative velocity with respect to the observer. Deprojecting the spectral distribution observed along a line of sight into a spatial distribution thus requires information on the velocity field. Oftentimes, it is assumed that gas is in perfect circular rotation around the Galactic centre and that the rotation curve is known, that is, the rotation velocity as a function of Galacto-centric distance. Deviations exist on a variety of scales, either due to peculiar motions, for example locally initiated by winds of massive stars or supernova explosions, or due to larger structures, for example spiral arms and the Galactic bar. Finally, distance ambiguities arise because even for perfectly circular rotation, the projection on the line of sight can produce two possible distance solutions for a given Doppler shift. Along certain directions, that is, towards or away from the Galactic centre, this ambiguity results in a complete lack of distance information for purely circular rotation models.

Historically, the first systematic surveys of emission lines (Muller & Oort 1951; Ewen & Purcell 1951; Wilson et al. 1970) were quickly followed by the first reconstructions of the gaseous Galaxy, albeit in 2D (longitude and radius). For HI, the pioneering work by Westerhout (1957) and Oort et al. (1958) already showed some structure that can be interpreted as spiral arms. New gas line surveys at a higher resolution, higher sensitivity, and with less noise were followed by renewed efforts at a reconstruction of the three-dimensional (3D) density of atomic and molecular hydrogen. These studies adopted different strategies for dealing with near-far distance ambiguities: If data within the solar circle and towards the centre and anticentre are avoided, then, under the assumption of perfectly circular or elliptic rotation, there is no distance ambiguity, so some studies (Levine et al. 2006a,b) limited themselves to data outside the solar circle.

Assuming a Gaussian profile for the vertical distribution of gas in the inner Galaxy independent of position in the disc approximately results in a latitudinal line emission profile at a given velocity that consists of the sum of two Gaussians. The amplitude of both contributions is proportional to the gas density at both (near and far) distances. Thus fitting a sum of two Gaussians to the latitudinal profiles offers some distance resolution. This method is known as the double-Gaussian method and variants of this have been adopted by a number of authors (Clemens et al. 1988; Nakanishi & Sofue 2006; Pohl et al. 2008). Some issues due to spread-out low-intensity regions, which cannot be reasonably decomposed into few Gaussians, can be avoided by focussing only on the highest-intensity emission (Koo et al. 2017). Alternatively, under the assumption of an axisymmetric gas distribution, a deprojection into Galactocentric rings is possible (Marasco et al. 2017). Finally, one can circumvent the aforementioned reconstruction issues altogether by forward-modelling a parametric model of the gas distribution (Jóhannesson et al. 2018); however, this necessarily limits the complexity of the model and hence the fidelity in reproducing the data.

What all these approaches have in common is that different lines of sight are treated separately. However, to a certain degree, the processes determining the spatial distribution of gas guarantee that the gas density is smooth. For the reconstruction, the assumption of spatial smoothness provides some regularisation that can counteract the above-mentioned ambiguities. This idea is also at the core of a reconstruction of the Galactic velocity field by Tchernyshyov & Peek (2017). There, it is introduced as a cost term penalising the difference of velocity field values between neighbouring voxels in an angular direction. Operationally, such smoothness can more generally be enforced as a prior on the correlation structure of gas densities and velocities, encoded, for instance, as a covariance matrix. These priors can be readily provided in the context of a Bayesian reconstruction. In such a reconstruction, we are not seeking out one best-fit model for the gas density; instead we attempt to reconstruct the posterior probability of the gas distribution, that is, the likelihood conditioned on the data in addition to prior knowledge.

We have recently presented the first Bayesian reconstruction of Galactic gas densities in 3D, making use of methodology from information field theory (IFT, Enßlin et al. 2009; Enßlin 2019). Our 3D maps of Galactic HI (Mertsch & Phan 2023) and H₂ (Mertsch & Vittino 2021) are based on the HI4PI survey (HI4PI Collaboration 2016) and the Dame et al. (2001) survey of CO. At a spatial resolution of 1 / 16 kpc = 62.5 pc, our maps are the most detailed 3D HI and H2 maps of the Galaxy to date. They show a structure on a variety of scales that had not been resolved previously, and they agree with averages obtained on larger scales, for example the density profiles averaged over Galacto-centric rings. However, there are also some obvious artefacts, both in the 3D maps and in the column density maps inferred from our reconstructions. These are mainly due to the following two issues: the lack of (local) resolution and the assumption of a rigid velocity field. The lack of local resolution is a consequence of the adopted Cartesian grid which is not well-suited to account for small- and intermediate-sized structures close by as the nearest grid points cover large portions of the sky when viewed from our vantage point in the Galaxy. The assumption of a fixed velocity field also clearly constitutes an oversimplification.

For this work, we improved the reconstruction in a number of ways. Most importantly, we relaxed the constraints on the assumed velocity field by allowing for deviations from the assumed circular rotation pattern. The additional degrees of freedom can – to a certain degree – be constrained by independent information on the velocity field, for instance from Galactic masers and young stellar objects (YSOs). In order to obtain a matching ensemble of densities and velocities, we unified our previously separate reconstructions (Mertsch & Vittino 2021; Mertsch & Phan 2023) into one common inference scheme. We also adopted a non-Cartesian grid that is much better suited for astronomical data, offering the spatial resolution in nearby regions demanded by the data. Finally, we also treated the radiative transport of the oftentimes optically thick emission, thus significantly improving the fidelity of our maps.

The remainder of the paper is structured as follows. We begin by laying out our methodology in Sect. 2, which consists of our datasets (Sect. 2.1), modelling of radiation transport (Sect. 2.2), numerical discretisation (Sect. 2.3), modelling of the interstellar medium (Sects. 2.4–2.6), and our inference scheme (Sect. 2.7). In Sect. 3, we present and discuss our reconstruction. We start by presenting the reconstructed gas densities (Sect. 3.1) and velocities (Sect. 3.2) before comparing them to previous reconstructions thereof (Sect. 3.3) and 3D maps of interstellar dust (Sect. 3.4). We conclude with a summary and outlook in Sect. 4.

2 Methodology

In order to infer various Galactic parameters – most prominently the 3D gas densities of HI and CO, and the Galactic velocity field – we need to set up a model for these quantities and infer them from data with an appropriate algorithm. We describe the datasets used in Sect. 2.1, with the primary data being gas line surveys which measure the intensity of emission lines as a function of relative velocity (due to Doppler shift). In Sect. 2.2, we describe our modelling of the transport of radiation from its source somewhere in the galaxy to the observer, including absorption effects in both HI and CO emission. We describe our discretisation scheme in Sect. 2.3. In Sects. 2.4 – 2.6, we describe the modelling of our gas densities, emission line-widths and velocities in terms of Gaussian processes. This forward-modelling approach allows us to draw samples from our (prior) model by generating specific realisations of said Gaussian processes. From these realisations, we can compute synthetic data which can be compared to the actually measured data. The inverse problem – finding the specific realisation(s) that led to the actually measured data – is formulated as a Bayesian inference problem in Sect. 2.7.

2.1 Datasets

Previous approaches to reconstructing the distribution of gases in the Galaxy were usually restricted to reconstructing a single gas constituent by using a single dataset and a fixed velocity model in their reconstruction scheme (e.g. Kulkarni et al. 1982; Burton & te Lintel Hekkert 1986; Clemens et al. 1988; Nakanishi & Sofue 2003, 2006; Levine et al. 2006a; Kalberla & Dedes 2008; Pohl et al. 2008; Koo et al. 2017; Mertsch & Vittino 2021; Mertsch & Phan 2023, among others). A notable exception is Tchernyshyov & Peek (2017), using data from HI and CO radio line emission surveys, a dust reconstruction, and Galactic Masers associated with high-mass star-forming regions in a single framework. We attempt to improve on that and reconstruct multiple coherent ensembles of mainly the HI and CO gas densities and a common gas velocity field. To this end, we also combine multiple datasets. Our two biggest datasets are HI and CO line spectra described in Sects. 2.1.1 and 2.1.2. We also use data directly relating to the Galactic velocity field. These are measurement of the Galactic velocity curve (Sect. 2.1.3), and measurements of Galactic masers and young stellar clusters (Sects. 2.1.4 and 2.1.5).

2.1.1 HI4PI survey

The main dataset for reconstructing the HI gas density is the HI4PI full-sky survey (HI4PI Collaboration 2016) mapping the 21-cm emission of atomic hydrogen in the Galaxy. It is a merged dataset of two individual surveys of comparable sensitivity and resolution: The Effelsberg–Bonn HI Survey (EBHIS) (Kerp et al. 2011; Winkel et al. 2016) and the third revision of the Galactic All-Sky Survey (GASS) (McClure-Griffiths et al. 2009; Kalberla et al. 2010; Kalberla & Haud 2015). The EBHIS data was collected using the Effelsberg 100 m telescope between 2008 and 2013. The GASS observations were performed using the Parkes 64 m telescope between 2005 and 2006. The survey provides high-sensitivity spectra with exceptional resolution in velocity (frequency) and direction. The data is published in the form of the brightness temperature T_B of the HI emission as a function of direction (HEALPix¹ pixel number) and local standard of rest (LSR) velocity v^LSR. In order to limit the computational expense, we restrict our dataset to 257 points within |v_LSR| < 320 km s⁻¹ for every direction. Additionally, we degrade the survey from its original N_side = 1024 resolution to N_side = 64 by averaging, so that it matches with the numerical grid we will be using in our reconstruction. We also remove three known extragalactic features from the dataset as listed in Table 1, most importantly the Magellanic Clouds. The measurement uncertainty of every raw data point is estimated using $${\sigma }^{\mathrm{H}\mathrm{I}}\left(T_{\mathrm{B}}\right)=0.09\mathrm{K}+0.025T_{\mathrm{B}},$$(1)

based on the analysis in Winkel et al. (2016). We also compute the corresponding velocity in the laboratory frame for every data point by reverting the transformation to the LSR. The formula used for this is $$v^{\mathrm{l}\mathrm{a}\mathrm{b}}=v^{\mathrm{L}\mathrm{S}\mathrm{R}}-\hat{r}\cdot {\mathit{v}}_{\odot }^{\mathrm{L}\mathrm{S}\mathrm{R}},$$(2)

with the peculiar velocity of the Sun in the LSR frame given by ${\mathit{v}}_{\odot }^{\mathrm{L}\mathrm{S}\mathrm{R}}\approx (10.27,15.32,7.74)\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$ in Galactic coordinates² and the unit radial vector $\hat{r}$ in direction of measurement. Projections of the data over velocity and latitude are shown in Figs. 1 and 2 respectively.

2.1.2 CO surveys

For reconstructing the CO-density, we use primarily observations of the 1 → 0 rotational transition of CO. For this, we combine two datasets: The CO survey compilation presented in Dame et al. (2001) that maps CO emission mainly at low Galactic latitudes and the newer northern sky survey presented in Dame & Thaddeus (2022) that maps the entire northern sky. The surveys were conducted with a 1.2 m telescope at the Center for Astrophysics I Harvard & Smithsonian (CfA) in Cambridge, Massachusetts and a similar telescope on Cerro Tololo in Chile. While the first dataset covers basically all areas, where significant emission was reported, the second survey extends far above the Galactic plane and finds many smaller clouds that have not previously been catalogued. The data is published in the form of brightness temperature T_B of the CO-emission as a function of longitude l, latitude b and LSR-velocity v^LSR. Due to the uneven sampling of the surveys on the (already distorted) l-b-grid, we have a strongly varying angular density of data points when projecting the surveys on the 2-sphere using the equal-area HEALPix pixelisation. Considering this, we use the moment-masked version of both surveys and project them on an N_side = 64 sphere, keeping track of the number of data points contributing to every value in the projection. In overlap regions, we use the newer data set. We estimate the measurement uncertainty of every raw data point uniformly as σ^CO(T_B) = 0.18 K as suggested by the authors. After the projection, the uncertainty estimate will no longer be uniform but depend on the number of raw data values contributing to one projected data value. As for the HI dataset, we calculate the laboratory-frame velocities by reverting the transformation to the LSR. Projections of the combined data over velocity and latitude can be seen in Figs. 3 and 4 respectively.

2.1.3 The Galactic rotation curve

For modelling the large-scale circular motion of matter around the Galactic centre, we do not use the rotation curve by Reid et al. (2019) as a prior since we want to use the Masers as data (see Sect. 2.1.5) on the velocity field directly and incorporating them into the rotation curve prior would count them double. Instead, we merge two recent analyses of the Galactic velocity curve:

For the inner part of the galaxy (r_gal ⪅ 5.5 kpc), we use the rotation curve by Sofue (2021), derived by applying the tangent point method to CO-spectra from the FUGIN project (Umemoto et al. 2017). At larger radii (r_gal ⪆ 5.5 kpc), this dataset is seamlessly merged with the rotation curve by Zhou et al. (2023), derived by combining and analysing photometric and astrometric data from the APOGEE (Majewski et al. 2017; Abdurro’uf et al. 2022), LAMOST (Luo et al. 2012; Deng et al. 2012) and 2MASS (Skrutskie et al. 2006) surveys in addition to the Gaia EDR3 (Gaia Collaboration 2021) data.

For merging, we take all but the last four data points from the ‘dR100 pc’ table from Sofue (2021) and all but the first data point from Table 4 from Zhou et al. (2023). We linearly interpolate 1000 values between 0 kpc and 35 kpc which we then smooth with a Gaussian kernel with a width of 1 kpc. These values are used in the following to evaluate the smoothed velocity curve v_circ,0(r_gal) by linear interpolation. The resulting velocity curve together with the individual data points can be seen in Fig. 5 as an inset-plot in the bottom right corner. The flattening towards r_gal → 0 originates from the linear interpolation and the finite value of the leftmost data point. This may seem unphysical but is expected to be of minor importance since the affected volume is small and the reconstructed velocity can differ from this prior.

Fig. 1 Nside = 128 projection of the HI4PI data on the two-sphere in a Mollweide projection on a logarithmic scale.

Fig. 2 Longitude-velocity diagram of the HI4PI data (integrated over latitude) corresponding to Fig. 1.

Fig. 3 Nside = 128 projection of the CO data on the two-sphere in a Mollweide projection on a logarithmic scale. Some notable structures are indicated by arrows.

Fig. 4 Longitude-velocity diagram of the CO data (integrated over latitude) corresponding to Fig. 3.

Fig. 5 Positions of masers (circles) and YSO clusters (diamonds) shown in a 2D projection of the Galactic disc. The marker colour shows the corresponding vLSR measurement. The black streaks indicate the 1σ uncertainties on the parallax (distance) measurements. The Sun is located at (0, 0) and the centre of the Galaxy at (8.178, 0). The background shows vLSR computed from an interpolated rotation curve as described in Sect. 2.1.3. The inset shows the underlying rotation curve vcirc,0 and the data used for its creation (cf. Sect. 2.1.3).

2.1.4 YSO cluster

Since Galactic gas clouds will not move on perfectly circular orbits around the Galactic centre, it is important to gather additional information on their peculiar velocities. In general, stars do not necessarily move in unison with the gas clouds they were once formed in. While stars can be reasonably modelled as collision-less particles, the gas cannot. Thus, their dynamical time evolution can be expected to differ quite significantly. This difference will be more severe, the more time has passed since the formation of the stars. For young star clusters that have recently been formed in the centre of dense gas clouds, we expect the difference to be small and that the mean velocity of these young clusters traces the gas velocities reasonably well. We use the recently compiled catalogue of near-by stellar clusters by Hunt & Reffert (2023). They analyse Gaia DR3 (Gaia Collaboration 2023) data by applying the HBDSCAN clustering method (Campello et al. 2015) in a modified version (Hunt & Reffert 2021). We perform the quality cuts as suggested in the appendix of Hunt & Reffert (2023) and select clusters with a median stellar age smaller than 20 Myr and at least six stars with a measured radial velocity. This is a compromise between selecting reliable data and retaining enough clusters (235) to still be of use. For every selected cluster, we use the measured parallax, the radial heliocentric velocity, the proper motion in right ascension and declination, and corresponding uncertainties. A top-down view with this data can be seen in Fig. 5 with stars indicating the positions of YSO clusters and their fill colour indicating the measured v_LSR.

2.1.5 Galactic masers

Additional information about the velocity field, especially for larger distances, can be obtained from measurements of Galactic masers. In conditions with enhanced excitation, molecules in gas clouds can undergo a population inversion which leads to strongly amplified emission lines. Since this emission clearly stands out from the background, relative velocity and parallax of this emission can be accurately measured. We use data of 199 Maser sources compiled in Reid et al. (2019), which is based mainly on results from the BeSSeL survey (Brunthaler et al. 2011) and VERA project (Kobayashi et al. 2005). We remove 31 Masers with a fractional parallax uncertainty exceeding 20% from the dataset because we cannot expect to adequately sample their distances without an excessive number of prior samples. The transformation of velocities to the LSR frame is reverted again using Eq. (2). The data of these masers will simply be appended to the YSO-cluster data. A top-down view with this data can be seen in Fig. 5 with circles indicating the positions of Galactic Masers and their fill colour indicating the measured v_LSR.

2.2 Radiation transport

In order to model the measured line spectra, we have to specify the measurement process. For both gas surveys, this will involve a modelling of the transport of radiation from the origin of its emission anywhere in the Galaxy to the position of the Sun. While the effects of absorption are especially important for the often optically thick emission lines of the first rotational transition of CO, there are also directions where it makes a sizeable difference for HI (Pohl et al. 2022). An optically thick, saturated emission line effectively hides the true amount of emitting gas from sight. There could be just enough gas to saturate the line, or arbitrarily more. In the reconstruction, we want our samples to reflect this non-linearity. Therefore, we have to include it in our modelling.

The evolution of specific intensity I_ν at frequency ν along a line of sight is described by the radiation transport equation (e.g. Kulkarni & Heiles 1988; Draine 2011), $$\mathrm{d}{I}_{\nu }=-I_{\nu }{\kappa }_{\nu }\mathrm{d}s+j_{\nu }\mathrm{d}s,$$(3)

where κ_ν is the attenuation coefficient at frequency ν, j_ν is the emissivity at frequency ν, and s denotes the line-of-sight distance from the observer. The first term describes the net change of the propagating radiation due to absorption and stimulated emission. The second term describes the increase of intensity due to spontaneous emission. We neglect any additional sources of radiation here for simplicity, in particular the contribution by continuum emission. The emissivity of the transition for randomly oriented emitting molecules or atoms is given as $$j_{\nu }=\frac{1}{4\pi }{n}_{\mathrm{u}}{A}_{\mathrm{u}\mathrm{l}}h{\nu }_{\mathrm{u}\mathrm{l}}{\varphi }_{\nu },$$(4)

and the attenuation coefficient as $${\kappa }_{\nu }=n_1\frac{g_{\mathrm{u}}}{g_1}\frac{A_{\mathrm{u}\mathrm{l}}}{8\pi }{\lambda }_{\mathrm{u}\mathrm{l}}^2(1-\frac{n_{\mathrm{u}}}{n_1}\frac{g_{\mathrm{u}}}{g_{\mathrm{l}}}){\varphi }_{\nu },$$(5)

with the multiplicities of upper (u) and lower (1) energy levels³ g_u = 3 and g₁ = 1, their number densities n_u and n₁, the decay rate A_ul, the Planck constant h, the frequency (and wavelength) ν_ul (and λ_ul), and the line profile ϕ_ν (Draine 2011). It is important to note that the attenuation coefficient, and therefore the process of radiation transport, depends on the number density and velocity of the gaseous emitters – quantities that we aim to infer.

Integrating Eq. (3) (while neglecting all extragalactic contributions) and expressing it in terms of the antenna temperature, $$T_{\mathrm{A}}(\nu )=\frac{c^2}{2k_B{\nu }^2}{I}_{\nu },$$(6)

yields the following solution to the transport equation: $$T_{\mathrm{A}}(\nu )={\int }_0^{\mathrm{\infty }}{S}_{\nu }^{\mathrm{\prime }}\exp (-{\int }_0^s{\kappa }_{\nu }\left(s^{\mathrm{\prime }}\right)\mathrm{d}{s}^{\mathrm{\prime }}){\kappa }_{\nu }(s)\mathrm{d}s.$$(7)

Here, we have used the source function $S_{\nu }^{\mathrm{\prime }}=\frac{j_{\nu }}{{\kappa }_{\nu }}\frac{c^2}{2k_B{\nu }^2}$, the speed of light in vacuum c and the Boltzmann constant k_B. What is left now, is to compute the attenuation coefficient and the source function for both HI and CO.

For HI, we have a very large spin (excitation) temperature compared to the transition energy, that is $k_B{T}_{\text{spin\hspace{0.17em}}}\gg h{\nu }_{\mathrm{u}\mathrm{l}}^{\mathrm{H}\mathrm{I}}$. Therefore, the upper energy level is always populated and we can approximate $n_{\mathrm{u}}^{\mathrm{H}\mathrm{I}}/n_1^{\mathrm{H}\mathrm{I}}\approx g_{\mathrm{u}}/g_1=3$. The attenuation coefficient at the relative velocity v (corresponding to the frequency ν by Doppler-shift) then takes the (simple) form $${\kappa }_v^{\mathrm{H}\mathrm{I}}\approx 679.168{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\left(\frac{\mathrm{K}}{T_{\text{spin\hspace{0.17em}}}}\right)\left(\frac{n^{\mathrm{H}\mathrm{I}}}{{\mathrm{c}\mathrm{m}}^{-3}}\right)\left(\frac{\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}{{\sigma }_v^{\mathrm{H}\mathrm{I}}}\right)e^{-v^2/2{\left({\sigma }_v^{\mathrm{H}\mathrm{I}}\right)}^2},$$(8)

with the intrinsic line-width ${\sigma }_v^{\mathrm{H}\mathrm{I}}$ of the emitted radiation under the assumption of a Maxwellian line profile. The source function can be calculated to be $$S_v^{\mathrm{\prime }\mathrm{H}\mathrm{I}}=T_{\text{spin}}.$$(9)

While we have to acknowledge that T_spin is not a uniform constant across the whole Galaxy, we have to balance physical accuracy with spatial resolution. We choose not to model T_spin as a full 3D field, because we expect that information about this quantity is limited to very few regions in our datasets. We will instead fix the HI spin temperature at T_spin = 200 K, as suggested in Pohl et al. (2022) as being able to well reproduce HI spectra towards the Galactic centre region. This aims at improving our modelling of the dense and cold regions where (self-)absorption is expected to be the most important.

For CO there are infinitely many rotational energy levels but only the lowest two contribute to the transition that concerns us. Thus, we have to compute the fraction of molecules in the J = 0 and J = 1 levels to infer the total CO number density. The population ratio is given by $$\frac{n_{\mathrm{u}}^{\mathrm{C}\mathrm{O}}}{n_1^{\mathrm{C}\mathrm{O}}}=\frac{n^{\mathrm{C}\mathrm{O}}(J=1)}{n^{\mathrm{C}\mathrm{O}}(J=0)}=3\text{ߙ}\exp (-\frac{2B_0}{k_B{T}_{\mathrm{e}\mathrm{x}\mathrm{c}}})\approx 3\text{ߙ}\exp (-\frac{5.53\mathrm{K}}{T_{\mathrm{e}\mathrm{x}\mathrm{c}}}),$$(10)

with the rotation constant B₀ and the excitation temperature T_exc. The relation to the total CO density can be approximately expressed, for instance for the lower energy level, as $$n_1^{\mathrm{C}\mathrm{O}}=\frac{n^{\mathrm{C}\mathrm{O}}}{Z}\approx \frac{n^{\mathrm{C}\mathrm{O}}}{\sqrt{1+{\left(\frac{T_{\text{exc\hspace{0.17em}}}}{2.77\mathrm{K}}\right)}^2}},$$(11)

by using the approximation to the partition function Z from Draine (2011). This leaves us with the attenuation coefficient $$\begin{array}{rl}{\kappa }_v^{\mathrm{C}\mathrm{O}}\approx & 1.75\times {10}^6{\mathrm{k}\mathrm{p}\mathrm{c}}^{-1}\left(\frac{1}{\sqrt{1+{\left(\frac{T_{\mathrm{e}\mathrm{x}\mathrm{c}}}{2.77\mathrm{K}}\right)}^2}}\right)(1-\exp (-\frac{5.53\mathrm{K}}{T_{\mathrm{e}\mathrm{x}\mathrm{c}}}))\\ & \cdot \left(\frac{n^{\mathrm{C}\mathrm{O}}}{{\mathrm{c}\mathrm{m}}^{-3}}\right)\left(\frac{\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}}{{\sigma }_v^{\mathrm{C}\mathrm{O}}}\right)e^{-v^2/2{\left({\sigma }_v^{\mathrm{C}\mathrm{O}}\right)}^2}.\end{array}$$(12)

The source function for CO can then be calculated to be $$S_v^{\mathrm{\prime }\mathrm{C}\mathrm{O}}\approx 5.53\mathrm{K}\left(\frac{1}{\exp \left(\frac{5.53\mathrm{K}}{T_{\text{exc\hspace{0.17em}}}}\right)-1}\right)\approx 5.55\mathrm{K},$$(13)

for an excitation temperature of T_exc = 8 K that we will keep fixed at this value, again aiming at improving our description of the cold cores of molecular clouds. The actual excitation temperature is expected to be spatially varying, with higher values towards the Galactic centre (cf. Fixsen et al. 1999, ~40 K) and smaller values in outer regions (cf. Burgh et al. 2007, ~4 K for the diffuse to translucent regime).

2.3 Numerical discretisation

For computational modelling of the radiation transport, and later also the gas densities, we have to represent all quantities of interest on a discrete grid. This grid should suit our purpose and data. Since the gas line surveys are given as functions of direction times relative velocity with the latter corresponding to some radial distance, a direction times radius grid is a good choice. This will also speed up the necessary line-of-sight integration for performing the radiation transport since it will reduce to a relatively simple sum along the radial grid dimension. Together with the condition that the individual voxels shall be as undistorted as possible, that is their separations in all directions should be similar, leaves us with logarithmically spaced grid points in radial direction⁴. This coincidentally leaves us with an increased resolution nearby which is welcome since structures within a couple of hundred parsecs cover large angular scales. More precisely, we discretise our 3D volume into HEALPix spheres at uniform distances in logarithmic radius log(r), denoted in the following as a HEALPix × log(r)-grid. We use a resolution of N_side = 64 for the HEALPix spheres, which equals 49 152 grid points distributed on each sphere. In the radial direction, we use 388 grid points, for a total of slightly over 19 · 10⁶ grid points.

For numerically evaluating Eq. (7) on this grid, we approximate the integrals by sums over radial grid points. We number the points of the outer integral by i and the points of the inner integral by j < i. The case j = i is treated separately by calculating the average over the upper integration bound. The increment ds is replaced by the cell size Δs_i. This yields the approximate discretised formula $$\begin{array}{rl}{T}_{\mathrm{A}}(v)=& \sum _i{S}_v^{\mathrm{\prime }}\exp (-\sum _{j<i}{\kappa }_v\left(s_j\right)\mathrm{\Delta }{s}_j)\\ & \cdot \left(\frac{1}{{\kappa }_v\left(s_i\right)\mathrm{\Delta }{s}_i}(1-\exp (-{\kappa }_v\left(s_i\right)\mathrm{\Delta }{s}_i))\right){\kappa }_v\left(s_i\right)\mathrm{\Delta }{s}_i,\end{array}$$(14)

where the factors κ_v(s_i, v)Δs_i cancel, but are left here for completeness. The large expression in brackets results from averaging the contribution of the bin with i = j, in which the emitting radiation will partly self-absorb. In this form the calculation is computationally simple to perform since, using our radial numerical grid, all the sums along the line of sight become simple sums along an array axis.

2.4 Gas densities

A key requirement for our model of gas densities is that it should include spatial correlations between neighbouring lines of sight, thus enforcing a spatially coherent density field. Additionally, we want to encode our prior knowledge regarding the shape of the Galaxy into our model. Therefore, we model the number densities as $$n^i(\mathit{r})=f_{\mathrm{r}}^i\left(r_{\mathrm{g}\mathrm{a}\mathrm{l}}\right)\cdot f_{\mathrm{z}}^i(z)\cdot a_n^i\cdot \exp (s_n^i\cdot g_n^i(\mathit{r})),$$(15)

for i ∈ {HI, CO} respectively. Here, $g_n^i$ are 3D Gaussian processes, $f_r^i$ and $f_z^i$ are one-dimensional Gaussian processes, and $a_n^i$ and $s_n^i$ are scalar parameters. The (spatial) arguments are r = (x, y, z)^T for the position vector in Galactic coordinates, and r_gal for the polar Galacto-centric radius.

The 3D homogeneous and isotropic Gaussian processes gⁱ are drawn according to normalised Matérn-1/2 covariance functions (or kernels)⁵. On small scales the power spectrum of this kernel approximately matches the Kolmogorov spectrum, that is a –11/3 power-law slope (Kolmogorov 1941). For generating these Gaussian processes on our HEALPix × log(r)-grid, we use an algorithm called Iterative Charted Refinement (Edenhofer et al. 2022) which ensures a linear scaling of required computations in the number of grid points. Inferring the covariance kernel with this method proved to be numerically unstable, so we fix the kernels with different length scales for HI (300 pc) and CO (75 pc), which have led to consistent results. These scales determine the preferred size of coherent structures. We work under the assumption that the Gaussian process is homogeneous. This needs to checked a posteriori and we find the chosen kernel size to be in good agreement with this hypothesis. We leave a more in-depth exploration of the kernel space to future work. It should be noted, that this still allows for significant amounts of small-scale fluctuations and coherent gas structures are not strictly limited in size by the chosen length scale.

These Gaussian processes are multiplied by the scalars $s_n^i$, following a Gaussian prior distribution, that specify the global scale of fluctuations and exponentiated to give a log-normal random field. This ensures the necessary positivity of the number density and also allows for the observed large dynamical range in density. The factors $a_n^i$, also following a Gaussian prior distribution, act as a global offset to $g_n^i$, thereby specifying the log-average values of the gas densities.

The assumption of a global homogeneous random field is rather crude since the gases are localised inside a disc with finite extent in radius and height. Ignoring this can lead to an overestimation of distances, and would allow for significant amounts of gas almost arbitrarily far away from the Galactic disc. To increase the similarity of our model with the Galactic disc, we multiply it by two one-dimensional log-normal processes $f_r^i$ (radial profile) and $f_z^i$ (vertical profile), depending on the polar Galactic radius r_gal and the Galactic height z respectively. This effectively amounts to modelling a random field with a non-homogeneous correlation structure. The radial profiles are modelled as $$f_{\mathrm{r}}^i\left(r_{\mathrm{g}\mathrm{a}\mathrm{l}}\right)=1\cdot \exp (0.2\cdot g_{f_{\mathrm{r}}}^i\left(r_{\mathrm{g}\mathrm{a}\mathrm{l}}\right)),$$(16)

with the kernels for $g_{f_r}^i$ being normalised Gaussian kernels with widths of 4 kpc for HI and 2 kpc for CO. The width of these kernels is chosen so that only the overall large-scale gas distributions are represented, but not individual small-scale features. We chose the rather small scale of fluctuations (0.2) as to not overestimate the large-scale radial symmetry of the Galactic disc. The exponential function again ensures positivity of the gas densities.

The vertical profiles are modelled as $$f_{\mathrm{Z}}^i(z)=\exp ([K_{f_{\mathrm{z}}}^i\ast (-\frac{|\cdot |}{z_{\mathrm{h}}^i})](z))\cdot \exp (0.1\cdot g_{f_{\mathrm{z}}}^i(z)),$$(17)

with the Kernels $K_{f_z}^i$ for $g_{f_z}^i$ being normalised Gaussian kernels with widths of 100 pc for HI and 50 pc for CO. The first exponential acts as a prior to these profile functions, with an exponential decay in height |z| as motivated from previous studies (e.g. Kalberla & Dedes 2008; Gaensler et al. 2008; Mertsch & Phan 2023). The convolution in the exponent ensures a smoothness matching that of the Gaussian process. The corresponding length scales ($z_{\mathrm{h}}^{\mathrm{H}\mathrm{I}}=500\mathrm{p}\mathrm{c},z_{\mathrm{h}}^{\mathrm{C}\mathrm{O}}=200\mathrm{p}\mathrm{c}$) are chosen deliberately large as to not overestimate the confinement of the gases to the disc.

For the solar radius R_⊙, which is needed for the calculation of the Galacto-centric radius r_gal, we take as prior the value R_⊙ = 8.178(35) kpc from GRAVITY Collaboration (2019).

2.5 Emission line widths

The line widths ${\sigma }_v^i$ with i ∈ {HI, CO} of the radiation emitted by the gases at some relative velocity v is dictated by the intrinsic gas temperature and the turbulent local fluid motions (Draine 2011). We will not attempt to model these effects separately, but instead model the phenomenological line-width directly as $${\sigma }_v^i(\mathit{r})={\sigma }_{v,0}^i+a_{\sigma }^i\cdot \exp (s_{\sigma }^i\cdot g_{\sigma }^i(\mathit{r})),$$(18)

for i ∈ {HI, CO}, respectively. Here, $g_{\sigma }^i$ are 3D Gaussian processes with the same covariance kernels as their densitycounterparts $g_n^i$, and $a_{\sigma }^i$ and $s_{\sigma }^i$ are scalar parameters. We add a constant scalar value ${\sigma }_{v,0}^i$, effectively setting a minimum allowed value, which corresponds to the coarsest data-spacing in velocity-space of their corresponding main datasets, described in Sects. 2.1.1 and 2.1.2. In particular, we fix ${\sigma }_{v,0}^{\mathrm{H}\mathrm{I}}=2.5\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$ and ${\sigma }_{v,0}^{\mathrm{C}\mathrm{O}}=1.3\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$. This is motivated mainly by numerical considerations, improving the stability of our algorithm by removing extremely narrow, high-intensity peaks. We do not multiply these fields by profile functions, like we did for the gas number densities, as we do not expect them to feature such a pronounced large-scale structure. Re-using the same correlation kernel as for the gas densities is partly numerically motivated since intermediate variables in enforcing the desired kernel can be re-used, leading to a reduced overall memory consumption. This choice also naturally leads to variations in the line-width on the same spatial scale as the cloud sizes.

2.6 Velocity model

For deprojecting the gas line surveys, we also need to model the line-of-sight Galactic gas velocities relative to the Sun. Our velocity model will consist of three parts: a circular motion part from a Galactic rotation curve; a 3D Gaussian process contribution; and the peculiar motion of the Sun relative to purely circular motion. In particular, we model the relative velocity of gases at position r to the Sun as $$\mathrm{\Delta }\mathit{v}(\mathit{r})=\underset{{\mathit{v}}_{\mathrm{g}\mathrm{a}\mathrm{s}}}{\underbrace{{\mathit{v}}_{\mathrm{c}\mathrm{i}\mathrm{r}\mathrm{c}}(\mathit{r})+\left(\begin{array}{c}{s}_{\mathrm{x}}^v\cdot g_{\mathrm{x}}^v(\mathit{r})\\ s_{\mathrm{y}}^v\cdot g_{\mathrm{y}}^v(\mathit{r})\\ s_{\mathrm{z}}^v\cdot g_{\mathrm{z}}^v(\mathit{r})\end{array}\right)}}-\underset{{\mathit{v}}_{\odot }}{\underbrace{\left(\begin{array}{c}{U}_{\odot }\\ V_{\odot }+v_{\text{circ\hspace{0.17em}}}(R_{\odot })\\ W_{\odot }\end{array}\right)}},$$(19)

with the contribution from purely circular motion v_circ, its amplitude v_circ, three 3D Gaussian processes $g_x^v,g_y^v$ and $g_z^v$, three scalar parameters $s_x^v,s_y^v$ and $s_z^v$, and the solar velocity relative to purely circular motion (U_⊙, V_⊙, W_⊙)^T. Together, the first two terms form the total gas velocity v_gas as seen from a non-rotating, stationary observer and the last term forms the total velocity of the Sun v_⊙ in the same frame.

The circular motion contribution is split as $$v_{\text{circ}}(\mathit{r})=v_{\text{circ},0}\left(r_{\text{gal}}\right)+g_{\text{circ}}^v\left(r_{\text{gal}}\right)\cdot 0.5\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1},$$

into a constant contribution v_circ,0 as described in Sect. 2.1.3 and deviations thereof drawn from a normalised one-dimensional Gaussian process with a Gaussian kernel with width of 1 kpc. The average amplitude of these deviations is chosen to be optimistically small, with further deviations having to be captured by the 3D Gaussian processes.

These Gaussian processes $g_x^v,g_y^v$ and $g_z^v$ have the same covariance kernel as the Gaussian process for the CO-number density, that is a normalised Matérn-1/2 covariance kernel with a length scale of 75 pc. This ensures that they can capture velocity fluctuations on scales of CO-clouds. There is a notable difference though, namely that they are not exponentiated and so the resulting fluctuation field is normal and not log-normal. We multiply these fields by the scalar parameters $s_x^v,s_y^v$ and $s_z^v$ to set the global scale of average fluctuations. Importantly, these fields are able to produce deviations from a purely circular motion around the Galactic centre, which is a very important ingredient, especially for modelling local gas.

From Eq. (19), we can compute the line-of-sight component by projection onto radial basis vectors as $$v_{\mathrm{l}\mathrm{o}\mathrm{s}}(\mathit{r})=\hat{r}\cdot \mathrm{\Delta }\mathit{v}(\mathit{r}),$$(20)

and use these in Eqs. (12) and (8) when computing the radiation transport.

We can also calculate the projection on the plane of sky to obtain the astronomical proper motion μ. In particular, we compute $${\mu }_b=\frac{\mathrm{d}b}{\mathrm{d}t}=\frac{\cos (b)}{r}(\mathrm{\Delta }{v}_z-\tan (b)[\mathrm{\Delta }{v}_x\cos (l)+\mathrm{\Delta }{v}_y\sin (l)]),$$(21) $${\mu }_l\cos (b)=\frac{\mathrm{d}l}{\mathrm{d}t}\cos (b)=\frac{1}{r}(\mathrm{\Delta }{v}_y\cos (l)-\mathrm{\Delta }{v}_x\sin (l)),$$(22)

with the longitude l, latitude b, and distance from the Sun r. The latter can be computed via $$r[\mathrm{p}\mathrm{c}]=(\tan (p[\mathrm{arcsec}]){)}^{-1},$$(23)

where p are the parallaxes obtained from measurements (cf. Sects. 2.1.4 and 2.1.5). To convert this proper motion from Galactic coordinates to equatorial coordinates, we apply the transformation matrix (cf. Poleski 2013) to obtain $$(\genfrac{}{}{0ex}{}{{\mu }_{\alpha }\cos (\delta )}{{\mu }_{\delta }})=\frac{1}{\cos (b)}\left(\begin{array}{c}{C}_1-C_2\\ C_2{C}_1\end{array}\right)(\genfrac{}{}{0ex}{}{{\mu }_l\cos (b)}{{\mu }_b}),$$(24)

with α being the right ascension, δ being the declination, and $$C_1=\sin \left({\delta }_G\right)\cos (\delta )-\cos \left({\delta }_G\right)\sin (\delta )\cos (\alpha -{\alpha }_G),$$(25) $$C_2=\cos \left({\delta }_G\right)\sin (\alpha -{\alpha }_G),$$(26)

being the matrix coefficients with the constants α_G ≈ 192.86° and δ_G ≈ 27.13°. These can be directly compared to the measured proper and radial motions described in Sects. 2.1.4 and 2.1.5.

2.7 Bayesian formulation

So far, we have set up a model of the gas densities, velocities, and line-widths in terms of Gaussian processes and described the transport and measurement process of the gas line emission. Chaining these individual models yields a forward model going from a physical state to the expected measurement data. Our goal here is to invert this forward model and learn about the true physical state that led to the actually measured data. Attempting to do so, we cannot expect to actually find the one true physical state of gases and velocities in our universe due to a number of different reasons. For one, our model is only an approximate description of reality, with many strong assumptions built in. Also, the measured data is not absolute (it comes with significant measurement uncertainties) and can be described by infinitely many different physical states (our problem is underdetermined). For some lines of sight, there is even a total lack of direct information as, for instance, large parts of the sky have not been observed in CO-emission. Thus, we should strive to quantify not only what possible physical states describe the measured data, but also their probability, thereby estimating the uncertainty of our reconstruction.

In a Bayesian framework, this translates to the following: ee aim to infer the posterior probability, $$\begin{array}{rl}P(n^{\mathrm{H}\mathrm{I}},n^{\mathrm{C}\mathrm{O}},{\mathit{v}}_{\mathrm{g}\mathrm{a}\mathrm{s}},\dots \mid d)=& \frac{1}{P(d)}\cdot P(d\mid n^{\mathrm{H}\mathrm{I}},n^{\mathrm{C}\mathrm{O}},{\mathit{v}}_{\mathrm{g}\mathrm{a}\mathrm{s}},\dots )\\ & \cdot P(n^{\mathrm{H}\mathrm{I}},n^{\mathrm{C}\mathrm{O}},{\mathit{v}}_{\mathrm{g}\mathrm{a}\mathrm{s}},\dots ),\end{array}$$(27)

of multiple Galactic parameters, most prominently among them the gas densities n^HI and n^CO, and the Galactic velocity field v_gas. It is constrained by some general measurement data d and implicitly also conditional to our modelling. This problem is of very high dimensionality, as – in principle – every grid point in our representation of the gas densities and velocities is a free parameter. Therefore, instead of attempting to infer the complete posterior probability, which is utterly unfeasible, we instead aim to obtain a representative set of samples thereof. The extremely high dimensionality rules out typical sampling methods, such as MCMC algorithms. We instead use an inference method called Geometric Variational Inference (geoVI; Frank et al. 2021), which approximates the posterior by a multivariate normal distribution in a coordinate system induced by the Fisher information metric. This algorithm is a generalisation of Metric Gaussian Variational Inference (MGVI; Knollmüller & Enßlin 2019), the algorithm used in our previous works on Galactic gas tomography (Mertsch & Vittino 2021; Mertsch & Phan 2023), and is expected to greatly improve our representation of nonGaussianities in the actual posterior. The algorithm works in an iterative manner: first, it is drawing samples around an expansion point according to the local curvature of the posterior probability, approximated by use of the Fisher information metric. These samples now serve as a stochastic approximation to the true posterior. Secondly, the expansion point is updated by minimising the ‘distance’ between the stochastically approximated posterior and the true posterior by means of the Kullback-Leibler-divergence (Kullback & Leibler 1951). This iteration is repeated until convergence, which happens when the estimated local curvature and the expansion point are self-consistent. The result is a collection of samples, $${(n_i^{\mathrm{H}\mathrm{I}},n_i^{\mathrm{C}\mathrm{O}},{\sigma }_{v,i}^{\mathrm{H}\mathrm{I}},{\sigma }_{v,i}^{\mathrm{C}\mathrm{O}},\mathrm{\Delta }{\mathit{v}}_i,U_{\odot ,i},V_{\odot ,i},W_{\odot ,i},R_{\odot ,i},p_i)}_{i=1\dots N},$$(28)

that approximate the posterior probability and implicitly contain the correlations between all the model parameters. Our priors are summarised in Table 2. The exact choices for the a and s parameters were guided by results of reconstruction attempts during testing the algorithm. The result of the reconstruction proved to be very insensitive to their exact choice, indicating that they are very tightly constrained by the data.

The inference algorithm is implemented in the publicly available python code-package NIFTY8 (Selig et al. 2013; Steininger et al. 2019; Arras et al. 2019; Edenhofer et al. 2024a), and has recently undergone a complete re-write. It is now based on JAX (Bradbury et al. 2018), which enables harnessing the power of auto-differentiation methods and GPU-accelerated computations. This makes evaluation of the likelihood and its derivatives quick and simple, thereby giving access to, for example, Fisher matrices without having to resort to finite difference methods. We ran our reconstruction on a single node with two NVIDIA Volta 100 GPU’s on the RWTH High Performance Computing cluster for four weeks. We used eight posterior samples for our reconstruction, half of which are antithetical mirror samples of the others.

3 Results and discussion

Our main results consist of eight posterior samples (cf. Eq. (28)), each containing all the reconstructed quantities. All data products are available via zenodo⁶. The gas densities are in units of cm⁻³; velocities and line-widths in km s⁻¹; the solar radius R_⊙ in kpc; and the parallaxes in mas. Especially the line-widths may strongly depend on the voxel volume as they contain the effects of any non-resolved gas motions in addition to small-scale turbulence and temperature effects, so they should only be interpreted with a great deal of caution.

A general feature of our reconstruction is the irregular resolution, a consequence of the HEALPix × log(r)-grid. The innermost (outermost) pixels have a radial extent of 0.82 pc (495 pc). At a radius of 8 kpc, the radial extent of our pixels reaches 132 pc. The angular resolution is approximately 55′ and independent of the radius.

Even though we primarily derived the distributions of HI and CO, we will present the latter after transforming it to H₂ via the commonly used linear relation between the H₂ column density N_H2 and the velocity-integrated CO brightness temperature. Using this, the local H₂ density can be expressed as a differential version of the measurement function along a line of sight $$n^{{\mathrm{H}}_2}(\mathit{s})=X_{\mathrm{C}\mathrm{O}}\cdot \int \mathrm{d}v\frac{\mathrm{d}{T}_{\mathrm{A}}}{\mathrm{d}s},$$(29)

with X_CO = 2 × 10²⁰ K⁻¹ (km/s)⁻¹ cm⁻² being the customarily adopted conversion factor (Bolatto et al. 2013)⁷. As before, the synthetic temperature data is given by Eq. (7), and computed via its discretised version (Eq. (14)). We also define a total Hydrogen density as $$n^{\text{tot\hspace{0.17em}}}=n^{\mathrm{H}\mathrm{I}}+2\cdot n^{{\mathrm{H}}_2}.$$(30)

Derived quantities thereof, such as projections like the surface mass densities Σ, are indexed accordingly.

3.1 Reconstructed gas densities

In Fig. 6, we show the sample average of the total gas surface density of the galactic disc in a top-down view. An animated GIF flipping through the individual reconstruction samples and interactive 3D visualisations of our reconstructed total gas density are available online⁸, including a density scatter plot, an isosurface rendering, and a full volumetric rendering (more computationally demanding). The Galactic centre is in the middle of the figure at approximately (x, y) = (8, 0), while the Sun is at the coordinate origin.

The gas is not solely concentrated in a few localised spots, but distributed diffusely throughout the Galactic disc with local accumulations on all scales. This is also true for the individual reconstructed samples, not only for the sample average. Additionally, as a consequence of our irregular grid, the gas distribution lacks small-scale structure far away from the Sun’s position.

We have overlaid the pattern of spiral arms derived by Reid et al. (2019) using the galactic Masers also mentioned in Sect. 2.1.5. While we can identify some structures in our gas maps consistent with these patterns, especially in the Perseus-, Sgr-Car-, and the inner part of the Sct-Cen-arm, there are also many regions traced by the spiral pattern that are almost completely devoid of gas. We do not expect our model to be able to deduce definitely, whether the Milky Way interstellar gas exhibits a pronounced spiral arm pattern, since the data we use cannot be expected to be conclusive in this matter. Additionally, the expectations obtained by analysing the line spectra (in particular longitude-velocity-diagrams), can be misleading, as large-scale velocity perturbations can create illusory spiral arm imprints, as showcased in Peek et al. (2022). We expect our model to be partly susceptible to this as well, with elevated uncertainty estimations in regions of enhanced distance ambiguity. Estimations of the number of spiral arms also differ strongly (e.g. two in Churchwell et al. 2009, four in Reid et al. 2019), and no clear winner has yet emerged (for a recent review on this topic, see e.g. Sellwood & Masters 2022). Even the processes governing the generation and maintenance of spiral patterns in galaxies are generally not sufficiently understood yet.

We compare our azimuthally averaged surface mass densities of HI and H₂ with previous studies thereof in Figs. 7 and 8. The HI gas density is slightly peaked in the Galactic centre but falls off very quickly within r_gal ≈ 1 kpc. Then the density steadily increases until r_gal ≈ 12 kpc where it reaches its maximum, roughly consistent with all the other data points. For larger radii though, our reconstructed HI gas density quickly decreases, leading to an overall more compact galaxy.

Part of this can be explained by our prior rotation curve, which falls off for radii larger than approximately 7 kpc which naturally leads to a less extended gas distribution. While an exhaustive discussion regarding the validity of these rotation curves is beyond the scope of this paper, it can be noted that the declining outer part of the rotation curve by Zhou et al. (2023) that we use, seems to be largely consistent with many other recent studies, such as Eilers et al. (2019), Sylos Labini et al. (2023), Jiao et al. (2023), Ou et al. (2024), and Wang et al. (2023). Comparing our rotation curve with a flat rotation curve, as used by many previous studies, we expect this to shift gas on scales ≲2 kpc for the radii considered. We believe the major reason for this strong decline in gas density to be an issue with the radial profile function: This function strongly favours radially symmetric gas distributions, which poorly approximates the true structure – especially in a primarily circular-rotation model. Previous studies (discussed in e.g. Levine et al. 2006b) show a strong asymmetry between the northern and southern parts of the HI disc in both the radial and vertical large-scale distribution. A more compact Galaxy, where gas is strongly confined within some maximum radius, can, although considered implausible to such an extent, explain the data and resolve the asymmetries. As discussed later (in Sect. 3.2), the reconstructed velocity rotation curve also steepens, which is (consistently) assisting this effort to move gas inwards. A more restrictive rotation curve and radial profile will generally lead to a worse data fit overall and introduces artefacts at large radii, particularly along lines of sight towards the Galactic centre and anticentre. Ultimately, we deem this steep falloff in density physically implausible, yet it can account for the data.

Other data points also show a decrease in density for radii r_gal ≥ 15 kpc, but at a much slower rate. The reconstruction of Mertsch & Phan (2023) even exhibits a slight increase in surface mass density for r_gal ≥ 20 kpc, which may be a similar-but-opposite artefact stemming from a prior that enhances the surface mass density at large radii. Remarkably, the variance in the overall radial surface mass distribution between the individual samples is very small for most of the galactic disc. This is not necessarily surprising since summary statistics like radial profiles tend to average out differences between the individual samples.

The H₂ gas density on the other hand is much more strongly peaked towards the Galactic centre. It drops quickly and reaches a plateau between r_gal ≈ 2 kpc and r_gal ≈ 6 kpc. The bright spot at x ≈ 14.5 kpc towards the Galactic centre that we believe to be an unphysical reconstruction artefact leads to the narrow peak at r_gal ≈ 6.5 kpc. For larger radii, the H₂ density quickly and steadily falls off, slightly stronger than the other data-points but in overall qualitative agreement.

In Figs. 9 and 10, we show projections of our HI and H₂ maps along the Cartesian axes of Galactic coordinates. We show the average of the posterior samples μ, their standard deviation σ and the resulting signal-to-noise ratio μ/σ. These plots allow distinguishing which features in Fig. 6 stem from which gas component. The very dense concentration of gas in the Galactic centre region is mainly in the form of H₂. This is also true for the bright spot at x ≈ 14.5 kpc towards the Galactic centre mentioned earlier. In this direction, there is next to no direct distance information since the prior relative velocity is constant for all distances. The fact that this uncertainty is not correctly reflected in our samples requires further investigation. It is possible that this is a problem of convergence into a local minimum, where the algorithm fails to properly approximate the true posterior distribution. Apart from this feature, there is a paucity of HI gas for very large heliocentric radii (> 15 kpc).

The second obviously suspicious feature is the large gas accumulation at approximately (x = 20, y = 2.5) kpc in the HI maps. We suspect that the origin of this artefact lies in the radial profile functions (cf. Eq. (16)). The radius at which this large gas accumulation sits, coincides in Galacto-centric radius with the maximum of the inferred radial profiling function. This may be a hint that this way of profiling imposes a stronger prior on the gas distribution than intended.

In the anticentre-direction, we can see a break in the large-scale HI distribution. For negative y, there is a lot of gas at larger radii and smaller z than at positive y. This feature was discussed and analysed by Levine et al. (2006b), who argued, on the grounds of large-scale smoothness, that this is an unphysical artefact. It was shown that it can be corrected by assuming an inward motion of the gas towards the anticentre direction, for example due to elliptic gas orbits. A similar fix would be a change in the x-component of the Sun’s peculiar velocity (w.r.t the LSR), which the authors argue to be unphysical since the best-fit value, reconciling the radial surface mass densities, depends on radius indicating that this is a global, rather than a local effect.

Another property, which is highlighted in these plots, is that the uncertainties in the HI maps are generally much smaller than those of the H₂ maps. The main reason for this are the sizeable measurement uncertainties on the H₂ data. The error-bars of the HI dataset can be basically neglected for our purposes, as the resulting uncertainties are much smaller than those due to for example distance ambiguities. For the H₂ dataset, these model-induced uncertainties add to a much higher detector noise, incomplete sky-coverage, and uneven, sometimes very sparse sampling in spectral direction. Measurement surveys have mostly been conducted and published in angular and spatial regions where emission is present or expected. Regions of no-measurement are not treated as zero-measurements in our inference scheme, so they contribute to increased uncertainties. This effect is much more important for large, coherent structures that we can actually resolve (both in data and reconstruction). This leads to nearby structures being much more tightly constrained than those further away, because they occupy many lines of sight (cf. Fig. 9, centre-right). Since our reconstruction algorithm takes into account correlations in three dimensions, the information gain from neighbouring, correlated lines of sight proves to be very important.

In Fig. 9 one can also already see that the HI disc is significantly warped and flaring, with generally lower gas densities in the central regions of our galaxy. The H₂ disc is much more concentrated towards the inner portions of the Galaxy and mostly flat. In Figs. 11 and 12 we analyse this further and present the average height (z₀) and full-width-at-half-maximum (FWHM) of the z-profile of the reconstructed total gas density as a function of x and y. Regarding the Galactic warp, Fig. 11 reveals a significant displacement of the galactic disc from z = 0 beyond r_gal ≈ 10 kpc. This feature has been found in many previous studies (e.g. Burton & te Lintel Hekkert 1986; Nakanishi & Sofue 2016; Mertsch & Phan 2023) and is again confirmed here. The debate about its origin is not settled yet, with possible solutions being for example a recent encounter of the Milky Way with a satellite galaxy (Poggio et al. 2020) or a tilted, dynamical dark matter halo (Han et al. 2023).

Regarding the Galactic flaring, Fig. 12 shows that the FWHM increases with Galacto-centric radius in all directions. Given this approximate radial symmetry, we compute the azimuthally averaged FWHM and compare our results with previous studies in Figs. 13 and 14. For HI, we find a significant flaring of the disc up to r_gal ≈ 17 kpc. We find a slightly higher value for the FWHM than previous studies for regions inside the solar circle and a slightly lower value for regions outside. This may be partly due to the z-profile (cf. Eq. (17)) having no r-dependence, thereby slightly favouring gas realisations with a flatter FWHM-profile. For larger radii, our FWHM decreases again though there is almost no gas in our reconstruction at these large radii so this is likely not very meaningful. For H₂, we have significant amounts of gas only inside the solar circle, where our FWHM-profile is almost completely flat. This is consistent with previous reconstructions by Mertsch & Vittino (2021), but in conflict with all the other data points. Since there is a lot of scatter in the derived FWHM of previous works, we do not judge this to be a significant issue. In the Galactic centre, the FWHM reaches approximately 100 pc, which is comparable to the extent of a single voxel at these radii. This showcases a significant problem in our H₂ reconstruction, namely that it is strongly resolution-limited at these radii.

Fig. 6 Top-down view of the total reconstructed Hydrogen density Σtot with the spiral model of Reid et al. (2019) overlaid. The position of the Sun is at (0, 0). Positive x points towards the Galactic centre. A movie is available online.

Fig. 7 Azimuthally averaged surface mass density ΣHI of HI gas as a function of the Galacto-centric radius rgal. This work is indicated by a blue band showing the 1σ region of our individual samples. We compare our results with the results of earlier analyses by Wouterloot et al. (1990), Kalberla & Dedes (2008), Marasco et al. (2017), Mertsch & Phan (2023), and Nakanishi & Sofue (2016).

Fig. 8 Azimuthally averaged surface mass density ${\mathrm{\Sigma }}^{{\mathrm{H}}_2}$ of H2 gas as a function of the Galacto-centric radius rgal. This work is indicated by a red band showing the 1σ region of our individual samples. We compare our results with the results of earlier analyses by Bronfman et al. (1988), Grabelsky et al. (1987), Digel (1991), Mertsch & Vittino (2021), and Nakanishi & Sofue (2006).

Fig. 9 Projections of our 3D HI distribution along the x, y, and z axis. The left column shows the average surface mass density, the middle column shows the standard deviation thereof, and the right column shows the signal-to-noise ratio, i.e. the ratio of the first two columns. Note the different ranges on the colour-axes.

3.2 Reconstructed velocities

We present the azimuthally and z-averaged rotational velocity (i.e. the realised rotation curve of our reconstructed samples) in Fig. 15. Our reconstructed rotational velocity is enhanced over the prior for radii smaller than approximately 14 kpc, after which it briefly decreases before reverting to the prior. The maximum amplitude of this deviation is approximately 10 km s⁻¹ at r_gal ≈ 7 kpc. This increase in rotational velocity is thus favoured for essentially all regions that contain significant amounts of gas in our reconstruction. It is not clear how much this conflicts with rotation curves derived from stellar populations as they usually come with large, sometimes neglected systematic uncertainties. Additionally, they often probe the rotational velocity of stellar populations, which could slightly differ from gas velocities even on large scales. The aforementioned dip in rotational velocity at r_gal ≈ 16 kpc increases the relative velocity for gas clouds in outwards pointing directions to values normally appearing only at larger radii. Here, an interplay between the velocity curve and radial profile function occurs. Gas clouds are effectively moved inwards by the modified velocity curve (compared to where one would naively map them using a fixed velocity field), which is in agreement with the decreasing radial density profile at large radii.

For the peculiar motion of the Sun with respect to a corotating observer, we find U_⊙ = 1.590 ± 0.014 km s⁻¹, V_⊙ = 12.682 ± 0.029 km s⁻¹ and W_⊙ = 21.505 ± 0.047 km s⁻¹. The values for U_⊙ and W_⊙ differ significantly from the common values of the local standard of rest (${\mathit{v}}_{\odot }^{\mathrm{L}\mathrm{S}\mathrm{R}}$ ≈ (10.27, 15.32, 7.74) km s⁻¹). The deviation in U_⊙ is likely linked to the effect discussed in Levine et al. (2006b) and already mentioned above, that is, large-scale gas asymmetries (possibly sourced by deviations from circular gas orbits) being reduced by an outward motion of the LSR. Our value for V_⊙ is much closer to the literature value. Derivations thereof are less constraining and subject to systematic errors introduced by assuming a certain rotation curve for the Galaxy. We therefore deem this value to be consistent with common estimates. Our derived value for W_⊙ is probably dominated by local high-velocity-clouds that were not removed from the dataset (see e.g. Westmeier 2018). The inference algorithm attempts to fit this data by reconstructing extremely high velocities (magnitude ≥100 km s⁻¹) at high latitudes within a few hundred parsecs. Since the prior on deviations in the 3D velocity field strongly suppresses such solutions, a higher W_⊙ provides a middle ground. We started a reconstruction for a few iterations with fixed standard values for (U_⊙, V_⊙, W_⊙) and find that virtually nothing changes, as the impact of W_⊙ is mostly limited to local, high-latitude clouds and simply gets absorbed into the random field.

Depictions of how accurately the direct velocity data is reproduced by the inferred velocity field can be seen in Figs. 16 and 17, which show the mean of the reconstruction of the line-of-sight velocity together with the tracer objects. Akin to Fig. 5 in Tchernyshyov & Peek (2017), we display the data points with an inner segment showing the data value and an outer segment showing the reconstructed value at the tracer position. The background displays the velocity evaluated at a height $z_0^{\mathrm{\prime }}$, defined as the z-coordinate corresponding to maximum gas density at each x and y position, smoothed with a Gaussian filter of 50 pc width. This creates a continuous representation of the densest gas structures. The background velocity field shows a rich small-scale structure with many local features describing peculiar gas motions on scales of tens to hundreds of parsecs. Most of our tracer objects’ data values are reasonably well reproduced, but there are also many mismatches. The reduced χ² value for the radial velocities is 11.6, dominated by a few very large values. Removing the 5% (19 in absolute numbers) velocity tracers with the worst fit, this drops to 5.2, indicating a moderately good fit. The method by Tchernyshyov & Peek (2017) achieves an overall better data fit (4.0 unregularised, 1.3 regularised) using considerably less data and a distinct methodology. It is plausible that not all Masers and YSOs serve as accurate tracers of the bulk gas flow at their respective locations, as they may either capture motion on scales smaller than our resolution or exhibit (partial) decoupling from the large-scale dynamics. Additionally, deviations may arise from the assumed prior correlation structure for both the gas densities and velocity fields, which imposes a smoothness scale that may be inconsistent with some of our tracer objects. We expect that a more flexible velocity model, especially on large (kpc) scales would improve this fit.

Fig. 10 Projections of our 3D H2 distribution along the x, y, and z axis. Left column shows the average surface mass density, middle column shows the standard deviation thereof, and right column shows the signal-to-noise ratio, i.e. the ratio of the first two columns. Note the different ranges on the colour-axes.

Fig. 11 Average height $z_0=⟨z{⟩}_{n^{\text{tot\hspace{0.17em}}}}$ of the galactic gas densities in our reconstruction as a function of x and y. Overlaid is the spiral pattern fitted by Reid et al. (2019).

Fig. 12 FWHM of the z profile of the galactic gas densities in our reconstruction as a function of x and y. Overlaid is the spiral pattern fitted by Reid et al. (2019).

Fig. 13 Azimuthally averaged FWHM of the mean reconstructed HI gas density. We compare our results with previous works by Wouterloot et al. (1990), Kalberla & Dedes (2008), Marasco et al. (2017), Mertsch & Phan (2023), and Nakanishi & Sofue (2016).

Fig. 14 Azimuthally averaged FWHM of the mean reconstructed H2 gas density. We compare our results with previous works by Digel (1991), Malhotra (1994), Nakanishi & Sofue (2006), Mertsch & Vittino (2021), Bronfman et al. (1988), and Grabelsky et al. (1987).

Fig. 15 Inferred rotational velocity curve of the galactic gas as a function of rgal. The green line shows the prior mean and the black dots show the data used to create it (cf. Sect. 2.1.3).

Fig. 16 Mean LSR velocity at height $z=z_0^{\mathrm{\prime }}$, defined as the z coordinate corresponding to the maximum gas density at each x and y position, smoothed with a Gaussian filter with a 50 pc width. The markers are akin to Fig. 5 but consist of two parts: The inner part shows the measurement value and the outer part shows the mean value of the reconstructed velocity fields at the position of the tracer object. Markers with a thick grey edge indicate tracers whose radial velocity component is badly reproduced by the reconstructed velocity field.

Fig. 17 Similar to Fig. 16, but limited to within 1.25 kpc.

3.3 Comparison to previous 3D gas reconstructions

We can also directly compare our reconstructed gas maps with previous reconstructions of HI in Mertsch & Phan (2023) and H₂ in Mertsch & Vittino (2021). We show side-to-side comparisons of the projected total surface mass density in Fig. 18. One obvious difference is the smaller radial extent of our reconstruction, courtesy of the radial profiling function as already discussed above. Another striking feature is that structures in our reconstruction tend to be more smoothed out. This is for two related reasons: First, the variable velocity field in our model results in infinitely many solutions for gas placement along a single line of sight (but they are not equally likely), while in the reconstructions of Mertsch & Vittino (2021) and Mertsch & Phan (2023) a fixed velocity model was assumed that only allows for a few distinct allowed distances along each line of sight. This leads to every individual sample in our reconstruction having a generally more diffuse gas distribution. Secondly, also due to our flexible velocity model, the placement of gas differs more strongly from sample to sample. This leads to much more diverse individual samples where whole clouds and cloud complexes shift in distance, thereby producing even more smoothness in the radial direction when an average is taken over all samples. Also, the different (prior) velocity profile leads to a different arrangement of cloud positions, which has the most impact in regions of large velocity gradients, that is towards the inner Galaxy.

In Fig. 19, we present the sample-average integrated sky projection of the total Hydrogen number column density, again in comparison to Mertsch & Vittino (2021) and Mertsch & Phan (2023). The older reconstructions show a strong discontinuity at l = 0 (central vertical line in the sky projection). Also, some individual structures are very faint or entirely missing. This is mainly due to the assumed constant velocity model offering a finite amount of possible positions towards these lines of sight where gas can be placed with the correct relative velocity. In this work however, the variable velocity field can in principle accommodate all velocities, automatically striking a balance between deviating from the prior and forming coherent structures at sensible distances in a Bayesian sense. This leads to a much better agreement of our reconstruction with the measurement data at the cost of introducing more degrees of freedom to our model.

A key feature of our HEALPix × log(r)-grid is the very high resolution nearby, which allows us to resolve local structures with high fidelity. This is of high interest, since nearby structures cover large angular scales. Therefore, these local structures will benefit the most from our modelling using correlated 3D fields. We illustrate the difference, again comparing to Mertsch & Vittino (2021) and Mertsch & Phan (2023), on a regular grid in Fig. 20, showing the local total surface mass density. They are very different morphologically, and it seems almost impossible to match structures between these two, partly also because of the stark difference in local resolution. In our map, structures such as Vela ((x, y) ≈ (0, 0.8) kpc or (Serpens-)Aquila Rift ((x, y) ≈ (0.4, 0.2) kpc) are clearly identifiable. Other structures, such as the Radcliffe wave (Alves et al. 2020) are partly visible, and many other structures are misplaced with respect to estimates derived from stellar distances (see next section).

3.4 Comparison to local 3D interstellar dust reconstructions

Locally, the quality of our reconstruction can be illustrated by comparing it to the distribution of other, better-known local quantities. One such quantity that is believed to be strongly correlated to gas densities is interstellar dust. Comparing the amount of reddening of stars along a line of sight allows inferring the amount of dust in-between. Doing this for the full sky yields a 3D map of local interstellar dust that is comparatively accurate since distances to stars are well constrained by parallax measurements.

We show a comparison with the dust reconstruction by Edenhofer et al. (2024b) in Fig. 21. It can immediately be seen that the overall morphology of the gas and dust distributions bear some similarity. The valley-shaped structure tilted to the upper right of Fig. 21 and some traces of a local bubble can be recognised. The dense cloud at the bottom of the figure, Vela, is also correctly placed. A couple other features are not reproduced but appear to be displaced in the radial direction. Given that we have only included sparse and implicit distance information with often large uncertainties, we can not hope for a much more accurate result than this. This comparison highlights that there lies great potential for future reconstructions of the ISM to use dust information for fixing the distances to gaseous clouds in the ISM. Previous works, for instance by Tchernyshyov & Peek (2017), Zucker et al. (2021), or Soler et al. (2023), have already successfully exploited this relation.

In Fig. 22, we aim to analyse the correlation between dust and gas further. To this end, we compute the dust extinction A_V and the total Hydrogen column density N_H up to every grid point within 1250 pc. Plotting these in a logarithmic histogram, we can visualise the nearly linear correlation between dust and gas. The large spread around a clear central line indicates that there are many gas clouds in our reconstruction that are radially displaced with respect to dust clouds, but they are in good agreement overall. Since we are comparing integrated quantities, this local misplacement does not ruin the linear correlation. We find the best fit value of $N_{\mathrm{H}}/A_{\mathrm{V}}\approx {2.37}_{-1.09}^{+2.01}\times {10}^{21}{\mathrm{c}\mathrm{m}}^{-2}{\mathrm{m}\mathrm{a}\mathrm{g}}^{-1}$. This value obtained for the dust-to-gas ratio is compatible with commonly found estimates (Bohlin et al. 1978; Rachford et al. 2009).

To estimate at what spatial scales our reconstruction agrees with the dust extinction maps, we integrate the dust and reconstructed gas densities up to some variable upper integration bound. Then, we compute the Pearson and Spearman correlation coefficients as a function of this upper integration bound. The Pearson correlation coefficient tests for a linear relationship between the two quantities, while the Spearman correlation coefficient tests for a monotonous relationship. We also compute weighted versions of the correlation coefficient taking into account our uncertainty estimate of the reconstruction. The result is shown in Fig. 23. We find that for radii below approximately 400 pc, both correlation coefficients are comparatively small. For larger radii, especially the Spearman coefficient reaches very large values, indicating a very good monotonous relationship between integrated dust and gas densities. The Pearson coefficient is a bit smaller, indicating an imperfect linear relationship. We suspect that this is mainly due to misplaced gas which completely overshadow potential non-linearities in the true relationship between these quantities. The uncertainty-weighted versions of the correlation coefficients result in a worse agreement for small r and a smoother growth of the correlation coefficients for large r. This may again indicate that gas clouds (or the absence thereof) are slightly misplaced with an overly confident error estimate locally. On the other hand, the dust reconstructions should not be understood as the absolute truth. For one, the dust map is probably not devoid of artefacts and can suffer from the effects of incomplete data in some regions. The other reason is that the assumption that gas and dust should be connected by the same linear relationship globally rests on a shaky ground. Typical deviations from this could be of the order of 10–20 per cent (Tricco et al. 2017), sometimes even larger. We conclude that our gas maps should prove useful for applications that require accuracy only on scales of 100s of parsecs.

Fig. 18 Side-by-side comparison of mean surface mass density Σtot. Left: the ‘BEG03’ model from Mertsch & Vittino (2021) and Mertsch & Phan (2023). Right: this work.

Fig. 19 Side-by-side comparison of mean total Hydrogen number column density NH. Left: the ‘BEG03’ model from Mertsch & Vittino (2021) and Mertsch & Phan (2023). Right: this work.

Fig. 20 Side-by-side comparison of local mean surface mass density Σtot. Left: the ‘BEG03’ model from Mertsch & Vittino (2021) and Mertsch & Phan (2023). Right: this work.

Fig. 21 Local total Hydrogen density Σtot in comparison with the projected dust extinction by Edenhofer et al. (2024b). The coloured background shows the dust distribution, with no values available in the outer black regions. The white contours show the distribution of total gas density.

4 Summary and prospects

We have presented new 3D reconstructions, mainly of the interstellar HI and H₂ (traced by CO) distributions, and their velocity field. These reconstructions are based on deprojecting gas line surveys by the HI4PI Collaboration (2016) as well as Dame et al. (2001) and Dame & Thaddeus (2022) with a variable Galactic velocity field and emission line widths. The underlying model is based on Gaussian processes with a fixed, isotropic, and homogeneous correlation structure. We included modelling of the effects of absorption under the assumption of constant excitation temperatures for HI and CO. We used a Bayesian variational inference method by Frank et al. (2021) to obtain a set of eight reconstruction samples (cf. Eq. (28)). This allowed us to implicitly capture the uncertainties of the reconstructed quantities and distinguish between real structures and noise artefacts. Additionally, our numerical discretisation scheme features an increased resolution near our position in the Galaxy. All data products are available via zenodo⁹.

While our gas maps are spatially coherent and successfully capture noise uncertainties, they are not free from artefacts and model insufficiencies. We showed that our reconstructions are particularly vulnerable to such problems for heliocentric radii larger than 15 kpc, that is, on the far side of the Galaxy. On scales larger than approximately 400 pc, we find a good agreement with local maps of interstellar dust. We investigated the large-scale radial structure, as well as warping and flaring of the Galactic disc, and we confirmed previous findings in these matters. We also present reconstructions of the 3D Galactic velocity field corresponding to our gas distributions. Due to the sparse data availability of objects tracing the Galactic velocity field, this naturally comes with large uncertainties. This leads us to believe that our gas reconstructions can be highly useful for applications that require a good data fit and can tolerate gas misplacements on these scales.

In the future, we plan to improve on this in multiple ways: First, our assumed correlation structure for the Gaussian processes is a major source of model uncertainty. The resulting reconstructions depend on the particular choice of the correlation kernels and since not all of them can be measured directly, we should treat them as free parameters with suitable priors in our reconstruction. It is unclear to what degree estimations from nearby Galaxies and clouds can be applied to the whole of the Milky Way. Studies of the Small and Large Magellanic Clouds suggest gas power spectrum slopes roughly between −3.0 and −4.0, depending on the position and scale (Szotkowski et al. 2019). Here, we kept our correlation kernels fixed due to numerical limitations of the used algorithm. We plan to alleviate this restriction and reconstruct the correlation kernels as well.

Second, as Fig. 21 shows, we can learn a lot about the 3D distribution of interstellar gas by comparison to dust. Treating external 3D dust maps as data in our reconstruction, we can, not only improve our gas distributions locally, but also far away as this should help us with distinguishing local and far structures. Additionally, we would gain valuable insight into the structure of our velocity field locally. By fixing the correlation structure of various fields with local data, this knowledge would then also transfer to regions outside the availability of 3D dust maps. One could also consider a joint reconstruction of gas (this work) and dust (e.g. Edenhofer et al. 2024b) in a common inference scheme. However, these two models are quite large by themselves, and combining them would be approximately twice as expensive computationally. Additionally, while there is a lot of information about the distances to gas clouds contained within the stellar reddening data, we do not expect the dust reconstruction to profit much from the information contained within the gas spectra.

Third, the assumed radial profiling function introduces significant issues and should be replaced with an improved large-scale prior on gas density. This adjustment should help suppress artefacts at large radii – particularly towards the Galactic centre and anticentre – while still yielding a physically realistic model of the Galaxy.

Fourth, our current velocity model lacks sufficient flexibility on small scales to account for clouds with velocities that deviate strongly from the prior rotation curve. On large scales, it also fails to capture potential deviations from purely circular motion. An even more flexible velocity model on all scales is expected to cure some problems of our current approach.

Finally, we anticipate that increasing the resolution of our reconstruction would improve our results by increasing the available information contained within correlations between neighbouring lines of sight. Moreover, a higher resolution allows for a higher fidelity of the reconstruction and facilitates studying smaller structures in the ISM.

To conclude, despite various shortcomings, our reconstructions present a clear improvement over previous works. We anticipate that they are useful for estimating uncertainties in modelling, for example, diffuse gamma-ray emission from cosmic-ray interactions with interstellar gas. They also present yet another step forward in unravelling the (spiral) structure of the Milky Way.

Fig. 22 2D logarithmic histogram of grid points within 1250 pc comparing the dust reconstruction by Edenhofer et al. (2024b) with this work. For every grid point, we compute the dust extinction AV and the total Hydrogen column density NH up to and including this point. We also show the result of a fit with a linear function (without constant term) and 1σRMS (root-mean-square error) neighbouring lines.

Fig. 23 Pearson and Spearman correlation coefficients between our reconstructed gas density and the dust densities reconstructed in Edenhofer et al. (2024b) as a function of the upper radial integration bound. Weighted versions of the correlation coefficients, taking into account our uncertainty estimate of the reconstruction, are also shown.

Movie

Movie 1 associated with Fig. 6 Access here

Data availability

The movie associated to Fig. 6 is available at https://www.aanda.org

Acknowledgements

Some of the results in this paper have been derived using the healpy and HEALPix package. Some figures in this publication have been created using matplotlib (Hunter 2007). Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number 495252601. This research was funded in part, by the Austrian Science Fund (FWF) [I 5925-N]. The authors gratefully acknowledge the computing time provided to them at the NHR Center NHR4CES at RWTH Aachen University (project number p0020479). This is funded by the Federal Ministry of Education and Research, and the state governments participating on the basis of the resolutions of the GWK for national high performance computing at universities (https://www.nhr-verein.de/unsere-partner). Gordian Edenhofer acknowledges that support for this work was provided by the German Academic Scholarship Foundation in the form of a PhD scholarship (‘Promotionsstipendium der Studienstiftung des Deutschen Volkes’). Philipp Frank acknowledges funding through the German Federal Ministry of Education and Research for the project ErUM-IFT: Informationsfeldtheorie für Experimente an Großforschungsanlagen (Förderkennzeichen: 05D23EO1). V.H.M.P. acknowledges support from the Initiative Physique des Infinis (IPI), a research training program of the Idex SUPER at Sorbonne Université.

References

All Tables

All Figures

	Fig. 1 Nside = 128 projection of the HI4PI data on the two-sphere in a Mollweide projection on a logarithmic scale.
In the text

	Fig. 2 Longitude-velocity diagram of the HI4PI data (integrated over latitude) corresponding to Fig. 1.
In the text

	Fig. 3 Nside = 128 projection of the CO data on the two-sphere in a Mollweide projection on a logarithmic scale. Some notable structures are indicated by arrows.
In the text

	Fig. 4 Longitude-velocity diagram of the CO data (integrated over latitude) corresponding to Fig. 3.
In the text

Fig. 5 Positions of masers (circles) and YSO clusters (diamonds) shown in a 2D projection of the Galactic disc. The marker colour shows the corresponding vLSR measurement. The black streaks indicate the 1σ uncertainties on the parallax (distance) measurements. The Sun is located at (0, 0) and the centre of the Galaxy at (8.178, 0). The background shows vLSR computed from an interpolated rotation curve as described in Sect. 2.1.3. The inset shows the underlying rotation curve vcirc,0 and the data used for its creation (cf. Sect. 2.1.3).

In the text

	Fig. 6 Top-down view of the total reconstructed Hydrogen density Σtot with the spiral model of Reid et al. (2019) overlaid. The position of the Sun is at (0, 0). Positive x points towards the Galactic centre. A movie is available online.
In the text

	Fig. 7 Azimuthally averaged surface mass density ΣHI of HI gas as a function of the Galacto-centric radius rgal. This work is indicated by a blue band showing the 1σ region of our individual samples. We compare our results with the results of earlier analyses by Wouterloot et al. (1990), Kalberla & Dedes (2008), Marasco et al. (2017), Mertsch & Phan (2023), and Nakanishi & Sofue (2016).
In the text

Fig. 8 Azimuthally averaged surface mass density ${\mathrm{\Sigma }}^{{\mathrm{H}}_2}$ of H2 gas as a function of the Galacto-centric radius rgal. This work is indicated by a red band showing the 1σ region of our individual samples. We compare our results with the results of earlier analyses by Bronfman et al. (1988), Grabelsky et al. (1987), Digel (1991), Mertsch & Vittino (2021), and Nakanishi & Sofue (2006).

In the text

	Fig. 9 Projections of our 3D HI distribution along the x, y, and z axis. The left column shows the average surface mass density, the middle column shows the standard deviation thereof, and the right column shows the signal-to-noise ratio, i.e. the ratio of the first two columns. Note the different ranges on the colour-axes.
In the text

	Fig. 10 Projections of our 3D H2 distribution along the x, y, and z axis. Left column shows the average surface mass density, middle column shows the standard deviation thereof, and right column shows the signal-to-noise ratio, i.e. the ratio of the first two columns. Note the different ranges on the colour-axes.
In the text

	Fig. 11 Average height $z_0=⟨z{⟩}_{n^{\text{tot\hspace{0.17em}}}}$ of the galactic gas densities in our reconstruction as a function of x and y. Overlaid is the spiral pattern fitted by Reid et al. (2019).
In the text

	Fig. 12 FWHM of the z profile of the galactic gas densities in our reconstruction as a function of x and y. Overlaid is the spiral pattern fitted by Reid et al. (2019).
In the text

	Fig. 13 Azimuthally averaged FWHM of the mean reconstructed HI gas density. We compare our results with previous works by Wouterloot et al. (1990), Kalberla & Dedes (2008), Marasco et al. (2017), Mertsch & Phan (2023), and Nakanishi & Sofue (2016).
In the text

	Fig. 14 Azimuthally averaged FWHM of the mean reconstructed H2 gas density. We compare our results with previous works by Digel (1991), Malhotra (1994), Nakanishi & Sofue (2006), Mertsch & Vittino (2021), Bronfman et al. (1988), and Grabelsky et al. (1987).
In the text

	Fig. 15 Inferred rotational velocity curve of the galactic gas as a function of rgal. The green line shows the prior mean and the black dots show the data used to create it (cf. Sect. 2.1.3).
In the text

Fig. 16 Mean LSR velocity at height $z=z_0^{\mathrm{\prime }}$, defined as the z coordinate corresponding to the maximum gas density at each x and y position, smoothed with a Gaussian filter with a 50 pc width. The markers are akin to Fig. 5 but consist of two parts: The inner part shows the measurement value and the outer part shows the mean value of the reconstructed velocity fields at the position of the tracer object. Markers with a thick grey edge indicate tracers whose radial velocity component is badly reproduced by the reconstructed velocity field.

In the text

	Fig. 17 Similar to Fig. 16, but limited to within 1.25 kpc.
In the text

	Fig. 18 Side-by-side comparison of mean surface mass density Σtot. Left: the ‘BEG03’ model from Mertsch & Vittino (2021) and Mertsch & Phan (2023). Right: this work.
In the text

	Fig. 19 Side-by-side comparison of mean total Hydrogen number column density NH. Left: the ‘BEG03’ model from Mertsch & Vittino (2021) and Mertsch & Phan (2023). Right: this work.
In the text

	Fig. 20 Side-by-side comparison of local mean surface mass density Σtot. Left: the ‘BEG03’ model from Mertsch & Vittino (2021) and Mertsch & Phan (2023). Right: this work.
In the text

	Fig. 21 Local total Hydrogen density Σtot in comparison with the projected dust extinction by Edenhofer et al. (2024b). The coloured background shows the dust distribution, with no values available in the outer black regions. The white contours show the distribution of total gas density.
In the text

Fig. 22 2D logarithmic histogram of grid points within 1250 pc comparing the dust reconstruction by Edenhofer et al. (2024b) with this work. For every grid point, we compute the dust extinction AV and the total Hydrogen column density NH up to and including this point. We also show the result of a fit with a linear function (without constant term) and 1σRMS (root-mean-square error) neighbouring lines.

In the text

	Fig. 23 Pearson and Spearman correlation coefficients between our reconstructed gas density and the dust densities reconstructed in Edenhofer et al. (2024b) as a function of the upper radial integration bound. Weighted versions of the correlation coefficients, taking into account our uncertainty estimate of the reconstruction, are also shown.
In the text