© The Authors 2024

1. Introduction

The tight correlations found between the cold molecular gas content of galaxies and other physical galaxy properties, such as stellar mass (M_*) and star formation rate (SFR) are powerful and widely used tools in the study of galaxy evolution. These quantities are inextricably linked because stars form in clouds of cold molecular gas. However, the exact physical processes governing the interplay of cloud collapse and feedback are not fully understood and may vary among galaxy populations (e.g., Sun et al. 2020). Scaling relations have the potential to reveal any systematic differences across regimes in M_*, SFR, and environment, but this requires data covering the associated parameter ranges. Due to observational challenges, molecular line data in low-mass, low-metallicity galaxies are lacking. As a consequence, scaling relations have so far been largely based on massive (M_* > 10¹⁰ M_⊙), metal-rich galaxies (see Saintonge & Catinella 2022 for a review).

Recent studies have made efforts to push observations of cold molecular gas to lower masses with statistically significant sample sizes (Cicone et al. 2017; Saintonge et al. 2017; Wylezalek et al. 2022), but there is no clear picture of the behavior of molecular gas scaling relations in the low-mass regime (M_* < 10⁹ M_⊙) yet. Particularly in the case of intensive properties (i.e., those that do not directly depend on the size of the system), such as the molecular gas fraction (f_mol = M_mol/M_*), the results reported in the literature differ significantly. Some recent studies find f_mol to be anticorrelated with stellar mass following a simple power-law across the entire mass range studied (e.g., Hunt et al. 2020), while others report a shallow or flat relation at lower masses (e.g., Calette et al. 2018; Saintonge et al. 2017) with a downturn above a certain stellar mass. This feature is sometimes attributed to the larger fraction of passive galaxies at higher masses (Jiang et al. 2015), but it appears even in samples selected to exclude passive galaxies (e.g., Saintonge et al. 2017; Calette et al. 2018). The situation is similar for the molecular gas depletion time (τ_mol = M_mol/SFR), which is a measure of how rapidly galaxies convert their molecular gas into stars. Some authors report that τ_mol is independent of M_* for main-sequence galaxies (e.g., Accurso et al. 2017; Boselli et al. 2014), but others find τ_mol and M_* to correlate, at least above M_* = 10⁹ M_⊙ (e.g., Hunt et al. 2020; Saintonge et al. 2017). A lack of relation between τ_mol and M_* would indicate that star formation proceeds at a universal rate in galaxies, while a correlation has been argued to be a consequence of global galaxy properties affecting processes on cloud scales (Sun et al. 2020; Saintonge & Catinella 2022).

These differences hint at the uncertainties associated with measuring the molecular gas content of galaxies and using it to constrain scaling relations. One such uncertainty is due to the fact that the H₂ molecule, which makes up the majority of the cold molecular gas, has no permanent dipole. As a consequence, it is not readily excited to emit radiation at the densities and temperatures found in molecular clouds, providing no direct tracer for the majority of the gas. Instead, emission arising from the J = 1 to J = 0 rotational transition of carbon monoxide (CO), the second most abundant molecule, is frequently used as a proxy for the molecular gas mass. However, the abundance of CO in the molecular gas depends on the local metallicity and radiation field; the conversion factor between observed CO(1−0) luminosity and inferred molecular gas mass varies with the physical properties of the interstellar medium (ISM). Many empirical and theoretical prescriptions for this conversion factor (${\alpha }_{\text{CO}}\equiv M_{\text{mol}}/L_{\text{CO}(\text{1-0})}^{\prime }$) have been proposed (Accurso et al. 2017; Wolfire et al. 2010; Glover & Low 2011; Bolatto et al. 2013), but there remains a large uncertainty in choosing a specific model, which goes some way to explain the different scaling relations reported in the literature.

In low-mass, low-metallicity galaxies, shielding is less effective so that very little CO can form, increasing the fraction of so-called dark molecular gas (Wolfire et al. 2010). This makes it even more difficult to obtain reliable estimates of the molecular gas mass of such galaxies, and studies targeting galaxies independently of their gas fraction generally have low CO detection rates in this regime. This leads to samples either including only the most gas-rich low-mass objects, or a large fraction of non-detections depending on the survey strategy. The former has the potential to bias the sample toward a particular, non-representative class of objects, while the latter necessitates the use of analysis methods suitable to include upper limits derived from non-detections. These differences in turn influence the scaling relations inferred from the data.

The environment in which galaxies evolve may present an additional confounding factor affecting their scaling relations. Galaxies residing in large clusters have been found to adhere to different scaling relations than their isolated counterparts (e.g., Brown et al. 2021; Zabel et al. 2022). Most galaxies however reside in smaller groups, rather than the rare large clusters. The effect of this less extreme environment is likely to be more subtle, with lower gas temperatures in the halo and fewer close interactions between members. Nonetheless, examples of the group environment affecting a member galaxy have been found in both observations and simulations, raising the question whether the group environment systematically affects galaxy evolution (e.g., Lee et al. 2022; Zhou et al. 2022). Tidal forces and feedback mechanisms can contribute to the transfer of gas from the ISM of the group members to the intra-group region while simultaneously affecting its physical properties, as observed in the Stephan’s quintet (Appleton et al. 2013, 2023) and in the M81 triplet (de Blok et al. 2018). If this were the case in most galaxy groups, the environment would need to be accounted for as a factor in determining statistical scaling relations. Single-dish CO surveys in the past decade have increased the available sample size to a level where it becomes feasible to study populations separately based on their environment, but this possibility has not yet been exploited to its full extent.

In this study, we present improved constraints on the scaling relations of cold molecular gas in local star-forming galaxies down to M_* = 10^8.5 M_⊙ based on a reanalysis of the APEX low-redshift legacy survey for molecular gas (ALLSMOG; Cicone et al. 2017) and extended CO Legacy Database for GASS (xCOLDGASS; Saintonge et al. 2017) surveys, which include rigorous upper limits for their non-detections. Additional homogeneous SFR, M_*, metallicity, redshift, and environment data are available from catalogs compiled for the SDSS parent sample, from which both of these surveys were originally drawn. Cicone et al. (2017) demonstrated that including the low-M_* sample from ALLSMOG can significantly affect the slopes of scaling relations, but did not take the extra steps to convert CO luminosities into molecular gas masses, and so did not specifically delve into the f_mol vs. M_* or τ_mol vs. M_* scaling relations that we aim to explore here. In addition, we exploit the size of the sample and the availability of the necessary ancillary data products to investigate the impact of the group environment on molecular gas scaling relations which, to our knowledge, has never been done. This is necessary groundwork to inform future studies that may want to control for the environment variable when investigating other effects, such as various modes of feedback, on the cold gas scaling relations.

This paper is structured as follows. We describe the sample and its underlying selection criteria in Sect. 2 and the treatment of censored data in Sect. 3. The resulting scaling relations are presented and compared to previous studies of similar samples in Sect. 4. In Sect. 5 we explore the effect of the group environment. Finally, we present our conclusions and discuss their implications in Sect. 6. Throughout this paper we assume a flat ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹ for consistency with Saintonge et al. (2017) when deriving distance dependent quantities.

2. Sample selection and methods

We rely on observations of the CO(1−0) and CO(2−1) lines to infer the cold molecular gas properties used to constrain scaling relations. From the available archival CO data we chose our sample to:

1. be representative of the local star-forming galaxy population with 10^8.5 M_⊙ < M_* < 10^10.5 M_⊙;

2. have available data on galaxy environment;

3. have homogeneous ancillary data such as star formation rates (SFR) and stellar masses.

The ALLSMOG (Cicone et al. 2017) and xCOLDGASS (Saintonge et al. 2017) surveys are drawn from the Sloan Digital Sky Survey (SDSS) DR7 (Abazajian et al. 2009) which ensures criteria (2) and (3) are met, while their survey strategies are well suited to achieve criterion (1), as is outlined in the following. Section 2.4 further describes how we select star-forming galaxies based on their star formation rates relative to the main sequence in the local universe. Both surveys also target galaxies with available HI 21 cm observations, which enables us to include the scaling relations of cool neutral hydrogen gas in our analysis.

2.1. xCOLDGASS

xCOLDGASS (Saintonge et al. 2017) is comprised of 532 local galaxies with M_* > 10⁹ M_⊙ and redshift 0.01 < z < 0.05 observed in CO(1−0) with the IRAM 30-m telescope. Of these, 366 were observed in the original COLDGASS survey and have M_* > 10¹⁰ M_⊙ and redshift 0.025 < z < 0.05, while the remaining 166 were targeted in an effort to extend the survey into the range of 10⁹ M_⊙ < M_* < 10¹⁰ M_⊙ and have 0.01 < z < 0.02. Targets were selected randomly from the SDSS parent sample based only on the mass and redshift criteria (see Saintonge et al. 2011, for details), resulting in a sample representative of local galaxies within the mass range. xCOLDGASS did not specifically target star forming galaxies so we apply an additional selection criterion for this study, which we describe in Sect. 2.3.

2.2. ALLSMOG

ALLSMOG (Cicone et al. 2017, see also Bothwell et al. 2014 for a previous incomplete release) contains a total of 97 local star-forming galaxies, of which 88 were observed in CO(2−1) with the APEX telescope and 9 more in CO(1−0) and CO(2−1) with the IRAM 30-m telescope. Their targets were selected from the SDSS in the same randomized manner as employed in xCOLDGASS with 10^8.5 M_⊙ < M_* < 10¹⁰ M_⊙, and 0.01 < z < 0.03. Additionally, they select targets whose gas-phase metallicity, as initially reported in the MPA-JHU catalog computed according to the Tremonti et al. (2004) calibration, is 12 + log(O/H) > 8.5 for targets observed with APEX and 12 + log(O/H) > 8.3 for observations with the more sensitive IRAM 30-m telescope. However, in this work, metallicities are recalculated homogeneously for all of the sample for which this computation is possible (see details in Sect. 2.3) using the more recent calibration of the [OIII]–[NII] strong line index by Curti et al. (2017, hereafter O3N2 C17), so the aforementioned limits do not directly translate. Since the metallicity limit applied to ALLSMOG could potentially bias our sample, we apply a weighting in our analysis to account for this effect, as discussed in more detail in Sect. 3.

2.3. BPT classification

In order to select the targets for which the gas-phase metallicity (needed to calculate the α_CO conversion factor) can be directly computed using strong line calibration methods, we employ the Baldwin, Phillips & Terlevich (BPT) diagram method (Baldwin et al. 1981). The BPT classification indicates the dominant excitation mechanism in the region covered by the optical spectra of the galaxies in the sample. The left panel of Fig. 1 shows the BPT diagram of the combined ALLSMOG and xCOLDGASS sample of 629 galaxies. We calculated the emission line ratios reported in Fig. 1 consistently for the entire sample using a full fit of the optical spectra performed with the GANDALF code (Sarzi et al. 2006), which fits stellar continuum and emission line templates simultaneously.

Fig. 1. BPT classification of the sample. Left: BPT diagram of the full sample. The dashed line is the demarcation between pure SF galaxies and AGN or composite galaxies according to Kauffmann et al. (2003). The solid line is the theoretical limit for the ionization state produced by star formation according to Kewley et al. (2001), splitting the sample further into AGN and composite systems. Right: distribution of stellar masses in the SF, AGN, composite and low-S/N subsamples.

There are 170 galaxies with an amplitude-to-noise ratio S/N < 3 in at least one of the relevant emission lines (Hα, Hβ, [OIII], and [NII]) which are not classified in the BPT diagram, and so are not shown in the left panel of Fig. 1. We designate these as “low-S/N” and show their stellar mass distribution in the right panel of Fig. 1. The remaining galaxies with S/N > 3 in all lines are then classified based on their position in the BPT diagram, following Kauffmann et al. (2003) and Kewley et al. (2001). The sample is thus split into star forming (SF) galaxies, active galactic nuclei (AGN), and composite systems, where the SF subsample covers the entire range of 8.5 < log(M_*/M_⊙) < 11.5, while the AGN subsample has a much higher average stellar mass and only covers 9.5 < log(M_*/M_⊙) < 11.5, as shown in the right panel of Fig. 1.

The optical spectra used in this study sample a 3″-diameter region, which is the size of the SDSS fiber. Since our targets are low-redshift galaxies, they are quite extended in the sky and the fiber will only probe the very central region of the galaxy. The BPT classification is therefore indicative of the dominant excitation mechanism in this central region only. It is possible that we classify galaxies as AGN dominated or composite systems even though globally they may be mainly star forming. However, if a galaxy is star forming in the central region, star formation is likely to be the dominant excitation mechanism globally as well, since AGN are generally located centrally. Thus, we can be reasonably certain that galaxies classified as star-forming in the BPT diagram are not affected by current AGN activity. This is important for the metallicity calculations outlined in Sect. 2.5, which are based on emission line ratios that may be affected by the presence of an AGN. Indeed, the purpose of the BPT classification discussed here is only to select the galaxies for which we can reliably apply strong line calibration recipes to compute the gas phase metallicity.

2.4. “Star-forming” classification

Figure 2 shows the SFR of the 625 ALLSMOG and xCOLDGASS galaxies for which we have a reliable SFR measurement plotted against their stellar mass. Both quantities are obtained from the MPA-JHU catalog¹ for SDSS DR8, thus ensuring a homogeneous methodology across the entire sample. Also shown are two empirical relations for the star-forming main sequence of galaxies (SFMS) reported by Renzini & Peng (2015) and Saintonge & Catinella (2022) with an added 0.4 dex scatter.

Fig. 2. SFR vs. stellar mass for the 625 ALLSMOG + xCOLDGASS galaxies that have a reliable SFR in the MPA-JHU catalog. The color scale corresponds to the gas phase metallicity, which was calculated from the O3N2 line index (Curti et al. 2017) for galaxies classified as SF or composite in the BPT diagram and according to the mass-metallicity relation (Curti et al. 2020) for those classified AGN and low-S/N. The blue dash-dotted line shows the empirical relation for the star forming main sequence found by Renzini & Peng (2015) including a 0.4 dex scatter (blue shaded). The green dashed line shows the star forming main sequence following Saintonge & Catinella (2022).

As we aim to study specifically star-forming galaxies in the local universe, we need to apply an additional selection to this sample. Applying a uniform selection across the entire sample also ensures that we homogenize the ALLSMOG (which specifically targeted star-forming galaxies) and xCOLDGASS (which did not) data. We chose a criterion based on the offset from the main sequence of star-forming galaxies following the definition of Renzini & Peng (2015). We calculate the offset from the main sequence as

$$\begin{array}{c}\hfill \mathrm{\Delta }(\mathrm{MS})={\mathrm{sSFR}}_{\mathrm{measured}}/{\mathrm{sSFR}}_{\mathrm{MS}},\end{array}$$(1)

where sSFR_measured = SFR_measured/M_{*, measured} is the measured specific star formation rate and sSFR_MS is the specific star formation rate predicted for the object by the Renzini & Peng (2015) relation based on its stellar mass. In our analysis of molecular gas scaling relations, we exclude galaxies with log Δ(MS) <  − 0.4 dex, to avoid contaminating our sample with quenched passive galaxies.

For the remainder of this work we refer to the sample selected as outlined above simply as the “star-forming sample”. When using a sample based on different selection criteria this is stated explicitly. The star-forming sample contains 353 galaxies, 283 of which are detected in at least one of the CO(1−0) and CO(2−1) emission lines and 70 of which are non-detections in CO.

Starburst galaxies, which lie above the main sequence in the M_*–SFR plane, are included in this selection. The possible additional dependence of α_CO on the dynamic state of the gas is properly taken into account when calculating molecular gas masses for the starburst galaxies in our sample, as described in Sect. 2.5.

There are a number of different approaches for selecting star-forming galaxies, which will generally result in different slopes and shapes of the main sequence relation (see Pearson et al. 2023 for an overview). In Fig. 2 we show both a linear main sequence relations (Renzini & Peng 2015) and one with a high-mass turnover (Saintonge & Catinella 2022). This illustrates that while the composition of our final sample will depend on the main sequence definition we adopt here, the key differences arise in the high-mass regime (M > 10^10.5 M_⊙). The low-mass regime, on which we focus in our analysis, is unlikely to be affected significantly by this choice. We test this assumption by comparing between the results obtained for the Renzini & Peng (2015) and Saintonge & Catinella (2022) main sequence samples respectively for key scaling relations later on. We further explore how our choice of selection criteria affects our results by comparing it to other common approaches in Appendix A.

2.5. Metallicity and CO-to-H₂ conversion factor

Obtaining a galaxy’s molecular gas mass from its integrated CO(1−0) line luminosity is not straightforward. Most of the molecular gas mass is contained in H₂, which we do not observe directly, while CO only makes up a small fraction of the total mass. More importantly, the abundance of CO depends on the local gas chemistry and radiation field, so that a conversion factor obtained from a local cloud in the Milky Way need not apply to another galaxy. The main parameter governing CO abundance is thought to be the gas-phase metallicity, with both observations and theoretical predictions pointing at a lower CO abundance in low-metallicity galaxies (e.g., Wolfire et al. 2010; Madden et al. 2020). The ratio of CO to H₂ gas changes with metallicity in this picture, necessitating a metallicity dependent α_CO conversion factor.

The color scale in Fig. 2 corresponds to the gas-phase metallicity of our sample in units of 12 + log(O/H), which, for galaxies that were classified as star-forming or composite in the BPT diagram, was calculated from measured fluxes of the Balmer α and β lines and the [NII]λ6583 and [OIII]λ5007 forbidden lines following Curti et al. (2017). Galaxies classified as AGN instead have their metallicity assigned according to the prediction from the mass–metallicity-relation as defined by Curti et al. (2020), because the strong-line diagnostics presented in Curti et al. (2017) are only calibrated for star-forming galaxies. Most AGN in our sample lie above M_* = 10¹⁰ M_⊙, where the mass-metallicity relation flattens, so that we do not expect our assumptions on their metallicity to add significant additional uncertainty to our molecular gas scaling relations.

To determine the α_CO value we use the prescription presented in Accurso et al. (2017):

$$\begin{array}{cc}\hfill log{\alpha }_{\mathrm{CO}}=& 14.752-1.623·[12+log(\mathrm{O}/\mathrm{H})]\hfill \\ \hfill & +0.062·log\mathrm{\Delta }(\mathrm{MS}).\hfill \end{array}$$(2)

Here, the gas phase metallicity in units of 12 + log(O/H) is calculated according to the calibration by Pettini & Pagel (2004) using the same line fluxes as outlined above (hereafter O3N2 PP04), in order to be fully consistent with Accurso et al. (2017). For AGN, we use the mass-metallicity relation defined by Kewley & Ellison (2008), which uses the PP04 calibration. Equation (2) contains a secondary dependence on the local radiation field through the parameter Δ(MS) = SFR_measured/SFR_MS, the distance in SFR off the star forming main sequence (SFMS). In this instance, SFR_MS corresponds to the redshift dependent analytic definition of the main sequence by Whitaker et al. (2012). Both metallicity and log Δ(MS) measurements come with associated uncertainties, that we propagate to estimate the uncertainty in the resulting conversion factor. The median uncertainty on the metallicity values, derived from the uncertainties of the line ratio measurements, is 0.025 in units of 12 + log(O/H) (O3N2 PP04). For log Δ(MS) the median error is 0.16 dex.

Accurso et al. (2017) give the parameter ranges in which their empirical α_CO is constrained as 7.9 < 12 + log(O/H) < 8.8 in metallicity (O3N2 PP04) and −0.8 < log Δ(MS) < 1.3 Our sample spans 7.94 < 12 + log(O/H) < 8.88 (O3N2 PP04) and −0.49 < log Δ(MS) < 1.32, so the majority of our galaxies lie in the well constrained regime. The full distribution of our star-forming sample in metallicity, log Δ(MS), and the resulting α_CO is shown in Fig. 3. For those galaxies with 12 + log(O/H) > 8.8 we follow the suggestions of Accurso et al. (2017) and calculate α_CO using 12 + log(O/H) = 8.8. For starburst galaxies with log Δ(MS) > 1.3, Accurso et al. (2017) recommend a lower value of α_CO, because dynamical effects dominate over metallicity in setting the α_CO in this regime. We adopt a value of α_CO = (1.7 ± 0.5) (K km s⁻¹ pc²)⁻¹ based on a recent study of local ultra luminous infrared galaxies (ULIRGs) by Montoya Arroyave et al. (2023), which is consistent with recent estimates of α_CO values in other local starbursts. This lower value is used to account for the effect of turbulence and high velocity dispersion in starbursts lowering the optical depth of the CO emission and so boosting the luminosity of CO line emission. The conversion factors thus obtained for the galaxies in our star-forming sample range from 1.7 (K km s⁻¹ pc²)⁻¹ (2.75 (K km s⁻¹ pc²)⁻¹ for non starbursts) to 82.2 (K km s⁻¹ pc²)⁻¹ (see Fig. 3). All values quoted here include a factor of 1.36 to account for the presence of helium in the molecular gas, based on the helium mass fraction measured in the Milky Way (Asplund et al. 2009).

Fig. 3. Distribution of our sample in gas phase metallicity (left), distance off the main sequence (center), and αCO (right). Shaded areas indicate the parameter space where Accurso et al. (2017) suggest adapting their prescription for αCO.

Not all galaxies in our sample are observed in the CO(1−0) line. For the ALLSMOG galaxies with only CO(2−1) observations, we assume $r_{21}\equiv L_{\text{CO}(\text{2-1})}^{\prime }/L_{\text{CO}(\text{1-0})}^{\prime }=0.7\pm 0.2$, following Leroy et al. (2022) to obtain CO(1−0) luminosities. We use aperture corrected values for all CO(1−0) and CO(2−1) luminosities following the corrections applied in the original ALLSMOG (Cicone et al. 2017) and xCOLDGASS (Saintonge et al. 2017) surveys. Finally, we calculate the molecular gas masses for our sample as

$$\begin{array}{c}\hfill M_{\mathrm{mol}}={\alpha }_{\mathrm{CO}}(12+log(\mathrm{O}/\mathrm{H}),\mathrm{\Delta }(\mathrm{MS}))·L_{\mathrm{CO}}^{\prime },\end{array}$$(3)

where $L_{\text{CO}}^{\prime }$ is the aperture corrected CO(1−0) line luminosity and we are using the α_CO prescription from Eq. (2). To ensure the robustness of our results regarding the chosen prescription for α_CO, we test two additional models. The results are described in Appendix B.

2.6. Environment classification

We use a galaxy group catalog² by Tinker (2021) to classify the environment of each galaxy in our sample. The catalog is constructed using a halo-based group finder algorithm which assigns galaxies to halos based on halo mass and distance in redshift space. Galaxies are classified either as satellites or centrals in the catalog based on these criteria. For a more detailed description of the algorithm see Tinker (2021) and references therein. We add a third class, isolated galaxies, which are assigned to a group in which they are the only member.

In our sample of 353 star-forming local galaxies we identify 38 centrals, 42 satellites, and 247 isolated galaxies. There are also 26 galaxies which are not contained in the group catalog, because it is not based on the original SDSS data but on the NYU value added galaxy catalog (Blanton et al. 2005). The groups in our sample have halo masses of 11.0 < log(M_h/M_⊙) < 14.9. Galaxy groups are usually defined by a halo mass of 13 < log(M_h/M_⊙) < 14 as opposed to clusters with log(M_h/M_⊙) > 14. In our sample only 5 galaxies reside in halos with log(M_h/M_⊙) > 14, so we probe mainly group environments and isolated objects.

3. Statistical methods

Scaling relations are frequently expressed as power-law functions of the form

$$\begin{array}{c}\hfill Y=A{X}^b,\end{array}$$(4)

making it convenient to fit the data in log-log space where the power-law is transformed into a linear function with slope b and intercept log(A). To fit this model to our data we use a Bayesian linear regression method implemented in the IDL procedure LINMIX_ERR (Kelly 2007). This procedure is capable of accounting for measurement errors in both the independent and dependent variables as well as left-censored data (i.e., upper limits) in the independent variable, making it ideally suited for our purposes.

To account for possible biases due to the somewhat skewed distribution of our sample in metallicity at low masses (see Sect. 2), we apply a weighting to our fitting procedure. To obtain the appropriate weights we compare the metallicity distribution of our sample to that of the sample of all star-forming galaxies selected from the MPA-JHU catalog (using the same criteria we applied to our sample). The weights are calculated for bins of 0.2 dex in stellar mass and 0.1 dex in gas-phase metallicity. It should be noted that this method is somewhat limited by the range of metallicities in our sample, since we cannot use weights to account for the lowest metallicities, which our sample does not cover. However, our sample covers a range of metallicities that contains at least 90%+ of galaxies in the MPA-JHU sample in each mass bin, giving us reasonable certainty that the weighting procedure can account for most of the potential bias in metallicity. We discuss this weighting procedure in more detail in Appendix C.

The LINMIX_ERR procedure calculates the linear correlation coefficient ρ_corr between the variables in addition to returning a best-fit linear function. As an additional measure for statistical independence in our data we adopt an estimator for Kendall’s tau statistic proposed by Oakes (2008) that takes into account the presence of upper limits. To compute this estimator we analyze pairs of bivariate data (X_i, Y_i) and (X_j, Y_j), where i ≠ j. Such a pair is said to be concordant if (X_i − X_j)⋅(Y_i − Y_j) > 0 and discordant if (X_i − X_j)⋅(Y_i − Y_j) < 0. The estimator is then computed as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{Kendall}}=\frac{C-D}{C+D},\end{array}$$(5)

where C is the number of definite concordances among the n(n − 1)/2 pairs of data and D is the number of definite discordances. Due to the presence of upper limits, some pairs cannot be determined to be discordant or concordant. This fact is taken into account by renormalizing the original Kendall’s tau (for details see Oakes 2008). The resulting measure can range in value from −1 to 1 with values close to zero indicating statistical independence of the two variables.

To test for nonlinear trends, we additionally bin the data and employ survival analysis in the form of the nonparametric Kaplan-Meier estimator (Kaplan & Meier 1958) to calculate the median for each bin. We follow the method outlined in Lee & Wang (2003), which is well suited to retrieve a median even for bins with detection fractions as low as 20% (e.g., Calette et al. 2018), albeit with large uncertainties. Consequently, we are able to estimate median cold gas masses and derived quantities even for the low-mass galaxies, where our data are dominated by CO non-detections.

We choose bin sizes such that each bin is populated by roughly the same number of galaxies to ensure that the calculated median is not dominated by random fluctuations due to low number statistics. As a consequence, the bins at both ends of the mass distribution (Fig. 1) are considerably larger than at the centre.

4. Gas scaling relations

Having established that our sample of 353 galaxies is representative of the local star-forming population with 10^8.5 M_⊙ < M_* < 10^10.5 M_⊙, we use it to derive scaling relations of M_mol, f_mol, and τ_mol for this population. We apply the methods described in Sect. 3 to ensure that non-detections are accounted for in our scaling relations and analyze their impact by comparing the results to the literature.

4.1. Molecular gas mass

We begin by investigating the relationships between the molecular gas mass, the star formation rate and the stellar mass. The former scaling relation is known as the integrated Kennicutt–Schmidt (KS) law (Schmidt 1959; Kennicutt 1989) while the latter is often referred to as the molecular gas main sequence (MGMS) and both are frequently used in galaxy evolution studies, allowing us to compare our results with previous work (e.g., Chown et al. 2022; Baker et al. 2022).

The molecular gas mass in our sample is strongly correlated with SFR and M_* as can be seen in Fig. 4. This is expected, since all these properties are extensive (i.e., they scale with overall galaxy mass). More massive galaxies will typically contain larger amounts of cold gas which leads to the formation of more stars. While more massive galaxies are also more often passive quenched systems, with smaller cold gas reservoirs, we exclude passive galaxies from our sample explicitly to remove their effect on the scaling relations.

Fig. 4. Scaling relations of the molecular gas mass. Left: molecular gas mass vs. SFR for the star-forming sample. For comparison we show the integrated inverse Schmidt–Kennicutt relation (blue dash-dotted) and corresponding intrinsic scatter reported by Chown et al. (2022). Right: molecular gas mass vs. M*, with the best-fit MGMS from Chown et al. (2022) (blue dash-dotted) and Baker et al. (2022) with only detections (purple dotted) and detections and upper limits (green dashed). The color scale corresponds to the adopted CO-to-H2 conversion factor. In both panels the solid gray line shows the best fit from the Bayesian linear regression to the unbinned data with the shaded area indicating the 95% confidence band. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25−75 percentile interval (in y). Open symbols with downward arrows denote upper limits. The navy area of the color bar indicates the specific value of the conversion factor used for extreme starbursts as described in Sect. 2.5.

For both the integrated SK-relation and the MGMS, the best-fit slope for our sample from the Bayesian linear regression is close to unity (see Table 1 and Fig. 4). We note that the integrated SK relation is more commonly reported with M_mol on the x and SFR on the y-axis, but we chose the assignment of the axes according to how the variables are treated in the fitting procedure. For comparison, we show empirical integrated SK and MGMS relations from recent studies (Chown et al. 2022; Baker et al. 2022) in Fig. 4. This illustrates that our results are consistent with previous works. However, it should be noted that both Chown et al. (2022) and Baker et al. (2022) use a constant CO-to-H₂ conversion factor and their samples consist of galaxies with log(M_*/M_⊙) > 9. We test how the assumption of a constant as opposed to a metallicity dependent conversion factor affects our results in Appendix B.

4.2. Molecular gas fraction

Given the expected strong scaling between extensive galaxy properties in our sample, the next step is to look at the corresponding intensive properties such as the gas fraction (f_mol = M_mol/M_*). Dividing by stellar mass removes the primary scaling and provides insight into the physical processes behind the formation of cold gas and stars, and their dependence on galaxy stellar mass. In Fig. 5, we plot f_mol against M_* for our sample and find a flat trend with no significant correlation as quantified by the linear correlation coefficient from the fitting procedure and the adopted estimator for Kendall’s rank correlation coefficient (see Table 1). The binned medians shown in Fig. 5 are also consistent with the fitted constant molecular gas fraction; The 25−75 percentile intervals derived from survival analysis overlap with the linear regression for all bins.

Fig. 5. Molecular gas mass fraction vs. stellar mass for the star-forming sample. The solid gray line shows the best fit from the Bayesian linear regression to the unbinned data with the shaded area indicating the 95% confidence band. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25−75 percentile interval (in y). Open symbols with downward arrows denote upper limits. The purple diamond shows the upper limit derived from spectral stacking of the 22 ALLSMOG CO(2−1) non-detections with masses below M* = 109 M⊙. The blue dash-dotted line shows the scaling relations reported by Calette et al. (2018) based on their sample of late-type galaxies.

While the Kaplan-Meier estimator is capable of retrieving a median even for the lowest mass bin, which is dominated by non-detections, the corresponding uncertainties become large. We tested the validity of the KM medians for the low-mass regime by stacking the spectra of the 22 ALLSMOG non-detections of CO(2−1) with masses below M_* = 10⁹ M_⊙ contained in our sample. Stacking did not result in a clear detection, but we computed an upper limit on the typical molecular gas fraction in this mass range using the same procedure as for individual sources following Cicone et al. (2017). This upper limit is shown as a purple diamond in Fig. 5 and is consistent with the KM median value for the lowest mass bin. The stacked spectrum and the procedure used to obtained it are shown in Appendix D. We note that since the stacking produces only an upper limit it does not rule out a downturn of the f_mol versus M_* relation below M_* ∼ 10⁹ M_⊙.

We have further tested the robustness of the KM medians against different binning choices (specifically, using 4, 5 and 6 bins instead of 9) and we confirm that each bin size delivers the same flat trend as seen here. The bin size shown here is the one that maximizes the number of bins while still including enough objects in each bin to enable the survival analysis.

We apply the same survival analysis method to the entire star-forming sample to determine the median gas fraction around which the local star-forming population is distributed, and find a median value of $log{f}_{\mathrm{mol}}=-0.{99}_{-0.19}^{+0.22}$ ($log{f}_{\mathrm{mol}}=-1.{01}_{-0.23}^{+0.24}$ when considering only galaxies with M_* < 10^10.5 M_⊙). Here the quoted errors correspond to the 25th and 75th percentiles.

A flat or nearly flat relation between f_mol and M_* in the stellar mass range 10⁹ M_⊙ < M_* < 10¹⁰ M_⊙ is well supported by previous findings (e.g., Saintonge et al. 2017; Calette et al. 2018). The picture is somewhat less clear at high masses (M_* ⪆ 10¹⁰ M_⊙), where reported scaling relations are sensitive to sample selection, and in the very low-mass regime (M_* < 10⁹ M_⊙). For the total xCOLDGASS sample, Saintonge et al. (2017) reported a mild negative trend between f_mol and M_* for main sequence galaxies with 10⁹ M_⊙ < M_* < 10^10.5 M_⊙ and a sharp drop-off at higher masses. The fact that we do not observe this drop off is a direct consequence of how our star-forming sample is selected. Saintonge et al. (2017) adopt a criterion based on the distance from the main sequence as defined in Saintonge et al. (2016) and updated in Saintonge & Catinella (2022) which flattens at high masses, as opposed to the linear main sequence derived in Renzini & Peng (2015) that we adopt here. We illustrate this Fig. 6, where we plot the same f_mol vs. M_* scaling relation for a sample selected using the offset from the main sequence as defined by Saintonge & Catinella (2022) and recover the same drop-off beyond M_* ∼ 10^10.5 M_⊙ they report, while the scaling relations at the low-mass end do not change significantly. We discuss the impact of the choice of main sequence definition in more detail in Appendix A, where we also show that the inclusion of starburst galaxies in our sample does not significantly change the molecular gas scaling relations we obtain.

Fig. 6. Same scaling relation as in Fig. 5 but for the star-forming plus starburst sample defined using the main sequence definition from Saintonge & Catinella (2022).

Calette et al. (2018) adopt a morphology-based criterion to divide their sample into late and early types, so that neither of their subsamples is directly comparable to our star-forming sample. They obtain a double power-law scaling relation from fitting to the means of binned data calculated using Kaplan-Meier estimators. At low masses (M_* ≪ 1.74 × 10⁹ M_⊙) their best fit f_mol vs. M_* slope for both late and early types agrees with what we find here, but at higher masses it is significantly steeper (see Fig. 5). This is not surprising since at low masses, star-forming galaxies dominate the population, while at higher masses, the morphology based selection is likely to include a significant fraction of passive quenched galaxies. When we adopt a main-sequence definition with a high mass turnover, we find a similar steepening trend as reported by Calette et al. (2018) as can be seen in Fig. 6.

In the low-mass regime (M_* < 10⁹ M_⊙), observations of molecular gas are much more challenging, and scaling relations are less well constrained as a consequence. One recent study (Hunt et al. 2020) reports a single power-law scaling relation with a negative slope for f_mol across a wide mass range of stellar masses (10⁷ M_⊙ < M_* < 10¹¹ M_⊙). This stands in contrast to our results, which indicate that the flat trend observed in the range 10⁹ M_⊙ < M_* < 10¹⁰ M_⊙ extends down to M_* = 10^8.5 M_⊙. However, while our sample includes non-detections as upper limits, the data analyzed in Hunt et al. (2020) is limited to CO detections, which is likely to introduce a bias towards extraordinarily gas-rich objects at the low-mass end (10⁷ M_⊙ < M_* < 10⁹ M_⊙). This could explain the increasing molecular gas fraction they observe at these masses.

The color scale in Fig. 5 additionally illustrates that the scatter in the molecular gas fraction is correlated with specific star formation rate (sSFR). Galaxies found above the best-fit relation tend to have a higher sSFR, while those below tend to be less star forming. A positive correlation between f_mol and sSFR is also reported in the literature (Saintonge et al. 2017, 2011) and expected, as this relation can be seen as a variation of the SK-law (see Saintonge & Catinella 2022 for a detailed review).

The physical origin of the scatter in f_mol is not clear however. The environment in which galaxies evolve could be a factor in this, as it is a parameter independent of galaxy mass. We explore this possibility in Sect. 5.

4.3. Molecular gas depletion time

Next to the overall availability of molecular gas, the rate of star formation also depends on the efficiency with which clouds convert gas mass into stars. While it is difficult to determine the true efficiency (i.e., how much of the available molecular gas will ultimately be converted into stars), we can measure the timescale of star formation from molecular gas, also called the molecular depletion time (τ_mol = M_mol/SFR), which is the reciprocal of the star formation efficiency (SFE = SFR/M_mol).

Figure 7 shows τ_mol plotted against stellar mass for our star-forming sample. Unlike f_mol, τ_mol shows some residual dependence on stellar mass despite being an intensive galaxy property. While this correlation is much weaker than the M_mol vs. M_* and M_mol vs. SFR relations (see Table 1), its best-fit slope is non-zero at the 3σ level. This is corroborated by the Kaplan-Meier estimated medians and 25–75 percentile intervals, which overlap with the fit to the unbinned data at all masses. Such a mass scaling indicates that higher mass galaxies consume their molecular gas reservoir through star formation less rapidly than their low-mass counterparts, despite their similar f_mol values.

Fig. 7. Molecular gas depletion time τmol = Mmol/SFR vs. stellar mass. The solid gray line shows the best fit from the Bayesian linear regression to the unbinned data with the shaded area indicating the 95% confidence band. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25−75 percentile interval (in y). Open symbols with downward arrows denote upper limits.

There seems to be no clear consensus in the literature concerning a possible correlation between molecular gas depletion time and stellar mass. Some authors report a mild dependence of τ_mol on stellar mass in the mass range we investigate here (Boselli et al. 2014; Saintonge et al. 2017; Hunt et al. 2020), while others see no statistically significant correlation (Accurso et al. 2017). A possible explanation for the dependence of τ_mol on stellar mass are properties of molecular clouds being affected by global galaxy properties (Sun et al. 2020; Rosolowsky et al. 2021) resulting in galaxies following different Kennicutt–Schmidt relations (Ellison et al. 2021). In particular the increased presence of stellar bulges in high mass galaxies has been linked to higher depletion times (Saintonge et al. 2011; Davis et al. 2014) and it has been suggested that their gravitational potential could stabilize the gas disk against fragmentation and collapse into star-forming clumps (Martig et al. 2009). This could explain how specific star formation rate decreases with increasing stellar mass, despite the flat molecular gas fraction we find here.

4.4. Neutral atomic gas mass

Since almost all galaxies in our sample have constraints on their neutral atomic gas mass from archival HI 21 cm data, we can compare the mass scaling of the neutral atomic gas fraction f_HI = M_HI/M_* to that of the molecular gas fraction f_mol. Figure 8 shows f_HI plotted against stellar mass with a power-law fit from Bayesian linear regression and binned medians from survival analysis analogous to the f_mol plots above. We find that, contrary to f_mol, f_HI is anticorrelated with stellar mass, and that a single power-law is a good fit to our sample (see Table 1).

Fig. 8. Neutral atomic gas mass fraction vs. stellar mass for the star-forming sample. The solid gray line shows the best fit from the Bayesian linear regression to the unbinned data with the shaded area indicating the 95% confidence band. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25−75 percentile interval (in y). Open symbols with downward arrows denote upper limits.

The anticorrelation of f_HI and M_* we find generally aligns with results from previous studies of similar samples (Saintonge et al. 2017; Calette et al. 2018; Hunt et al. 2020). Some authors report a “break” or transition in the relation at a certain mass with a steeper relation above that characteristic mass (e.g., Hunt et al. 2020; Calette et al. 2018, although the latter note that both a single and a double power law fit the data well). We find no evidence for such a transition mass in our data, the binned medians agree well with a single-slope power law.

5. Role of the (group) environment

In order to investigate environmental effects on the cold gas content of galaxies, we use information from the SDSS galaxy group catalog by Tinker (2021). This allows us to split our sample into galaxies that reside in groups or clusters and those that are isolated. The non-isolated galaxies are further split into centrals and satellites as described in Sect. 2.6.

5.1. Molecular gas

The top panel in Fig. 9 shows the molecular gas fraction plotted against the stellar mass but with the sample split into the three environment classes described above. The Bayesian linear regression performed on each subsample yield best-fit slopes that are consistent with each other and with the complete star forming sample within their respective 1σ errors. Details on the best-fit parameters are listed in Table 2. There is no significant correlation between M_* and f_mol in any of the subsamples, nor is there a significant offset between the samples.

Fig. 9. Molecular gas scaling relations by environment. Top: molecular gas mass fraction vs. stellar mass by environment. Bottom: molecular gas depletion time vs. stellar mass by environment. The environment is classified as central (red), satellite (orange), and isolated (blue). The best fits from the Bayesian linear regression to the unbinned data are indicated by the blue-black-dashed (isolated), red-black-dotted (centrals), and orange-black-dash-dotted (satellites) lines respectively. The shaded areas indicate the 95% confidence bands. Open symbols with downward arrows denote upper limits. We note that for a few galaxies there was no environment information available from the group catalog.

The largest discrepancy is between isolated and satellite galaxies. For the isolated galaxies, the best-fit slope is steeper than for the full sample, and this is the only subsample for which it is not within 3σ of zero. The satellite subsample on the other hand has the shallowest best-fit slope. However, taking into account the uncertainties on the best fit parameters, both subsamples are consistent with the best fit-relation obtained for the full star-forming sample.

The bottom panel in Fig. 9 shows the molecular gas depletion time versus stellar mass plot with the sample divided into the same sub-categories of isolated, satellite, and central galaxies. Once again the best-fit parameters for the isolated, satellite, and central subsamples remain consistent with those obtained for the entire star-forming sample to within the estimated uncertainties, and no significant variations are observed as a function of environment.

5.2. Neutral atomic gas

Neutral atomic hydrogen resides at much lower densities and larger radii in galaxies than its molecular counterpart, and should thus be more readily influenced by the environment. The environmental impact on HI are well documented in the case of cluster galaxies, which have been found to be HI deficient (e.g., Haynes et al. 1984), but the picture for group galaxies is less clear (Appleton et al. 2013, 2023; de Blok et al. 2018).

In Fig. 10 we show the neutral atomic gas fraction f_HI plotted against stellar mass for our three environment based subsamples. Much like in the case of molecular gas, the three separate Bayesian linear regression fits remain consistent with the results for the entire star-forming sample.

Fig. 10. Neutral atomic gas fraction vs. stellar mass for central galaxies (red), satellites (orange), and isolated galaxies (blue).

6. Discussion and conclusions

In our reanalysis of the combined ALLSMOG and xCOLDGASS star-forming galaxies, we include CO non-detections as 3σ upper limits. We apply a metallicity and radiation field dependent α_CO conversion factor that accounts for the dominant effect of turbulence in starburst galaxies (Accurso et al. 2017). Our results indicate that the molecular gas fraction does not correlate with stellar mass for star-forming galaxies across the mass range 10^8.5 M_⊙ < M_* < 10^10.5 M_⊙. For M_* > 10^10.5 M_⊙, the f_mol vs. M_* relation becomes highly sensitive to sample selection criteria, as evidenced by the comparison of our results with previous studies (Saintonge et al. 2017; Calette et al. 2018) and our own tests (Appendix A). We adopt a criterion based on the offset from the star-forming main sequence (following the definition of Renzini & Peng 2015) in this work, and find the molecular gas fractions of the resulting sample are consistent with aforementioned flat trend for M_* > 10^10.5 M_⊙. Using a main sequence definition with a high-mass turnover instead (Saintonge & Catinella 2022) results in a sample with a sharp drop-off in the median molecular gas fraction above M_* ∼ 10^10.5 M_⊙. From this evidence we conclude that, in the local star-forming population in the mass range of 10^8.5 M_⊙ < M_* < 10^10.5 M_⊙, the abundance of cold molecular gas does not depend on stellar mass. The molecular gas depletion time τ_mol on the other hand shows a mild positive dependence on stellar mass, indicating that more massive galaxies are less efficient at forming stars than their low-mass counterparts. This is consistent with previous reports on gas consumption timescales (e.g., Boselli et al. 2014; Saintonge et al. 2017; Hunt et al. 2020) and has been linked to the presence of stellar bulges (Saintonge et al. 2011; Davis et al. 2014) which could stabilize the gas against collapse (Martig et al. 2009).

In contrast to our results, an anticorrelation of M_* and f_mol has been reported in previous work on low-mass galaxies (e.g., Hunt et al. 2020). This discrepancy between scaling relations in the literature is likely influenced by differences in methodology, particularly sample selection and assumptions on the α_CO. Results differ in particular between studies based only on CO detections (e.g., Hunt et al. 2020), and those taking non-detections into account as upper limits (e.g., Saintonge et al. 2017; Calette et al. 2018) as is the case in this work. Examples where authors report scaling relations derived for the same sample with both approaches show that this can lead to significant differences. Baker et al. (2022) report a ∼20% steeper best-fit slope for the M_mol vs. M_* relation in their sample when fitting only detections compared to when taking non-detections into account.

The inclusion of upper limits requires statistical methods suited to extract information from censored data. We employ three such techniques here, namely Bayesian linear regression (Kelly 2007), survival analysis in the form of the Kaplan-Meier estimator (Kaplan & Meier 1958), and an estimator for Kendall’s tau statistic (Oakes 2008). For our sample, all three methods yield results that are consistent with one another. We additionally test spectral stacking for galaxies with M_* < 10⁹ M_⊙ where our data are dominated by non-detections (see Appendix D), and find the results to be consistent with the other methods. A detailed comparison between commonly employed approaches can be found in Saintonge & Catinella (2022). Finally, we note that while the inclusion of non-detections is important to avoid biasing the sample and recovering the flat trend we observe, we need to increase the number of CO detections at M_* < 10⁹ M_⊙ significantly to place stronger constraints on the f_mol vs. M_* scaling relations in the low-mass regime. This would also allow to distinguish between a flat f_mol vs. M_* relation and a downturn in the low-mass regime, which cannot be done with upper limits from non-detections. Such constraints would require deeper observations of their molecular gas content in the same unbiased manner that the original ALLSMOG and xCOLDGASS surveys employed but with a much more sensitive telescope, such as the Atacama Large Aperture Submillimeter Telescope (AtLAST, see e.g. Mroczkowski et al. 2023, and in prep.; Ramasawmy et al. 2022; Klaassen et al. 2020).

Another important difference between studies is the choice of α_CO. The prescription we adopt for this work from Accurso et al. (2017) attempts to take into account variations of the α_CO on gas-phase metallicity, local radiation field and dynamical state of the molecular gas. Other studies account only for the dependence on metallicity but vary in functional form and effective slope of the dependence, especially at low metallicities (Hunt et al. 2020; Glover & Low 2011; Wolfire et al. 2010). Baker et al. (2022) forego the metallicity dependence altogether and apply a constant α_CO value based on data from the Milky Way. To assess the effect of the adopted conversion factor, we repeated our analysis using the prescription recommended by Hunt et al. (2020), which has a shallow slope compared to most other works, and a constant value corresponding to what is observed in the Milky Way (Bolatto et al. 2013). We find that a metallicity-dependent conversion factor is necessary to recover the flat f_mol scaling we report, but that the exact functional form can vary without fundamentally changing our results.

In low-mass, low-metallicity galaxies, where much of the molecular gas may be “CO-dark” (e.g., Madden et al. 2020), the CO(1−0) line may not be a good tracer of the overall molecular gas mass at all. Some models suggest that the connection between CO(1−0) emission and molecular gas mass becomes increasingly complicated in the low-metallicity regime, depending strongly on the “clumpiness” of the ISM in addition to metallicity and radiation field (Ramambason et al. 2024). Other tracers, such as [CII]158 μm and rotational lines of [CI], have been suggested as more effective in recovering molecular gas masses in low-metallicity galaxies (Ramambason et al. 2024; Papadopoulos et al. 2021; Bisbas et al. 2021), and could be used instead of relying on increasingly complex prescriptions for α_CO. However, these tracers are relatively more challenging to observe, in particular [CI], which can be observed from ground in local galaxies in both its transitions (at rest-frame frequencies of 492 GHz and 809 GHz), but these observations require a sensitive telescope placed on a high dry site with optimal submillimeter weather conditions, such as AtLAST.

While we find the molecular gas fraction to be independent of stellar mass in local star-forming galaxies, there is significant scatter of 0.31 dex around this relation. This scatter correlates with sSFR and is consistent with the scatter in SFR observed around the main sequence. We investigate the galaxy (group) environment as a possible factor in the origin of this scatter, by analyzing the differences between central, satellite, and isolated galaxies in our sample, classified using the SDSS group catalog by Tinker (2021). Our results indicate that the f_mol vs. M_*, τ_mol vs. M_*, and f_HI vs. M_* scaling relations do not differ significantly between the three subsamples. From this we conclude that the group environment is not likely to be the primary driver of the scatter around the constant molecular gas fraction we find in our sample. It is possible that a more nuanced separation of galaxies by environment, using additional parameters such as group size or separation between members, could reveal environmental effects that our classification is not sensitive to, but we do not have the sample size to perform such an analysis here.

Acknowledgments

We would like to thank the anonymous referee for their insightful comments. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 951815 (AtLAST). S. Shen acknowledges support from the European High Performance Computing Joint Undertaking (EuroHPC JU) and the Research Council of Norway through the funding of the SPACE Centre of Excellence (grant agreement No. 101093441). C. Cicone, P. Severgnini, and C. Vignali acknowledge a financial contribution from the Bando Ricerca Fondamentale INAF 2022 Large Grant, “Dual and binary supermassive black holes in the multi-messenger era: from galaxy mergers to gravitational waves”. This publication is based on data acquired with the Atacama Pathfinder Experiment (APEX) under programme ID 192.A-0359. APEX is a collaboration between the Max-Planck-Institut für Radioastronomie, the European Southern Observatory, and the Onsala Space Observatory. The calibrated and reduced APEX spectra of the ALLSMOG survey are publicly available through the ESO Phase 3 archive (ALLSMOG Phase 3 data release; http://archive.eso.org/wdb/wdb/adp/phase3_spectral/form?collection_name=ALLSMOG).

References

Appendix A: Impact of the star-forming sample selection

We explore the effect of the chosen criterion for classifying galaxies as star-forming on the resulting scaling relations by repeating the analysis with different criteria. First, we test how the inclusion of starbursts in our sample impacts the f_mol vs. M_* scaling relation. To this end, we change the criterion of log Δ(MS) >  − 0.4 dex to |log Δ(MS)| < 0.4 dex while maintaining the definition of the main sequence following Renzini & Peng (2015). This ensures that galaxies with specific star formation rates high above the main sequence are excluded from the sample, in addition to the passive quenched galaxies. The f_mol vs. M_* relation for this altered sample is shown in the top panel of Fig. A.1. Comparing this to the main sequence and starburst sample from Fig. 5 makes it clear that the inclusion of starbursts does not significantly change the overall scaling relation. This is unsurprising insofar as the different physics of starbursts are, to some degree, taken into account in the α_CO prescription we adopt in our analysis, as discussed in Sect. 2.5. We perform the same test for the curved main sequence definition from Saintonge & Catinella (2022) and find the same result (see bottom panel of Fig. A.1 compared to Fig. 6).

Fig. A.1. Molecular gas fraction fmol = Mmol/SFR vs. stellar mass for main sequence galaxies according to the definitions of Renzini & Peng (2015) (top) and Saintonge & Catinella (2022) (bottom), respectively. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. Open symbols with downward arrows denote upper limits. The associated error bars show the extent of the bin (in x) and the 25–75 percentile interval (in y)

Second, we test approaches to sample selection which do not rely on an analytic main sequence definition. One such approach is to select galaxies that have log(sSFR/[yr⁻¹]) > −11 regardless of stellar mass. The results, shown in Fig. A.2, are similar to those obtained when assuming the Renzini & Peng (2015) main sequence (Fig. 5). Another frequently used selection method (for example in Baker et al. 2022) relies on the BPT classification (see Sect. 2.3). Using the BPT diagram, we exclude AGN and use only galaxies classified as SF or composite as our sample. The results, shown in Fig. A.3, are again close to those of the linear main sequence and constant sSFR cutoff samples.

Fig. A.2. Molecular gas fraction fmol = Mmol/SFR vs. stellar mass for the sSFR selected sample. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25-75 percentile interval (in y). Open symbols with downward arrows denote upper limits.

Fig. A.3. Molecular gas fraction fmol = Mmol/SFR vs. stellar mass for the BPT-selected SF plus composite sample. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25-75 percentile interval (in y). The associated error bars show the extent of the bin (in x) and the 25-75 percentile interval (in y). Open symbols with downward arrows denote upper limits.

We note that, while the choice of how to demarcate the star-forming population clearly affects the scaling relation between M_* and f_mol at high masses (above M_* ∼ 10^10.5 M_⊙), there is no significant effect on the low-mass regime. Our predictions for the molecular gas scaling relations in the very low-mass regime M_* < 10⁹ M_⊙ should therefore be robust against differences in sample selection.

Appendix B: Choice of CO-to-H₂ conversion factor

To test the influence of the chosen prescription for the CO-to-H₂ conversion factor α_CO on our results, we re-derive the f_mol vs. M_* scaling relation using two additional prescriptions. The first is the simple assumption of a constant α_CO with a Milky Way derived value of 4.35 (K km s⁻¹ pc²)⁻¹ (Bolatto et al. 2013). The second is a prescription proposed by Hunt et al. (2020) based on a sample of local low-mass (M_* ∼ 10⁷ − 10¹¹ M_⊙) star-forming galaxies. It follows a simple power law scaling with gas-phase metallicity for low-metallicity galaxies while keeping α_CO constant in the high-metallicity regime:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{CO}}=\{\begin{array}{cc}{\alpha }_{\mathrm{CO},\odot }{(Z/Z_{\odot })}^{-1.55}\hfill & \mathrm{if}Z<Z_{\odot }\hfill \\ {\alpha }_{\mathrm{CO},\odot }\hfill & \mathrm{if}Z\ge Z_{\odot }.\hfill \end{array}\end{array}$$(B.1)

Here α_CO, ⊙ = 4.35 (K km s⁻¹ pc²)⁻¹ is the same Milky Way derived conversion factor as above and Z_⊙ = 8.69 is the solar metallicity in units of 12 + log(O/H). This prescription is functionally similar to previous suggestions for a metallicity dependent α_CO (Glover & Low 2011; Wolfire et al. 2010) but has a somewhat shallower slope than most (see also Bolatto et al. 2013 for a comparison between prescriptions).

Adopting the metallicity dependent prescription from Hunt et al. (2020) decreases the correlation between M_* and f_mol and boosts the correlation found between M_* and τ_mol somewhat, but the overall picture remains unchanged. The best-fit slopes for both relations are consistent with those obtained using the Accurso et al. (2017) factor, to within 3σ. Adopting a constant conversion factor instead results in a positive slope in the f_mol vs. M_* relation, as it systematically underestimates the molecular gas masses for low-mass, low-metallicity galaxies. For the same reason it results in a steeper slope for the τ_mol vs. M_* relation.

Appendix C: Metallicity weighted fit

At low masses (M_* < 10⁹ M_⊙), our sample is drawn entirely from the ALLSMOG survey. Since the original survey selected galaxies above a certain gas phase metallicity (see Sect. 2 for details), we need to ensure our sample is not biased towards metal-rich galaxies in this mass range. We do so by comparing the metallicity distribution of our sample against that of all star-forming galaxies from the full MPA-JHU catalog (selected using the same criterion of Δ(MS) > −0.4 dex we outline in Sect. 2). For this comparison we bin both samples in stellar mass (with 0.2 dex bin size) and metallicity (0.1 dex bin size). The results for the first three mass bins are shown in Fig. C.1. These examples illustrate that the distribution of gas-phase metallicity in our sample is indeed somewhat skewed towards high metallicities at low stellar masses. In the 8.7 < log M_*/M_⊙ < 8.9 bin we achieve the worst coverage, but our sample range still covers just over 90% of the MPA-JHU distribution.

Fig. C.1. Gas-phase metallicity distribution of the full MPA-JHU star-forming sample vs. our sample in the three lowest stellar mass bins. The histograms are normalized to the real distribution in our sample for each bin.

To alleviate this bias, we calculate weights for each mass-metallicity bin by dividing the normalized histogram of the complete MPA-JHU sample by that of our sample. During the fitting procedure, each galaxy in the mass range 10^8.5 M_⊙ < M_* < 10^10.5 M_⊙ where our sample shows signs of metallicity bias, is weighted by the appropriate factor for its bin.

Appendix D: Stacking of low-mass non-detections

In the 8.5 < log M_*/M_⊙ < 9.0 mass range, most of the data in our sample are upper limits. While these upper limits are crucial in recovering unbiased scaling relations in the low-mass regime, they do not provide as strong a constraint as detections would. Here, we test the possibility of improving the constraints derived from non-detections by spectral stacking.

Our sample contains 22 non-detections with stellar masses below 10⁹ M_⊙ that were observed in the CO(2-1) emission line with APEX as part of the ALLSMOG survey. The individual spectra of these non-detections are re-binned to a spectral resolution of 0.2 km/s, because the native spectral resolution is not consistent between different observations. Spectra are then co-added, weighted by their inverse variance (calculated over two spectral windows of [-900,-300] km/s and [300,900] km/s relative to the rest frequency of the CO(2-1) line inferred from the galaxies’ redshifts), and re-binned again on a 10 km/s grid. The stacked spectrum shows no clear detection of the CO(2-1) line (see Fig. D.1). Nonetheless, we can use the rms noise of the stacked spectrum to calculate an upper limit on the average CO(2-1) line flux of the galaxies in the bin following the procedure outlined in Cicone et al. (2017). The resulting upper limit is about a factor of five lower than those calculated for the individual galaxies, since the rms noise scales with the inverse square root of the number of spectra in the stack when co-adding white noise.

Fig. D.1. Stacked spectrum of 22 CO(2-1) non-detections with M* < 109 M⊙ from the ALLSMOG survey. The title lists the rms noise and the derived upper limit (UL) on the line flux.

From the upper limit on the CO(2-1) flux, we calculate an upper limit on the mean molecular gas mass of the galaxies in the stack, following the same procedure as for the individual galaxies, as outlined in Section 2.5. For the stellar mass, luminosity distance, and redshift needed in the calculations, we use the mean of all galaxies in the stack. The α_CO conversion factor varies between ∼3 (K km s⁻¹ pc²)⁻¹ and ∼20 (K km s⁻¹ pc²)⁻¹ within the bin (with one notable outlier at 30 (K km s⁻¹ pc²)⁻¹) and shows no trend with stellar mass. We adopt a α_CO = 20 (K km s⁻¹ pc²)⁻¹ in our computation, which corresponds to the higher end of the typical values in the bin, ensuring that the result is a meaningful upper limit. The resulting upper limit is above, and therefore consistent with, the KM median of the lowest mass bin (see Fig. 5).

All Tables

All Figures

Fig. 1. BPT classification of the sample. Left: BPT diagram of the full sample. The dashed line is the demarcation between pure SF galaxies and AGN or composite galaxies according to Kauffmann et al. (2003). The solid line is the theoretical limit for the ionization state produced by star formation according to Kewley et al. (2001), splitting the sample further into AGN and composite systems. Right: distribution of stellar masses in the SF, AGN, composite and low-S/N subsamples.

In the text

Fig. 2. SFR vs. stellar mass for the 625 ALLSMOG + xCOLDGASS galaxies that have a reliable SFR in the MPA-JHU catalog. The color scale corresponds to the gas phase metallicity, which was calculated from the O3N2 line index (Curti et al. 2017) for galaxies classified as SF or composite in the BPT diagram and according to the mass-metallicity relation (Curti et al. 2020) for those classified AGN and low-S/N. The blue dash-dotted line shows the empirical relation for the star forming main sequence found by Renzini & Peng (2015) including a 0.4 dex scatter (blue shaded). The green dashed line shows the star forming main sequence following Saintonge & Catinella (2022).

In the text

	Fig. 3. Distribution of our sample in gas phase metallicity (left), distance off the main sequence (center), and αCO (right). Shaded areas indicate the parameter space where Accurso et al. (2017) suggest adapting their prescription for αCO.
In the text

Fig. 4. Scaling relations of the molecular gas mass. Left: molecular gas mass vs. SFR for the star-forming sample. For comparison we show the integrated inverse Schmidt–Kennicutt relation (blue dash-dotted) and corresponding intrinsic scatter reported by Chown et al. (2022). Right: molecular gas mass vs. M*, with the best-fit MGMS from Chown et al. (2022) (blue dash-dotted) and Baker et al. (2022) with only detections (purple dotted) and detections and upper limits (green dashed). The color scale corresponds to the adopted CO-to-H2 conversion factor. In both panels the solid gray line shows the best fit from the Bayesian linear regression to the unbinned data with the shaded area indicating the 95% confidence band. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25−75 percentile interval (in y). Open symbols with downward arrows denote upper limits. The navy area of the color bar indicates the specific value of the conversion factor used for extreme starbursts as described in Sect. 2.5.

In the text

Fig. 5. Molecular gas mass fraction vs. stellar mass for the star-forming sample. The solid gray line shows the best fit from the Bayesian linear regression to the unbinned data with the shaded area indicating the 95% confidence band. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25−75 percentile interval (in y). Open symbols with downward arrows denote upper limits. The purple diamond shows the upper limit derived from spectral stacking of the 22 ALLSMOG CO(2−1) non-detections with masses below M* = 109 M⊙. The blue dash-dotted line shows the scaling relations reported by Calette et al. (2018) based on their sample of late-type galaxies.

In the text

	Fig. 6. Same scaling relation as in Fig. 5 but for the star-forming plus starburst sample defined using the main sequence definition from Saintonge & Catinella (2022).
In the text

Fig. 7. Molecular gas depletion time τmol = Mmol/SFR vs. stellar mass. The solid gray line shows the best fit from the Bayesian linear regression to the unbinned data with the shaded area indicating the 95% confidence band. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25−75 percentile interval (in y). Open symbols with downward arrows denote upper limits.

In the text

Fig. 8. Neutral atomic gas mass fraction vs. stellar mass for the star-forming sample. The solid gray line shows the best fit from the Bayesian linear regression to the unbinned data with the shaded area indicating the 95% confidence band. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25−75 percentile interval (in y). Open symbols with downward arrows denote upper limits.

In the text

Fig. 9. Molecular gas scaling relations by environment. Top: molecular gas mass fraction vs. stellar mass by environment. Bottom: molecular gas depletion time vs. stellar mass by environment. The environment is classified as central (red), satellite (orange), and isolated (blue). The best fits from the Bayesian linear regression to the unbinned data are indicated by the blue-black-dashed (isolated), red-black-dotted (centrals), and orange-black-dash-dotted (satellites) lines respectively. The shaded areas indicate the 95% confidence bands. Open symbols with downward arrows denote upper limits. We note that for a few galaxies there was no environment information available from the group catalog.

In the text

	Fig. 10. Neutral atomic gas fraction vs. stellar mass for central galaxies (red), satellites (orange), and isolated galaxies (blue).
In the text

Fig. A.1. Molecular gas fraction fmol = Mmol/SFR vs. stellar mass for main sequence galaxies according to the definitions of Renzini & Peng (2015) (top) and Saintonge & Catinella (2022) (bottom), respectively. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. Open symbols with downward arrows denote upper limits. The associated error bars show the extent of the bin (in x) and the 25–75 percentile interval (in y)

In the text

	Fig. A.2. Molecular gas fraction fmol = Mmol/SFR vs. stellar mass for the sSFR selected sample. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25-75 percentile interval (in y). Open symbols with downward arrows denote upper limits.
In the text

Fig. A.3. Molecular gas fraction fmol = Mmol/SFR vs. stellar mass for the BPT-selected SF plus composite sample. The black diamonds correspond to the binned medians computed using the Kaplan-Meier estimator. The associated error bars show the extent of the bin (in x) and the 25-75 percentile interval (in y). The associated error bars show the extent of the bin (in x) and the 25-75 percentile interval (in y). Open symbols with downward arrows denote upper limits.

In the text

	Fig. C.1. Gas-phase metallicity distribution of the full MPA-JHU star-forming sample vs. our sample in the three lowest stellar mass bins. The histograms are normalized to the real distribution in our sample for each bin.
In the text

	Fig. D.1. Stacked spectrum of 22 CO(2-1) non-detections with M* < 109 M⊙ from the ALLSMOG survey. The title lists the rms noise and the derived upper limit (UL) on the line flux.
In the text