© ESO, 2017

1. Introduction

For over 60 years, the properties of the neutral gas in the Galaxy have been studied extensively by several authors. The pioneering work by Westerhout (1957) allowed imaging of the spatial distribution of the H I in the Milky Way for the first time. In Westerhout’s work, the transformation from the observed space to the physical space was achieved by assuming an a-priori model for the Galactic velocity field, that is, for the rotation curve. Later on, a similar analysis was repeated several times, with relatively minor modifications to the original procedure, by using the most up-to-date observations for both the atomic (Kulkarni et al. 1982; Burton & te Lintel Hekkert 1986; Nakanishi & Sofue 2003; Kalberla & Dedes 2008) and the molecular hydrogen (Bronfman et al. 1988; Nakanishi & Sofue 2006). These works revealed that the gas distribution in the Milky Way is relatively complex: the H I disc is flared and warped (Levine et al. 2006b) and a multi-arm spiral structure is visible in both the atomic and the molecular hydrogen (e.g. Dame et al. 1986; Levine et al. 2006a). Insights into the dark matter distribution in the disc can be obtained when the mapping is coupled to a self-consistent mass model (Kalberla et al. 2007).

Despite its many successes, this approach has a number of shortcomings. Perhaps the most severe of them was pointed out by Burton (1974), who noticed that the derived gas density distribution is dramatically sensitive to small changes in the assumed velocity field, as velocity perturbations by only a few km s^-1 can mimic variations in density up to a factor of 100. An additional complication is introduced by the so-called “near-far problem”, which occurs in the Galactic region inside the solar circle (hereafter, the “inner” Galaxy) where the mapping is not unique: for each sight-line, two distinct distances – one in front and one beyond the tangent point – can be associated to the same line-of-sight velocity. This produces ambiguity in the reconstruction of the density field, and even though clever solutions have been proposed (e.g. Nakanishi & Sofue 2003), the detailed gas distribution in the inner Galaxy remains uncertain. Finally, not all these studies show a direct comparison between the synthetic line-profiles – produced by the derived gas density distribution moving in the modelled velocity field – and the data, which is key to assessing the goodness of the model.

Since any prejudice on the rotation curve can significantly affect the derived gas distribution, it is mandatory to determine the former independently from the latter. The rotation curve in the inner Galaxy is classically derived by using the so-called tangent point method (Kwee et al. 1954). This method assumes that the gas has pure circular motions, and it is based on extracting the terminal-velocity from the line profiles and assigning it to gas located at the tangent-point radius R_tp = R_⊙sin(l), where the near-far ambiguity vanishes. However, the identification of the terminal-velocity in the profiles is not straightforward, because gradients in gas density and velocity dispersion along the line-of-sight have a large impact on the tail of the profiles. Hence, the terminal velocity is often derived by assuming that the velocity dispersion and the density of the gas are constant for a wide range of line-of-sight velocities (Shane & Bieger-Smith 1966; Celnik et al. 1979; Rohlfs & Kreitschmann 1987; Malhotra 1995). This leads to a certain degree of circularity in the method, as an assumption on the gas density profile must be made in order to derive the terminal velocities, and so the rotation curve, which thereafter can be used to image the spatial distribution of the gas itself.

In this paper, we present a novel approach to deriving the kinematics and the distribution of the atomic and molecular hydrogen in the inner Galaxy via the modelling of the observed emission-line profiles. Our modelling is based on the following assumptions: the Galactic disc is axisymmetric and the gas moves on co-planar and concentric circular orbits, without coherent radial or vertical motions. While these assumptions may appear crude, they allow us to simultaneously fit both the gas kinematics (rotation velocity and velocity dispersion) and the gas distribution (midplane density and scale height) of our model to the data. This approach circumvents most of the aforementioned issues, and provides a simple model of the inner Galactic disc that is in remarkable agreement with the observations. We note that the modelling of the gas kinematics in nearby galaxies is classically done with the so-called “tilted ring” technique (e.g. Rogstad et al. 1974; Begeman 1987), which is based on the same assumption adopted here¹.

This paper is structured as follows. In Sect. 2 we describe, in detail, our models and present the method adopted to fit them to the data. In Sects. 3 and 4 we show our results on the kinematics and distribution of the H I and H₂ as derived by our modelling. We discuss our results in Sect. 5 and present our conclusions in Sect. 6.

In this work, we use two coordinate systems. The first is the Galactic coordinate system (l,b,v), where l and b are the Galactic longitude and latitude and v is the line-of-sight velocity in the local standard of rest. The second is a cylindrical coordinate system (R,φ,z), where R and z are the distances perpendicular to and along the rotation axis of the Galaxy and φ is the azimuthal coordinate, assumed positive in the direction of the Galaxy rotation. We use the Galactic constants R_⊙ = 8.3 kpc and v_⊙ = 240 km s^-1 (McMillan 2011). In Appendix C we show how our results change if we assume a different value for v_⊙.

2. Method

To model the atomic and the molecular hydrogen components of the Milky Way, we make use of the Leiden-Argentine-Bonn (LAB) all-sky 21-cm survey (Kalberla et al. 2005) and of the ¹²CO (J = 1−0 rotational transition line) survey of Dame et al. (2001). We use a hanning-smoothed version of the LAB data, with an angular (velocity) resolution of ~ 0.6° (~ 2 km s^-1) and voxel size of 0.5° × 0.5° × 2 km s^-1. The CO data of Dame et al. (2001) have instead angular (velocity) resolution of ~ 0.15° (~ 1.3 km s^-1) and voxel size of 0.125° × 0.125° × 1.3 km s^-1. The different angular and spectral resolutions may complicate the interpretation of our results, thus we prefer to bring both datasets to a common resolution. In order to increase the signal-to-noise ratio and reduce the number of points to fit, we smooth both datacubes to an angular resolution² of 1° × 1°. The CO datacube was also smoothed to a velocity resolution of 2 km s^-1, and re-binned to match the voxel size of the LAB data. For the remainder of the paper we use these “degraded” versions of the datacubes.

We focus on the regions where the emission lines are produced by the inner (R<R_⊙) Galaxy: 0°<l< 90° and v> 0, or first quadrant (QI), and at 270°<l< 360° and v< 0, or fourth quadrant (QIV). Because of the Galaxy rotation, QI appears to be receding from an observer placed at the location of the Sun, whereas QIV appears to be approaching. We model QI and QIV separately.

We assume that the Galactic disc can be decomposed into a series of concentric and co-planar rings, centred at the Galactic centre. Following Olling (1995), the vertical distribution of the gas is approximated as a Gaussian. The gas in each ring is described by four parameters: the rotation velocity v_φ, the velocity dispersion σ, the midplane volume density n₀ and the scale height h_s, which is the half width at half maximum (HWHM) of the Gaussian distribution (so that gives the gas density n as a function of z for a given ring). We assume that the gas kinematics does not change with z (but see Sect. 5.7). These four parameters are then fit to the data: our fitting strategy consists of building a synthetic observation of our model (Sect. 2.3) and comparing it with the data (Sect. 2.4).

We model the inner disc by using 180 rings for each Galactic quadrant, each ring is at distance R_⊙ | sin(l) | from the Galactic centre: this gives us two points per resolution element per quadrant. We do not model the warp in the Galactic disc as its amplitude is negligible in the regions considered (Levine et al. 2006b).

In the following, we describe how we build synthetic datacubes for the H I emission of our model, and how we compare them to the observations. An equivalent procedure is used for the molecular gas, with a few differences discussed in Sect. 4.

2.1. Modelling the gas on the midplane

We first focus on the H I located on the midplane (b = 0°) and fit the parameters v_φ, σ and n₀ of each ring to the LAB data.

Our iterative method, similar to that adopted by Kregel & van der Kruit (2004) to study the H I kinematics in nearby edge-on galaxies, is sketched in Fig. 1. We consider an observer located on the midplane at R = R_⊙ and we are interested in modelling the H I brightness temperature profile T_B(v) at a given Galactic longitude l and at latitude b = 0°. This profile is built from the contribution of all rings with 0 ≤ i ≤ i_max along the line-of-sight, where R₀ ≡ R_⊙ and R_imax ≡ R_⊙sin(l) = R_tp. We note that rings at radii smaller than R_tp do not contribute to this profile. For each new l we fix v_φ,i, σ_i and n_0,i of all rings with i<i_max to the values determined by the previous iterations, we evaluate a synthetic brightness temperature profile T_B(v) (see Sect. 2.3) and fit it to the observed profile by varying only the parameters of the ring at i = i_max. Details of the fitting procedure are given in Sect. 2.4. We apply this method to determine the values of the parameters for all rings recursively, starting from R_⊙ (l = 90° for QI, or l = 270° for QIV) down to the Galactic centre (l = 0° or 360°).

The parameters of the first ring in the iteration cannot be determined and must be assumed: we fix them to v_φ,0 = v_⊙, σ₀ = 8 km s^-1 and n_0,0 = 0.5 cm^-3. We tested that a factor of two variation in n_0,0 or σ₀ has virtually no impact on our results. This happens because the ring number density is the largest around R_⊙, thus in the solar neighbourhood the model parameters stabilise quickly regardless of the values chosen for the outermost ring.

Our approach differs substantially from the classical tangent point method, where a) each line profile is treated independently, and b) the emission at the terminal-velocity is always associated to gas at the tangent point. In fact, under the assumptions of axisymmetry and circular rotation, (a) and (b) are redundant. Our iterative fitting procedure intrinsically accounts for the fact that line profiles at different longitudes are not independent, but instead originate from a common underlying disc. Also, for any sight-line we do not force the gas at the tangent point radius (i.e. at ring i = i_max) to contribute to the highest-velocity region of the line profile, because the fit to the innermost ring is unconstrained. We argue that, within the framework of axisymmetry, ours is the best possible approach for modelling the data without making any a-priori assumption on the trends of the various parameters.

2.2. Modelling the gas above the midplane

Once the parameters for the midplane are determined, we proceed to compute the scale height of the gas in each ring. We focus on gas at b ≠ 0 and adopt the same iterative method of Sect. 2.1. The difference here is that, instead of modelling a single line-profile for each longitude l, we model the latitude-velocity (b-v) slice T_B(b,v) and we fit it to the LAB data by adjusting the scale height of the tangent-point ring.

The main caveat of this procedure lies in the assumption that the gas is in cylindrical rotation, that is, the rotation velocity of the gas in each ring does not change with z. However, Marasco & Fraternali (2011) showed that a layer of slow-rotating extra-planar H I is present in the Milky Way. Therefore, we decide to confine our fit to latitudes | b | ≤ 5°, which allows us to minimise contaminations from this slow-rotating component while maintaining enough data points to study the vertical distribution of the gas at most radii (see Sect. 5.7 for further details). Profiles at | b | > 5° are relevant only for modelling regions close to R_⊙: for instance, for a disc half-thickness of 100 pc, a sight-line of b = 6° would already “pierce” through the disc at R = 7.3 kpc and the line profile would contain little information on the gas located at smaller radii. As contamination from extra-planar H I should be significant close to R_⊙, we prefer not to model the high-latitude profiles in the data. We refer the reader to Marasco & Fraternali (2011) for a model of the Galactic H I emission at high latitudes.

2.3. Simulating a line profile

In our model, the brightness temperature profile T_B(v) in a given direction of the sky (l,b) is determined in two steps. The first consists in deriving the H I column density profile N_HI(v). Gas at a generic distance d with respect to an observer placed at the Sun location (see Fig. 1) has the following cylindrical coordinates: and its line-of-sight velocity is given by (1)We note that Eq. (1) is valid for circular motions only. This gas will contribute to the total column density in (l,b) by an amount ΔN_HI = n(R,z) Δd, being Δd a generic step along the line-of-sight. We assume that this column density is spread in the velocity domain by following a Gaussian distribution with mean equal to v(d) and standard deviation equal to σ(d), which is the velocity dispersion of the gas at radius R(d). In practice, we consider a series of steps along the line of sight 0 <d_i<d_max, where d_max is the distance at which the line of sight intercepts the solar circle (see Fig. 1), and evaluate the H I column density profile as (2)where v, σ , and n are derived by linear interpolation of the various rings and Δv is the channel separation in the datacube (2 km s^-1). We assume that Δd_i is either 10 pc or the distance along the line-of-sight between two consecutive rings, whichever is the smallest.

As a second step, we must convert N_HI(v) to a brightness temperature T_B(v). This is a delicate step, which we begin to tackle by first reiterating some basic equations for the H I radiative transfer. In the absence of background sources, the solution to the equation of radiative transfer for an extended (i.e. filling entirely the telescope beam), homogeneous and isothermal layer of H I at a given line-of-sight velocity v is (3)where τ(v) is the H I optical depth and T_S is the H I excitation (or spin) temperature. Following Binney & Merrifield (1998, page 473), the H I column density between two generic line-of-sight velocities v₁ and v₂ can be expressed in terms of the gas spin temperature and the optical depth as (4)where c = 1.823 × 10¹⁸ cm^-2 K^-1 km^-1 s. If we differentiate Eq. (4) we can easily re-arrange for the optical depth: (5)where the quantity N_HI(v,v + Δv) / Δv has units of cm^-2 km^-1 s and is in practice given by Eq. (2). In the optical thin limit τ ≪ 1 (i.e. low column density and/or high spin temperature), Eq. (3) gives T_B ≃ T_Sτ = N_HI/c, so the measured brightness temperature is proportional to the intrinsic column density per unit velocity. In the optical thick limit τ ≫ 1, instead, T_B ≃ T_S. Thus Eq. (3) implies that the brightness temperature cannot exceed the gas spin temperature.

The interstellar medium (ISM), however, is not isothermal and homogenous, but instead consists of an ensemble of clouds, shells, and filaments with different sizes, temperatures, and densities, roughly in pressure equilibrium (e.g. Wolfire et al. 2003). Radiative transfer in such a complex medium cannot readily be described by Eqs. (3) and (5). In general, there is no unique relation between column density and brightness temperature: given any N_HI(v), the observed brightness temperature profile will depend on the physical properties of the clouds that make up the ISM, on their distribution along the line of sight, and on the observing setup (e.g. Dickey & Lockman 1990, see also Appendix A.1 of this work). However, thanks to the many theoretical (McKee & Ostriker 1977; Wolfire et al. 2003) and observational (Gibson et al. 2000; Heiles & Troland 2003; Dickey et al. 2003; Strasser & Taylor 2004) studies on the Galactic ISM carried out in the last decades, we have formed a good understanding of the physical properties of the neutral gas. Using results from these studies, in Appendix A we build a model for the multiphase Galactic ISM and compute the exact N_HI(v)−T_B relation for an observing setup analogous to that of the LAB survey. We show that such a relation can still be approximately described by Eqs. (3) and (5), using an “effective” spin temperature of ~ 152 K, the peak brightness temperature measured in the LAB survey. This value is two to three times larger than the typical T_S determined in studies of H I absorptions against background continuum sources. For instance, Dickey et al. (2003) and Heiles & Troland (2003) found median T_S of 65 K and 50 K, respectively. These studies, however, attempt to characterise only the coldest and densest phase of the H I, which is the only phase responsible for the absorptions but contributes only in part to the overall 21-cm emission. Our aim is to have the simplest possible prescription for the conversion between N_HI(v) and T_B(v) that best describes the 21-cm LAB data, thus our effective spin temperature of 152 K is a better motivated choice and we adopt it as our fiducial T_S for the purpose of computing Eqs. (3) and (5). We caution the reader that this value must be interpreted neither as the typical T_S of the Galactic H I nor as representative for one particular phase of the neutral gas. It is, in fact, only a practical choice. In Sect. 5.6, we discuss how our results change for different T_S.

In summary, the H I line profile is derived by first evaluating Eq. (2), then Eq. (5), and finally Eq. (3).

2.4. Fitting technique

In order to fit the ring parameters x = (v_φ,σ,n₀,h_s) to the data, we must define a way to compare the model H I line profile, , with the observed profile . This comparison is performed via a Monte Carlo Markov Chain (MCMC) algorithm. MCMC samples the parameter space in a way proportional to the probability of the model given the observations , which is proportional to the product of the likelihood and the prior (Bayes theorem). We define the likelihood as: (6)where ℛ(x) is the residual between the model and the data, ϵ is the error bar in the data (assumed to be constant throughout the whole dataset) and the sum is extended to all pixels in the data velocity domain, which extends from − 400 km s^-1 to + 400 km s^-1. We note that the choice of x that minimises ℛ is also the one that maximises the likelihood. The value of ϵ is important as it affects the width of the posterior probability distribution, which gives the uncertainty on the fit parameters. In principle ϵ should be the rms-noise of the LAB survey but, given that our model is a simple axisymmetric approximation of the Galactic disc, the requirement that it fit the data to the noise level is not realistic and would lead to a severe underestimation of the uncertainty associated to our results. Therefore we define an ad-hoc error ϵ^∗ as the ratio between the minimum value of ℛ(x) and the number of pixels used to calculate it, and for each sight-line we use either ϵ^∗ or the LAB rms-noise, whichever is the largest.

As fitting boundary conditions, we require that 0 <v_φ< 300 km s^-1, n₀> 0, h_s> 0. In addition, we require that , that is, the gas total velocity dispersion be at least equal to the thermal velocity dispersion, this latter value being 1.5 km s^-1 for our fiducial model (T_S = 152 K). We set our prior to be 1 in the allowed parameter ranges, and 0 otherwise.

In practice, when fitting the midplane emission (Sect. 2.1), we use 2000 iterations of the chain to evaluate ϵ^∗ and further 8000 iterations to sample the full posterior probability distribution for the parameters v_φ, σ and n₀. We use the mean and the standard deviation of these distributions as representative for the central value and the error bar of our parameters. The emission above the midplane (Sect. 2.2) is fit with a similar procedure, but we use 100 iterations for ϵ^∗ and a further 400 to sample the posterior probability, given that the midplane parameters are already set and we only fit for h_s.

In Appendix B we test our fitting technique on mock H I observations of an idealised disc, and we demonstrate that it can recover the system input parameters with very good accuracy.

3. Results: atomic hydrogen

3.1. Midplane: one-component model

Figure 2 shows the longitude-velocity (l−v) diagrams at latitude b = 0° predicted by our model (leftmost panels) for the receding and the approaching sides of the disc, and compares them with the observations (rightmost panels). The agreement with the LAB H I data is surprisingly good in most locations. In particular we are able to reproduce the brighter H I features visible in the data, which often occur at velocities different from the terminal velocities (like the bright H I spot visible around l = 350°,v = −10 km s^-1). However, the tail of the H I emission is more extended in the data than in our model, which suggests that the model underestimates the velocity dispersion of the gas close to the tangent points. This is not caused by our choice of the spin temperature (152 K), as the optically thin model (not shown here) appears identical to that shown in the left panels of Fig. 2. One possibility is that, as the H I brightness temperatures in the data span over three orders of magnitudes, our fitting technique (Sect. 2.4) is biased towards the brightest part of the profiles and is less sensitive to the faint H I emission. To test this, we re-performed our fit by substituting T_B with log (T_B) in Eq. (6). This significantly compresses the brightness temperature dynamical range. We found that, in this way, the H I tails are reproduced significantly better, but at the expense of the fitting of the brightest regions.

We conclude that our base model cannot reproduce the brightest and faintest regions of the H I Galactic disc simultaneously, and consider the inclusion of a second gas component in our disc model.

3.2. Midplane: two-component model

We now assume that the H I in the Galactic disc can be described by two components that share the same rotational velocity but have different densities and velocity dispersions. This increases the number of free parameters for each ring by two. Both components are fit simultaneously to the data by using the procedure described in Sect. 2.4, but we substitute T_B with log (T_B) in Eq. (6) to account for the large brightness temperature dynamical range.

The l−v diagrams for the best-fit model with two-components are shown in the middle panels of Fig. 2. We have assumed that both components have T_S of 152 K, but this choice does not affect the appearance of the l−v diagrams. The differences with respect to the single component model are minimal, except that now the tail of the H I emission is more extended and the data are reproduced remarkably well at all locations and in both quadrants. This can be better appreciated in Fig. 3, which shows three representative H I line profiles at l = 20°,40°, and 60° for the data and for both our models. Clearly, both models are in good agreement with the data, but the two-component model reproduces the tails of the profiles significantly better. Hereafter, we refer to this two-component, best-fit disc model with T_S = 152 K as our fiducial model.

Figure 4 shows the rotation curve, the velocity dispersion and midplane density profiles for our fiducial model. Profiles for the approaching and the receding sides of the Galaxy are shown separately, along with an average value representing a full axisymmetric model. The panel on top in Fig. 4 shows the best-fit values for rotation velocity (points + error-bars) as derived via our MCMC routine. A smooth rotation curve (solid thin lines) is derived by smoothing these points with a Gaussian kernel with 0.2 kpc of FWHM. In the smoothing procedure, each point has been given a weight inversely proportional to its error bar to ensure that the most uncertain values do not affect our final profiles. From this panel, it is clear that the H I rotation curves in QI and QIV are very similar to one another. Such similarity in v_φ points towards a high degree of axisymmetry in the inner Galactic disc. Differences in velocity larger than ~ 10 km s^-1 occur only at R< 2.5 kpc: at such small radii, we expect the Galactic bar to have a non-negligible impact on the kinematics of both gas (e.g. Fux 1999) and stars (e.g. Portail et al. 2017). Further discussion on the effect of the bar can be found in Sect. 5.4.

The central and bottom panels of Fig. 4 show the σ and n₀ profiles separately for the two components, and for both QI and QIV. Here, we show only the smoothed profiles. There is a clear distinction between the two components: the first one (component A) has low velocity dispersion (~ 8 km s^-1) and high density (~ 0.3−0.5 cm^-3), while the other (component B) has high σ (15−20 km s^-1) and low n₀ (~ 0.04 cm^-3). As mentioned, the high-σ component is the one needed to reproduce the faint tail of the line profiles. The contribution of this component to the H I mass within the solar circle is 15% in QI and 20% in QIV. The average velocity dispersion and density profiles (⟨σ⟩ and ⟨n₀⟩) are determined as The average σ is very large in the Galactic centre, where it reaches values of 25−30 km s^-1, but it decreases quickly with R, dropping down to ~ 10 km s^-1 around R = 1.5 kpc. For larger R, ⟨σ⟩ slowly decreases down to values of ~ 8 km s^-1 at R = R_⊙. Similar trends for σ(R) have also been derived in external galaxies (Tamburro et al. 2009). The averaged midplane density peaks in the Galactic centre but then shows a large depression at 0.5 <R< 2.5 kpc. At larger R, ⟨n₀⟩ oscillates around a typical value of 0.4 cm^-3. The density profiles of QI and QIV show a series of peaks and depressions, and are discussed in more depth in Sect. 5.5.

We discuss in more depth the nature of the second component in Sect. 5.3. Here, we stress that the averaged radial profiles of this new model are remarkably similar to those determined for our single component model (not shown here). The rotation velocities, in particular, are practically the same at all radii. The averaged velocity dispersions alone show some minor differences, which are however never larger than ~ 2 km s^-1 and typically below 1 km s^-1. This not only suggests that the models presented are indeed the best possible fit to the data, but also indicates the overall robustness of our fitting procedure.

3.3. Scale height and column density

Our fiducial model is the one that best reproduces the Galactic H I distribution and kinematics at b = 0°. We can now study how the gas distributes above the midplane using the method described in Sect. 2.2. Since we believe that the two H I components do not represent physically distinct phases of the Galaxy’s ISM – as discussed in Sect. 3.2 – we do not treat them separately and fit them both with a single scale height.

Figure 5 shows the l−v diagrams taken above and below the midplane at b = ± 1°, ± 2°, ± 4° for our fiducial case and for the LAB data. In general, the model follows nicely the trend shown by the data, where the high-velocity, low-level H I emission visible close to the galactic centre (l< 20° in QI and l> 340° in QIV) vanishes at | b | ≥ 2° because H I above the midplane is lacking in the innermost regions of the disc. In the model, l−v diagrams at opposite latitudes are identical by construction and are shown in the left-most panels of Fig. 5. This is not the case for the data, where l−v diagrams at opposite latitudes show some differences (in particular at b = ± 2°) that are likely caused by the presence of the Galactic bar and by the vertical displacement of the Sun with respect to the Galaxy midplane (e.g. Humphreys & Larsen 1995). Our model does not take these features into account and predicts l−v distributions that are in between those shown by the data. We also note that the model does not reproduce the data at b ≠ 0° as accurately as on the midplane. This is because we are assuming that the H I rotation velocity and velocity dispersion do not change with z, and attempt to fit a simple Gaussian for the vertical gas distribution of each ring. We discuss this further in Sect. 5.7.

Figure 6 shows the resulting scale height and H I surface density profiles for our fiducial model. The H I scale height increases inside-out, moving from ~ 100 pc for R< 3 kpc to ~ 200 pc at larger radii. Also in this case the regions QI and QIV show very similar trends, pointing towards a high degree of symmetry in the inner disc. Unfortunately, a precise value for the scale height at R> 7 kpc is difficult to derive, as the H I profiles at l> 60° (QI) and l< 300° (QIV) change very little with b. Therefore, for the purpose of determining the surface density at radii larger than 7 kpc, we fixed the scale height to a value of 217 pc, derived by averaging the approaching and receding h_s at R = 7 kpc. The bottom panel of Fig. 6 compares the H I surface density Σ_HI predicted by our model (thick solid line) with a selection of those presented in the literature (dashed lines). In general, the Σ_HI of our model is in very good agreement with that determined by Binney & Merrifield (1998), while it lies above that determined by Nakanishi & Sofue (2003) and well below that predicted by (Kalberla & Dedes 2008). In addition, we predict the H I depression at 0.5 <R< 2.5 kpc to be deeper than all previous determinations, and a steep rise in the inner few hundred parsecs.

4. Results: molecular hydrogen

As is often done in the literature, we assume that most molecular hydrogen is traced by the CO (J = 1−0) emission line, whose source of excitation is primarily collisions with H₂ molecules, and use the CO data of Dame et al. (2001) to model the Galactic H₂ inside the solar circle. We caution that there exists a fraction of molecular hydrogen that is not coupled with CO, and therefore remains “invisible” in our data (the so-called “CO-dark” H₂, Langer et al. 2014; Tang et al. 2016). What we model here is instead the “CO-bright” phase of the molecular hydrogen.

The modelling of the Galactic CO follows a prescription similar to that adopted for the H I, with the following differences. First, we use only a single gas component. Second, we do not consider the effect of the optical thickness. Even though molecular clouds are optically thick in the CO (J = 1−0) rotational line, it is commonly assumed that their ensemble is not; clouds do not shade each other because they are well separated in both space and velocity (e.g. Gordon & Burton 1976; Binney & Merrifield 1998). Thus a typical CO line-profile is produced by the sum of the profiles of each individual cloud. Assuming that the population of molecular clouds does not vary within the inner Galaxy, the integral of the CO line profile between two line-of-sight velocities v₁ and v₂ can be considered proportional to the H₂ column density: (9)where T_B is the brightness temperature of the CO line and X_CO is the CO-to-H₂ mass conversion factor (or simply X_CO factor). Following Bolatto et al. (2013), we adopt X_CO = 2.0 × 10²⁰ cm^-2 K^-1 km^-1 s, with an uncertainty of ± 30%. Equation (9) can be inverted to derive T_B(v), and the quantity N_H2(v) can be evaluated with a formula analogous to Eq. (2). We assume v_φ = v_⊙, σ = 4 km s^-1 and n₀ = 0.5 cm^-3 for the first ring of our model (see Gordon & Burton 1976). Finally, in our fit we do not impose a lower limit on the velocity dispersion because of the low kinetic temperature and large molecular weight of the CO.

Figure 7 shows the l−v diagrams of the midplane CO emission for our model (on the left) and for the data (on the right). Comparing the H I and the CO data reveals some important differences. As molecular gas is organised mostly in clumps which are distributed along spiral arms, the CO emission appears to be more patchy (in the l−v space) than the H I emission, despite the fact that both datacubes have been smoothed to the same resolution. Clearly, this patchiness cannot be properly reproduced by our model, which assumes a smooth and axisymmetric distribution of molecular gas. Apart from this difference, the model reproduces quite well the pattern shown by the data. This is especially evident in QI, where several features that are not located in the proximity of the tangent velocity – and therefore are not directly fit – are well predicted by the model. The most evident of these features originates around l = 0°,v = 0 km s^-1 and crosses the diagram diagonally up to l ~ 30°,v = 120 km s^-1, representing part of the Galactic “molecular ring” (e.g. Clemens et al. 1988). As we can reproduce reasonably well both the peaks and the tails of the observed CO profiles, there is no need to include a second component in the line-fitting.

The rotation curve, the velocity dispersion and midplane density profiles for the Galactic H₂ are shown in Fig. 8. As already noticed for the H I, the molecular gas of the Milky Way also shows a remarkable degree of symmetry between the approaching and the receding sides of the disc. For both v_φ and σ, the main differences between the two sides arise at R< 1.5 kpc, where non-circular motions induced by the bar play an important role. The H₂ velocity dispersion peaks at the Galaxy centre with a value of 15−20 km s^-1. Similarly to the atomic hydrogen, σ drops quickly within the inner 1 kpc down to a value of ~ 4 km s^-1, and it slowly decreases at larger radii. Despite the molecular gas being concentrated in spiral arms, the midplane H₂ densities in QI and QIV follow a similar global trend. The exception is that, in a few locations (R = 2.7,3.5,5.2 kpc), peaks in QI correspond to depressions in QIV (or viceversa). The atomic gas shows similar features (Sect. 3.2), though their locations do not perfectly overlap with those of the H₂.

We now model the CO distribution in the direction perpendicular to the midplane. Our procedure is analogous to that used in Sect. 3.3 for the H I, with two differences. First, we limit the fit to the latitudes | b | < 3°, because a) the signal-to-noise of the data decreases significantly at higher latitudes; and b) in the regions of our interest, the l−v diagrams at | b | > 3° are similar to those at | b | = 3° and do not seem to contain additional information on the vertical distribution of the CO. Second, we impose h_s< 200 pc as a prior. Our fitting technique does not converge well if we leave the CO scale height completely free to vary, and it is reasonable to assume that the disc of molecular gas is thinner than the H I disc in our Galaxy.

Figure 9 compares some representative l−v diagrams extracted from our best-fit model at different latitudes with those of the data. We note that the noise in the data increases significantly with | b |. As expected, the main difference between the model and the data is that the latter is more clumpy and shows more structures, suggesting that the molecular gas is not diffusely distributed in the disc as the model assumes. However, the model reproduces quite well the observed drop in CO emission as a function of latitude, indicating that the scale heights that we are determining are reliable.

The molecular hydrogen scale height and surface density profiles are shown in Fig. 10. The scale heights in the two quadrants show a very similar trend with R. The average h_s declines slowly with R, varying from ~ 100 pc in the Galactic centre to ~ 50 pc in the solar neighbourhood. As for the atomic hydrogen, in this case we also prefer to fix the scale height at R> 7 kpc to the value of 52 pc, determined by averaging the approaching and the receding h_s obtained at R = 7 kpc. In a recent review, Heyer & Dame (2015) presented a compilation of H₂ scale-height profiles derived in the last 30 years via different methods applied to different datasets. In all cases, the H₂ HWHM scale-height remains approximately flat in the region 3 < R < 8 kpc, with mean and standard deviation of 50 ± 11 pc, compatible with the values determined here (see also Table 1). The H₂ surface density profile indicates that most gas is concentrated either close to the Galactic centre or in the region 2.5 < R < 5.0 kpc, while a large depression occurs at 1 <R< 2.3 kpc which was already noticed by Sanders et al. (1984). The average Σ_H2 is characterised by three distinct peaks: the first around R = 3.3 kpc, the second at R = 4.5 kpc and the third at R = 6.4 kpc. The first peak is however a by-product of our averaging, as the CO distribution in QI and QIV in the region 2 < R < 4 kpc is opposite in nature. The other two peaks are instead “genuine”, being present in both the approaching and the receding regions – although the third is especially prominent in QIV. Our H₂ surface density profile is compatible with that determined by Bronfman et al. (1988; dashed line in Fig. 10), but shows more structures: peaks are higher and depressions are deeper. The agreement with the profile derived by Nakanishi & Sofue (2006) is however less good, especially at R < 2.5 kpc. We note that the Bronfman et al. and the Sofue et al.’s profiles do not agree with each other at R < 3 kpc.

5. Discussion

5.1. Comparison between atomic and molecular gas

The results obtained in this work are summarised in Fig. 11, which shows the radial profiles for the rotation velocity (panel a), surface density (panel b), velocity dispersion (panel c), midplane density (panel d), and half width at half maximum scale height (panel e) for both the atomic and the CO-bright molecular hydrogen of the Milky Way. As in the previous sections, each profile is derived by averaging the approaching and the receding parts, which are smoothed to the spatial resolution of 0.2 kpc, while Eqs. (7) and (8) are used to compute the average H I velocity dispersion and midplane density. We show one point per resolution element.

The uncertainty σ_tot on each quantity, which is represented by the error bars in Fig. 11, is derived by summing two terms in quadrature. The first term is σ_asym, which accounts for the asymmetries between the two regions QI and QIV. Following Swaters (1999), we estimate σ_asym as the difference between the (smoothed) approaching and receding profiles divided by 4. For the H I velocity dispersion and midplane density, σ_asym is computed by subtracting – instead of adding up – the two terms in the square brackets of Eqs. (7) and (8), and dividing the result by 4. The second term is σ_fit, which takes into account the MCMC error-bars given by the standard deviations of the posterior probability distributions as described in Sect. 2.4. For a given parameter at a given R, σ_fit is computed as the median MCMC error-bar within 0.2 kpc from that radius. Error-bars from QI and QIV are both included in the computation of σ_fit. An additional 30% uncertainty on the H₂ midplane and column density is summed in quadrature to account for the uncertainty in the X_CO factor (Bolatto et al. 2013). Finally, a generous ad-hoc 40 pc error is assumed for both the H I and H₂ scale heights at R> 7 kpc, where our modelling is poor and h_s is fixed to a constant value.

The similarities between the atomic and molecular gas in terms of kinematics and distribution are striking. The H I and H₂ rotation curves remarkably overlap with each other. They both rise quickly, peak to ~ 260 km s^-1 at 0.6 kpc from the centre, decline down to 210 km s^-1 at R = 2.5, and then gently rise again. The rotation velocity flattens around R = 6.5 kpc to a value that depends on the choice of v_⊙, which here is 240 km s^-1 (see Appendix C for the rotation curves re-scaled to v_⊙ = 220 km s^-1). The main difference between the two components is that that the rotation traced by the H₂ shows several “bumps” that pulls the velocity downwards with respect to the values determined for the H I. Our interpretation is that, given that the molecular gas a) is less diffuse and has a more “patchy” distribution with respect to the H I, and b) is likely to be more concentrated around the Galaxy spiral arms, there is no guarantee that some CO emission will always occur at the tangent point for all lines-of-sight. When this does not happen, our fitting algorithm interprets the lack of emission in terms of a drop in v_φ. These bumps would therefore be artificial and caused by the fact that H₂ is not homogeneously distributed within the disc. However, the magnitude of these velocity drops is very small, typically less than 10 km s^-1, suggesting that the molecular gas in the Galactic disc can be treated as a diffuse and continuous component, at least as a first approximation.

The H I and H₂ surface density profiles are characterised by large error bars, due to the propagation of the uncertainties on the scale heights and midplane densities. In spite of that, they show several similarities: they both peak around the Galactic centre, have a deep depression in the region 0.5 <R< 2.5 kpc and then rise again at larger distances. We find all these correlations to be very significant, especially because they have been derived by a blind fit to two totally independent datasets. The surface density in the Galactic centre is dominated by the molecular gas, but its precise value is highly uncertain: we estimate 12 ± 2 M_⊙ pc^-2 for the H I and 115 ± 40 M_⊙ pc^-2 for the H₂. The H I density profile at R> 3 kpc is fairly flat, between 3 and 7 M_⊙ pc^-2, whereas the H₂ density declines with radius and abruptly drops down to ~ 1 M_⊙ pc^-2 at R> 7 kpc. Apart from the general trend, both surface density profiles show two clear peaks at R = 4.5 kpc and R = 6.4 kpc. In Sect. 5.5 we discuss two possible origins for these features: one is that they are due to actual ring-like overdensities in the gas distribution, and the other is that they arise because, at those particular locations, the line-of-sight intersects the regions of spiral arms, where the gas density is above the mean. Interestingly, there is little trace of these features in the velocity dispersion profiles or in the H₂ scale height profile, while the H I scale height peaks at slightly lower radii with respect to the density enhancements. An additional peak can be seen around R ≃ 3 kpc but, as mentioned in Sect. 4, at 2 <R< 4 kpc, the H₂ distributions in QI and QIV are anti-correlated thus the correspondence of the average H₂ profile with the H I profile in this regions could be, to some extent, coincidental.

Table 1 shows the typical properties of the Galactic atomic and molecular hydrogen averaged in the region 3 <R< 8.3 kpc, where fluctuations in the various quantities are limited, along with their total masses within the solar circle. For each model parameter we show the mean value computed in the region considered, the standard deviation around the mean, and the mean uncertainty – that is, the mean error bar shown in Fig. 11 – in that region. Values for the rotation velocity depend on the choice of v_⊙ (see Appendix C), and therefore are not shown. We note that, given that the molecular gas density is a decreasing function of R, its mean value is not truly representative for the density in the region considered. The total masses within the solar circle are computed as the mean and the standard deviation of 10⁴ different realisations, each derived by randomly extracting the H I and H₂ surface densities profiles from the values (points with error bar) shown in panel (b) of Fig. 11³. Our fiducial value for the H₂ mass within the solar circle is in broad agreement with that of Heyer & Dame (2015). We stress that these masses, even when corrected for the helium content, do not represent a complete census of the gas content in the inner disc of the Milky Way. Diffuse ionised gas contributes around one third of the mass of the neutral atomic component (Reynolds 1993). Also, studies of [C II] emission-lines have revealed the presence of a CO-dark H₂ component that may encompass a significant fraction (20−50%) of the molecular gas in the Galaxy (Langer et al. 2014). Clearly, these components are missed by our analysis.

In order to derive the gas properties in the solar neighbourhood, it is best to extrapolate from measurements at lower radii, as the uncertainties associated to the measurements at R = R_⊙ are large. The extrapolation is computed as follows. First, for each parameter radial profile, we produce 10⁴ realisations by randomly extracting the parameter values in the region R> 7.0 kpc from the points with error bars shown in Fig. 11. Each realisation is then fit with a straight line, and the value of the fit extrapolated at R = 8.3 kpc is recorded. In Table 2 we show, for the various parameters, the mean and the standard deviations of these extrapolated values. Clearly, the properties of the H I in the solar neighbourhood are similar to those determined at lower radii, while the molecular gas is less abundant at R = R_⊙ with respect to the rest of the inner disc.

5.2. An independent estimate for the local H I column density

Recently, McKee et al. (2015) re-analysed archival data in order to provide a census of surface densities for different baryonic components in the solar neighbourhood. Using previous results by Heiles (1976) and Dickey & Lockman (1990), and adopting a correction for the optical depth based on Strasser & Taylor (2004), they derive an H I surface density in the solar neighbourhood Σ_{HI, ⊙} of 7.8 ± 1.2 M_⊙ pc^-2 (from their Table 2, after removing the 40% correction for the He and heavier elements), which is in tension with our determination of 4.5 ± 0.7 M_⊙ pc^-2. They also quote a fiducial value for the local H₂ surface density of 0.7 ± 0.2 M_⊙, which is compatible with our estimate given our large error-bar.

One might indeed doubt the values presented in our Table 2, given that our modelling technique is less robust in the proximity of the solar circle. For this reason, we now provide a model-independent estimate of Σ_{HI, ⊙} using an approach similar to that adopted by Heiles (1976), which is based on a simple geometrical relation between the column density perpendicular to the disc and that measured towards a generic line of sight.

First, we mask in the LAB datacube the H I emission produced by high and intermediate-velocity clouds, which can potentially contaminate our computation. This masking follows the procedure described by Marasco & Fraternali (2011) in their Sect. 2.1, and is based on a 35 km s^-1 cut in “deviation” velocity (see Wakker 1991). As discussed in Marasco & Fraternali (2011), the masking is very conservative and indeed we find that it has only a marginal impact on our results, as by neglecting it, the surface density increases by only ~ 0.2 M_⊙ pc^-2. Next, we produce an all-sky H I column density map by integrating in velocity all line profiles in the (masked) LAB data. From this map, at any given latitude b we compute a median column density by considering many sight-lines: consistently with our previous analysis, we focus separately on 0°<l< 90° (region QI) and on 270°<l< 360° (region QIV), for which we compute two distinct . Under the assumption of a constant gas scale-height within ~ 1 kpc from the Sun, Σ_{HI, ⊙} can be computed as (10)where is the correction for the H I optical depth () and Σ_HI,LB(b) is a sight-line-dependent correction due to the Local Bubble, a severely gas-depleted region that surrounds the solar neighbourhood (e.g. Cox & Reynolds 1987). Here we adopt , in agreement with our findings (see Sect. 5.6) and with those of Strasser & Taylor (2004). Given the complex shape of the Local Bubble, the rightmost term in Eq. (10) is difficult to determine and we ignore it for a moment.

In Fig. 12 we show how Σ_{HI, ⊙} varies as a function of | b |, for | b | > 10°, and for regions QI and QIV separately. Ideally, Σ_{HI, ⊙} should stay constant at all latitudes. Instead, there is a clear discontinuity around | b | ~ 45°, above which the mean column density (shown by green circles in Fig. 12) drops by a factor of ~ 2: this can be readily explained as being due to the “hourglass” shape of the Local Bubble (Lallement et al. 2003) that – extending preferentially in the direction perpendicular to the midplane – drastically contaminates the column density measurements at high latitudes. At lower latitudes, however, the contamination should be largely reduced as the Bubble “wall” is much closer to the Sun. For instance, assuming a distance to Bubble wall of 60 pc (Lallement et al. 2003), we find an overall correction to Σ_{HI, ⊙} of about 0.2 M_⊙ pc^-2 at b = 10° using the mean H I disc parameters of Table 1. In the region 10°< | b | < 45° the mean Σ_{HI, ⊙} seems approximately flat: here we derive a H I surface density of 5.3 ± 0.5 ± 0.2 M_⊙ pc^-2, where the first term is the mean in the region considered, followed by the standard deviation and by the mean uncertainty computed as the difference between QI and QIV divided by 4. We find a slightly lower and more precise value (5.0 ± 0.3 ± 0.1 M_⊙ pc^-2) if we also exclude data at 28°< | b | <°40, where QI and QIV show a significantly different trend.

These estimates are compatible with – but slightly larger than – the surface density reported in Table 2. While this suggests that the value of Σ_{HI, ⊙} shown in Table 2 may be slightly underestimated, in Sect. 5.7 we argue that the whole H I surface density profile may be overestimated because of contamination from extra-planar gas. In light of these considerations, we still take Σ_{HI, ⊙} = 4.5 ± 0.7 as our fiducial value.

5.3. On the nature of the second H I component

It is hard to tell whether or not the additional low-density, high-σ component used in Sect. 3.2 to fit the Galactic H I data is spurious in nature, that is, simply serves the purpose of absorbing the insufficiencies of our model, or has true physical meaning. On the one hand, the two components are not related to the classical “cold” and “warm” neutral ISM (e.g. Wolfire et al. 2003), as their velocity dispersion is too high to be caused by pure thermal broadening. On the other hand, double-component H I profiles are observed also in external galaxies: Ianjamasimanana et al. (2012) studied stacked profiles in The H I Nearby Galaxy Survey (THINGS; Walter et al. 2008) and found that they can be decomposed into a narrow and a broad component, with velocity dispersions analogous to those found here.

One possibility is that the second, high-σ component is required to account for the various non-axisymmetric features along the line of sight, like streaming motions in the spiral arms or stochastic substructures. These features locally modify the gas density and velocity dispersion and might induce deviations from the pure differential rotation. All these effects contaminate the H I line profiles in a way that is difficult to predict, but that can be globally taken into account by including a second, “artificial” H I component with a different density and velocity dispersion.

Another possibility is that diffuse star formation in the disc injects turbulence in a fraction of the H I medium. In this case, the high-σ component would represent a genuine, turbulent phase of the Galactic ISM. The kinetic energy associated to our second component is 6.3 × 10⁵³ erg, assuming a constant velocity dispersion of ~ 15 km s^-1 and accounting for the He content. Assuming a turbulence dissipation timescale of 10⁷ yr (Tamburro et al. 2009), an energy injection rate of 2 × 10³⁹ erg s^-1 would be required to maintain the gas in this state. This rate can be easily achieved by supernova feedback: for a Galactic supernova rate of two events per century (Diehl et al. 2006), and assuming the canonical 1 × 10⁵¹ erg per event, a feedback efficiency of only 0.3% would suffice to match this requirement. Thus this second scenario is energetically plausible.

A third possibility is that this high-σ component is the midplane counterpart of the high-latitude Galactic extra-planar H I. Extra-planar gas is indeed more turbulent than the gas in the disc and, overall, contributes to only a small fraction (5−10%) of the total H I mass (Fraternali et al. 2001; Oosterloo et al. 2007; Marasco & Fraternali 2011). If the extra-planar H I is produced by the Galactic fountain mechanism (Marasco et al. 2012), the neutral phase of the returning fountain flow could be observable on the midplane as a low-density, high-σ component⁴. In this scenario, as in that previously discussed, it is feedback from star formation that ultimately produces the second component. A possible argument against this scenario is that, if the fountain clouds that return to the disc have substantially lower rotational speed with respect to the midplane gas, their emission should be systematically shifted to velocities lower than the terminal velocity, and therefore their contribution to the H I profiles should be limited to the brightest regions of the l−v diagram and not appear in its high-velocity, low-T_B tail.

5.4. The effect of the bar

Perhaps the main weakness of our model is that it is based on the assumption that the gas follows pure circular orbits, which would only be correct if the gravitational potential was perfectly axisymmetric. It is well established, however, that the Milky Way is a barred galaxy (e.g. Dwek et al. 1995; Binney et al. 1997). The first attempt to characterise the gas distribution and kinematics in the central regions of the Galaxy was made by Liszt & Burton (1980), who proposed a scenario where the gas in the innermost 2 kpc moves along elliptical paths in a ringlike structure which is tilted with respect to disc midplane (see also Ferrière et al. 2007). Later on, numerical simulations of gas flows in a barred potential have revealed that, in the presence of a bar, the gas kinematics and distribution become asymmetric and non-stationary (Fux 1999; Rodriguez-Fernandez & Combes 2008; Sormani et al. 2015). As a consequence, the exact shape of the H I and CO longitude-velocity diagrams in the innermost regions of the Galaxy – say 340°<l< 20° – will depend on the details of the bar shape and pattern speed, on its orientation with respect to our view-point, and on the observing epoch.

In general, the rotation curve at R ≲ 2.5 kpc shown in Fig. 11 should not be taken as representative for the circular velocity in that region. The circular velocity should rise more gently than what can be derived from a simple tangent-point analysis of the l−v diagrams: this was first recognised by Binney et al. (1991), and considered recently by dynamical models of the inner Galactic disc based on the density and kinematics of stars (Portail et al. 2017; see also Fig. 17). Unfortunately, the only “correct” approach to study the gas dynamics in the innermost few kpc of the Galaxy is to follow the gas flow via hydrodynamical simulations and adjust the potential parameters in order to reproduce the details of the l−v diagrams, like in Sormani & Magorrian (2015). This is, however, beyond the scope of the current work.

5.5. The effect of the spiral arms

As already mentioned, the atomic and the molecular gas distributions show two prominent peaks at R = 4.5, 6.4 kpc. Remarkably, these peaks can be seen in both the receding and the approaching regions of the disc (Figs. 6 and 10), and they correspond to distinctive “bumps” visible in the l−v diagrams around l = 30°,50° (QI) and l = 310°,330° (QIV). We noted that a third gas overdensity can be seen around R = 3 kpc: this is visible in both QI and QIV in H I, but it is separated in two distinct peaks in H₂. The fact that at least two of these overdensities are visible in both quadrants is significant. The simplest explanation for these features would be the presence of two distinct ring-like gas overdensities within the Galactic disc at radii of 4.5 and 6.4 kpc. Ring-like features are observed in the neutral gas of nearby galaxies, such as NGC 3184 or NGC 2841 (Walter et al. 2008).

Another possibility is that these features arise when the line-of-sight intersects the regions of spiral arms, where the gas density is above the mean. Unlike the previous scenario, in this case the density enhancement would be local, thus the derived column density would not be representative for the mean column density at that radius. Burton (1971) already noticed the “bumps” in the H I l−v diagram of QI, and interpreted them in the framework of the density-wave theory (Lin et al. 1969) as streaming motions induced by the Galactic spiral arms. According to Burton (1971), streaming motions alone are sufficient to explain the irregularities in the l−v diagram of QI, even in the presence of a uniform gas density distribution (see his model V). Clearly, our model does not include self-consistently non-circular motions induced by the spiral arms. However, we can check whether or not the sight-lines that are responsible for the peaks in the density distributions intersect the Galactic spiral arms at the locations where the overdensities occur.

In Fig. 13 we show a face-on view of the inner disc of the Milky Way with the four-arm spiral pattern of Steiman-Cameron et al. (2010) overlaid on top. The spiral arm model is based on all-sky intensity maps of [C II] and [N II] lines, which trace the cool phase of the Galactic ISM. The dashed red and blue lines in Fig. 13 represent the sight-lines (in QI and QIV, respectively) for which the tangent point radius (R_⊙sin(l)) is 3,4.5 or 6.4 kpc. Clearly, the gas overdensity at 6.4 kpc can be explained as the l = 50° sight-line in QI intercepts the Sagittarius arm, while the l = 310° sight-line in QIV intercepts the Scutum-Crux arm. The overdensity at 4.5 kpc appears to be produced as the sight-lines in QI and QIV intercept the Scutum-Crux and the Norma arm, respectively, although in this case, the correspondence between spiral arms and tangent points is less precise. Finally, the feature around 3 kpc is likely related to the long-known 3-kpc arm (van Woerden et al. 1957; Dame & Thaddeus 2008) (not shown in Fig. 13), or, more generally, to the region where the Perseus and the Scutum-Crux arms connect to the edge of the Galactic Bar.

5.6. The effect of the spin temperature

We now discuss the effect of different H I spin temperatures on our results. In Fig. 14 we compare the properties of the H I disc derived by our fiducial model (T_S = 152 K, shaded region) with those predicted by our best-fit optically thin (T_S = 8000 K, dot-dashed line) and our best-fit optically thick (T_S = 80 K, dashed line) cases. Clearly, differences concern mainly the gas density. In the optically thin regime, n₀ is smaller because the H I emission is not partially self-absorbed like in the fiducial case, so less gas is required to produce the same flux. The opposite happens in the optically thick regime. We note that the density profiles predicted by these three regimes are not consistent with one another. The H I velocity dispersions predicted by the three regimes are virtually identical, despite the fact that, at higher temperature, H I has a larger thermal broadening, corresponding to a higher velocity dispersion floor (see Sect. 2.4).

Since the scale heights of the three models are almost the same at all radii, the difference in the midplane density implies a difference in surface density and total mass. The model with T_S = 80 K reaches column densities up to ~ 10 M_⊙ pc^-2, similar to those found by Kalberla & Dedes (2008, see also Fig. 6, but we exclude that such a low T_S is representative for the typical spin temperature of the Galactic H I, as we discuss in Appendix A. The H I masses within the solar circle predicted by these models are 6.4 × 10⁸M_⊙ for the optically thin regime and 1.4 × 10⁹M_⊙ for the optically thick one, which corresponds, respectively, to 30% less and 52% more than the mass predicted by our fiducial model (9.1 × 10⁸M_⊙). Interestingly, a 30% discrepancy in H I mass between an opaque and an optically thin regime has also been found by Strasser & Taylor (2004) in our Galaxy and by Braun et al. (2009) in M31. We note that different spin temperatures have little impact on the H I surface density in the solar neighbourhood, which sits robustly around 5 M_⊙ pc^-2; nevertheless we reiterate that this value is quite uncertain given that the scale height at R> 7 kpc cannot be accurately determined.

5.7. The effect of the extra-planar H I

The H I scale height and column density have been derived by assuming that the gas kinematics above and below the midplane are the same as those at z = 0. However, Marasco & Fraternali (2011) found that the Milky Way is surrounded by a faint (~ 5−10% of the total H I mass) layer of slow-rotating, inflowing extra-planar H I, analogous to that found in external galaxies like NGC 2403 (Fraternali et al. 2002), NGC 891 (Oosterloo et al. 2007) and NGC 3198 (Gentile et al. 2013). This medium extends up to a few kpc above the disc and is thought to be produced by the Galactic fountain mechanism (Marasco et al. 2012). As we test in Appendix B, neglecting this additional component leads to overestimation of the H I scale height and, as a consequence, the H I column density. A way to tackle this problem is to fit for the scale height by including in the model an additional component that represents the Galactic extra-planar H I. The parameters of this component are taken from Marasco et al. (2013). This approach has been tested on mock H I observations of a system consisting of a disc plus an extra-planar component, and leads to the correct recovery of the input disc scale height (see Appendix B for the details).

In Fig. 15 we compare the scale height and surface density profile of our fiducial model, which does not include the extra-planar H I, with those derived by including this additional component. Differences are small, and arise only at R> 3−4 kpc. The surface density predicted by this new model is slightly lower than, but still consistent with, that predicted by our fiducial case. The disc H I mass within the inner Galaxy decreases by about 12% (from 9.0 to 7.9 × 10⁸M_⊙) when including the extra-planar gas, which is less than the difference that we would derive by neglecting the H I opacity in our calculations (Sect. 5.6). We conclude that the presence of the extra-planar H I has only a modest impact on our results. By analogy with Tables 1 and 2, in Tables 3 and 4 we report, respectively, the typical gas properties averaged in the region 3 <R< 8.3 kpc and those extrapolated to R = R_⊙ for this new model that includes the extra-planar gas.

5.8. Comparison with external galaxies

How do the gas properties of our Galaxy compare with those determined in nearby spirals? It is not trivial to find in the literature a proper comparison sample, which should contain spirals with a Milky-Way-like morphology and with spatially resolved H I and H₂ data. Here, we have focused on the sample of Leroy et al. (2008), who have derived the atomic and molecular hydrogen surface density profiles in 23 nearby galaxies using H I data from THINGS (Walter et al. 2008), CO (J = 1−0) data from the BIMA Survey of Nearby Galaxies (BIMA SONG; Helfer et al. 2003), and CO (J = 2−1) data from the HERA CO-Line Extragalactic Survey (HERACLES; Leroy et al. 2009). Within the sample of Leroy et al., we have selected those systems that are classified as SAB or SB according to RC3 (de Vaucouleurs et al. 1991): NGC 3184, NGC 3198, NGC 3351, NGC 3521, NGC 3627 and NGC 6946. In addition, we have included in our sample the H I surface density profiles of NGC 925 and NGC 2903, both SAB systems, from de Blok et al. (2008). Finally, we have complemented our sample with the H I and H₂ density profiles of the SAB galaxy NGC 5236 (M83)⁵ from Crosthwaite et al. (2002). All systems have numerical Hubble T values ranging from 3 to 7, equivalent to b-to-d stages in the classical de Vaucouleurs classification.

In Fig. 16 we compare the H I and H₂ surface density profiles of the systems in our sample with those of the Milky Way. Overall, the H I radial distribution of our Galaxy does not differ from that of the other spirals, especially at R> 3 kpc. With the exception of NGC 925, all galaxies show a depression in their inner (R ≲ 3 kpc) H I content, analogous to – but less prominent than – that of the Milky Way. Differences are more evident when it comes to the molecular gas content, as the Milky Way appears to be one of the most H₂-poor galaxies in the sample, along with NGC 3198 and NGC 3351 (shown by solid lines in Fig. 16). Interestingly, these two galaxies are the only SB systems in our sample, and they both have a H₂ distribution analogous to that of the Milky Way, including a depression around R ~ 2 kpc similar to that of our Galaxy. Both NGC 3198 and NGC 3351 are however less massive than the Milky Way, having v_flat (i.e. the rotational speed in the region where the rotation curve is flat) of 150 km s^-1 (de Blok et al. 2008) and 196 km s^-1 (Leroy et al. 2008) respectively⁶. We note that virtually all galaxies show some features in their H₂ surface density profile around R = 2−3 kpc.

Even though our galaxy sample is small, our analysis suggests that the way the cold gas is distributed in our Galaxy is not untypical for a barred spiral, and in particular for an SB system. There is general consensus that gas flow in a barred potential can produce a depression in the gas surface density in the regions of the bar (e.g. Athanassoula 1992). However, the details of the gas distribution within and around the bar depend on the physical condition of the bar, in particular on its pattern speed and exact shape (e.g. Sormani et al. 2015). We argue that variations in the bar shape and speed can justify the variety in surface density profiles shown by the systems in our sample.

In Fig. 17 we compare the H I rotation curve of the Milky Way with those of the other galaxies in our sample. The rotation curves of the external galaxies are taken from the Spitzer Photometry & Accurate Rotation Curves (SPARC) database (Lelli et al. 2016; NGC 2903, NGC 3198, NGC 3521, NGC 6946), from de Blok et al. (2008; NGC 925, NGC 3627) and from Crosthwaite et al. (2002; NGC 5236). NGC 3184 and NGC 3351 are part of the THINGS sample but have not been studied by de Blok et al., given their low inclination. Despite their variety in dynamical range, all rotation curves in our sample of external galaxies (except for NGC 6946) show a shallow rise in the inner parts followed by a flattening in the outer regions. The Galactic rotation curve follows instead a different trend: its steep rise followed by wiggles and a rapid decline in the innermost ~ 2.5 kpc strikes as a unique feature. As discussed in Sect. 5.4, the Galactic rotation curve at R ≲ 2.5 kpc is strongly contaminated by non-circular motions triggered by the bar and it is hardly representative for the circular velocity in that region. It is hard to compare the impact of non-circular motions on the derived rotation curve of the Milky Way with that on the rotation curves of external galaxies, given the completely different technique adopted to extract the rotation velocities. It is quite possible that the ‘true’ rotation curve of our Galaxy, once cleaned from non-circular motions, may appear similar to some of those in our sample (for instance, like that of NGC 3521).

6. Conclusions

We have modelled the kinematics and the distribution of the atomic and CO-bright molecular hydrogen in the inner (R<R_⊙) disc of the Milky Way. We have assumed that the Galactic disc can be decomposed into a series of concentric and co-planar rings, and that the gas in each ring can be fully described by four parameters: the rotation velocity, the velocity dispersion, the midplane density and the scale height. The model is axisymmetric, but we have modelled the receding (QI) and the approaching (QIV) sides of the Galaxy separately. Non-circular motions are not included in our model.

The parameters of each ring are fitted to the H I data of the LAB survey (Kalberla et al. 2005) and to the CO data of Dame et al. (2001). Our fitting strategy consisted in simulating the various H I and CO line profiles in the regions of interest and comparing them with the observed profiles. The fit proceeds iteratively from the outermost sight-lines (l = 90° for QI or l = 270° for QIV) to the innermost one (l = 0°), so that the parameters of each ring depend on those determined for the other rings at larger radii. With our method, for a given ring a) all parameters can be fit simultaneously to the data without assuming a-priori a radial profile for the gas density or velocity dispersion, b) the terminal velocity in the line profile is not related a-priori to the gas at the tangent point radius. To our knowledge, this is the first time that such a novel approach has been used to model both the H I and the H₂ in the Galaxy.

Our results can be summarised as follows:

• Our disc model reproduces very well the l−v diagrams for both the H I and the CO emissions in the inner Galaxy up to a few degrees above the midplane. The agreement with the H I is more evident as the atomic gas is distributed more smoothly within the disc.

• The approaching and the receding sides of the inner Galaxy show very similar rotation, velocity dispersion and scale height profiles, suggesting a large degree of axisymmetry in the inner disc. Discrepancies between the two sides occur mainly at R < 2.5 kpc and are most likely induced by the Galactic bar.

• The H I and H₂ rotation curves overlap remarkably well with each other. They rise in the centre and reach a flat part at R> 6.5 kpc (for v_⊙ = 240 km s^-1). Some bumps visible in the molecular gas rotation are likely caused by the patchy distribution of molecular clouds in the disc.

• The H I (H₂) velocity dispersion peaks at the Galactic centre with a value of ~ 30 km s^-1 (~ 20 km s^-1), decreases rapidly in the innermost 1.5 kpc (1 kpc) and remains constant at larger radii, with a typical value of 8.9 km s^-1 (4.4 km s^-1). In the solar neighbourhood, we derive a value of ~ 7.8 km s^-1 (~ 3.7 km s^-1).

• The H I (H₂) scale height is fairly flat for R > 3 kpc, with a Gaussian HWHM of 202 ± 28 pc (64 ± 12 pc).

• The H I and H₂ surface density profiles show several common kpc-scale features; they both peak at Galactocentric radii of 4.5 and 6.4 kpc, and show a deep depression at 0.5 <R< 2.5 kpc. The peaks represent either genuine, axisymmetric ring-like overdensities within the disc, or local features due to the line-of-sights intercepting the Galactic spiral arms. The depression is most likely caused by the Galactic bar. The profiles are in good agreement with those determined by Binney & Merrifield (1998), but are in tension with other determinations, in particular with the H I surface density derived by Kalberla & Dedes (2008).

• We estimate an H I (H₂) surface density in the solar neighbourhood of 4.5 ± 0.7 M_⊙ pc^-2 (1.3 ± 0.7 M_⊙ pc^-2).

• When compared with nearby barred spiral galaxies, the gas distribution in the Milky Way is not exceptional. In particular, it seems analogous to that of other SB galaxies like NGC 3198 and NGC 3351.

Additionally, we have built a simple model for the Galactic ISM and demonstrated that a spin temperature of ~ 150 K is optimal to describe the typical H I opacity in the LAB dataset. Neglecting the H I opacity leads to underestimation of the H I mass within the inner Galaxy by ~ 30%.

We have found that the H I emission is better described by assuming two components, the first with high density (~ 0.4 cm^-3 on the midplane) and low-velocity dispersion (8 km s^-1), and the second with low density (~ 0.04 cm^-3 on the midplane) and high-velocity dispersion (15−20 km s^-1). We argue that the second component originates either from non-circular motions induced by non-axisymmetric features in the disc (like spiral arms), or from star formation feedback. In the former scenario, this component would serve the sole purpose of absorbing insufficiencies in our model, and therefore would be artificial.

We have stressed that the gas flow in a barred galaxy such as the Milky Way results in an asymmetric, non-stationary gas distribution and kinematics in the regions of the bar. For this reason, we discourage the reader to trust our rotation velocity measurements in the innermost ~ 2.5 kpc.

---

Acknowledgments

A.M. would like to thank Renzo Sancisi for critically reading this manuscript and for providing constructive comments, and Erwin de Blok for providing the rotation curves and the surface density profiles of THINGS galaxies. J.M.v.d.H. acknowledges support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 291531. The authors thank an anonymous referee for insightful comments.

References

Appendix A: The N_HI(v)−T_B relation in the Galaxy for a realistic ISM

In this appendix we study the relation between the intrinsic H I column density per unit velocity, N_HI(v), and the observed brightness temperature, T_B, for a light-cone that intersects the Galactic disc. The medium swept up by this light-cone is modelled as an ensemble of clouds with properties that are typical for the interstellar medium of the Milky Way. The brightness temperature profile is derived by assuming spatial and velocity resolutions analogous to those of the LAB survey.

In Appendix A.1 we describe the radiative transfer equations for a generic ensemble of H I clouds with different sizes, temperatures and densities. Next, in Appendices A.2 and A.3 we set up a simple model for the Galactic ISM that uses a limited number of free parameters to fully describe the cloud internal properties, overall distribution and kinematics for any generic sight-line. Given a set of free parameters, H I line profiles can be derived as shown in Appendix A.4. Finally, the N_HI(v)−T_B relation as resulting from our modelling is discussed in Appendix A.5.

Appendix A.1: Radiative transfer in a multiphase clumpy medium

As we discussed in Sect. 2.3, Eqs. (3) and (5) give the solution to the equation of radiative transfer and the optical depth for an homogenous and isothermal H I layer that uniformly fills the telescope beam.

We now consider the case of two overlapping layers of H I, one in the foreground (layer A) and one in the background (layer B), both filling the telescope beam. Given that the emission from layer B is partially absorbed by layer A, the T_B is not simply given by the sum of the two components. Instead we have (A.1)It is trivial to generalise Eq. (A.1) to a number N of overlapping layers. We consider the layers ordered along the line of sight such that layer i + 1 is located behind layer i (where 1 <i<N), then (A.2)We consider now the case of a spherical, homogeneous and isothermal cloud of H I with radius R_cl and H I mass M_cl. We assume that the sphere is unresolved, that is, R_cl<R_beam, where R_beam is the equivalent radius of the telescope beam (in physical units). Then brightness temperature will be given by the right-hand side of Eq. (3) diluted by a factor Ω = (R_cl/R_beam)². In the case of a collection of N clouds that do not shade each other along the line of sight, T_B will be (A.3)where τ_i(v), the optical depth of the ith cloud, is still given by Eq. (5), N_HI being the cloud surface density (thus ).

Finally, by using both Eqs. (A.2) and (A.3), it is possible to derive the T_B for an ensemble of N clouds that partially shade each other: (A.4)where F_i(v) is a correction term due to foreground absorption from clouds in front of cloud i (at velocity v), and is given by (A.5)where the first sum is extended to all m “shaded” sectors in which cloud i is partitioned by the foreground clouds (see example below), and δ_k is the area of sector k divided by the area of cloud i (so that ). We note that F ≤ 1 by construction.

The foreground correction depends on how the clouds are distributed in space. As an example, we compute F for the background cloud in a system of three clouds that are partially shading each other and have the same velocity (Fig. A.1). The cloud in the background (cloud 2) is divided into sectors A, B, C and D by the two foreground clouds. Sector A is unshaded, sectors B and C are shaded by clouds 1 and 0, respectively, while sector D is shaded by clouds 1 and 0 simultaneously. Thus the foreground correction for cloud 2 will be given by and δ_A + δ_B + δ_C + δ_D = 1.

We note that, for a generic distribution of clouds, computing the various δ_k terms in Eq. (A.5) is not trivial. In the calculations presented here, we found it more practical to randomly sample each cloud i with a number η of points, and study which foreground clouds shade the various regions of cloud i point-by-point. We adopt η = 250, which we found to be the best compromise between accuracy and computational speed.

Equation (A.4) can be used to predict any H I brightness temperature profile if a model for the cloud distribution in mass, radius, temperature, location and velocity is provided. We discuss such a model below.

Appendix A.2: A simple model for the neutral ISM

The model that we adopt is a simplified version of that of McKee & Ostriker (1977, hereafter MO77). We assume that the H I is made by a collection of homogeneous and isothermal clouds, in pressure equilibrium with each other at the ISM pressure of P_ISM/k_B = 3000 K cm^-3 (Wolfire et al. 2003). These clouds follow distributions in temperature, size, mass and velocity that we discuss below. To simplify the calculation, even though these clouds are distributed in a 3D space, we assume that they effectively behave as discs of constant surface density, and have no physical depth.

Each cloud consists of two components: an outer envelope of warm material with optically thin properties, and a cold optically thick core. In the MO77 model, the typical temperatures of these components are 8000 K for the warm envelope and 80 K for the cold core. However, Heiles & Troland (2003, hereafter HT03) studied in detail the interstellar H I emission and absorption lines in the solar neighbourhood using a Gaussian decomposition technique, finding a broad range of temperatures for both the warm and the cold neutral medium. The temperature distributions found by HT03 are shown in Fig. A.2. We note that spin temperatures can be measured only for the coldest (absorbing) H I clouds, whereas only an upper limit to the H I kinetic temperature T_K can be inferred for warmer clouds from the broadening of their emission profile. Spin and kinetic temperatures are virtually identical for T_S ≲ 1000 K, but T_K ≥ T_S at higher T_S (Liszt 2001). However, for T_S> 1000 K and for the column density range that we investigate here, the H I optical depth is sufficiently small that even overestimating T_S by thousands of K produces little variation in the line profiles. Therefore, in the following we make no distinction between T_S and T_K and assume that clouds in our model follow the temperature distributions shown in Fig. A.2.

In the MO77 model, two additional envelopes of warm and hot ionised gas surround each cloud, and the mass ratio between the different components is set by the detailed balance between the cloud physical properties and an external ionising UV and X-ray background. We do not include such a level of detail in our model. Instead, we assume that the mass fraction of the warm neutral gas, f_w = m_w/ (m_w + m_c), is constant for each cloud. HT03 found that the average warm-to-total H I ratio in the solar neighbourhood is ~ 0.6, remarking that this is likely a lower limit. The size of the warm clouds follows a power law distribution with a slope of − 3.6, according to the observations of Miville-Deschênes et al. (2003), who pointed out that this is in agreement with the Kolmogorov prediction for an incompressive turbulent fluid. Consistently with the MO77 model, we assume that the minimum and maximum radii for the warm clouds are R_min = 2.1 pc (below which clouds are easily disrupted by supernova shockwaves) and R_max = 10.8 pc (above which clouds are gravitationally unstable), respectively. Assuming a uniform size distribution in the same range has very little impact on our results. Given the steepness of the size distribution, the value of R_max does not affect our results; we attempted to use R_max = 100 pc and found no variations in the N_HI(v)−T_B relation for our “fiducial” model (see below). Our results are somewhat more sensitive to the minimum cloud size. For instance, using R_min = 1 pc gives a slightly more “optically thin” N_HI(v)−T_B relation. It is difficult to define precise boundaries for the cloud sizes, as they depend on the local condition of the ISM. Most importantly, we show below that the main driver of the N_HI(v)−T_B relation is the warm gas fraction, while all the other variables seem to play a secondary role. For these reasons, we decide to fix the cloud size boundaries to the values proposed by MO77.

Given that the temperature, pressure and size for a given warm cloud are set, then its mass m_w will be and its optical depth is . The choice of f_w sets the cold core mass m_c and then its size and optical depth. We verified that, with these choices for the temperature and size distributions, virtually all warm (cold) components have τ ≪ 1 (τ> 1).

Appendix A.3: Cloud distribution and kinematics

In order to derive a brightness temperature profile via Eq. (A.4), we must set up a virtual observation for our model of neutral interstellar medium described above. For this purpose we consider a synthetic light-cone that, radiating from an observer placed within the Galactic disc at R = R_⊙, intersects the ISM of the inner Galaxy at latitude b = 0. For simplicity, we consider a Galaxy model with a constant rotation velocity v_φ(R) = v_⊙ and a constant H I density n_HI.

At a chosen longitude l, the maximum line-of-sight distance within the inner disc is d_max(l) = 2 R_⊙cos(l) (so 0 <d<d_max) and, at this distance, the total volume encompassed by our synthetic light-cone will be , where α is the FWHM angular resolution of our virtual observation. We distribute our warm and cold clouds randomly within our light-cone until its mean density reaches the value of n_HI. The physical parameters of the clouds are extracted randomly from the distributions discussed in the previous section (see also Appendix A.4). A cloud at a given distance d along the light-cone will have line-of-sight velocity given by Eq. (1). For simplicity, we treat the cold core-warm envelope as two separate discs that share the same centre and velocity.

Then we need another step. If all H I in a cloud was at velocity v₀ given by Eq. (1), then the line profile would be a delta function. In reality, thermal and turbulent motions broaden the gas velocity distribution around v₀. In our calculation, what matters is how the H I optical depth of a given cloud is distributed in velocity, for which we assume (A.6)where τ₀ is the integrated optical depth of the cloud (proportional to the total surface density divided by T_S), is a normal distribution with standard deviation equal to σ, σ_th is the thermal velocity dispersion of the gas (equal to , where k_B is the Boltzmann constant and m_p is the proton mass) and σ_turb is an additional turbulent component. We have verified that a different prescription for turbulence, based on cloud-to-cloud relative motions rather than on the internal motions of clouds, gives results similar to the current implementation.

Appendix A.4: Building synthetic H I profiles

We briefly summarise the steps required to compute the H I brightness temperature profile for a given Galactic longitude l, and therefore to study the N_HI(v)−T_B relation.

First, we distribute H I clouds within a light-cone as follows:

• Given the longitude and the beam FWHM of the virtualobservation, we define the geometry of the light-cone withinwhich clouds will be randomly distributed.

• A warm neutral cloud is created at a random position within the light-cone. The cloud size and temperature are extracted randomly from the distributions discussed in Appendix A.2. Given that P_ISM is fixed, the cloud mass, surface density and optical depth are set univocally. The location of the cloud sets its line-of-sight velocity (given a flat rotation curve). The cloud temperature, along with the chosen internal turbulence σ_turb, sets the cloud velocity dispersion.

• A cold neutral cloud is created at the same location, with a temperature extracted randomly from the distribution discussed in Appendix A.2. The choice of f_w, the warm-to-total H I ratio, sets univocally all the other cloud parameters.

• The previous two steps are repeated until the H I volume density within the light-cone reaches the desired value, n_HI. The properties of each cloud are stored. The total number of clouds varies depending on the model parameters and on the longitude (for instance, the “fiducial” model shown in Appendix A.5 at l = 70° uses ~ 10⁴ clouds).

At this stage, an array N_HI(v_k) is built by using the surface densities, the line-of-sight velocities and the velocity dispersions of all clouds in the light-cone. We focus on the velocity range − 30 <v_k< 300 km s^-1 and adopt a velocity bin (v_{k + 1}−v_k) of 1 km s^-1.

Next, we compute the brightness temperature profile as the sum of the contributions of each individual cloud, as shown by Eq. (A.4). Each cloud i, however, contributes only to a limited portion of the profile, and specifically only around its line-of-sight velocity v_i. Therefore, for each cloud, we focus on the velocity range [v_i−3σ_i,v_i + 3σ_i], σ_i being the cloud velocity dispersion. The following steps must be repeated for all clouds in the light-cone and, for a given cloud, for all velocity bins v_k in the appropriate velocity range:

• We first identify a subsample of clouds that can potentiallycontribute to the foreground absorption. These clouds areselected as those that a) are located between cloud i and the observer, and b) have line-of-sight velocity that differs from v_k by less than three times their velocity dispersion.

• Cloud i is then randomly sampled by 250 points, which are used to compute F_i(v) via Eq. (A.5) (see Appendix A.1 for further details). The terms τ_j(v_k) in Eq. (A.5) are given by Eq. (A.6). We note that, in computing the overlaps between clouds, projection effects must be taken into account.

• The contribution of cloud i to the T_B at velocity v_k can be now computed. We note that, in this context, the quantity Ω_i in Eq. (A.4) represents the ratio between the angular size of cloud i and the beam FWHM.

Once all clouds have been processed, the resulting T_B profile is smoothed to a velocity resolution of 2 km s^-1.

As inputs, our model requires the warm-to-total H I ratio f_w, the H I volume density n_HI, the turbulent component of the velocity dispersion σ_turb and the beam FWHM. We fix the latter to 0.6°, as here we compare our synthetic H I profiles with those of the (hanning-smoothed) LAB survey at full (0.6°) angular resolution. We note that, for this value of the beam, the largest light-cone intersecting the inner Galaxy – which is obtained at l = 0° – has base radius (farthest away from the observer) of 87 pc, about half of the H I scale height at R = R_⊙. This implies that, in our model, we can neglect the vertical gradient in gas density.

Appendix A.5: Results

We now show how varying the parameters of our model affects the shape of the H I line profiles and the relation between the intrinsic H I column density per unit velocity and the brightness temperature. We define a “fiducial” model with parameters f_w = 0.6, n_HI = 0.5 cm^-3, σ_turb = 5 km s^-1 and study how the line profile and the N_HI(v)−T_B relation change by varying one parameter at a time. We focus on the latitude l = 70°, which gives a sufficiently large light-cone (5.7 kpc in length) to ensure a robust cloud statistics: in the tests below, the number of clouds ranges from 5000 to 22 000.

Figure A.3 shows the effect of varying the model parameters (f_w,n_HI,σ_turb) on the H I line profiles (top panels) and on the N_HI(v)−T_B relation (bottom panel). Remarkably, in all cases, the latter falls in between the predictions for two isothermal layers of H I at the spin temperature of 8000 K and 80 K, two values often quoted as representative for the temperatures of the warm and the cold phases of the atomic hydrogen. As shown by panels (a) and (d), increasing f_w boosts both the amplitude of the line profile and the brightness temperature at a given N_HI(v). We remark that 0.6 is likely to be a lower limit to the warm-to-total H I ratio (HT03). For f_w = 1, the N_HI(v)−T_B relation follows closely that of a single-phase warm (~ 10⁴ K) homogeneous medium. Although we do not aim to fit the observed line profile, we note that the models with f_w = 0.6 or 0.8 are those that better represent the LAB data. The gas properties are likely to vary as a function of the distance from the Galactic centre, thus fluctuations in the value of f_w across the disc are expected. However, the model with f_w = 1 not only fits the data poorly, but also largely exceeds the maximum brightness temperature measured in the LAB survey (152 K). This suggests that a cloud model for the Galactic ISM with no cold cores is inconsistent with the observations.

As shown by panels (b) and (e) of Fig. A.3, increasing or decreasing n_HI produces large variations in the amplitude of the line profile but, interestingly, does not affect dramatically the N_HI(v)−T_B relation. Increasing the cloud internal turbulence has the predictable effect of smearing the line profile, but has no impact on the N_HI(v)−T_B relation, as shown by panels (c) and (f). We note that the model with σ_turb = 5 km s^-1 is the one that better reproduces the tail in the observed line profile. This is consistent with our results based on fitting the LAB data (Sect. 3), as the inferred H I velocity dispersion never drops below 5 km s^-1 (even for the low-σ component alone in our two-component model, see Fig. 4).

One may ask whether or not the N_HI(v)−T_B relation depends significantly on the sight-line chosen. Unfortunately, the computational time increases dramatically by decreasing l (i.e. by increasing the light-cone length) and this makes the study of different sight-lines unpractical. Reducing α, the light-cone aperture, decreases the cone volume and allows us to study different sight-lines, but the resulting model cannot be considered as representative of the LAB data. In Fig. A.4 we compare the N_HI(v)−T_B relations derived at the longitudes varying between 10° and 70° for a model that has our fiducial parameters, but with α = 0.2°. Overall, the effect of varying the sight-line is not dramatic, as the scatter in the N_HI(v)−T_B relation visible in Fig. A.4 is comparable to that produced by varying n_HI (see panel (e) in Fig. A.3).

The above tests indicate that the warm-to-total H I ratio (f_w) is the main driver for the N_HI(v)−T_B relation, whereas the ISM density, the cloud turbulence and the adopted sight-line play only a secondary role. From our analysis, it seems the values of f_w between 0.6 and 0.8 are more consistent with the observations. In Fig. A.5 we show that Eq. (3) provides an excellent fit to the N_HI(v)−T_B relation (l = 70°, α = 0.6°) with f_w = 0.6 and 0.8. The resulting best-fit spin temperatures are 104 K and 176 K respectively, or 140 K as a mean. Remarkably, this spin temperature range brackets the maximum brightness temperature measured in the LAB data (152 K). Given that T_B ≤ T_S, in our calculation we decided to fix the H I spin temperature to the value of 152 K.

Appendix B: Testing the method

We test here the goodness of our modelling method (Sect. 2). We perform a number of tests based on building mock H I datacubes from axisymmetric, idealised gas distributions, and then running our fitting technique in an attempt to retrieve the original input parameters.

As a first test, we consider the simple case of a H I disc where all radial profiles for the various parameters are flat. We set v_φ = 240 km s^-1, h_s = 0.15 kpc, n_0,A = 0.45 cm^-2 cm^-3, n_0,B = 0.05 cm^-3, σ_0,B = 8 km s^-1 and σ_0,B = 13 km s^-1, where A and B represents two components with different densities and velocity dispersions but with the same rotation velocity and scale height (see Sect. 3.2). We build a mock H I observation of this model and we run our fitting method on the derived mockcube in order to verify whether we are able to recover all the input parameters. The leftmost panel of Fig. B.1 shows that this is indeed the case, as the fit values (points with error bars) match the original input values (solid lines) very well. In particular, the scale height is recovered with almost no uncertainties, whereas n₀ shows some scatter for the densest component.

In our second test, we consider a mock system built by two components: the same H I disc of the previous test, and an additional component representing a layer of slow rotating, inflowing extra-planar H I (see Sect. 5.7). We model the latter with the following density profile (from Oosterloo et al. 2007): (B.1)where (B.2)We assume R_g = 1.08 kpc, h₀ = 175 pc, h_R = 10.01 kpc, γ = 6.30, δ = 1.89 (see Marasco et al. 2013), and n₀ = 5 × 10^-5 cm^-3, so that the total H I mass of this component is ~ 3 × 10⁸M_⊙ (Marasco & Fraternali 2011). As for the kinematics, we assume that this component rotates slower with increasing height from the midplane, with a vertical lag of 15 km s^-1 kpc^-1 with respect to the H I in the disc, and that it is inflowing onto the disc with a speed of 30 km s^-1 and 20 km s^-1 in the radial and vertical direction, respectively (Marasco & Fraternali 2011). We build a mock observation of this model and we attempt to recover disc scale height without making any correction for the presence of the additional component. Our results, shown in the middle panel of Fig. B.1, suggest that in the presence of an extra-planar H I component, the disc scale height at R> 3 kpc is slightly overestimated with respect to the input value. The clear peak visible around R = 4 kpc (or l ~ 30°) is due to the fact that, at this longitude, the line profiles at | b | > 0° are largely contaminated by extra-planar H I emission. This can produce a spurious increase in H I surface density and total mass.

A possible way to correct for this effect is to include in our model the additional extra-planar component that we used to produce the mock observations in the first place, and proceed normally to derive the scale height of the H I disc. This procedure allows us to recover the input scale height profile much better, as shown in the rightmost panel of Fig. B.1. We use this procedure on the LAB data in Sect. 5.7.

Appendix C: Re-scaling the model to a different v_⊙

It is straightforward to re-scale our models to a different value of v_⊙, the rotation velocity of the solar circle, as it affects solely the rotation curve. Once the parameters of the model are set, any choice of v_⊙ would produce exactly the same mock datacube if the model rotation curve v_φ(R) is re-scaled as follows: (C.1)where v_φ′(R) is the new rotation curve and .

As an example, Fig. C.1 illustrates how the H I and H₂ Galactic rotation curve vary by using v_⊙ = 220 km s^-1 rather than v_⊙ = 240 km s^-1, the value used in the main text. We note that, after their computation, the rotation curves are smoothed to a resolution of 0.2 kpc and so the rotation velocity at R = R_⊙ may slightly differ from v_⊙. Interestingly, with this lower value of v_⊙ the rotation curve flattens at smaller radii (~ 4 kpc instead of 6.5 kpc). We reiterate that, for gas density, velocity dispersion and scale height profiles, such as those in Fig. 11, the model with this new rotation curve fits the data as well as the previous one.

All Tables

All Figures