xGASS: The connection between angular momentum, mass and atomic gas fraction in nearby galaxies

We assume that the baryonic component of disc galaxies can be fully represented by a combination of stars, atomic gas and molecular gas. We use molecular gas masses from xCOLD GASS, (this includes a helium contribution). Atomic gas is comprised of atomic hydrogen (Hi, taken from xGASS) and helium (assumed to be 35% of the Hi). We also assume that all the hydrogen and inferred helium are located in the disc. Using bulge-to-disc decompositions, the stellar mass of the disc component can be separated (for details see Cook et al. 2019 and Hardwick et al. 2022). Therefore, the baryonic mass of the disc is given as follows:

where $M_{\star,D}$ is the stellar mass of the disc, $M_{HI}$ is the mass of atomic hydrogen gas and $M_{\rm{mol}}$ is the mass of molecular gas.

Despite all our galaxies having Hi and stellar masses, only 45% of our sample has molecular gas masses from xCOLD GASS. Therefore, the remaining 363 galaxies in our sample need their molecular gas mass estimated. We use the SFR surface density ($\Sigma_{\rm{SFR}}=\rm{SFR}/\pi R_{e}^{2}$) vs. molecular mass surface density ($\Sigma_{\rm{mol}}=M_{\rm{mol}}/\pi R_{e}^{2}$) relation (i.e., the inverted Kennicutt-Schmidt law, Kennicutt 1989; Schmidt 1959) to estimate molecular gas for the galaxies not observed or with only upper limits in xCOLD GASS. Here, $R_{e}$ is the r-band half-light radius of the total galaxy. Figure 1 shows this relation for the subset of our sample observed with xCOLD GASS.

For this sub-sample, the SFR and molecular gas surface densities are related by the following formula;

which was determined by fitting a power-law (blue line) in hyper-fit¹¹1This same code has been used to determine all the best fits in this paper. (Robotham & Obreschkow, 2015) to the galaxies with H₂ detections. Using equation 2, we obtain molecular gas masses from SFRs for the galaxies where a CO detection is not available. For two of the galaxies in our sample, we do not have SFRs or H₂ detections (GASS123007 and GASS38198). For these, we determine the molecular gas mass by assuming that their Hi-to-H₂ mass ratio is the average value observed in our sample (i.e., 20%). As H₂ contributes very little to both $M_{\rm{bar}}$ and $j_{\rm{bar}}$, these assumptions do not affect our results.

To calculate the baryonic specific AM of the disc ($j_{\rm{bar,D}}$) we take the mass-weighted average of the specific AM of each component; stars, Hi and molecular gas, as follows:

where $j_{\star,D}$, $j_{\rm{mol}}$ and $j_{HI}$ are the specific AM of the disc stars, molecular gas and Hi, respectively. We use the values of $j_{\star,D}$ that are available from Paper 1.

We assume that $j_{\rm{mol}}$ is equal to $j_{\star,D}$. This is because it has been shown in works such as Bigiel & Blitz (2012) and Smith et al. (2016) that the surface density profile of the H₂ follows closely that of the stars. The specific AM is primarily set by the surface mass density distribution (assuming that the rotation curves are similar for both components), so it follows that the specific AM of the H₂ is the same as the stars.

To determine $j_{HI}$ from spatially unresolved data, we proceed as follows. Firstly, we assume the Hi surface density profile ($\Sigma_{HI}$) to be universal for all galaxies. This was shown to be a good approximation for the sample in Wang et al. (2016) when radii were normalised by the Hi diameter. To determine the Hi diameter we use the intrinsically tight Hi mass to size relation from Wang et al. (2016):

where $D_{HI}$ is the diameter of the Hi disc measured at a surface density of $\Sigma_{HI}=1\,\rm{M}_{\odot}\,pc^{-2}$. We tested the reliability of this approach by also allowing for a variation in diameter consistent with 1$\sigma$ above and below the relation. Accounting for the scatter in the Hi mass-diameter relation changes $j_{HI}$ by a mean difference of $\pm 0.06$ dex. The $\Sigma_{HI}$ profiles from Wang et al. (2016) are only given out to $1.3R_{HI}$, (where $R_{HI}=0.5D_{HI}$), which for the majority of our galaxies (92% of the sample) is less than the extent used to estimate $j_{\star}$ (i.e., $10R_{e}$). For consistency with all other quantities $j_{HI}$ is also calculated within $10R_{e}$. To do so, we linearly extrapolate the universal $\Sigma_{HI}$ profile (in log-log space). This extrapolation slightly increases the value of $j_{HI}$ (less than 0.03 dex). Conversely, for the 46 galaxies, where $10R_{e}<1.3R_{HI}$, $j_{HI}$ is reduced by an average of 0.02 dex up to a maximum of 0.15 dex, when compared to $j_{HI}$ calculated within $1.3R_{HI}$.

Lastly, we assume that the velocity profile is constant for all radii and given by the width of the Hi emission line. This is the same assumption made in Paper 1 for calculating stellar specific AM. This assumption was tested in Paper 1 by comparing constant velocity profiles to the template rotation curves in Catinella et al. (2006), where we showed that the most extreme case of low-mass dwarf galaxies (i.e. $M_{\star}\approx 10^{9}M_{\odot}$), had a maximum systematic offset in j of $\sim 0.08$ dex, while high-mass galaxies had less discrepancy ($<0.04$ dex). Therefore, this assumption only introduces minor systematic effects which do not impact our result (for more details see Paper 1). We emphasise that, despite using a constant velocity profile, our approach does not rely on the approximation of $j\propto V_{max}R_{e}$ (Romanowsky & Fall, 2012), and instead estimates j by directly integrating the AM profile obtained from the surface density profiles of each baryonic component. We found this to be a more accurate way to determine AM, as the Romanowsky & Fall (2012) approximation can overestimate a galaxies AM by up to 25% at low masses (e.g., Obreschkow & Glazebrook, 2014; Bouché et al., 2021).

Once all of these assumptions are combined, we obtain $j_{HI}$. When comparing $j_{\star,D}$ and $j_{HI}$, on average $j_{HI}$ is roughly twice the value of $j_{\star,D}$, although there is a considerable spread ($\sigma=0.22$ dex). Note that this factor $\sim$2 was also found in things (Obreschkow & Glazebrook, 2014) using resolved Hi kinematic maps.

We reiterate that although we have outlined many assumptions in this subsection, these choices are not likely to affect our results. Although we assume that stars and gas share the same rotation curve (by using Hi line widths to calculate AM for both), this does not impact our results. Stellar rotation curves rise more slowly than cold gas ones due to asymmetric drift, so our estimates of $j_{\star,D}$ may be slightly overestimated. However, $j_{\rm{bar}}$ would not change considerably if $j_{\star,D}$ was altered slightly. In addition, assumptions made when calculating $j_{\star,D}$ were robustly tested in Paper 1. As mentioned above, the assumptions involving $j_{HI}$ have been tested, and even if these assumptions are altered, they do not affect the result significantly. In addition, the molecular mass is only a small component of our galaxies, such that excluding molecular material, only alters $M_{\rm{bar}}$ by 0.02 dex and $j_{\rm{bar}}$ by 0.01 dex on average. In summary, the typical uncertainty on our estimate of $j_{\rm{bar}}$ is 0.14 dex, which is dominated by the assumptions made on $j_{HI}$. This is comparable to the statistical error on $M_{\rm{bar}}$, i.e., 0.11 dex.

As a first step, in Figure 2 we show the baryonic Fall relation for galaxies in xGASS.

The left panel shows the baryonic mass and specific AM for just the disc component ($M_{bar,D}$ and $j_{bar,D}$), while on the right we show the total (i.e., disc and bulge component; $M_{bar,T}$ and $j_{bar,T}$). Although the majority of the paper will focus on just the disc components, we also show the right panel to compare with previous works. The best-fit linear relation to our data is shown as a thin black line, with 1$\sigma$ scatter shown by the dashed lines. The blue line shows the baryonic Fall relation found by Murugeshan et al. (2020) and the magenta line is from Mancera Piña et al. (2021a), which both agree well with the best-fitting relation to our work. This is not a trivial result given the difference in samples and techniques used to determine this relation. Therefore, the agreement with previous work gives us confidence that our approach is solid and we can move forward in our analysis.

As for previous works, we suggest caution when interpreting the slopes of the baryonic Fall relation at face value, as our sample is not selected by baryonic mass. xGASS is selected to only include galaxies with stellar masses greater than $10^{9}\,\rm{M}_{\odot}$ with an oversampling at the high stellar mass end to increase statistics (as shown in Paper 1 this distribution is conserved for our sub-sample). As such, the distribution in baryonic mass of our sample is not representative, with a deficiency of galaxies at lower baryonic masses. Below a baryonic mass of approximately $10^{9.8}\,\rm{M}_{\odot}$, at fixed baryonic mass, our sample will preferentially miss galaxies at high gas fractions, whose stellar mass lies below the cut of $10^{9}\,\rm{M}_{\odot}$. In addition, we know from the results of works such as Mancera Piña et al. (2021a, see also Paper 1) that galaxies with higher gas fractions also have higher $j_{\rm{bar}}$. Therefore, we would expect the Fall relation to be shallower for a sample that is uniformly distributed in baryonic mass. This echoes the conclusions of Paper 1, which found that the sample used to determine a Fall relation will affect the slope obtained.

O16 developed a formalism to look at the connection between specific baryonic AM, baryonic mass and atomic gas. This resulted in a dimensionless effective stability parameter q;

where $j_{\rm{bar,D}}$ is the baryonic specific AM of the disc component, $\sigma_{HI}$ is the one-dimensional velocity dispersion of Hi, G is the gravitational constant and $M_{\rm{bar,D}}$ is the baryonic mass of the disc component. We use $M_{\rm{bar,D}}$ and $j_{\rm{bar,D}}$ as defined in sections 3.1 and 3.2 respectively. $\sigma_{HI}$ is assumed to be a constant 10  km s^-1, which is set by the temperature of the warm neutral medium (following the assumption of O16). For a flat exponential disc, the q parameter relates to the expected fraction of atomic gas via;

Here the fraction of atomic gas is defined as the fraction of disc baryonic mass in neutral atomic hydrogen and helium;

This stability model assumes that any atomic gas that is unstable (following the Toomre 1964 stability parameter Q) will collapse to form stars. As a consequence, if $f_{\rm{atm}}>\rm{min}\{1,2.5\ q^{1.12}\}$, then the galaxy is over-saturated in Hi and is bound to form stars quickly (i.e., within a giant molecular clouds free-fall time). However, in the opposite case ($f_{\rm{atm}}<\rm{min}\{1,2.5\ q^{1.12}\}$) a galaxy will remain stable, with no self-regulating processes to alter its atomic gas fraction. Therefore, the O16 relation (equation 6) can be thought of as an upper limit of the amount of atomic gas that a galaxy can hold (for further details of the models assumptions see O16).

Figure 3 shows the relationship between the q stability parameter and $f_{atm}$. We compare our sample (grey) to the data presented in O16 (coloured points) and their model (black line). The binned median of our data (thick blue line) shows good agreement with the theoretical relation shown in black, however, our data show a larger scatter. For our data the median 16th to 84th percentile is 0.52 dex, whereas it is 0.27 dex for the original O16 data. This result does not depend on the method used to determine $M_{\rm{mol}}$, as $f_{atm}$ only weakly depends on the molecular gas mass, and is mostly constrained by the Hi mass. The significant scatter observed for xGASS is an important result, as previous work have used the small observed scatter in the $f_{atm}-q$ relation to support a scenario in which $q$ is the primary galaxy property physically connected to gas fraction (e.g., Obreschkow et al., 2016; Lutz et al., 2018; Murugeshan et al., 2019). Instead, for our sample, the $f_{atm}-q$ relation shows a spread comparable to that of the well-known relation between gas fraction and $NUV-r$ colour (median 16th to 84th percentile 0.49 dex).

As mentioned above, the O16 model predicts that galaxies will not lie above the relation for long, as their Hi gas will be unstable and collapse to form stars. O16 expect most galaxies to lie within 40% of the analytical relation given the approximate uncertainties of the model, with only 6 out of their 105 galaxies lying above this point. In our data, above this limit, we have roughly double the fraction (61/559) than found in O16. To understand if this is simply due to issues with the data, we analysed the 61 outlier galaxies via the Tully-Fisher relation (Tully & Fisher, 1977) and found that only 11 of them appear to have underestimated rotational velocities. Therefore, as this issue does not affect the rest of the population scattering above the theoretical line, it suggests that the larger scatter compared to what is seen in O16 is not due to data quality but intrinsic to our sample.

Previous works have shown that the down-ward scatter of the $f_{atm}-q$ relation can be due to environmental effects removing the gas from the disc without affecting the stability parameter (e.g. Li et al., 2020; Cortese et al., 2021). If we divide our sample into central (N=432) and satellites (N=122) according to the classification from the Yang et al. (2007) group catalogue, we find that even in our sample, at fixed q, satellites show a lower $f_{atm}$ (0.2 dex on average) than central galaxies. This is intriguing, as xGASS includes less extreme density environments than the Virgo cluster analysed in Li et al. (2020). However, even if we focus on central galaxies alone, the scatter in the $f_{atm}-q$ relation remains substantial, with a median 16th to 84th percentile of 0.42 dex. This means that the large spread in $f_{atm}$ at fixed q does not only trace environmentally-driven gas stripping.

To further understand the physical origin for this large scatter, we tested many potential galaxy parameters that could be correlated with the vertical offset from the theoretical $f_{atm}$ - q relation, defined as follows;

Namely, we investigated: the bulge-to-total mass ratio (B/T), Sersic index, Hi depletion time ($M_{HI}$/SFR), baryonic mass ($M_{\rm{bar,D}}$), specific SFR and the offset of galaxies from the star-forming main sequence ($\Delta$MS, see Section 2). These parameters were chosen as they either related to morphology (which has been shown to relate to the scatter of the Fall relation; Romanowsky & Fall 2012; Cortese et al. 2016; Hardwick et al. 2022) or are an aspect of the model that we wanted to make sure was being encapsulated properly. We use the Spearman rank correlation coefficients ($\rho_{s}$) to quantify the correlation between each quantity and the offset from the theoretical line, as shown in Table 1 in descending order of correlation. We chose to use Spearman rank to measure the correlation as it does not make any assumptions on the type of correlation between the two quantities.

The strongest correlation with $\Delta f_{atm}$ is $\Delta$MS, which is known to be strongly linked to a galaxy’s atomic gas reservoir (e.g. Catinella et al., 2018; Janowiecki et al., 2020; Saintonge & Catinella, 2022). There is also a non-negligible correlation between $\Delta f_{atm}$ and $M_{\rm{bar,D}}$. This implies that the $f_{atm}-q$ plane may not provide the best representation for the link between AM and gas content in our data, and that xGASS galaxies may be better described by a model with different dependencies on $f_{atm}$ and $M_{\rm{bar,D}}$. Therefore, the next step is to use a more empirical approach and look for the best way parameterise the correlation between $j_{\rm{bar},D}$, $M_{\rm{bar},D}$ and $f_{atm}$.

By fitting a plane to specific AM, mass and gas fraction, we aim to quantify what functional form better describes our data and how this differs from the O16 model. We chose to fit two planes; one for the baryonic parameters ($j_{\rm{bar,D}}$, $M_{\rm{bar,D}}$ and $f_{atm}$) and one for stellar parameters ($j_{\star,D}$, $M_{\star,D}$ and $f_{\star,atm}=1.35M_{HI}/M_{\star,D}$). Although the majority of this work is focused on baryonic parameters, we also chose to investigate this parameter space using just stellar properties as, a) the stellar parameters have been more extensively studied in the literature (e.g., Mancera Piña et al., 2021b; Hardwick et al., 2022) and b) they allow us to check if any of our conclusions are affected by the assumptions made to estimate the baryonic AM. We fit these 3-dimensional planes using the following formula;

where $i$ is either disc baryons or stars ($\star$). The best-fitting parameters to equation 9 are shown in Table 2.

For comparison, when the O16 relation is converted into this formalism, the free parameters of equation 9 would be $\alpha=1$, $\beta=1/1.12\approx 0.89$ and $\gamma=\log_{10}(G/\,[\rm{kpc}\,($ km s^-1$)^{2}\,\rm{M}_{\odot}^{-1}])-(1/1.12)\log_{10}(2.5)-\log_{10}(\sigma/[$ km s^-1$])\approx-6.72$ . Table 2 shows that, as hinted in the previous section, the best-fit parameters to our data are significantly different from those predicted by the O16 model. In particular, there is a slightly weaker dependence on the mass ($\alpha$) term, and a significantly weaker dependence on the gas fraction ($\beta$) term, which is approximately half that expected by the O16 model. The empirical fit also has a smaller scatter (0.15 dex vertical scatter in $j_{\rm{bar}}$, compared to 0.23 dex vertical scatter observed for our data in the O16 projection).

We obtain very similar results to that of Mancera Piña et al. (2021b) who, as mentioned in Section 1, also fit a 3D plane to baryonic AM, mass and atomic gas fraction. The coefficients of our planes are consistent within 3$\sigma$, despite their work focusing on a small sample (N = 157) of gas-rich disc galaxies. The only major difference between the results of this work and that of Mancera Piña et al. (2021b) was that they are very confident in their error analysis and quoted intrinsic scatter of their plane (which did not include measurement errors). Conversely, as we believe that we are not able to estimate proper errors for each individual galaxy in our sample, we chose to give $\sigma$ which includes both measurement and intrinsic scatter of the plane. If we assume that our error estimates are correct, we find enticing evidence for nearly zero intrinsic scatter in our plane, i.e., a plane much tighter than what found in Mancera Piña et al. (2021b). However, as mentioned above, these errors are only approximations and we advice caution in over-interpreting this result.

Figure 4 shows the stellar (panel a) and baryonic (panel b) relations in the mass - AM plane. Projections of the best-fitting 3D relations from Table 2 are shown for fixed gas fractions as coloured lines. To easily compare, points are colour coded in the same way (i.e., by their gas fractions). The 2D best-fit Fall relation is shown as a thin black line (with the vertical 1$\sigma$ scatter as black dashed lines).

For comparison, Figure 5 shows the same baryonic Fall relation presented in 4b, but this time the lines of constant gas fraction correspond to the prediction from the O16 model. The comparison between Figure 4b and 5 helps visualising the differences between the O16 model and the best-fit empirical relation. Specifically, as already mentioned, the empirical relation has a weaker dependence on atomic gas fraction,which can be seen by the larger separation of the lines of constant gas fraction in Figure 5. The weaker dependence on baryonic mass is also evident, as the lines at fixed gas fraction have a shallower slope than the O16 projection.

A possible issue with the results presented so far is that xGASS was selected to have a nearly flat stellar distribution (Catinella et al., 2018; Hardwick et al., 2022), and the oversampling of high mass galaxies may affect our best-fitting results. To test this, we also fit the empirical relation including weights for each galaxy. Following Catinella et al. (2018), these weights were calculated to recover the $M_{\star}$ distribution of a volume-limited sample using the local stellar mass function as defined in Baldry et al. (2012). We find that, even after including the weights, our best-fit parameters are within error of the values presented in Table 2. Therefore, for simplicity, we only present the values obtained without the weights.

Interestingly, the scatter of the empirical relation still shows a weak correlation with $\Delta$MS, with a Spearman rank correlation coefficient of $\rho_{s}=0.25\pm 0.05$. Thus, it is worth exploring whether such residual dependence needs to be included in our fit. We trialled fitting a 4D plane, with $\Delta$MS as a fourth parameter and found that the dependence on $f_{atm}$ and $M_{\rm{bar},D}$ was largely unchanged when compared to the 3D plane, with a very small dependence ($-0.08\pm 0.02$) given to the $\Delta$MS term. In addition to this, the standard deviation between the 4D and 3D fits was unchanged. We conclude that any third-order dependence on $\Delta$MS is negligible or at least within the scatter of the 3D fit presented here.

The choice to use a planar fit for our empirical relation was motivated by the aim of keeping a similar approach to the one adopted by O16, but we also tested for non-linear trends in the residuals of the empirical fit, to see if this choice was sensible. There is a slight curvature in the residuals, when binned in $f_{atm}$, as is shown in Figure 6. This is a potential indication that the 3D planar fit may be too restrictive. However, the curvature is very weak and would likely not change the fit considerably if fit with a curved surface. Therefore, while we cannot exclude a curvature in the plane, our sample size is too small to justify the use of a 3D curved surface to parametrise the $j_{\rm{bar}}$, $M_{\rm{bar}}$, $f_{atm}$ plane.