© ESO 2020

1. Introduction

The molecular gas in galaxies is the place of birth of stars. Gravitational collapse and fragmentation of molecular clouds leads to the formation of stars (Toomre 1964; McKee & Ostriker 2007; Ostriker et al. 2010; Hopkins et al. 2011; Kim et al. 2011; Kennicutt & Evans 2012; Krumholz 2014; Agertz et al. 2015; Hayward & Hopkins 2017; Krumholz et al. 2017). In galaxies, molecular gas serves as the immediate precursor to star formation. In contrast, atomic gas plays the role of long-term supply for the star formation. The molecular gas in spiral galaxies dominates the inner parts of the interstellar medium (ISM), whereas the atomic gas dominates the outer regions. As the formation of stars must go through the molecular phase, this phase of the ISM is of immense importance in relation to galaxy formation and evolution (Myers et al. 1986; Shu et al. 1987; Scoville & Good 1989; Saintonge et al. 2013; Agertz & Kravtsov 2015). Not only that, but the abundance and the distribution of this component can also lead to a better understanding of the processes and conditions by which the gas is converted into stars. For example, it is well known that in spiral galaxies, the molecular gas surface density directly correlates to the star formation rate density (Kennicutt 1998; Bigiel et al. 2008; Leroy et al. 2009, 2013; Onodera et al. 2010; Schruba et al. 2010; Feldmann et al. 2011; Saintonge et al. 2011; Kruijssen & Longmore 2014; Kreckel et al. 2018); however, the same phenomenon in dwarf galaxies is not apparent (Schruba et al. 2012). Hence, investigating the structure and the distribution of the molecular gas in galaxies could lead to a better understanding of the interplay between gas and star formation.

Estimating the three-dimensional distribution of the molecular and the atomic gas in external galaxies is difficult because of line-of-sight projection effects (Patra 2018, 2019). In the Milky Way, however, it is possible to estimate the three-dimensional distribution of the molecular gas by locating the positions of the molecular clouds and remapping them into a common three-dimensional grid (Grabelsky et al. 1987; Solomon et al. 1987; Wouterloot et al. 1990; Rathborne et al. 2009; Roman-Duval et al. 2010, 2016). However, the estimated distance to the clouds is very often found to be ambiguous as the kinematic distances introduce significant fractional errors (e.g. Sanders et al. 1984; Scoville et al. 1987; Bronfman et al. 1988).

To overcome these errors, there have been several remarkable studies in recent years to explore the molecular gas in external galaxies. These studies have mostly focused on the following two aspects of the molecular discs: firstly, the global distribution of the molecular gas and its connection to the star formation (see, e.g. Leroy et al. 2009; Schruba et al. 2012; Utomo et al. 2018; Sun et al. 2018; Schinnerer et al. 2019); and secondly, the phases of the molecular gas (Stark 1984; Scoville 1990; Combes & Becquaert 1997; Caldú-Primo et al. 2013; Mogotsi et al. 2016; Caldú-Primo & Schruba 2016; Roman-Duval et al. 2016). The studies of the spectral properties in connection to the phases in the molecular discs reveal that the molecular gas in galaxies might exist in the following two phases: one with a low kinetic temperature and hence a narrow spectral width (∼7 km s⁻¹); the other with a larger spectral width, which is comparable to that of the atomic gas (∼12 km s⁻¹) (Wilson & Scoville 1990; Wilson & Walker 1994; Caldú-Primo et al. 2013; Mogotsi et al. 2016). This high spectral width component is much more diffused and not easily detectable in the individual spectrum of localised CO emission, though it has a significant share in the total molecular mass (see, e.g. Goldsmith et al. 2008; Liszt et al. 2010; Pety et al. 2013; Roman-Duval et al. 2016). It should be mentioned here that these velocity widths refer to a much higher kinetic temperature than what is expected in the molecular clouds and hence, these widths indicate the turbulent motions of the clouds rather than their pure thermal widths.

Interestingly, the diffuse component of the molecular disc found to have a velocity dispersion similar to that of the atomic gas (H I) indicates that they might occupy the same dynamical phase space with similar scale heights (Caldú-Primo et al. 2013; Mogotsi et al. 2016; Caldú-Primo & Schruba 2016). This demands a theoretical understanding of the origin and the sustenance of this diffuse molecular component in environments where it is much harder to protect the molecular gas via shelf-shielding. An explanation for the existence of the diffuse molecular gas arises in the context of the presence of the molecular gas in the arm and the inter-arm regions of spiral galaxies. Observationally, it has been found that the giant molecular clouds (GMCs) are detected in higher numbers in the arm regions than the inter-arm regions in spiral galaxies (Dame et al. 1986; Hunter et al. 1997; Sawada et al. 2004; Nakanishi & Sofue 2006). To explain this, two different propositions are suggested. Firstly, it is proposed that, during the passage of the spiral density wave, there is an accelerated phase transition that converts the atomic gas in the inter-arm region to molecular gas in the arm very rapidly. The molecular gas again gets converted into atomic and ionised gas when the density wave gets passed through the inter-arm region (Blitz & Shu 1980; Cohen et al. 1980; Scoville & Hersh 1979; Heyer & Dame 2015). In the second scenario, it is proposed that there is a little phase transition. Instead, in the inter-arm regions (or in the low-density environment), there already exist tiny molecular clouds which can self-shield themselves. The existence of the GMCs in the arm region is then attributed to the dynamically driven coagulation and fragmentation of these tiny molecular clouds in the arms and the inter-arms, respectively (Scoville & Hersh 1979; Vogel et al. 1988; Hartmann et al. 2001; Kawamura et al. 2009; Miura et al. 2012). This, in turn, generates large molecular clouds in the arm region by the merger of smaller clouds, and the more massive clouds get fragmented dynamically once the spiral density wave passes through (Koda et al. 2016). For the first scenario to be true, no CO emission is expected to be observed in the inter-arm regions, which the earliest observations with insufficient sensitivity confirmed (Engargiola et al. 2003; Heyer et al. 2004; Fukui et al. 2009; Vogel et al. 1988; Rand & Kulkarni 1990). However, recent sensitive observations have revealed a significant amount of CO in the inter-arm regions as well (Liszt et al. 2010; Koda et al. 2011; Pety et al. 2013).

The existence of the molecular gas at large heights from the midplane can be easily explained under the assumption of the second scenario. The molecular gas does not get destroyed at these low-density environments; instead, it is protected within small clouds. In that case, the observed molecular spectra in these regions is decided by the turbulent motions of these clouds, which is much higher than the turbulent motions of the GMCs. This, in turn, results in a two-component molecular disc with two different molecular velocity dispersions. Not only that, but as this component is made up of tiny molecular clouds, it is diffused in nature and would be easily resolved out in high-resolution interferometric observations.

Several recent observations, in fact, strongly suggest the existence of this diffuse low-density molecular gas, which has a vertical thickness comparable to that of the atomic gas. For example, Garcia-Burillo et al. (1992) found a thick molecular disc in the edge-on galaxy NGC 891, which extends up to ∼1.4 kpc from the midplane. The existence of a thick molecular disc in this galaxy was further strengthened by the sensitive observations of the cold dust by Bocchio et al. (2016) using Spitzer telescope. They detected cold dust emission in NGC 891 with a vertical scale height of ∼1.44 kpc (see also, Guelin et al. 1993; Rejkuba et al. 2009). On the other hand, for the Galaxy, Dame & Thaddeus (1994) used the 1.2 m telescope of the Harvard-Smithsonian Center for Astrophysics to conduct a very sensitive survey of diffuse molecular gas in the inner Galaxy beyond 100 pc from the midplane. They found a faint, thick molecular disc approximately three times as wide as that of the central dense CO layer and which is comparable to the width of the HI layer. Likewise, many other spectral studies have found a similar velocity dispersion in the molecular discs, as observed in the atomic discs. For example, Combes & Becquaert (1997) measured the velocity dispersions in two nearly face-on galaxies, NGC 628 and NGC 3938, to find a very similar velocity dispersion of the molecular and the atomic gas. Using the data from The HI Nearby Galaxy Survey (THINGS; Walter et al. 2008) and the HERA CO Line Extra-galactic Survey (HERACLES; Leroy et al. 2009), Caldú-Primo et al. (2013) stacked the line-of-sight of HI and CO spectra to produce high signal to noise ratio (S/N) stacked spectra. Using these high S/N spectra, they also observed similar velocity dispersions in both the molecular and the atomic discs (see, Caldú-Primo & Schruba 2016, also for a similar analysis in Andromeda galaxy). Later, Mogotsi et al. (2016) used the same data to study the HI and the CO spectra of bright individual regions to find a somewhat lower velocity dispersion in the molecular disc than what is observed in the atomic disc. This is consistent with the scenario that diffuse high dispersion molecular gas cannot be detected in a low S/N individual spectrum. Pety et al. (2013) studied the molecular gas in M 51 using both single-dish and interferometric observations and found that the single-dish observation recovers much more flux than the interferometric observation. They conclude that at least 50% of gas in the molecular disc of M 51 is in the diffused state, which is resolved out in the interferometric observation. A similar fraction of the diffuse molecular gas in other galaxies was also found in several other studies (see, e.g. Goldsmith et al. 2008; Liszt et al. 2010; Caldú-Primo et al. 2015; Roman-Duval et al. 2016).

With these observational inferences and theoretical propositions, it is worth investigating in detail the observational aspects due to the presence of two-component molecular discs in galaxies and, in particular, the impact of the diffused molecular component on the observed thickness and the spectral properties of the molecular discs. To do that, we modelled the galactic disc as a four-component system with a stellar disc, an atomic disc, and a two-component molecular disc in hydrostatic equilibrium under their mutual gravity in the external force field of the dark matter halo. Under this assumption, the density distributions of different baryonic components can be described by the joint Poisson’s-Boltzmann equation, which upon solving, would provide a detailed three-dimensional distribution of the molecular and atomic gas. Furthermore, these density distributions can be used to build dynamical models of the molecular and atomic discs, which again can be used to produce the column density maps and the spectral cubes. These simulated maps and cubes can then be studied to understand the observational signatures of a two-component molecular disc as compared to a single-component one. In this paper, we formulate the joint Poisson’s-Boltzmann equation for a four-component galactic system and solve it numerically to understand the distribution of the molecular gas and its observational consequences in a sample of eight nearby spiral galaxies from the HERACLES and THINGS survey.

2. Sample

We chose our sample galaxies from the parent sample of the HERACLES survey. As part of the HERACLES survey, 18 galaxies were observed in CO using the 30 m IRAM telescope. These galaxies were chosen to be a part of the THINGS (for available HI data) and Spitzer Infrared Nearby Galaxies Survey (SINGS; Kennicutt et al. 2003, for available stellar surface density data). In hydrostatic equilibrium, the gravity in the vertical direction is balanced by the pressure. Hence, the surface densities of different baryonic discs are essential inputs to the hydrostatic equation. Not only that, but the dark matter halo also provides a significant amount of gravity, which decides the vertical distribution of gas under hydrostatic equilibrium. Given these circumstances, the surface density of different baryonic discs and the mass model, which can be used to obtain the dark matter density distribution, for a galaxy are primary requirements in solving the hydrostatic equation. However, all of the 18 galaxies observed in the HERACLES survey did not have the mass models due to complications in extracting their rotation curves (e.g. low and high inclination, the presence of strong non-circular motion in HI data, etc.) (see de Blok et al. 2008, for more details). These constraints restrict ten galaxies for which hydrostatic modelling could not be performed. We hence solved the hydrostatic equation for the eight remaining galaxies for which all the necessary data are available. In Table 1, we list the basic properties of our sample galaxies. In Col. (1), we list the names of our sample galaxies, Col. (2) shows the distance to the galaxy in megaparsec, Cols. (3) and (4) represent the inclination and the position angle of the optical discs, respectively, and Col. (5) shows the optical radius, r₂₅, the 25th magnitude B-band isophote. Column (6) represents the radial extent of the molecular discs in our galaxies as obtained by Schruba et al. (2011, see their Fig. A.1).

3. Modelling the galactic discs

3.1. Setting up the hydrostatic equation

We assume the galactic disc to be a four-component system that consists of a stellar disc, an atomic disc, and a two-component molecular disc. The molecular disc considered here has the following two components: one with high velocity dispersion (thick disc) and diffuse gas that has large scale heights; the other with low velocity dispersion (thin disc), which is settled close to the midplane. All of these baryonic discs are considered to be in hydrostatic equilibrium under their mutual gravity in the external force field of the dark matter halo. For simplicity, we assume that these baryonic discs are concentric, coplanar, and azimuthally symmetric. Under these assumptions, the joint Poisson’s equation for an elemental mass can be written in cylindrical polar coordinate as

$$\begin{array}{c}\hfill \frac{1}{R}\frac{\partial }{\partial R}\left(R\frac{\partial {\mathrm{\Phi }}_{\mathrm{total}}}{\partial R}\right)+\frac{{\partial }^2{\mathrm{\Phi }}_{\mathrm{total}}}{\partial z^2}=4\pi G{\displaystyle (\sum _{i=1}^4{\rho }_i+{\rho }_{\mathrm{h}}),}\end{array}$$(1)

where Φ_total is the total gravitational potential due to all of the disc components and the dark matter halo. We note that ρ_i denotes the volume density, where i runs for stars, atomic gas, thin disc molecular gas, and thick disc molecular gas. Additionally, ρ_h denotes the mass density of the dark matter halo, which can be described by either an NFW profile (Navarro et al. 1997)

$$\begin{array}{c}\hfill {\rho }_{\mathrm{h}}(R)=\frac{{\rho }_0}{\frac{R}{R_{\mathrm{s}}}{(1+\frac{R}{R_{\mathrm{s}}})}^2}\end{array}$$(2)

or by an isothermal profile

$$\begin{array}{c}\hfill {\rho }_{\mathrm{h}}(R)=\frac{{\rho }_0}{1+{\left(\frac{R}{R_{\mathrm{s}}}\right)}^2},\end{array}$$(3)

where ρ₀ denotes the characteristic density, and R_s represents the scale radius. These two halo parameters completely describe a spherically symmetric dark matter halo.

For an individual component in hydrostatic equilibrium, the gradient in pressure (in the z direction) is balanced by the gradient in potential. Thus,

$$\begin{array}{c}\hfill \frac{\partial }{\partial z}\left({\rho }_i{⟨{\sigma }_z^2⟩}_i\right)+{\rho }_i\frac{\partial {\mathrm{\Phi }}_{\mathrm{total}}}{\partial z}=0,\end{array}$$(4)

where ⟨σ_z⟩_i is the vertical velocity dispersion of the ith component, which is an input parameter.

By eliminating Φ_total from Eqs. (1) and (4), we get

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {⟨{\sigma }_z^2⟩}_i\frac{\partial }{\partial z}\left(\frac{1}{{\rho }_i}\frac{\partial {\rho }_i}{\partial z}\right)=& -4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{HI}}+{\rho }_{\mathrm{H}2}+{\rho }_{\mathrm{H}2}^{\prime }+{\rho }_{\mathrm{h}})\hfill \\ \hfill & +\frac{1}{R}\frac{\partial }{\partial R}\left(R\frac{\partial {\mathrm{\Phi }}_{\mathrm{total}}}{\partial R}\right),\hfill \end{array}\end{array}$$(5)

where ρ_s, ρ_HI, ρ_H2, and ${\rho }_{\mathrm{H}2}\prime$ are the mass densities of stars, HI, molecular gas in the thin disc, and molecular gas in the thick disc, respectively.

The above equation can be further simplified using

$$\begin{array}{c}\hfill {\left(R\frac{\partial {\mathrm{\Phi }}_{\mathrm{total}}}{\partial R}\right)}_{R,z}={(v_{\mathrm{rot}}^2)}_{R,z}.\end{array}$$(6)

Assuming that the vertical gradient in (v_rot)_R, z is small, one can approximate (v_rot)_R, z by the observed rotation curve v_rot, which is a function of R alone. Thus Eq. (5) reduces to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {⟨{\sigma }_z^2⟩}_i\frac{\partial }{\partial z}\left(\frac{1}{{\rho }_i}\frac{\partial {\rho }_i}{\partial z}\right)=& -4\pi G({\rho }_{\mathrm{s}}+{\rho }_{\mathrm{HI}}+{\rho }_{\mathrm{H}2}+{\rho }_{\mathrm{H}2}^{\prime }+{\rho }_{\mathrm{h}})\hfill \\ \hfill & +\frac{1}{R}\frac{\partial }{\partial R}\left(v_{\mathrm{rot}}^2\right).\hfill \end{array}\end{array}$$(7)

Equation (7) represents four coupled, second-order ordinary partial differential equations in the variables ρ_s, ρ_HI, ρ_H2, and ${\rho }_{\mathrm{H}2}\prime$. The solutions of these equations at any radius R, provide the density distributions as a function of z. Thus solving these equations at all radii renders a detailed three-dimensional density distribution of different baryonic discs.

3.2. Input parameters

In vertical hydrostatic equilibrium, the gravity produced by the baryonic discs and the dark matter halo is balanced by the pressure. As the baryons in the galactic discs provide a significant amount of gravity, the surface density profiles of different baryonic components are one of the primary inputs to Eq. (7). We obtained the stellar surface density profiles for our sample galaxies as derived by Leroy et al. (2009) using the 3.6 μm data from the SINGS survey. We used the atomic and the molecular surface density profiles from Schruba et al. (2011). Schruba et al. (2011) derived the HI and the molecular surface density profiles using the THINGS and HERACLES survey data, respectively. The atomic and the molecular surface density profiles were determined by averaging the HI and the CO spectra, respectively, in the tilted rings of width 15″ as extracted by the process of the rotation curve fitting (see de Blok et al. 2008, for more details). To enhance the sensitivity of the molecular surface density profiles, the CO spectra within a tilted ring were co-added after shifting their central velocities to zero. These co-added and stacked spectra have a much higher S/N than the individual spectrum, which is crucial for the inclusion of the low-density thick disc molecular gas. In Fig. 1, we plotted the surface density profiles for different baryonic discs in our sample galaxies. As can be seen from the figure, for all our sample galaxies, the stellar surface density dominates the atomic or the molecular surface densities, though, the gas discs tend to extend to a larger radius than the stellar disc.

Fig. 1. Surface density profiles for our sample galaxies. Different panels indicate the surface densities for different galaxies, as mentioned in the top right corner of each panel. The cyan pluses, the red triangles, and the blue circles represent the stellar, the atomic, and the molecular surface density profiles, respectively. The stellar surface density profiles are taken from Leroy et al. (2009), whereas the atomic and the molecular surface density profiles are obtained from Schruba et al. (2011).

The vertical velocity dispersions of different baryonic components are the other important inputs to Eq. (7). They provide the necessary pressure support to balance gravity. As the vertical velocity dispersion of the stellar disc does not significantly influence the vertical scale height (defined as the half width at half maxima) of the gas discs (Banerjee et al. 2011), we calculated the stellar velocity dispersion analytically assuming the stellar disc to be an isothermal single component system in hydrostatic equilibrium (see Leroy et al. 2008, for more details). The velocity dispersions in the gaseous discs, on the other hand, can be estimated using spectroscopic observations with radio telescopes. For example, early low-resolution HI observations resulted in an HI velocity dispersion (σ_HI) between 6 and 13 km s⁻¹ in nearby spiral galaxies (Shostak & van der Kruit 1984; van der Kruit & Shostak 1984; Kamphuis & Sancisi 1993). In an extensive analysis of the second moments of the HI spectral cubes of THINGS galaxies, Tamburro et al. (2009) found a mean velocity dispersion of ∼10 km s⁻¹ at the optical radius, r₂₅. Later, Ianjamasimanana et al. (2012) used the same sample to stack line-of-sight HI spectra and found a σ_HI = 12.5 ± 3.5 km s⁻¹ (σ_HI = 10.9 ± 2.1 km s⁻¹ for galaxies with inclination less than 60°). Caldú-Primo et al. (2013) studied the σ_HI and σ_CO in a sample of 12 nearby galaxies from the THINGS survey and found a σ_HI/σ_CO = 1.0 ± 0.2 with a median of σ_HI = 11.9 ± 3.1 km s⁻¹. However, Mogotsi et al. (2016) used the same sample to study the high S/N individual spectra and found a σ_HI/σ_CO = 1.4 ± 0.2 with σ_HI = 11.7 ± 2.3 km s⁻¹. As can be seen from these studies, a σ_HI ∼ 12 km s⁻¹ is reasonable to assume in HI discs of our sample galaxies.

Stark (1984) observed molecular clouds in the Galaxy and found that the velocity dispersions of the low mass clouds are higher than that of the high mass clouds. The low mass clouds have a σ_CO ∼ 9.0 km s⁻¹, whereas the high mass clouds have a σ_CO ∼ 6.6 km s⁻¹, which is also reconfirmed by their subsequent observations (Stark & Lee 2005). In investigating the spectral cubes of the THINGS and the HERACLES survey data, Caldú-Primo et al. (2013) and Mogotsi et al. (2016) conclude that due to its diffuse and faint nature, the thick disc molecular gas could only be detected in high S/N stacked spectra. They found a σ_CO ∼ 7.3 ± 1.7 km s⁻¹ for the thin disc molecular gas. The σ_CO for the thick disc, on the other hand, is expected to be the same as the σ_HI.

Given these results, for our sample galaxies, we assume velocity dispersions of ∼7 km s⁻¹ and ∼12 km s⁻¹ in the thin and the thick molecular discs, respectively. The fraction of the total molecular gas in the thick disc is not well constrained by observations. Pety et al. (2013) estimated that at least ∼50% of the molecular gas in M 51 is in the thick disc. Here, in this analysis, we assume that the molecular gas is equally divided into the thin and the thick disc. However, we also explore the effect of a different assumed thick disc molecular gas fraction on our results.

The next input parameter is the rotation curve. The last term on the RHS of Eq. (7) represents the radial term, which might influence the solution of Eq. (7) considerably in the inner parts of a galaxy where the rotation curve is not flat. This term can be calculated using the knowledge of the observed rotation curve. In Fig. 2, we plotted the rotation curves of our sample galaxies as extracted by de Blok et al. (2008) using the THINGS survey data. As we use a second-order derivative of the rotation curve in Eq. (7), we parameterised it to prevent any possible divergence while solving the hydrostatic equation. We employed a commonly used Brandt’s profile to represent the rotation curves. The Brandt’s profile can be given as (Brandt 1960)

$$\begin{array}{c}\hfill v_{\mathrm{rot}}(R)=\frac{V_{\max }(R/R_{\max })}{{(1/3+2/3{\left(\frac{R}{R_{\max }}\right)}^n)}^{3/2n}},\end{array}$$(8)

Fig. 2. Rotation curves of our sample galaxies as derived using the HI data from the THINGS survey by de Blok et al. (2008). The rotation curves were fitted using a Brandt profile (Eq. (8)) for parametrisation. In each panel, the black circles with error bars represent the rotation curve, whereas the red solid line represents the Brandt profile fit. The best-fit parameters can be found in Table 2.

where V_max is the maximum rotational velocity and R_max is the radius at which this maximum velocity is attained. Additionally, n represents the power-law index, which decides how fast the rotation curve rises as a function of the radius. In Table 2, we show the best fit parameter values for the rotation curves of our sample galaxies.

The last input required to solve the hydrostatic equation is the density distribution of the dark matter halo. The dark matter halo provides a considerable amount of gravity, which influences the hydrostatic structure in the vertical direction, especially at the outer regions where the dark matter dominates the gravity. We used the dark matter distributions from the mass models of our sample galaxies as derived by de Blok et al. (2008). For our sample galaxies, neither an NFW profile nor an ISO profile describes all the galaxies systematically better than the other. Hence, for a particular galaxy, we selected a dark matter profile that describes the observed rotation curve better (see de Blok et al. 2008, for more details). In Table 3, we present the dark matter halo types and their profile parameters, which we used to solve Eq. (7) for our sample galaxies.

3.3. Solving the hydrostatic equilibrium equation

With these inputs, Eq. (7) can be solved to obtain detailed three-dimensional density distributions of different baryonic discs in our sample galaxies. A comprehensive description on how to solve Eq. (7) assuming the galactic disc to be a three-component system (stars + atomic gas + single-component molecular gas) can be found in Patra (2018, 2019). Here we used a similar approach, but with one extra component in the molecular disc. In this subsection, we briefly illustrate the scheme we used to solve a four-component galactic disc. Equation (7) refers to four second-order partial differential equations in stars, atomic gas, and two molecular gas components coupled through gravity-term in the RHS.

As these equations are second-order differential equations, one needs at least two initial conditions to solve them. We chose the following conditions to solve our equations:

$$\begin{array}{c}\hfill {\left({\rho }_i\right)}_{z=0}={\rho }_{i,0}\mathrm{and}{\left(\frac{\mathrm{d}{\rho }_i}{\mathrm{d}z}\right)}_{z=0}=0.\end{array}$$(9)

These two conditions signify that the volume densities of individual baryonic discs and their derivatives (along z) at the midplane are known. The second condition is a result of the force balance in the midplane (Spitzer 1942). Due to symmetry in the vertical direction, the gravitational potential at the midplane is expected to be maximum. This, in turn, would result in a density maxima. Hence, the first derivative of the densities, for all the baryonic discs, would be zero at z = 0. On the other hand, the first condition is not generic and assumes that the midplane density is known. However, it is not a directly measurable quantity, though, it can be estimated using the knowledge of the observed surface density, Σ_i. Since it is an unknown quantity, for an individual component, we first assumed a trial value, ρ_0, i, t, and solved Eq. (7) to obtain a trial solution, ρ_i, t(z). This trial solution was then integrated to calculate the trial surface density, Σ_i = ∫ρ_i(z)dz, and compared with the observed value, Σ_i, o, to update the trial midplane density in the next iteration. We repeated this exercise until we found a trial midplane density, which produces a surface density matched to the observed value with better than 0.1% accuracy.

Using the scheme mentioned above, Eq. (7) can be solved for a disc without any coupling, that is, in the absence of any other baryonic discs. However, as Eq. (7) refers to a coupled equation, all of the components should be solved simultaneously. We implemented this using an iterative approach similar to what has been used in many earlier studies (Banerjee & Jog 2008; Patra et al. 2014; Patra 2018, 2019). In the first iteration, we solved for individual components assuming no coupling (i.e. solved separately). In this case, the densities of the other components are assumed to be zero. For example, if we are solving for stars (ρ_s), then ρ_HI, ρ_H2, and ${\rho }_{\text{H2}}^{\prime }$ are considered to be zero. Thus at the end of the first iteration, we obtained the individual solutions, that is, ρ_s, 1(z), ρ_HI, 1(z), ρ_H2, 1(z), and ${\rho }_{\text{H2},\text{1}}^{\prime }(z)$ (the subscript “1” indicates the iteration number). In the next iteration, we solved individual components in a similar way, but this time in the presence of the other components, that is, using the density solutions from the previous iteration. For example, when we solved for stars in the second iteration, instead of considering ${\rho }_{\text{HI}}={\rho }_{\text{H2}}={\rho }_{\text{H2}}^{\prime }=0$, we used ρ_HI = ρ_HI, 1, ρ_H2 = ρ_H2, 1, and ${\rho }_{\text{H2}}^{\prime }={\rho }_{\text{H2},\text{1}}^{\prime }$ and produce ρ_s, 2. This ρ_s, 2 is a slightly better coupled-solution than ρ_s, 1. ρ_HI, 2, ρ_H2, 2, and ${\rho }_{\text{H2},\text{2}}^{\prime }$, which were also obtained in a similar way. We continued this iterative process until the density solutions for all the components converge with better than 0.1% accuracy. Thus, in this scheme, we iteratively approach towards a coupled solution starting from a single-disc one. We note that for our sample galaxies, the solutions quickly converge within a couple of iterations.

It can be seen from Fig. 1 that the molecular discs in our sample galaxies are shorter than the stellar or the atomic discs. Hence, we solved Eq. (7) up to a radius where the CO is detected with confidence (see Table 1). The central regions of spiral galaxies are very often found to host highly energetic activities (e.g. enhanced star formation, outflow, etc.), which might disturb the hydrostatic condition. The dark matter density of an NFW halo also rises very sharply in the central region. Due to these reasons, we did not solve Eq. (7) in a region, R <  1 kpc, to avoid any possible divergences. The atomic discs in our sample galaxies extend much further than the molecular discs. At these radii (beyond the molecular disc), the galactic discs reduce to two-component systems (stars + atomic gas), and subsequently Eq. (7) should be modified before solving. As Eq. (7) cannot be solved analytically, we solved it numerically using an eighth-order Runge-Kutta method as implemented in the python package scipy. We solved Eq. (7) in our sample galaxies at every 100 pc radial distance. As the mean resolution of the HERACLES survey is ∼500 pc, a resolution of 100 pc is expected to be good enough to capture the radial variation of the molecular surface density. However, in the vertical direction, one needs a much higher spatial resolution to capture the rapid variation of the molecular gas density, which varies in 10 s of the parsec-scale. To achieve this, we employed an adaptive resolution in the vertical direction and always find it to be better than a few parsecs for our sample galaxies. It should be noted that as we investigate the observational signatures of a two-component molecular disc against a single-component one, we solved the hydrostatic equation for our sample galaxies assuming their molecular discs to be both a single-component and a two-component system.

4. Results and discussion

In Fig. 3, we plotted sample solutions of the hydrostatic equilibrium equation for NGC 7331 at a radius of 8 kpc. The top panel shows the solutions for a three-component galactic disc, whereas the bottom panel shows the same for a four-component disc. Ideally, in the absence of coupling altogether, the solution of the hydrostatic equation for an isothermal disc should follow a sech² law (see, e.g. Bahcall 1984a,b). However, we find that due to coupling between the baryonic discs and the dark matter halo, the solution of Eq. (7) deviates from a sech² law and behaves more like a Gaussian. It can be seen from the figure that the stellar disc extends to a much larger height than the gaseous discs in both cases. As compared to a three-component disc, in a four-component disc, there is a lower amount of low-velocity dispersion molecular gas, settling close to the mid-plane. This effectively reduces the total gravity near the midplane, resulting in thicker atomic and molecular discs in a four-component system (bottom panel) than what is seen in a three-component system (top panel). The molecular gas in the thick disc can be seen to extend to much larger heights than the molecular gas in a thin disc. This height is often comparable to that of the atomic disc (magenta dotted and blue-dashed lines in the bottom panel), leading to a hypothesis that the thick disc molecular gas and the atomic gas are part of the same dynamical component. Here, we solved the hydrostatic equation assuming an equal amount of gas in the thin and thick disc of the two-component molecular disc. For comparison purposes, we also solved Eq. (7) for a thick disc molecular gas fraction of 0.7, and the respective solutions for the molecular gas are plotted in the bottom panel of Fig. 3 (thick dashed-dotted black and dashed magenta lines). As can be seen from the figure, an increased fraction of the thick disc molecular gas fraction leads to a thicker molecular disc, thus increasing its effective scale height.

Fig. 3. Sample solutions of Eq. (7) for NGC 7331 at R = 8 kpc. Top panel: density solutions for a three-component galactic disc having a single-component molecular disc. Bottom panel: same for a four-component disc. In this case, the molecular disc is assumed to be a two-component system with a thin and a thick disc. The red dotted and the solid blue lines in both of the panels represent the volume densities of the stellar and atomic discs, respectively. The black dashed line in the top panel represents the volume density of the thin disc molecular gas in a three-component disc. Bottom panel: the thin dashed-dotted black and dashed magenta lines represent the volume density of molecular gas in the thin and thick discs, respectively, for an assumed thick disc gas fraction of 0.5. Whereas, the thick dashed-dotted black and the dashed magenta lines represent the same for an assumed thick disc molecular gas fraction of 0.7. See the text for more details.

Density solutions can often provide meaningful insights into several physical processes. For example, by considering a hydrostatic equilibrium in 12 nearby star-forming disc galaxies, Bacchini et al. (2019a,b) show that the star formation law exhibits much less scatter when it is determined using volumetric gas density instead of gas surface density. Not only that, but using the density solutions and the observed velocity dispersion, one can estimate the pressure at the midplane, which is an important parameter when deciding if the gas participates in star formation or not (see, e.g. Usero et al. 2015; Bigiel et al. 2016; Helfer & Blitz 1997). However, many of these studies do not solve the hydrostatic equilibrium equation self-consistently to estimate the relevant physical quantities. For example, Bacchini et al. (2019a,b) calculated the volume density of the gas discs by considering a preassumed vertical profile for the stellar (exponential) and gas (Gaussian) distribution.

Similarly, Gallagher et al. (2018) estimated the midplane pressure in a sample of nearby, large galaxies by ignoring the dark matter halo and assuming a constant stellar scale height. They calculated the midplane pressure by summing the gas self-gravity and the weight of the gas in the potential of the stellar disc by assuming an exponential vertical profile. The midplane pressure can be given as

$$\begin{array}{c}\hfill P_{\mathrm{mid}}\approx \frac{\pi G{\mathrm{\Sigma }}_{\mathrm{gas}}^2}{2}+{\mathrm{\Sigma }}_{\mathrm{gas}}{\left(2G{\rho }_{\ast }\right)}^{1/2}{\sigma }_{\mathrm{gas},z},\end{array}$$(10)

where Σ_gas is the total gas surface density (atomic + molecular), ρ_* is the volumetric mass density of stars and dark matter together, and the σ_gas, z is the vertical velocity dispersion of the gas. Gallagher et al. (2018) ignored the dark matter contribution to ρ_*, arguing that the total mass of the dark matter within the stellar layer at any radius would be negligible. Thus, they estimated the stellar density at the midplane assuming a known scale height, h_*, given by (van der Kruit 1988)

$$\begin{array}{c}\hfill {\rho }_{\ast }\approx \frac{{\mathrm{\Sigma }}_{\ast }}{2h_{\ast }},\end{array}$$(11)

where Σ_* is the stellar surface density. They further assumed a constant scale height for the stellar disc, as computed by Leroy et al. (2008), using the knowledge of the flattening ratio of the stellar discs in galaxies (Kregel et al. 2002). The scale height can be given as l_*/h_* = 7.3 ± 2.2. Where l_* represents the exponential scale length of the stellar disc. It should be noted here that the observed flattening ratio might not capture the true flaring of the stellar disc as the observations are only sensitive to the high luminosity stellar component. In that case, this prescription can considerably underestimate the stellar scale height and hence overestimate the ρ_*.

To examine the impact of this assumed constant stellar scale height and the gravity of the dark matter halo, we solved Eq. (7) without dark matter and calculated the gas volume densities and the midplane pressure. In the top panel of Fig. 4, we show the effect of the dark matter halo on the estimated vertical density distribution of the molecular gas. For simplicity, here, we solved Eq. (7) assuming a single-component molecular disc; however, the conclusions also remain the same for a two-component molecular disc. We solved the hydrostatic equilibrium equation for two galaxies from our sample, that is, NGC 3521 (thick lines) and NGC 7331 (thin lines). These two galaxies represent two extreme ends of the dark matter dominance in our sample (see, e.g. Table 3). As can be seen from the top panel of the figure, the presence of the dark matter halo significantly alters the vertical density profiles. Ignoring the dark matter would result in an underestimation of the volume density close to the midplane and an overestimation of it at more considerable heights. The effect of the dark matter halo is expected to be more pronounced at outer radii, where the dark matter dominates the gravity. At these radii, it can play a significant role in regulating star formation by modulating the volume density and/or by providing the mass required to produce gravitational instabilities (see, e.g. Das et al. 2020).

Fig. 4. Top panel: comparison of the vertical molecular density solutions with and without dark matter halo. The thick lines represent the density solutions for NGC 3521, and the thin lines represent the same for NGC 7331. Bottom panel: effect of the dark matter halo on the midplane pressure. The thick and thin lines represent the same as the above panel. The solid red line represents the calculated midplane pressure, including the gravity of the dark matter halo, whereas, the blue dashed line represents the same by excluding the dark matter halo. The black dashed-dotted lines represent the calculated midplane pressure using the prescription of Gallagher et al. (2018). See the text for more details.

In the bottom panel of Fig. 4, we show the estimated midplane pressures for NGC 3521 (thick lines) and NGC 7331 (thin lines) using the following three prescriptions: firstly, by solving the hydrostatic equation self-consistently and including the dark matter halo (the solid red lines); secondly, by solving it self-consistently without the dark matter halo (the blue dashed lines); and thirdly, by using the same method used by Gallagher et al. (2018) (as given by Eq. (10)) (black dashed-dotted line). As can be seen from the figure, the formalism of Gallagher et al. (2018) significantly overestimates the midplane pressure. For NGC 3521 and NGC 7331, we find that this overestimation could be by a factor of a few (at the midplane) to as large as by ∼20 (at the outskirts). A substantial contribution to this discrepancy could be attributed to the fact that they overestimate the ρ_* significantly by underestimating the flaring of the stellar disc. They also did not solve the hydrostatic equation in a self-consistent manner by considering the gravitational coupling of the baryonic discs. Further, it can also be seen from the figure that the absence of the dark matter halo (hence reduced potential) underestimates the midplane pressure. For our two sample galaxies, the dark matter gravity can reduce the midplane pressure by ∼20% at an inner radius to ∼5−10% at outer radii.

We solved Eq. (7) self-consistently for all of our sample galaxies to obtain the vertical density distribution of different baryonic components. We further used these density solutions to asses the thickness of the gas discs by estimating the vertical scale heights, defined as the half width at half maxima of the vertical density distribution. In Fig. 5, we plotted the scale height profiles as a function of the radius. Observationally, for a two-component molecular disc, both the disc components are combined together to produce a single measurement. Hence, for these molecular discs, we estimated the total molecular density by combining the densities of the thin and thick disc. It can be seen from the figure that the atomic and molecular discs flare as a function of the radius. For NGC 7331, the atomic scale height is found to be a factor of ∼2 higher than the molecular scale height at all radii. Even so, the scale height of the thick molecular disc and the atomic disc are found to be the same due to their same assumed velocity dispersion. Nonetheless, it should be emphasised that the scale height alone is not a quantitative measure of the observable thickness of a baryonic disc. The measurable thickness also depends on the midplane density of a baryonic component, and hence, in spite of having the same scale heights, the atomic and thick molecular discs could have a different observable thickness.

Fig. 5. Atomic and molecular scale heights in our sample galaxies, assuming their molecular discs to be a two-component system. Each panel in the figure shows the scale heights for different galaxies, as mentioned in the top left corners of the respective panels. The solid red and magenta lines represent the scale heights for the thin and thick molecular discs, respectively, whereas the dashed blue lines show the atomic scale heights. The solid green lines represent the scale heights of the total molecular disc, that is, thin + thick. The cyan and black dashed lines show the exponential fits to the atomic and molecular scale heights. The respective scale lengths of the exponential fits are quoted in the top left (for atomic scale height) and the bottom right (for molecular scale height) corners of each panel in the units of kiloparsecs. For comparison purposes, the molecular scale heights assuming a single-component molecular disc in our sample galaxies are shown by the blue dashed-dotted lines.

It can also be seen from Fig. 5 that the scale height of a thin disc (solid red lines) and the single-component molecular disc (dashed-dotted blue lines) are very similar, though the scale height of the thick disc is much higher. This illustrates that the velocity dispersion is much more influential than the surface density in deciding the scale height. We find that the scale height of the thin and thick molecular discs in our sample galaxies vary between ∼20−200 pc and 100−500 pc, respectively, depending on the radius. In contrast, the total molecular scale height varies between ∼50 and 300 pc. This is higher than the scale height of the single-component molecular disc. To quantify the flaring of the atomic and the molecular discs, we fitted their scale height profiles using a rising exponential function. As can be seen, an exponential function describes the flaring of the scale height reasonably well in our sample galaxies. The respective scale lengths of the exponential scale heights are quoted in the top left (for atomic discs) and bottom right (for molecular discs) corners of each panel in Fig. 5 in the kiloparsecs.

The scale height profiles presented in Fig. 5 were produced by assuming an equal amount of molecular gas in the thin and the thick disc (Pety et al. 2013). However, the observational estimation of scale hight profiles is not well constrained, and a different fractional abundance is very well possible in galaxies (see, e.g. Goldsmith et al. 2008; Liszt et al. 2010). Hence, we explored the effect of this thick disc molecular gas fraction on the scale height of the gas discs by solving Eq. (7) in one of our sample galaxies, NGC 2841, by assuming a thick disc gas fraction of 0.3, 0.5, and 0.7. In Fig. 6, we plotted the atomic and molecular scale height profiles for this galaxy for a different thick disc molecular gas fraction. As can be seen, the atomic scale heights (solid lines) do not change considerably. However, the scale height of the molecular discs varies significantly if one assumes a different amount of molecular gas in the thin and thick disc. For NGC 2841, the molecular scale height has been found to vary by as much as a factor of ∼2 at the centre to ∼20% at the outer radii.

Fig. 6. Atomic and molecular scale height profiles for NGC 2841 for different assumed thick disc molecular gas fraction. We assumed three different thick disc molecular gas fractions of 0.3, 0.5, and 0.7. The solid lines represent the scale height of the atomic disc, whereas the dashed-dotted lines denote the scale height for the molecular disc. The scale height of the total molecular disc can change by as much as two times due to the variation of the thick disc molecular gas fraction from 0.3 to 0.7.

Next, we used the molecular scale height profiles (Fig. 5) to investigate the universality of the flaring of the molecular discs. To do that, we normalised the molecular scale heights of our sample galaxies to unity at a radius of 0.3r₂₅ by adopting a similar approach as used by Schruba et al. (2011). In Fig. 7, we plotted these normalised molecular scale height profiles. As can be seen from the figure, despite having very different surface density profiles and different dark matter halos, the flaring is considerably general and follows a fairly tight relation with the normalised radius, r₂₅. To quantify the flaring further, we averaged the normalised scale height values within a radial bin of 0.1r₂₅ and fitted this mean with an exponential function of the form h_scl = h₀exp(R/l_scl), where l_scl represents the scale length of the exponential flaring. We find a l_scl = (0.48 ± 0.01)r₂₅. For comparison purposes, in Fig. 7, we also plotted the exponential fit to the normalised scale height profile as obtained for single-component molecular discs in our sample galaxies (red dashed-dotted line). For single-component molecular discs, we find a flaring scale length of $l_{\text{scl}}^{\prime }=(0.46\pm 0.01)$r₂₅. This value differs marginally from the value obtained for two-component discs, which indicates that the nature of the flaring is very similar in both the single-component and the two-component molecular discs. We also examined the effect of a different assumed thick disc molecular gas fraction on this flaring scale length. To do that, we carried out the above exercise for thick disc molecular gas fractions of 0.3 and 0.7 as well. We find $l_{\text{scl}}^{\prime }=(0.47\pm 0.01)$r₂₅ and $l_{\text{scl}}^{\prime }=(0.48\pm 0.01)$r₂₅ for a thick disc fraction of 0.3 and 0.7, respectively, which are consistent with the results for the thick disc fraction of 0.5 within error bars. Hence, the nature of the flaring seems to be intrinsic and not profoundly dependent on the assumed thick disc gas fraction.

Fig. 7. Normalised scale height profiles of our sample galaxies for two-component molecular discs. The solid thin lines represent the scale height profiles normalised to 1 at 0.3r25. The blue squares with error bars and the thick blue dashed line indicate the mean normalised scale heights and an exponential fit to it, respectively. The thick red dashed-dotted line shows an exponential fit to the normalised scale height profiles assuming single-component molecular discs in our sample galaxies. The respective scale lengths of the exponential fits are quoted in the top left of the figure. See the text for more details.

However, the molecular scale height is not a directly measurable quantity. Instead, the molecular column density and the spectral cubes are what is observed with radio telescopes. Hence, building total column density maps and spectral cubes for our sample galaxies will be illuminating to investigate the difference between a single-component and a two-component molecular disc. To do that, using the density solutions of the molecular discs and the observed rotation curves, we built three-dimensional dynamical models for our sample galaxies. These models were then inclined to the observed inclinations and projected into the sky plane to generate the spectral cubes and the moment maps. These spectral cubes and the moment maps were further convolved with the telescope beam (13″ × 13″) to mimic the real observations. Thus, we produced model molecular column density maps and spectral cubes for our sample galaxies. To see how these models compare with observations, in Fig. 8, we overplotted the model molecular column density maps (in solid black contours) on top of the observed maps (colour scale) for two representative galaxies, NGC 2841 (top panel) and NGC 7331 (bottom panel). From the top panel of the figure, it can be seen that the observed molecular column density distribution in NGC 2841 is somewhat patchy and lacks a clear azimuthal symmetry. However, as our modelling explicitly assumes an azimuthal symmetry, the modelled molecular column density for NGC 2841 does not match the observation fully. On the other hand, for NGC 7331, the model matches the observation reasonably well (bottom panel). For comparison purposes, we also plotted the model intensity distribution for a three-component disc (single-component molecular disc) using broken contours in both of the panels. As can be seen from the figure, there is no practical difference in the intensity distributions of the molecular discs for a three-component and four-component galactic disc.

Fig. 8. Modelled and observed molecular column density maps for NGC 2841 (top panel) and NGC 7331 (bottom panel). The colour scales in both of the panels represent the CO intensity distribution as observed in the HERACLES survey, whereas the contours represent the modelled intensity distribution. The contour levels are (1, 1.3, 1.6, 1.9, ...) K km s−1 for NGC 2841 and (1.5, 4.4, 7.4, ...) K km s−1 for NGC 7331. The solid contours represent the modelled intensity distribution for a two-component molecular disc, whereas the dashed contours indicate the same for a single-component molecular disc. The vertical dashed lines in each panel enclose a region where the solutions cannot be trusted due to a possible violation of the hydrostatic condition.

For further comparison purposes, in Fig. 9, we plotted the column density maps for a single-component (top panel) and two-component molecular disc (bottom panel) for NGC 7331. As can be seen from the figure, qualitatively, the intensity distribution of the molecular gas looks very similar in both of the panels. This similarity in the intensity maps is observed for other galaxies in our sample as well. We plotted similar figures for all of our sample galaxies in Fig. A.1. For a better comparison of the column density maps, we took vertical cuts along the minor axis and extracted column density profiles. We estimated these column density profiles at different radii for both the simulated and observed molecular discs. In Fig. 10, we show these column density profiles at four different radii for NGC 7331. As can be seen, the simulated profiles (the thick, solid red and blue dashed lines) can be compared to the observed profiles nicely (the thin red and black solid lines), indicating fairly reasonable modelling of the molecular disc in NGC 7331. Other galaxies in our sample also exhibit a similar trend, which we show in Fig. A.2. From these figures, we note that the two simulated molecular discs (single-component and two-component) produce almost the same vertical profiles, even though they have considerably different scale heights (see Fig. 5).

Fig. 9. Modelled intensity maps for a single-component (top panel) and a two-component (bottom panel) molecular disc in NGC 7331. The colour scales in both of the panels indicate the observable CO intensity in the units of K km s−1. As can be seen, there is no apparent difference in the column density maps of the two discs. See the text for more details.

Fig. 10. Vertical column density profiles of the molecular discs in NGC 7331. The red and black thin solid lines are the observed vertical profiles for two different halves of the molecular disc. The thick blue dashed lines and the thick red solid lines represent the simulated vertical column density profiles for the single component and the two-component molecular discs, respectively. Different panels show the vertical profiles at different radii, as quoted in the top left corners of the respective panels. As can be seen, there is no detectable difference in the vertical profiles of the single-component and two-component molecular discs. See the text for more details.

This could happen due to a combined effect of low spatial resolution and a non-edge-on inclination. To check how spatial resolution can influence the observed column density profiles, we imaged the molecular discs of NGC 7331 (single-component and two-component) at different spatial resolutions. For this, we first inclined the molecular discs to 90°, since the maximum difference is expected at an edge-on orientation, and we convolved the resulting column density maps with different observing beams. These maps were then used to extract the column density profiles for comparison purposes. It should be emphasised here that the raw resolution we used is approximately a few parsecs, which is much higher than the adopted observing resolutions. In the top panel of Fig. 11, we show the corresponding vertical column density profiles at R = 4 kpc. The imaging resolutions are quoted in the top left corner of each sub-panel. As can be seen from the figure, for NGC 7331, at a high spatial resolution, the two molecular discs indeed produce different vertical column density profiles; however, the difference slowly washes away as one degrades the resolution. At the mean spatial resolution of the HERACLES survey, that is, ∼500 pc, it fails to produce any detectable difference between the column density profiles of the two discs. We observed the same trend for our other sample galaxies as well. We present their results in Fig. A.3 (first and third rows).

Fig. 11. Top panel: vertical column density profiles of the edge-on molecular discs in NGC 7331 (single-component and two-component) at a radius of 4 kpc for different imaging resolutions. The spatial resolutions are quoted in the top left corners of each sub-panel. Bottom panel: dependence of the normalised vertical column density profile on inclination. In both of the panels, the solid red and blue dashed lines represent the molecular column density profiles for a single-component and two-component molecular disc, respectively. See the text for more details.

We further investigated the effect of inclination on the column density profiles by inclining the molecular discs of NGC 7331 to several different inclinations. To decouple any influence from the imaging resolution, we used a spatial resolution of ∼5 pc for this purpose. In the bottom panel of Fig. 11, we plotted the normalised vertical column density profiles for discs with different inclinations (quoted in the top left corner of each sub-panel). As can be seen from the figure, at an inclination of 90°, the discs could be discriminated, in principle; however, a slight change in the inclination eliminates any difference promptly. The effect of inclination on vertical profiles for our other sample galaxies are presented in Fig. A.3 (second and fourth rows). These results indicate that the vertical profiles are highly sensitive to inclination, and it is tough to identify the presence of a two-component molecular disc if the galaxy is not an edge-on system.

Next, we explored the spectral properties of the discs to examine if they can be used to identify a two-component molecular disc in galaxies. Using the dynamical models of NGC 7331 (and other sample galaxies), we produced the spectral cubes for a single-component and two-component molecular disc for typical observing setups. We adopted the mean spatial resolution of the HERACLES survey, that is, ∼500 pc, and used a spectral resolution of 2.6 km s⁻¹. To simulate the real observing scenario, we added appropriate noise levels to the molecular spectral cubes. We note that the molecular surface densities of our sample galaxies were calculated by stacking the spectra within radial bins (Schruba et al. 2011). The stacking results in a much higher S/N on the final radial surface density profiles as compared to a moment map produced by adding the fluxes across the channels. In that sense, the noise in a spectral channel is much higher than the lowest surface density detected (in a radial bin) by stacking. Due to this reason, we did not simply add the noise as measured from the observed spectral cube. Instead, we added such an amount of noise to our simulated spectral cubes so that it produces an S/N of five at the farthest radial bin. A minimum S/N of five is necessary to determine the centroids of the spectra while stacking (de Blok et al. 2008). We emphasise that the choice of our noise levels does not affect our results in any significant way. Adopting the same technique, which led to the strong proposition for two-component molecular discs in galaxies, as used by Caldú-Primo et al. (2013), we stacked the molecular spectra in both of the discs. As the original signature of a two-component molecular disc would have a larger spectral width than a single-component disc, an artificial broadening due to blending should be minimised for a more explicit inspection. Hence, we excluded a region of 0.2R₂₅ around R = 0 from stacking to avoid significant blending due to the rotation of the galactic disc (Caldú-Primo et al. 2013). Even in doing so, there would be blending outside of this central region as well, which significantly depends on the disc inclination. To examine the effect of the inclination on the spectral widths, we produced the spectral cubes and stacked spectra for different molecular discs with different inclinations.

Thus, generated stacked spectra for a single-component (solid red lines) and a two-component molecular disc (blue dashed lines) in NGC 7331 are compared in Fig. 12. Different panels show the results for different inclinations, as quoted in the top left corner of each panel. As can be seen from the figure, a two-component disc systematically produces wider stacked spectra at all inclinations as compared to a single-component disc, though the difference narrows down as a function of inclination. However, we note that the difference in the widths (∼1−2 km s⁻¹) is not so extensive as to identify a two-component molecular disc just based on the absolute width of the stacked spectra. This result is consistent across all of our sample galaxies. We present the respective stacked spectra in Fig. A.4.

Fig. 12. Comparison of the widths of normalised stacked spectra for a single-component and two-component molecular disc at different observing inclinations. The solid red lines indicate the stacked spectra for a single-component molecular disc, whereas the blue dashed lines show the same for a two-component disc. Different panels represent the effect of different assumed inclinations. The inclinations are quoted in the top left corners of the respective panels. The stacked spectra of a two-component molecular disc were found to produce a systematically higher width than that of a single-component disc. See the text for more details.

Along with the spectral width, we investigated the spectral shapes as well to examine the spectral behaviour of a two-component molecular disc. For a two-component disc, it is expected that the stacked spectra have the response of both the thin and thick disc velocity dispersion and hence are best fitted by a double-Gaussian. Whereas, a single-Gaussian is presumed to be the best to describe the stacked spectra of a single-component disc. However, spectral blending (inclination dependent) complicates the stacked spectra, and hence, it is not trivial to predict the expected spectral shape a priori. To get a more in-depth insight, we fitted the stacked spectra of both of the discs with a single and double-Gaussian and study them as a function of inclination. In the top panel of Fig. 13, we show the fits for a single-component disc in NGC 7331. Ideally, in the absence of any blending, the spectra should be best fitted by a single-Gaussian with the same width as the velocity dispersion of the molecular disc (7 km s⁻¹). However, as can be seen from the figure, at high inclinations, there is a significant blending, and as a result, the stacked spectra are best fitted by a double-Gaussian. Surprisingly, the width of the two Gaussian components for a single-component disc at an inclination of ≃60° (rightmost panel on the top row) are found to be ∼7 km s⁻¹ and ∼12 km s⁻¹, which accurately mimics a two-component molecular disc. Hence, the spectral blending at high inclinations can induce the signatures of a two-component disc in the spectral cube of a single-component disc. However, at inclinations ≲40°, the blending is less, and the residual of a single-Gaussian fit to the stacked spectra does not show significant features. The spectral width is also recovered reasonably well by a single-Gaussian fit. On the other hand, the stacked spectra of a two-component molecular disc are always best fitted by double-Gaussians, irrespective of the disc inclination (bottom row in Fig. 13); however, the recovered spectral widths substantially deviate from the assumed velocity dispersions at high inclinations (≳40°). Nevertheless, at inclinations ≲40°, the recovery of the velocity dispersions is satisfactory. These results indicate that for a cleaner inspection of the presence of a two-component molecular disc in galaxies, the spectral cubes should be stacked for discs with an inclination ≲40°. We note that while building the spectral cubes, we assumed a spectral resolution of 2.6 km s⁻¹ and a spatial resolution of ∼500 pc. These kinds of resolutions are typical and not very coarse so as to consider these results non-typical. The widening of the stacked spectra indeed would depend on the observed spatial resolution and the rotation curve, however, and the parameters for NGC 7331 are typical to those of normal spiral galaxies. In fact, in many physical properties, NGC 7331 is so similar to the Galaxy that very often, it is called the “Milky Way’s twin”. For consistency purposes, we also carried out the same exercise for all of our sample galaxies (see Fig. A.5). Despite having a wide range of rotation curves and physical properties, they show the same signatures in their spectral shapes, as we find for NGC 7331. This confirms that the spectral shape of molecular discs with low inclination could be used to identify the presence of a two-component molecular disc in galaxies. We note that we assume an equal amount of molecular gas in the thin and thick disc to produce the spectral cubes. Nevertheless, due to the high S/N of the stacked spectra, our conclusion is not sensitive to the thick disc molecular gas fraction in a significant way. However, a detailed analysis is required to understand the dependence of the absolute molecular spectral widths on the thick disc molecular gas fraction in galaxies, which we plan to do next.

Fig. 13. Result of Gaussian fittings to the stacked spectra of the single component and two-component molecular disc. Large top panels: stacked spectra and their fits with single and double Gaussians for the single-component disc at different inclinations, whereas large bottom panels: same for a two-component molecular disc. Smaller panels in the top and bottom: residuals of the respective fits. The cyan crosses represent the stacked spectra, whereas the solid red lines and the blue dashed lines represent a single Gaussian and a double Gaussian fit to the spectra, respectively. The widths of the fitted Gaussian components (sigma) are quoted in the top left corners of each panel. The first number on the left (in red) is the sigma of a single Gaussian fit, whereas the blue numbers just below them are the widths for a double Gaussian fit. The assumed inclinations of the molecular discs are quoted at the top right corners of the respective panels.

5. Summary and conclusion

We consider the molecular discs in spiral galaxies to be a single component or a two-component (thin + thick) system in hydrostatic equilibrium under the mutual gravity of the baryonic discs in the external potential of the dark matter halo. We set up the respective joint Poisson’s-Boltzmann equation of hydrostatic equilibrium for the baryonic discs and solved them numerically in a sample of eight spiral galaxies from the HERACLES and THINGS survey to determine the three-dimensional density distribution of the molecular gas in these galaxies. We assumed a σ_CO = 7 and 12 km s⁻¹ for the thin and thick molecular discs, respectively. From the solutions of the hydrostatic equation, we find that the molecular gas in the thick disc extends to more significant heights than the molecular gas in a single-component disc or thin disc. We also find the thick disc molecular gas to have a very similar scale height to that of the atomic gas, driving the possibility of forming a common dynamical component together.

Exploiting the joint Poisson’s-Boltzmann equation, we evaluate the impact of the dark matter halo on deciding the vertical density distribution and the midplane pressure. In our two sample galaxies, NGC 3521 and NGC 7331, we find that ignoring the dark matter halo potential can significantly underestimate the vertical density close to the midplane and overestimate the same at more significant heights. We also find that the dark matter can considerably influence the midplane pressure as well. Ignoring dark matter gravity can reduce the midplane pressure by as much as ∼20% at the inner radii and ∼5−10% at the outer radii. A comparison of the midplane pressure with previous studies (e.g. Gallagher et al. 2018) reveals that solving the Poisson’s-Boltzmann equation self-consistently is crucial in estimating the midplane pressure correctly.

We further used the density solutions of the hydrostatic equation to estimate the vertical atomic and molecular scale heights in our sample galaxies. We find that the molecular scale height for a single-component molecular disc varies between 20 and 200 pc, whereas the molecular scale height is found to vary between 50 and 300 pc for a two-component disc. While solving the hydrostatic equilibrium equation, we primarily assumed that the molecular gas is equally distributed between the thin and thick discs. Even so, we also find that the vertical density profile and the molecular scale height is sensitive to the assumed thick disc molecular gas fraction and can change by a factor of approximately two for a change in the thick disc molecular gas fraction from 0.3 to 0.7. We further find that both the molecular discs flare as a function of the radius, and the normalised scale heights follow a tight exponential law. The scale length of this exponential flaring is found to be (0.48 ± 0.01)r₂₅ and (0.46 ± 0.01)r₂₅ for two-component and single-component molecular discs, respectively. However, as the molecular scale height is not a directly measurable quantity, we used the density solutions to build dynamical models, column density maps, and spectral cubes for our sample galaxies. We examined the column density maps of our sample galaxies in detail to find that a non-edge-on inclination and a low spatial resolution leads to no apparent difference in the column density distribution between the single component and the two-component molecular discs.

We study the column density maps of our sample galaxies at different assumed inclinations and spatial resolutions. Using the column density maps at 90° inclination, we demonstrate that at this edge-on orientation, with a high spatial resolution (≲100 pc), a single component molecular disc produces a considerably thinner observed disc as compared to a two-component molecular disc. This difference quickly disappears as the spatial resolution reaches ∼500 pc, which is typical of earlier large surveys. However, recent observational capabilities allow one to image the full disc of nearby galaxies with a much higher spatial resolution. For example, the PdBI Arcsecond Whirlpool Survey (PAWS) images the M 51 galaxy with a spatial resolution of ∼40 pc (Pety et al. 2013; Schinnerer et al. 2013). The Physics at High Angular resolution in Nearby GalaxieS with ALMA survey (PHANGS-ALMA) aims to image the molecular gas in 74 nearby galaxies with a spatial resolution of ∼45−120 pc (Sun et al. 2018; Utomo et al. 2018; Schinnerer et al. 2019; Chevance et al. 2020). The high-quality data from these surveys will be particularly useful to distinguish between a single-component and two-component molecular disc in galaxies. Besides the spatial resolution, we also find that the inclination is also very critical in imprinting any observable difference in the column density maps. A small departure from the edge-on orientation (∼5 ° −10°) washes away any difference between the observed vertical column density profiles of the two discs.

Furthermore, in using the dynamical models of our sample galaxies, we simulated the spectral cubes for a single-component and two-component molecular disc for different inclinations. Using a traditional approach, we stacked the molecular spectra in these spectral cubes to produce global stacked spectra. Single and double Gaussian fittings to these spectra reveal that at inclinations i ≳ 40°, spectral blending is significant and increases the spectral widths considerably. Not only that, but also at high inclinations, this broadening fictitiously produces a broad component into the stacked spectra, which exactly mimics the existence of a thick molecular disc in the spectral cube of a single-component molecular disc. However, at low inclinations (i ≲ 40°), the stacked spectra of a single and a two-component molecular disc are best fitted by a single and a double component Gaussian, respectively. Hence, we conclude that the spectral cubes of low inclination (i ≲ 40°) galaxies are ideal for examining the existence of a thick molecular disc in spiral galaxies.

In summary, we conclude that given the typical spatial resolution and inclination, it is tough to distinguish between a single component and two-component molecular disc just by inspecting the column density maps or the spectral cubes. However, at an edge-on orientation with a high spatial resolution, it is possible to identify if the observed molecular disc is produced by a single or two-component molecular disc. Whereas for a face-on galaxy (or low inclination, i ≲ 40°), the spectral cube can be used more effectively to identify the presence of a two-component molecular disc. The stacked spectrum of a two-component disc in these cases is expected to be better fitted by a double Gaussian than a single Gaussian.

Acknowledgments

NNP would like to thank the anonymous referee for producing a thoughtful review of the manuscript, which helped to improve the quality of the paper. This paper uses the data products extensively from two large surveys, that is, the THINGS and the HERACLES. NNP would like to acknowledge the respective team members for making the data publicly available. NNP would like to thank Dr. Samir Choudhuri and Dr. Gunjan Verma for their valuable comments and suggestions, which have helped to improve the quality and readability of this manuscript.

References

Appendix A: Additional figures

Fig. A.1. Modelled molecular discs of all of our sample galaxies (same as Fig. 9). First and third rows from the top: CO intensity maps for single-component molecular discs. Second and fourth rows from the top: same for two-component molecular discs. The names of the galaxies are quoted in the top right corner of the respective panels. The colour bars represent the simulated CO intensities in the units of K km s−1.

Fig. A.2. Vertical column density profiles of the molecular discs of all of our sample galaxies. The figure legends are the same as in Fig. 10.

Fig. A.3. Comparison of the vertical column density profiles for a single-component and a two-component molecular disc in our sample galaxies (same as Fig. 11). First and third rows: vertical column density profiles for different spatial resolutions, whereas second and fourth column: same for different inclinations. Different sub-panels represent different spatial resolutions or inclinations, as quoted in the top left of each sub-panel. The vertical profiles were extracted at a radius of 4 kpc for all of the galaxies except NGC 2976 for which it was computed at 1.5 kpc (as its molecular disc does not extend to a radius of 4 kpc). The solid red lines represent the vertical profiles for single-component discs, whereas the blue dashed lines represent the same for two-component molecular discs.

Fig. A.4. Comparison of the normalised stacked spectra of our sample galaxies (same as Fig. 12) for different inclinations. Different panels show the stacked spectra for different galaxies, as quoted in the top right corner of every panel. The solid red lines represent the stacked spectra for single-component molecular discs, whereas the blue dashed lines represent the same for two-component molecular discs. The inclinations are quoted in the top left corners of each sub-panel.

Fig. A.5. Same as Fig. 13, but for galaxies NGC 925 (top panel), NGC 2841 (middle panel), and NGC 2976 (bottom panel).

Fig. A.5. continued. For galaxies NGC 3198 (top panel), NGC 3521 (middle panel), and NGC 5055 (bottom panel).

Fig. A.5. continued. For galaxies NGC 6946 (top panel) and NGC 7331 (bottom panel).

All Tables

All Figures

Fig. 1. Surface density profiles for our sample galaxies. Different panels indicate the surface densities for different galaxies, as mentioned in the top right corner of each panel. The cyan pluses, the red triangles, and the blue circles represent the stellar, the atomic, and the molecular surface density profiles, respectively. The stellar surface density profiles are taken from Leroy et al. (2009), whereas the atomic and the molecular surface density profiles are obtained from Schruba et al. (2011).

In the text

Fig. 2. Rotation curves of our sample galaxies as derived using the HI data from the THINGS survey by de Blok et al. (2008). The rotation curves were fitted using a Brandt profile (Eq. (8)) for parametrisation. In each panel, the black circles with error bars represent the rotation curve, whereas the red solid line represents the Brandt profile fit. The best-fit parameters can be found in Table 2.

In the text

Fig. 3. Sample solutions of Eq. (7) for NGC 7331 at R = 8 kpc. Top panel: density solutions for a three-component galactic disc having a single-component molecular disc. Bottom panel: same for a four-component disc. In this case, the molecular disc is assumed to be a two-component system with a thin and a thick disc. The red dotted and the solid blue lines in both of the panels represent the volume densities of the stellar and atomic discs, respectively. The black dashed line in the top panel represents the volume density of the thin disc molecular gas in a three-component disc. Bottom panel: the thin dashed-dotted black and dashed magenta lines represent the volume density of molecular gas in the thin and thick discs, respectively, for an assumed thick disc gas fraction of 0.5. Whereas, the thick dashed-dotted black and the dashed magenta lines represent the same for an assumed thick disc molecular gas fraction of 0.7. See the text for more details.

In the text

Fig. 4. Top panel: comparison of the vertical molecular density solutions with and without dark matter halo. The thick lines represent the density solutions for NGC 3521, and the thin lines represent the same for NGC 7331. Bottom panel: effect of the dark matter halo on the midplane pressure. The thick and thin lines represent the same as the above panel. The solid red line represents the calculated midplane pressure, including the gravity of the dark matter halo, whereas, the blue dashed line represents the same by excluding the dark matter halo. The black dashed-dotted lines represent the calculated midplane pressure using the prescription of Gallagher et al. (2018). See the text for more details.

In the text

Fig. 5. Atomic and molecular scale heights in our sample galaxies, assuming their molecular discs to be a two-component system. Each panel in the figure shows the scale heights for different galaxies, as mentioned in the top left corners of the respective panels. The solid red and magenta lines represent the scale heights for the thin and thick molecular discs, respectively, whereas the dashed blue lines show the atomic scale heights. The solid green lines represent the scale heights of the total molecular disc, that is, thin + thick. The cyan and black dashed lines show the exponential fits to the atomic and molecular scale heights. The respective scale lengths of the exponential fits are quoted in the top left (for atomic scale height) and the bottom right (for molecular scale height) corners of each panel in the units of kiloparsecs. For comparison purposes, the molecular scale heights assuming a single-component molecular disc in our sample galaxies are shown by the blue dashed-dotted lines.

In the text

Fig. 6. Atomic and molecular scale height profiles for NGC 2841 for different assumed thick disc molecular gas fraction. We assumed three different thick disc molecular gas fractions of 0.3, 0.5, and 0.7. The solid lines represent the scale height of the atomic disc, whereas the dashed-dotted lines denote the scale height for the molecular disc. The scale height of the total molecular disc can change by as much as two times due to the variation of the thick disc molecular gas fraction from 0.3 to 0.7.

In the text

Fig. 7. Normalised scale height profiles of our sample galaxies for two-component molecular discs. The solid thin lines represent the scale height profiles normalised to 1 at 0.3r25. The blue squares with error bars and the thick blue dashed line indicate the mean normalised scale heights and an exponential fit to it, respectively. The thick red dashed-dotted line shows an exponential fit to the normalised scale height profiles assuming single-component molecular discs in our sample galaxies. The respective scale lengths of the exponential fits are quoted in the top left of the figure. See the text for more details.

In the text

Fig. 8. Modelled and observed molecular column density maps for NGC 2841 (top panel) and NGC 7331 (bottom panel). The colour scales in both of the panels represent the CO intensity distribution as observed in the HERACLES survey, whereas the contours represent the modelled intensity distribution. The contour levels are (1, 1.3, 1.6, 1.9, ...) K km s−1 for NGC 2841 and (1.5, 4.4, 7.4, ...) K km s−1 for NGC 7331. The solid contours represent the modelled intensity distribution for a two-component molecular disc, whereas the dashed contours indicate the same for a single-component molecular disc. The vertical dashed lines in each panel enclose a region where the solutions cannot be trusted due to a possible violation of the hydrostatic condition.

In the text

	Fig. 9. Modelled intensity maps for a single-component (top panel) and a two-component (bottom panel) molecular disc in NGC 7331. The colour scales in both of the panels indicate the observable CO intensity in the units of K km s−1. As can be seen, there is no apparent difference in the column density maps of the two discs. See the text for more details.
In the text

Fig. 10. Vertical column density profiles of the molecular discs in NGC 7331. The red and black thin solid lines are the observed vertical profiles for two different halves of the molecular disc. The thick blue dashed lines and the thick red solid lines represent the simulated vertical column density profiles for the single component and the two-component molecular discs, respectively. Different panels show the vertical profiles at different radii, as quoted in the top left corners of the respective panels. As can be seen, there is no detectable difference in the vertical profiles of the single-component and two-component molecular discs. See the text for more details.

In the text

Fig. 11. Top panel: vertical column density profiles of the edge-on molecular discs in NGC 7331 (single-component and two-component) at a radius of 4 kpc for different imaging resolutions. The spatial resolutions are quoted in the top left corners of each sub-panel. Bottom panel: dependence of the normalised vertical column density profile on inclination. In both of the panels, the solid red and blue dashed lines represent the molecular column density profiles for a single-component and two-component molecular disc, respectively. See the text for more details.

In the text

Fig. 12. Comparison of the widths of normalised stacked spectra for a single-component and two-component molecular disc at different observing inclinations. The solid red lines indicate the stacked spectra for a single-component molecular disc, whereas the blue dashed lines show the same for a two-component disc. Different panels represent the effect of different assumed inclinations. The inclinations are quoted in the top left corners of the respective panels. The stacked spectra of a two-component molecular disc were found to produce a systematically higher width than that of a single-component disc. See the text for more details.

In the text

Fig. 13. Result of Gaussian fittings to the stacked spectra of the single component and two-component molecular disc. Large top panels: stacked spectra and their fits with single and double Gaussians for the single-component disc at different inclinations, whereas large bottom panels: same for a two-component molecular disc. Smaller panels in the top and bottom: residuals of the respective fits. The cyan crosses represent the stacked spectra, whereas the solid red lines and the blue dashed lines represent a single Gaussian and a double Gaussian fit to the spectra, respectively. The widths of the fitted Gaussian components (sigma) are quoted in the top left corners of each panel. The first number on the left (in red) is the sigma of a single Gaussian fit, whereas the blue numbers just below them are the widths for a double Gaussian fit. The assumed inclinations of the molecular discs are quoted at the top right corners of the respective panels.

In the text

Fig. A.1. Modelled molecular discs of all of our sample galaxies (same as Fig. 9). First and third rows from the top: CO intensity maps for single-component molecular discs. Second and fourth rows from the top: same for two-component molecular discs. The names of the galaxies are quoted in the top right corner of the respective panels. The colour bars represent the simulated CO intensities in the units of K km s−1.

In the text

	Fig. A.2. Vertical column density profiles of the molecular discs of all of our sample galaxies. The figure legends are the same as in Fig. 10.
In the text

Fig. A.3. Comparison of the vertical column density profiles for a single-component and a two-component molecular disc in our sample galaxies (same as Fig. 11). First and third rows: vertical column density profiles for different spatial resolutions, whereas second and fourth column: same for different inclinations. Different sub-panels represent different spatial resolutions or inclinations, as quoted in the top left of each sub-panel. The vertical profiles were extracted at a radius of 4 kpc for all of the galaxies except NGC 2976 for which it was computed at 1.5 kpc (as its molecular disc does not extend to a radius of 4 kpc). The solid red lines represent the vertical profiles for single-component discs, whereas the blue dashed lines represent the same for two-component molecular discs.

In the text

Fig. A.4. Comparison of the normalised stacked spectra of our sample galaxies (same as Fig. 12) for different inclinations. Different panels show the stacked spectra for different galaxies, as quoted in the top right corner of every panel. The solid red lines represent the stacked spectra for single-component molecular discs, whereas the blue dashed lines represent the same for two-component molecular discs. The inclinations are quoted in the top left corners of each sub-panel.

In the text

	Fig. A.5. Same as Fig. 13, but for galaxies NGC 925 (top panel), NGC 2841 (middle panel), and NGC 2976 (bottom panel).
In the text

	Fig. A.5. continued. For galaxies NGC 3198 (top panel), NGC 3521 (middle panel), and NGC 5055 (bottom panel).
In the text

	Fig. A.5. continued. For galaxies NGC 6946 (top panel) and NGC 7331 (bottom panel).
In the text