Scaling relations and baryonic cycling in local star-forming galaxies: I. The sample

TL;DR

This study compiles a comprehensive local galaxy sample to analyze the relationships between metallicity, gas content, stellar mass, and star formation rate, revealing that these parameters are primarily governed by two underlying dimensions.

Contribution

The paper introduces the MAGMA sample, covering extensive parameter space, and demonstrates that galaxy properties can be effectively described by a two-dimensional hypersurface using PCA.

Findings

Galaxy relations are primarily two-dimensional.

A piecewise 3D PCA predicts metallicity with ~0.1 dex accuracy.

MAGMA serves as a benchmark for galaxy evolution models.

Abstract

Metallicity and gas content are intimately related in the baryonic exchange cycle of galaxies, and galaxy evolution scenarios can be constrained by quantifying this relation. To this end, we have compiled a sample of ~400 galaxies in the Local Universe, dubbed "MAGMA" (Metallicity And Gas for Mass Assembly), which covers an unprecedented range in parameter space, spanning more than 5 orders of magnitude in stellar mass (Mstar), star-formation rate (SFR), and gas mass (Mgas), and a factor of ~60 in metallicity [Z, 12+log(O/H)]. Stellar masses and SFRs have been recalculated for all the galaxies using IRAC, WISE and GALEX photometry, and 12+log(O/H) has been transformed, where necessary, to a common metallicity calibration. To assess the true dimensionality of the data, we have applied multi-dimensional principal component analyses (PCAs) to our sample. In confirmation of previous work, we find that even with the vast parameter space covered by MAGMA, the relations between Mstar, SFR, Z and Mgas (MHI+MH2) require only two dimensions to describe the hypersurface. To accommodate the curvature in the Mstar-Z relation, we have applied a piecewise 3D PCA that successfully predicts observed 12+log(O/H) to an accuracy of ~0.1dex. MAGMA is a representative sample of isolated star-forming galaxies in the Local Universe, and can be used as a benchmark for cosmological simulations and to calibrate evolutionary trends with redshift.

Tables

Table 1. Table 1: Logarithmic differences of M star subscript 𝑀 star M_{\rm star}  and SFR between the original description and adopted values a

Sample	log($M_{\mathrm{star}}$/M⊙)	log(SFR/M⊙ yr-1)	Number
	[old $-$ new]	[old $-$ new]
(1)	(2)	(3)	(4)
xGASS-CO	$-0.018\pm \mathrm{\hspace{0.17em}0.02}$	$-0.045\pm \mathrm{\hspace{0.17em}0.19}$	181
ALLSMOG	$-0.278\pm \mathrm{\hspace{0.17em}0.32}$	$-0.083\pm \mathrm{\hspace{0.17em}0.46}$	38
HRS	$-0.170\pm \mathrm{\hspace{0.17em}0.17}$	$-0.024\pm \mathrm{\hspace{0.17em}0.20}$	68
KINGFISH	$0.122\pm \mathrm{\hspace{0.17em}0.06}$	$0.069\pm \mathrm{\hspace{0.17em}0.11}$	35
NFGS	$0.005\pm \mathrm{\hspace{0.17em}0.16}$	$-0.214\pm \mathrm{\hspace{0.17em}0.57}$	26
BCD	$-0.048\pm \mathrm{\hspace{0.17em}0.21}$	$0.211\pm \mathrm{\hspace{0.17em}0.29}$	17
DGS	$0.062\pm \mathrm{\hspace{0.17em}0.05}$	$0.115\pm \mathrm{\hspace{0.17em}0.30}$	13
Virgo SFDs	$0.103\pm \mathrm{\hspace{0.17em}0.09}$	$0.078\pm \mathrm{\hspace{0.17em}0.62}$	7
Extra	$0.069\pm \mathrm{\hspace{0.17em}0.05}$	$0.107\pm \mathrm{\hspace{0.17em}0.12}$	7

Table 2. Table 2: 4D PCA results for MAGMA a

Method	PC4(1)	PC4(2)	PC4(3)	PC4(4)	PC4	PC4	PC3	PC1$+$PC2
	12$+$log(O/H)	log	log	log($x$)	std. dev.	proportion of variance
	(PP04N2)	($M_{\mathrm{star}}$/M⊙)	(SFR/M⊙ yr-1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
				$x=$$M_{\mathrm{HI}}$/M⊙
PCA	0.920	$-0.355$	0.164	$0.019$	0.127	0.010	0.051	0.94
PPCA	$0.95\pm 0.01$	$-0.29\pm 0.03$	$0.11\pm 0.03$	$-0.01\pm 0.01$	0.14	0.01
BSPCA	$0.92\pm 0.01$	$-0.36\pm 0.03$	$0.16\pm 0.03$	$0.02\pm 0.02$	0.13	0.01
				$x=$$M_{\mathrm{H2}}$/M⊙
PCA	0.913	$-0.378$	0.133	$0.079$	0.124	0.010	0.051	0.94
PPCA	$0.94\pm 0.01$	$-0.31\pm 0.03$	$0.09\pm 0.03$	$0.05\pm 0.02$	0.14	0.01
BSPCA	$0.91\pm 0.02$	$-0.38\pm 0.03$	$0.13\pm 0.03$	$0.08\pm 0.04$	0.12	0.01
				$x=$$M_{\mathrm{gas}}$/M⊙
PCA	0.917	$-0.366$	0.153	$0.048$	0.126	0.011	0.054	0.94
PPCA	$0.95\pm 0.01$	$-0.24\pm 0.02$	$0.13\pm 0.02$	$-0.09\pm 0.03$	0.14	0.02
BSPCA	$0.91\pm 0.02$	$-0.37\pm 0.03$	$0.16\pm 0.03$	$0.05\pm 0.03$	0.13	0.01
				$x=$$L_{\mathrm{CO}}^{\prime }$/$K\mathrm{km}{s}^{-1}{\mathrm{pc}}^2$
PCA	0.955	$-0.185$	0.181	$-0.147$	0.117	0.007	0.027	0.97
PPCA	$0.97\pm 0.01$	$-0.16\pm 0.04$	$0.13\pm 0.02$	$-0.13\pm 0.03$	0.13	0.01
BSPCA	$0.95\pm 0.01$	$-0.18\pm 0.04$	$0.18\pm 0.03$	$-0.15\pm 0.02$	0.12	0.01
Method	PC4(1)	PC4(2)	PC4(3)	PC4(4)	PC4	PC4	PC3	PC1$+$PC2
	12$+$log(O/H)	log	log	log	std. dev.	variance	variance	variance
	(KD02)	($M_{\mathrm{star}}$/M⊙)	(SFR/M⊙ yr-1)	($M_{\mathrm{H2}}$/M⊙)b
PCA	$0.86$	$-0.45$	$0.22$	$0.03$	0.153	0.014	0.048	0.95
PPCA	$0.92\pm 0.02$	$-0.37\pm 0.04$	$0.15\pm 0.04$	$0.02\pm 0.02$	0.17	0.02
BSPCA	$0.86\pm 0.03$	$-0.45\pm 0.04$	$0.22\pm 0.04$	$0.03\pm 0.03$	0.15	0.01

Table 3. Table 3: MAGMA sample means

Quantity	Meana	Std. dev.	Mean	Mean
			($M_{\mathrm{star}}$$\le$	($M_{\mathrm{star}}$$>$
			$M_{\mathrm{break}}$b)	$M_{\mathrm{break}}$)
(1)	(2)	(3)	(4)	(5)
log($M_{\mathrm{star}}$/M⊙)	$9.684$	$0.75$	$9.482$	$10.565$
log(SFR/M⊙ yr-1)	$-0.229$	$0.66$	$-0.375$	$0.409$
12$+$log(O/H)	$8.580$	$0.22$	$8.546$	$8.728$
log($M_{\mathrm{HI}}$/M⊙)	$9.279$	$0.70$	$9.166$	$9.774$
log($M_{\mathrm{H2}}$/M⊙)	$8.720$	$0.69$	$8.550$	$9.463$
log($M_{\mathrm{gas}}$/M⊙)	$9.438$	$0.64$	$9.313$	$9.985$
log($L_{\mathrm{CO}}^{\prime }$/	$7.950$	$0.96$	$7.728$	$8.918$
$K\mathrm{km}{s}^{-1}{\mathrm{pc}}^2$)

Table 4. Table 4: 3D PCA results for MAGMA a

Method	PC3(1)	PC3(2)	PC3(3)	PC3	PC3
	log($M_{\mathrm{star}}$/M⊙)	log(SFR/M⊙ yr-1)		std. dev.	variance
(1)	(2)	(3)	(4)	(5)	(6)
			12$+$log(O/H)
PCA	$0.346$	$-0.169$	$-0.923$	0.127	0.015
PPCA	$0.29\pm 0.02$	$-0.11\pm 0.03$	$-0.95\pm 0.01$	0.14	0.02
BSPCA	$0.34\pm 0.02$	$-0.17\pm 0.03$	$-0.92\pm 0.01$	0.13	0.02
			log($M_{\mathrm{HI}}$/M⊙)
PCA	0.700	$-0.704$	$-0.122$	0.264	0.047
PPCA	$0.69\pm 0.03$	$-0.71\pm 0.04$	$-0.12\pm 0.07$	0.29	0.06
BSPCA	$0.70\pm 0.03$	$-0.70\pm 0.05$	$-0.13\pm 0.08$	0.26	0.05
			log($M_{\mathrm{H2}}$/M⊙)
PCA	0.680	$-0.730$	$-0.061$	0.267	0.048
PPCA	$0.55\pm 0.20$	$-0.68\pm 0.19$	$0.31\pm 0.25$	0.29	0.06
BSPCA	$0.61\pm 0.18$	$-0.65\pm 0.22$	$0.33\pm 0.23$	0.26	0.05
			log($M_{\mathrm{gas}}$/M⊙)
PCA	$0.742$	$-0.613$	$-0.272$	0.263	0.049
PPCA	$0.65\pm 0.18$	$-0.59\pm 0.18$	$-0.35\pm 0.22$	0.29	0.06
BSPCA	$0.70\pm 0.09$	$-0.60\pm 0.16$	$-0.30\pm 0.18$	0.26	0.05
			log($L_{\mathrm{CO}}^{\prime }$/$K\mathrm{km}{s}^{-1}{\mathrm{pc}}^2$)
PCA	0.835	$-0.356$	$-0.420$	0.229	0.027
PPCA	$0.82\pm 0.04$	$-0.42\pm 0.11$	$-0.36\pm 0.08$	0.27	0.04
BSPCA	$0.83\pm 0.01$	$-0.36\pm 0.08$	$-0.42\pm 0.06$	0.23	0.03

Table 5. Table 5: 4D PCA results for the Bothwell et al. ( 2016b ) sample a

Method	PC4(1)	PC4(2)	PC4(3)	PC4(4)	PC4	PC4	PC3	PC1$+$PC2
	12$+$log(O/H)	log	log	log	std. dev.	proportion of variance
	(KD02)	($M_{\mathrm{star}}$/M⊙)	(SFR/M⊙ yr-1)	($M_{\mathrm{H2}}$/M⊙)b
Including the 17 high-$z$ galaxies and 8 BCDs (166)
PCA	$0.872$	$-0.408$	$-0.125$	$0.242$	0.169	0.014	0.042	0.943
PPCA	$0.905\pm 0.016$	$-0.381\pm 0.018$	$-0.069\pm 0.038$	$0.165\pm 0.044$	0.182	0.017
PPCA	$0.903\pm 0.017$	$-0.383\pm 0.017$	$-0.071\pm 0.042$	$0.167\pm 0.047$	0.182	0.017
BSPCA	$0.871\pm 0.028$	$-0.404\pm 0.027$	$-0.125\pm 0.042$	$0.239\pm 0.054$	0.167	0.014
BSPCA	$0.865\pm 0.033$	$-0.405\pm 0.025$	$-0.131\pm 0.048$	$0.249\pm 0.062$	0.168	0.014
With 8 BCDs but without the 17 high-$z$ galaxies (149)
PCA	$0.891$	$-0.378$	$-0.173$	$0.186$	0.163	0.023	0.074	0.90
PPCA	$0.908\pm 0.015$	$-0.355\pm 0.020$	$-0.155\pm 0.038$	$0.145\pm 0.038$	0.173	0.027
PPCA	$0.909\pm 0.015$	$-0.357\pm 0.019$	$-0.150\pm 0.039$	$0.142\pm 0.042$	0.173	0.027
BSPCA	$0.888\pm 0.033$	$-0.374\pm 0.039$	$-0.171\pm 0.046$	$0.184\pm 0.069$	0.159	0.023
BSPCA	$0.886\pm 0.032$	$-0.378\pm 0.039$	$-0.169\pm 0.042$	$0.188\pm 0.066$	0.160	0.023

Equations

$$\log (M_{\mathrm{star}})=1.050\log (\nu L_{\nu })+\mathrm{0.387.}$$

$$\begin{array}{c}5.2\lesssim \log ({\mathrm{M}}_{\mathrm{star}}/M_{\odot })\lesssim 11.2\\ 5.4<\log ({\mathrm{M}}_{\mathrm{gas}}/M_{\odot })<10.8\\ -4.6\lesssim \log (\mathrm{SFR}/{\mathrm{M}}_{\odot }\text{leavevmode}\mathrm{y}{\mathrm{r}}^{-1})\lesssim 1.6\\ 7.5\lesssim 12+\log (O/H)\lesssim 9.3\end{array}$$

$$o$$

$$o=\{\begin{array}{ll}(0.43\pm 0.03)m-(0.21\pm 0.04)s& \text{if log(Mstar/M⊙) ≤ 10.3}\\ (0.25\pm 0.11)m-(0.11\pm 0.05)s& \text{if log(Mstar/M⊙) > 10.3}\end{array}$$

Full text

1 Introduction

As long as star formation occurs in their gas reservoirs, galaxies evolve increasing their stellar mass ($M_{\rm star}$) and their metal content, depending on the relative efficiency of inflows/outflows, dynamical interactions, and environmental processes. In other words, at any time, $M_{\rm star}$ and metallicity ($Z$) reflect the combined effect of both the integrated history of star formation and the degree of interaction with the surrounding environment. Not surprisingly, the causal links between gas mass ($M_{\rm gas}$), star formation rate (SFR), $M_{\rm star}$, and $Z$, manifest in a number of observed correlations between these quantities, often referred to as scaling relations. Some among the most notable examples are: (i) the correlation between $M_{\rm star}$ and SFR (dubbed the “Main Sequence”, MS: e.g., Brinchmann et al. 2004; Noeske et al. 2007; Daddi et al. 2010; Elbaz et al. 2011; Renzini & Peng 2015); (ii) the correlation between $M_{\rm gas}$ and SFR (the “Schmidt-Kennicutt”, SK, relation; e.g., Schmidt 1959; Kennicutt 1998; Bigiel et al. 2008; Leroy et al. 2009); and (iii) the “mass-metallicity relation”, MZR, between $M_{\rm star}$ and $Z$ (e.g., Lequeux et al. 1979; Tremonti et al. 2004; Maiolino et al. 2008). In star-forming galaxies, $Z$ is typically measured by the abundance of oxygen, O/H, in the ionized gas, as it is the most abundant heavy element produced by massive stars.

These scaling relations among fundamental properties of galaxies are potentially insightful tools to explore demographics of galaxies and their evolution. In particular, the mutual correlations among physical properties in galaxies imply that the observed residuals from the main relations (in other words, their intrinsic scatters) could be correlated with other variables. Many studies have investigated such a notion, and this type of analysis has proved to be a powerful diagnostic, providing simple quantitative tests for analytical models and numerical simulations.

Only recently has it been possible to incorporate gas properties in studies of baryonic cycling, thanks to the growing number of available gas measurements (atomic and molecular), including: the Arecibo Legacy Fast ALFA Survey (ALFALFA, Haynes et al. 2011, 2018) the Galaxy Evolution Explorer (GALEX) Arecibo SDSS Survey (GASS, Catinella et al. 2010, 2018); the COLD-GASS survey (Saintonge et al. 2011a, 2017); the Nearby Field Galaxy Survey (NFGS, Jansen & Kannappan 2001; Wei et al. 2010; Stark et al. 2013); the Herschel Reference Survey (HRS, Boselli et al. 2010; Cortese et al. 2011; Boselli et al. 2014a); and the APEX Low-redshift Legacy Survey for MOlecular Gas (ALLSMOG, Bothwell et al. 2014; Cicone et al. 2017). These surveys have provided important new observations of Hi and CO, in order to derive H2 and compare gas content with other galaxy properties. Results suggest that the relation of atomic gas to $M_{\rm star}$ and SFR drives a galaxy’s position relative to the MS (e.g., Huang et al. 2012; Gavazzi et al. 2013; Saintonge et al. 2016), and that Hi gas fractions increase with decreasing $M_{\rm star}$ and stellar mass surface density, at least down to log($M_{\rm star}$/M⊙) = 9 (e.g., Cortese et al. 2011; Gavazzi et al. 2013; Catinella et al. 2018). Incorporating molecular gas H2 in the analysis suggests that the strongest correlations are between H2 content and SFR; in particular, molecular depletion time depends strongly on specific SFR (sSFR $\equiv$ SFR/$M_{\rm star}$) (e.g., Saintonge et al. 2011a, b; Boselli et al. 2014b; Hunt et al. 2015; Saintonge et al. 2017).

Important clues to baryonic cycling also come from systematic studies of the intrinsic scatter of the MZR, finding that a fundamental metallicity relation (FMR) exists between $M_{\rm star}$, $Z$ and SFR, that minimizes the scatter in the MZR (see e.g., Ellison et al. 2008; Mannucci et al. 2010). According to the FMR, galaxies lie on a tight, redshift-independent two-dimensional (2D) surface in 3D space defined by $M_{\rm star}$, $Z$ and SFR, where at a given $M_{\rm star}$, galaxies with higher SFR have systematically lower gas-phase $Z$ (see e.g., Hunt et al. 2012; Lara-López et al. 2013; Hunt et al. 2016a; Hashimoto et al. 2018; Cresci et al. 2018). Many theoretical models have investigated this finding, explaining it in terms of an equilibrium between metal-poor inflows and metal-enriched outflows (e.g., Davé et al. 2012; Dayal et al. 2013; Lilly et al. 2013; Graziani et al. 2017). Observational results suggest that the FMR may be more strongly expressed via the gas mass rather than via the SFR (see e.g., Bothwell et al. 2013; Brown et al. 2018). In this light, the FMR might be interpreted as a by-product of an underlying relationship between the scatter of the MZR and the gas content (e.g., Zahid et al. 2014). In particular, Bothwell et al. (2016b), with an analysis that included $M_{\rm star}$, SFR, O/H, and molecular gas mass, $M_{\rm H2}$, suggest that the true FMR exists between $M_{\rm star}$, O/H and $M_{\rm H2}$, which is linked to SFR via the SK star-formation law.

Virtually all previous studies of gas scaling relations in galaxies have focused on galaxies more massive than $10^{9}$ M⊙. In this paper, we extend previous studies to lower stellar masses, reporting the analysis of the mutual dependencies of physical properties in a sample of $\sim 400$ local galaxies, with simultaneous availability of $M_{\rm star}$, SFR, $M_{\rm HI}$, $M_{\rm H2}$ (thus also total gas, $M_{\rm gas}$), and O/H, spanning an unprecedented range in $M_{\rm star}$, from $\sim 10^{5}$ M⊙ to $3\times 10^{11}$ M⊙. In Sect. 2, we first describe the individual sub-samples, and then homogenize the stellar mass and SFR estimates by incorporating mid-IR (MIR) fluxes from the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) and photometry from GALEX (Morrissey et al. 2007). With updated principle component analysis (PCA) techniques, Sect. 3 explores the correlations in the four- and three-dimensional (4D, 3D) parameter spaces defined by $M_{\rm star}$, SFR, O/H, and $M_{\rm gas}$, together with the two separate gas components $M_{\rm HI}$ and $M_{\rm H2}$. There is a particular focus in Sect. 4 on the MZR scatter and the ramifications of including a significant population of low-mass galaxies in the sample.

2 Combined sample: MAGMA

We have compiled a sample of 392 local galaxies, with simultaneous availability of $M_{\rm star}$, SFR, gas masses (both atomic, $M_{\rm HI}$, and molecular, $M_{\rm H2}$, the latter obtained by measurements of CO luminosity, $L_{\rm CO}^{\prime}$) and metallicities [12$+$log(O/H)]. We assembled our sample by combining a variety of previous surveys at $z\sim 0$ with new observations of CO in low-mass galaxies. The details of the parent surveys such as metallicity calibration, stellar-mass, and SFR determinations are provided below. The following four selection criteria are adopted:

1. (1)

only galaxies with robust ($\gtrsim 3\sigma$) detections of $M_{\rm star}$, SFR, 12$+$log(O/H), $M_{\rm HI}$, and $L_{\rm CO}^{\prime}$ are considered; 2. (2)

galaxies were eliminated if they were thought to host active galactic nuclei (AGN) based on the BPT classifications111The Baldwin-Philips-Terlevich (BPT) diagram classification (Baldwin et al. 1981) relies on the emission-line properties of galaxies, based on the [Sii]/H$\alpha$ versus [Oiii]/H$\beta$ ratios. provided by the original surveys; 3. (3)

when Hi-deficiency measurements Hi–def were available (e.g., Boselli et al. 2009, 2014a), following Boselli et al. (2014b), only galaxies with Hi–def $\leq$0.4 were retained; 4. (4)

properties of galaxies in common among two or more parent surveys have been taken from the sample that provided more ancillary information (e.g., high quality spectra, resolved maps, uniform derivation of parameters, etc.).

Hi-deficiency is defined as the logarithm of the ratio of the observed Hi mass of a galaxy and the mean Hi mass expected for an isolated galaxy with the same optical size and morphological type (e.g., Haynes & Giovanelli 1984). The Hi-deficiency requirement is included to ensure that our sample is representative of isolated, field galaxies, not having undergone potential stripping effects from residence in a cluster. Because we require metallicity and gas measurements, we have dubbed our compiled sample MAGMA (Metallicity And Gas in Mass Assembly). The final MAGMA sample has been drawn from the following nine parent surveys/papers:

• –

**xGASS-CO: ** xGASS-CO is the overlap between the extended GALEX Arecibo SDSS Survey (xGASS: Catinella et al. 2018) and the extended CO Legacy Database for GASS (xCOLD GASS: Saintonge et al. 2017). xGASS222The full xGASS representative sample is available on the xGASS website, http://xgass.icrar.org in digital format. is a gas fraction-limited census of the Hi gas content of $\sim 1200$ local galaxies, spanning over 2 decades in stellar mass (M_{\rm star}$$\,=\,10^{9}-10^{11.5} M⊙). The xCOLD GASS survey333The full xCOLD GASS survey data products are available on the xCOLD GASS website http://www.star.ucl.ac.uk/xCOLDGASS/. contains IRAM-30m CO(1-0) measurements for 532 galaxies also spanning the entire SFR-$M_{\rm star}$ plane at M_{\rm star}$$>10^{9} M⊙. Stellar masses are from the MPA-JHU444http://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/ catalogue, where $M_{\rm star}$ is computed from a fit to the spectral energy distribution (SED) obtained using SDSS broad-band photometry (Brinchmann et al. 2004; Salim et al. 2007). SFRs are computed as described by Janowiecki et al. (2017) by combining NUV with mid-IR (MIR) fluxes from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010). When these are not available (the case for $\sim$70% of the xGASS sample), SFRs are determined using a “ladder” technique (Janowiecki et al. 2017; Saintonge et al. 2017). Data Release 7 (SDSS DR7, Abazajian et al. 2009), calibrated by Saintonge et al. (2017) to the [Nii]-based strong-line calibration by Pettini & Pagel (2004, PP04N2). In addition to omitting AGN and Seyferts (see above), we have also excluded galaxies in Saintonge et al. (2017) with an “undetermined” or “composite” classification, based on the BPT diagram; metallicities from PP04N2 for such galaxies tend to be highly uncertain. xGASS-CO, the overlap between xGASS and xCOLDGASS, includes 477 galaxies, with 221 non-AGN galaxies with robust CO detections. The subset of xGASS-CO that respects our selection criteria (i.e., with Hi and CO detections and not excluded for potentially uncertain O/H calibration) consists of 181 galaxies.

• –

**HRS: ** The Herschel Reference Sample (Boselli et al. 2010) is a $K$-band selected, volume-limited sample comprising 323 galaxies. HRS555A full description of the survey and the ancillary data can be found at https://hedam.lam.fr/HRS/. is a fairly complete description of the Local Universe galaxy population although underrepresented in low-mass galaxies (see Boselli et al. 2010). Stellar masses $M_{\rm star}$ and SFRs were obtained from Boselli et al. (2015), where $M_{\rm star}$ values were derived according to the precepts of Zibetti et al. (2009) using $i$-band luminosities and $g-i$ colors, and SFRs are the mean of the four methods investigated by Boselli et al. (2015). These include radio continuum at 20 cm, H\alpha$$+24 $\mu$m luminosities, FUV$+$24 $\mu$m luminosities, and H$\alpha$ luminosities corrected for extinction using the Balmer decrement666Boselli et al. (2015) do not publish the individual estimates, so we were unable to select the hybrid method based on 24 $\mu$m luminosities that would be more consistent with other samples discussed here.. 12$+$log(O/H) was taken from Hughes et al. (2013) based on the PP04N2 calibration, and gas quantities, $M_{\rm HI}$ and $L_{\rm CO}^{\prime}$ were taken from Boselli et al. (2014a). 86 HRS galaxies have Hi and CO detections, but 18 of these have Hi–def$>$0.4 (as given by Boselli et al. 2014a), so we are left with 68 HRS galaxies that satisfy our selection criteria.

• –

**ALLSMOG: ** The APEX Low-redshift Legacy Survey of MOlecular Gas (Bothwell et al. 2014; Cicone et al. 2017) comprises 88 nearby, star-forming galaxies with stellar masses in the range $10^{8.5}<M_{\rm star}/M_{\odot}<10^{10}$, and gas-phase metallicities 12$+$log(O/H) $>$ 8.4. ALLSMOG777The full ALLSMOG survey data products are available on the web page http://www.mrao.cam.ac.uk/ALLSMOG/. is entirely drawn from the MPA-JHU catalogue of spectral measurements and galactic parameters of SDSS DR7. Stellar mass and SFR values of ALLSMOG galaxies are taken from the MPA-JHU catalogue, and the SFR is based on the (aperture- and extinction-corrected) H$\alpha$ intrinsic line luminosity. We have used the PP04N2 O/H calibration given by Cicone et al. (2017). To convert the ALLSMOG CO(2–1) values from Cicone et al. (2017) to the lower-J CO(1–0) available for the remaining samples, we assume $R_{21}\,=\,0.8$ as they advocate. The subset of ALLSMOG that respects our selection criteria consists of 38 galaxies.

• –

**KINGFISH: ** The Key Insights on Nearby Galaxies: a Far-Infrared Survey with Herschel, KINGFISH888An overview of the scientific strategy for KINGFISH and the properties of the galaxy sample can be found on the web page https://www.ast.cam.ac.uk/research/kingfish. (Kennicutt et al. 2011), contains 61 galaxies with metallicity in the range $7.54\,\leq$ 12$+$log(O/H) $\leq\,8.77$ and stellar masses in the range $[2\times 10^{7}-1.4\times 10^{11}]\,$M⊙. Stellar masses and SFRs are taken from Hunt et al. (2019). The $M_{\rm star}$ values were computed from the SFR-corrected IRAC 3.6 $\mu$m luminosities according to the luminosity-dependent mass-to-light (M/L) ratio given by Wen et al. (2013), and are within $\sim$0.1 dex of those derived by comprehensive SED fitting (see Hunt et al. 2019). SFRs are inferred from the far-ultraviolet (FUV) luminosity combined with total-infrared (TIR) luminosity following Murphy et al. (2011). Atomic gas masses $M_{\rm HI}$ and CO measurements for $M_{\rm H2}$ are taken from Kennicutt et al. (2011), with refinements from Sandstrom et al. (2013) and Aniano et al. (2020). “Representative” metallicities evaluated at 0.4 times the optical radius Ropt from Moustakas et al. (2010) were converted from the Kobulnicky & Kewley (2004, KK04) to the PP04N2 calibration according to the transformations given by Kewley & Ellison (2008); more details are given in Hunt et al. (2016a) and Aniano et al. (2020). After omitting NGC 2841 and NGC 5055, because their metallicities exceeded the valid regime for the Kewley & Ellison (2008) transformations, the required data are available for 38 KINGFISH galaxies. Three of these have Hi–def$>$0.4 (given by Boselli et al. 2014a) so we ultimately select 35 galaxies from KINGFISH.

• –

**NFGS: ** The Nearby Field Galaxy Survey (Jansen et al. 2000; Kewley et al. 2005; Kannappan et al. 2009) consists of 196 galaxies spanning the entire Hubble sequence in morphological types, and a range in luminosities from low-mass dwarf galaxies to luminous massive systems. Stellar masses, $M_{\rm star}$, are given by Kannappan et al. (2013) and are based on NUV$+$*ugrizJHK$+$*IRAC 3.6 $\mu$m SEDs. We have taken (spatially) integrated SFR (based on H$\alpha$) and O/H values from Kewley et al. (2005), and transformed 12$+$log(O/H) from their Kewley & Dopita (2002, KD02) calibration to PP04N2 according to the formulations by Kewley & Ellison (2008). Stark et al. (2013) provide CO and $M_{\rm H2}$ measurements, and $M_{\rm HI}$ is tabulated by Wei et al. (2010) and Kannappan et al. (2013). After removing NGC 7077, that appears in the following dwarf sample, there are 26 galaxies that meet our selection criteria.

• –

**BCDs: ** The Blue Compact Dwarf galaxies (BCDs) have been observed and detected in 12CO($1-0$) with the IRAM 30m single dish (Hunt et al. 2015, 2017). They were selected primarily from the primordial helium sample of Izotov et al. (2007), known to have reliable metallicities 12$+$log(O/H) measured through the direct electron-temperature (Te) method. An additional, similar, set of BCDs has been detected in 12CO($1-0$) (Hunt et al. 2020, in prep.) with analogous selection criteria. For both sets of BCDs, $M_{\rm star}$ is derived as for KINGFISH galaxies, namely from IRAC 3.6 $\mu$m or WISE 3.4 $\mu$m luminosities, after correcting for free-free, line emission based on SFR, and dust continuum when possible (see also Hunt et al. 2015). This method has been shown to be consistent with full-SED derived $M_{\rm star}$ values to within $\la$0.1 dex (Hunt et al. 2019). For the galaxy in common with the NFGS, NGC 7077, the two $M_{\rm star}$ estimates are the same to within 0.07 dex. SFRs are based on the Calzetti et al. (2010) combination of H$\alpha$ and 24 $\mu$m  luminosities. Hi masses are given by Hunt et al. (2015) and Hunt et al. (2020, in prep.). As mentioned above, 12$+$log(O/H) is obtained from the direct Te method (for details see Hunt et al. 2016a). The subset of BCDs that respect our selection criteria with Hi and CO [12CO($1-0$)] detections comprises 17 galaxies with metallicities 12$+$log(O/H) ranging from 7.7 to 8.4; to our knowledge, this is the largest sample of low-metallicity dwarf galaxies in the Local Universe detected in CO.

• –

**DGS: ** The Dwarf Galaxy Survey, DGS999Information on the DGS sample, as well as data products, can be found on the website http://irfu.cea.fr/Pisp/diane.cormier/dgsweb/publi.html. (Madden et al. 2013, 2014), is a Herschel sample of 48 local metal-poor low-mass galaxies, with metallicities ranging from 12$+$log(O/H) = 7.14 to 8.43 and stellar masses from $3\times 10^{6}$ to $\sim 3\times 10^{10}\leavevmode\nobreak\ M_{\odot}$. The DGS sample was originally selected from several deep optical emission line and photometric surveys including the Hamburg/SAO Survey and the First and Second Byurakan Surveys (e.g., Ugryumov et al. 1999, 2003; Markarian & Stepanian 1983; Izotov et al. 1991). Although stellar masses are given by Madden et al. (2013) with corrected values in Madden et al. (2014), these are calculated according to Eskew et al. (2012) using the Spitzer/IRAC luminosities at 3.6 $\mu$m and 4.5 $\mu$m. Hunt et al. (2016a) gives $M_{\rm star}$ for these same galaxies by first subtracting nebular continuum and line emission, known to be important in metal-poor star-forming dwarf galaxies (e.g., Smith & Hancock 2009); a comparison shows that the $M_{\rm star}$ values in Madden et al. (2014) are, on average, 0.3 dex larger than those by Hunt et al. (2016a). Thus, in order to maximize consistency with the other samples considered here, like for the KINGFISH, BCDs, and the Virgo star-forming dwarfs (see below), we have used stellar masses based on WISE and/or IRAC 3.4-3.6 $\mu$m luminosities using the recipe by Wen et al. (2013) after subtracting off non-stellar emission estimated from the SFR (see also Hunt et al. 2012, 2015, 2019). Metallicities for the DGS are taken from De Vis et al. (2017), using their PP04N2 calibration. We have also recalculated the SFRs using H$\alpha$ and 24 $\mu$m luminosities as advocated by Calzetti et al. (2010) and reported in Hunt et al. (2016a). Of the 48 DGS galaxies discussed by Rémy-Ruyer et al. (2014), 7 are included in the BCD sample observed in CO by Hunt et al. (2015, 2017); 5 have CO detections from Cormier et al. (2014); and 9 elsewhere in the literature (Kobulnicky et al. 1995; Young et al. 1995; Greve et al. 1996; Walter et al. 2001; Leroy et al. 2005, 2006; Gratier et al. 2010; Schruba et al. 2012; Oey et al. 2017). However, one of these, UM 311, is a metal-poor Hii region within a larger galaxy (see Hunt et al. 2010). There is a discrepancy between the $M_{\rm star}$ values given by Madden et al. (2014) and Hunt et al. (2016a) of more than a factor of 100; this is roughly the difference between the larger UM 311 complex and the individual metal-poor Hii regions, illustrating that the gas content of the individual Hii regions is highly uncertain. We thus eliminate UM 311 from the DGS subset, and include the remaining 13 DGS galaxies that respect our selection criteria.

• –

**Virgo star-forming dwarfs: ** This sample of star-forming dwarf galaxies (SFDs) in Virgo is taken from a larger survey, namely the Herschel Virgo Cluster Survey, HeViCS101010Information on HeViCS and public data can be found in http://wiki.arcetri.astro.it/HeViCS/WebHome. (Davies et al. 2010). Here to supplement the low-mass portion of our combined sample, we have included the dwarf galaxies in HeViCS studied by Grossi et al. (2015, 2016). They were selected from the larger sample by requiring a dwarf morphology (e.g., Sm, Im, BCD) and detectable far-infrared (FIR) emission with Herschel. $M_{\rm star}$ was estimated according to the WISE 3.4 $\mu$m luminosities, and SFRs were calculated using H$\alpha$ luminosities and correcting for dust using WISE 22 $\mu$m emission as proposed by Wen et al. (2014). Both quantities are originally given with a Kroupa (2001) initial mass function (IMF), and have been corrected here to a Chabrier (2003) IMF according to Speagle et al. (2014). Metallicities, 12$+$log(O/H), were based on the SDSS spectroscopy and use the PP04N2 calibration reported by Grossi et al. (2016), with the exception of VCC 1686 for which O/H was derived using the mass-metallicity relation given in Hughes et al. (2013). Of 20 targets observed, 11 were detected in 12CO($1-0$) with the IRAM 30m (Grossi et al. 2016). Atomic hydrogen Hi is detected in all these (Grossi et al. 2016), but of the 11 galaxies with both Hi and CO detections, 4 have Hi–def$>$0.4 (from Boselli et al. 2014a); thus 7 Virgo SFDs satisfy our selection criteria.

• –

Extra single sources: This subset includes individual galaxies that are not included in any survey, but for which our required data of $M_{\rm star}$, SFR, O/H, CO, and Hi measurements exist. These include 7 low-metallicity galaxies: DDO 53 and DDO 70 (Sextans B) from Shi et al. (2016), NGC 3310 (Zhu et al. 2009), NGC 2537 (Gil de Paz et al. 2002), WLM (Elmegreen et al. 2013), Sextans A (Shi et al. 2015), NGC 2403 (Schruba et al. 2012). For these sources, as above for consistency, we used $M_{\rm star}$ and SFR from Hunt et al. (2016a). We were able to compare global $M_{\rm star}$ for one of these, WLM, and once reported to a common distance scale, our value of $M_{\rm star}$ agrees with that from Elmegreen et al. (2013) to within 0.1 dex. Metallicities for these objects are based on the direct Te method, and taken from Engelbracht et al. (2008), Marble et al. (2010), and Berg et al. (2012). In the following figures, the 7 galaxies from these additional sources are labeled as “Extra”.

2.1 Galaxy parameters, data comparison, and potential selection effects

Because of potential systematics that could perturb our results, we first “homogenized” the MAGMA sample by recalculating $M_{\rm star}$ and SFR in a uniform way. The following sections compare the newly-derived values with the original ones described above. We also analyze O/H, CO luminosities, and overall properties of the individual samples in order to assess any systematics that could affect the reliability of our results.

2.1.1 Stellar mass

To estimate $M_{\rm star}$, we use $3.4-3.6$ $\mu$m luminosities, together with a luminosity-dependent M/L. The photometry at 3.4 $\mu$m was acquired from the ALLWISE Source Catalogue (e.g., Wright et al. 2010), taking the photometry measured in an elliptical aperture or within a circular aperture of radius of 49$\aas@@fstack{\prime\prime}$5, whichever was larger. Galactic extinction was corrected for using the Schlafly & Finkbeiner (2011) $A_{B}$ values111111These were taken from the NASA/IPAC Extragalactic Database (NED), funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. and the reddening curve by Draine (2003). Luminosities were calculated from the apparent magnitudes according to the zero points of Jarrett et al. (2013). They were then converted to masses using the M/L ratio for $3.4-3.6$ $\mu$m luminosities given by Hunt et al. (2019), calibrated on the CIGALE (Boquien et al. 2019) stellar masses (see also Leroy et al. 2019):

[TABLE]

where $\nu L_{\nu}$ is the WISE (3.4 $\mu$m) or IRAC (3.6 $\mu$m) luminosities in units of erg s*-1*. Relative to the $M_{\rm star}$ values obtained with detailed SED fitting, Eqn. (1) gives a slightly lower scatter and offset relative to the formulation of Wen et al. (2013), and a negligible offset relative to constant M/L ratios at these wavelengths advocated by various groups (e.g., Eskew et al. 2012; Meidt et al. 2012, 2014; McGaugh & Schombert 2014, 2015).

When WISE photometry was unavailable, or had low signal-to-noise (mainly for the BCDs), we used IRAC 3.6 $\mu$m photometry from Engelbracht et al. (2008) or from Hunt et al. (2020, in prep.). We have assumed that IRAC 3.6 $\mu$m and WISE 3.4 $\mu$m monochromatic luminosities are identical to within uncertainties as discussed in detail by Hunt et al. (2016a) and Hunt et al. (2019) (see also Leroy et al. 2019). All stellar masses are calculated according to a Chabrier (2003) IMF (for more details, see Hunt et al. 2019).

Figure 1 (left panel) compares the original $M_{\rm star}$ values described in the preceding section and the new ones derived here. The mean differences (in log) are reported in Table 1. There are apparently no systematics among the different samples, given that the deviations reported in Table 1 are typically smaller or commensurate with their scatters. However, there is some tendency for the new WISE W1-derived $M_{\rm star}$ to be larger than the originals, typically derived from optical SED fitting. The best agreement is for xGASS-CO, but none of the samples, except for KINGFISH, shows a significant offset. Moreover, for virtually all the samples the scatter is good to within a factor of 2. This corresponds to an uncertainty of $\sim$0.3 dex, consistent with the overall uncertainty of mass-to-light ratios (e.g., Hunt et al. 2019). The KINGFISH galaxies show the same offset relative to SED fitting results found by Hunt et al. (2019). Here we use the best-fit CIGALE-calibrated power-law slope with luminosity, and there the Wen et al. (2013) power-law dependence was used; in any case, the scatter is small because the same photometry (from Dale et al. 2017) is adopted in both cases.

There are four ALLSMOG galaxies that show particularly high discrepancies relative to our homogenized estimates of $M_{\rm star}$: 2MASXJ1336+1552, CGCG 058$-$066, UGC 02004, and VIII Zw 039. The previous stellar masses are roughly an order of magnitude smaller than the new values. We have inspected the SDSS images of these, and they tend to be clumpy, with a series of brightness knots throughout their disks. In these cases, the stellar masses automatically estimated by SDSS tend to regard the clumps, rather than the galaxy as a whole. If these galaxies are eliminated from the comparison of the new homogenized values and the previous values for ALLSMOG, the mean difference (see Table 1) of old minus new log($M_{\rm star}$) becomes $-0.18\,\pm\,0.16$. These galaxies have been retained in our overall analysis.

2.1.2 SFR

Possibly the most difficult parameter is the SFR; the parent samples of MAGMA infer SFR originally using many different methods, ranging from extinction-corrected H$\alpha$ (e.g., ALLSMOG, NFGS: Cicone et al. 2017; Kannappan et al. 2013), to hybrid FUV$+$IR or H\alpha$$+IR (e.g., xGASS-CO, KINGFISH, BCDs, DGS, Virgo SFDs, “Extra”: Saintonge et al. 2017; Hunt et al. 2019, 2016a; Grossi et al. 2015) to the mean of different methods (e.g., HRS: Boselli et al. 2015). To calculate SFR for MAGMA, we have adopted the hybrid formulations of Leroy et al. (2019) based on linear combinations of GALEX and WISE luminosities, estimated for their sample of 15 750 galaxies within distances of $\sim$50 Mpc. Their expressions (see Table 7 in Leroy et al. 2019) are calibrated on the GALEX-SDSS-WISE Legacy Catalogue (GSWLC) by Salim et al. (2016, 2018), which were, in turn, obtained by integrated population synthesis modeling relying on CIGALE fits to $\sim$700 000 low-redshift galaxies. Thus, they are consistent with, and on the same scale, as our CIGALE-calibrated stellar masses. Here we have converted their Kroupa (2001) IMF to the Chabrier (2003) one used here, according to Speagle et al. (2014).

Leroy et al. (2019) give several “recipes” for SFR in hybrid combinations: we have preferentially used the expression with the smallest scatter, namely luminosities of GALEX FUV combined with WISE W4. The W4 luminosities are calculated in the same apertures as the W1 luminosities used for $M_{\rm star}$. For GALEX, we adopted the magnitudes in the Revised catalog of GALEX UV sources by Bianchi et al. (2017) that correspond to integrated values within elliptical apertures, and checked to make sure that the aperture size was commensurate with the WISE apertures. As for the $M_{\rm star}$ estimates, the GALEX and WISE luminosities are corrected for Galactic extinction using the Schlafly & Finkbeiner (2011) $A_{B}$ values and the reddening curve from Draine (2003). According to Leroy et al. (2019), the FUV$+$W4 formulation gives a mean scatter of $\sim$0.17 dex in log(SFR) for the $\sim$16 000 galaxies they analyzed. SFRs derived from FUV$+$W4 were available for 277 MAGMA galaxies ($\sim$71% of the sample), but if not, we adopted NUV$+$W4 (available for 66 galaxies, $\sim$17%), which gives a slightly higher scatter ($\sim$0.18 dex). Overall, these two formulations gave the lowest systematic uncertainties for SFR, and are available for $\sim$88% of the MAGMA galaxies. If GALEX was unavailable, we relied on W4 alone (for 32 galaxies, $\sim$9% of the sample), or otherwise on the hybrid recombination line (H$\alpha$) luminosity combined with 24 $\mu$m luminosities (for 10 galaxies) as prescribed by Calzetti et al. (2010) or on the original SFR value (7 galaxies: 1 ALLSMOG, 1 NFGS, 1 Virgo SFD, Sextans A, DDO 154, and regions of DDO 53 and DDO 73). All SFRs were converted to the Chabrier (2003) IMF.

Overall, as shown in Fig. 1 (right panel), the original SFRs and the new values agree reasonably well, with small mean differences, and always zero to within the scatter (see Table 1). The agreement with the original SFRs from xGASS-CO and HRS is particularly good, with virtually zero offsets and scatters of $\sim$0.2 dex. Both of these samples derived SFR using hybrid schemes, not unlike the ones reported by Leroy et al. (2019) used here. NFGS and the Virgo SFDs show the largest scatters, and for NFGS we attribute this to their use of H$\alpha$ luminosities only, corrected for extinction (see Kewley et al. 2005). The original SFRs for the Virgo SFDs were derived following Wen et al. (2014), but the scatter is dominated by the galaxies with the lowest SFRs (and $M_{\rm star}$); this effect for low-mass, low-metallicity dwarfs was also noted by Wen et al. (2014), so is not unexpected.

2.1.3 Metallicity

As mentioned above, all gas-phase metallicities in our combined sample are either direct Te methods or calibrated through the [Nii]-based PP04N2 Pettini & Pagel (2004) calibration. When the original O/H calibration is not PP04N2, we have converted it to PP04N2 according to Kewley & Ellison (2008). The PP04N2 calibration has been shown to be the most consistent with Te methods (see also Hunt et al. 2016a). Extinction corrections for this calibration are very small because the lines are very close in wavelength: $\lambda$Nii = 6549.86 Å, $\lambda$Nii = 6585.27 Å, H$\alpha$ = 6564.614 Å, so the extinction correction is negligible, for a Cardelli et al. (1989) extinction curve. Even for visual extinction $A_{\rm V}$ = 5 mag, the relative correction is of the order of 1%. This is well within the signal-to-noise of the H$\alpha$ flux itself.

Since metallicity is found to decline from galactic centers to the peripheries (e.g., Kewley et al. 2005; Pilyugin et al. 2014a, b; De Vis et al. 2019), such gradients represent a possible source of systematics; therefore it is worth examining their impact on our metallicity estimates. Some of the O/H for our sample are integrated (e.g., NFGS, HRS, xGASS-CO, ALLSMOG: Kewley et al. 2005; Boselli et al. 2013; Saintonge et al. 2017), some are “representative”, evaluated at 0.4 Ropt (e.g., KINGFISH: Moustakas et al. 2010), and some are nuclear (e.g., the dwarf samples: BCDs, DGS, Virgo SFDs, and “extra” sources). Kewley et al. (2005) have quantified the difference among the global metallicities and the ones measured in the nuclear regions for a sample of 101 star-forming galaxies selected from the NFGS. Independently of the galaxy type, they find that such difference amounts to $\sim 0.1$ dex, a value which is similar to the typical statistical uncertainty of metallicities of 0.1–0.2 dex. Metallicity gradients in late-type dwarf irregulars or BCDs are generally negligible (e.g., Croxall et al. 2009) or at most comparable to those in more massive spirals (e.g., Pilyugin et al. 2015). Thus, we conclude that metallicity gradients should not markedly affect our conclusions.

2.1.4 Molecular gas mass

Like metallicity, molecular gas mass is another delicate issue. Except for ALLSMOG and some KINGFISH galaxies, we use only CO(1–0) in order to avoid excitation issues; as mentioned before, to convert to convert the CO(2–1) values to CO(1–0), a ratio of $R_{21}\,=\,0.8$ was assumed (see also Leroy et al. 2009). Here, the molecular gas masses have been calculated from $L_{\rm CO}^{\prime}$, using the conversion $M_{\rm H2}$ = $L_{\rm CO}^{\prime}$ $\alpha_{\rm CO}$ (where $\alpha_{\rm CO}$ is the H2 mass-to-CO light conversion factor), and adopting a metallicity-dependent calibration, following Hunt et al. (2015). Specifically, for galaxies with $Z/Z_{\odot}<1$ (i.e., 12$+$log(O/H)$<$8.69, see Asplund et al. 2009), we applied $\alpha_{\rm CO}=\alpha_{\rm CO_{\odot}}(Z/Z_{\odot})^{-1.96}$; for metallicities $Z/Z_{\odot}\geq 1$ we used a constant solar value of $\alpha_{\rm CO}=\alpha_{\rm CO_{\odot}}\,=\,3.2\,$M⊙ $({\rm K\,km\,s^{-1}\,pc^{2}})^{-1}$.

As mentioned previously, for MAGMA O/H, we have adopted either Te or the PP04N2 metallicity calibrations. However, for the calculation of $\alpha_{\rm CO}$, we have also investigated the effect of adopting an alternative strong-line calibration, namely the formulation by Kewley & Dopita (2002, KD02). To emulate Bothwell et al. (2016a, b), we have also explored the $\alpha_{\rm CO}$ formulation from Wolfire et al. (2010) using KD02 metallicities. This assumes that the $\alpha_{\rm CO}$ varies exponentially with visual extinction, $A_{\rm V}$, with a weak metallicity dependence for $A_{\rm V}$ (see Bolatto et al. 2013, for more details). These results will be discussed in Sect. 4.1.

Aperture corrections for CO single-dish measurements are also a source of uncertainty, because of the beam size compared to the dimensions of the galaxy. Because of the relatively large distances of the xGASS-CO galaxies (0.025 $<z<$ 0.05 for M_{\rm star}$$\geq 10^{10} M⊙), the median correction to total CO flux applied by Saintonge et al. (2017) is fairly small, a factor of 1.17. Aperture corrections for the ALLSMOG galaxies (Cicone et al. 2017) are even smaller, corresponding to a median covering fraction of 0.98.

We have compared the CO luminosities $L_{\rm CO}^{\prime}$ of the six galaxies common to HRS and KINGFISH. This is an interesting comparison because the KINGFISH galaxies were mapped in CO(2–1) with the HERA CO Line Extragalactic Survey (HERACLES, Leroy et al. 2009), and the HRS measurements were mostly single-dish CO(1–0) observations with few maps (Boselli et al. 2014a). The mean difference between the two datasets is 0.04 dex (larger luminosities for the HRS sample), with a standard deviation of 0.24 dex. Thus there is no systematic difference in $L_{\rm CO}^{\prime}$ between these two samples, and, moreover, the spread is similar to the typical aperture corrections of $\sim$0.2 dex or less for xGASS-CO and ALLSMOG.

2.1.5 Overall properties

Finally, we compare median differences of each parent sample within MAGMA, relative to the sample as a whole. This is done in some detail in Appendix A where we show this comparison graphically. Our analysis shows that ultimately there are no significant systematic differences among the individual parent samples, despite the different criteria for their original selection. We therefore expect that MAGMA, as a whole, is representative of field galaxies in the Local Universe, and can be used to assess the gas scaling relations driving baryonic cycling.

The $M_{\rm star}$ distributions of our combined sample are shown in Fig. 3, while sSFR, $Z$ [measured in units of 12$+$log(O/H)] gas fractions $M_{\rm gas}$/($M_{\rm gas}$ + $M_{\rm star}$) distributions are shown in Fig. 3. The combined MAGMA sample covers the following unprecedented ranges in parameter space, spanning more than 5 orders of magnitude in $M_{\rm star}$, SFR, and $M_{\rm gas}$, and almost 2 orders of magnitude in metallicity121212Here and elsewhere throughout this paper, “log” means decimal logarithm unless otherwise noted.:

[TABLE]

To demonstrate the general applicability of results obtained for MAGMA to the general (field) galaxy population in the Local Universe, we have included in Figs. 3 and 3 the parameter distributions for the SDSS-DR7 catalogue consisting of $\sim$79000 galaxies from Mannucci et al. (2010); hereafter we refer to this sample as SDSS10. For a consistent comparison with MAGMA, like for the samples described above, we have transformed the original O/H calibration from Maiolino et al. (2008) based on KD02 to PP04N2 according to the formulation of Kewley & Ellison (2008) (for more details, see also Hunt et al. 2016a); according to Kewley & Ellison (2008), this transformation has an accuracy of $\sim$0.05 dex. The distributions shown in Figs. 3 and 3 have been renormalized to the number of galaxies in the MAGMA sample, to be able to compare the number distributions directly. The SDSS10 has a relatively narrow spread in O/H; there are 5 MAGMA galaxies ($\sim 1$%) beyond the highest metallicities in SDSS10, and all are from the xGASS-CO sample, corresponding also to some of the most massive galaxies. The MAGMA $M_{\rm star}$ median is 0.6 dex lower than for the SDSS10, and the sSFR median (Fig. 3) is $\sim$0.4 dex higher, illustrating that the MAGMA sample contains more low-mass galaxies than SDSS10. Interestingly, there are only 25 SDSS10 galaxies (0.03%) with log($M_{\rm star}$/M⊙)$>$11.5; thus because of the normalization in Fig. 3 they do not appear.

Both MAGMA and SDSS10 contain a large percentage of massive galaxies relative to local volume-limited samples such as the Local Volume Legacy (e.g., LVL, Dale et al. 2009; Kennicutt et al. 2009) or galaxy-stellar mass functions (GSMF). For the $z\sim 0$ GSMF determined by Baldry et al. (2012) we would expect $\sim$0.1% of the galaxies to be more massive than the break mass, M${}_{*}\,=\,5\times 10^{10}$; $\sim$14% of the galaxies in SDSS10 and $\sim$5% of those in MAGMA are more massive than this. The preponderance of massive galaxies in these two samples, relative to a volume-limited one, is due to flux limits, and the necessity of ensuring spectroscopic measurements (for SDSS10) and CO measurements (for MAGMA). In any case, as shown in Figs. 3 and 3, the parameter coverage of MAGMA does not deviate significantly from SDSS10 at high $M_{\rm star}$ and O/H, but substantially extends the parameter space to lower $M_{\rm star}$ and O/H values.

Given that we required that CO be detected, even at low metallicity, there could be a chance that the MAGMA sample is dominated by starbursts, i.e., galaxies with SFRs significantly above the main sequence. We examine this possibility in Fig. 4, where we compare the distributions for MAGMA galaxies of SFR in different bins of $M_{\rm star}$ with the SDSS10 data from Mannucci et al. (2010) as above. Also shown is the main sequence of star formation given by fitting LVL and KINGFISH by Hunt et al. (2019), here approximated by a Gaussian distribution with a width of 0.3 dex (see also Renzini & Peng 2015). Except for possibly the lowest-mass bins, log($M_{\rm star}$/M⊙) $\leq$ 9, Fig. 4 demonstrates that the MAGMA sample is well approximated by main-sequence SFR distributions. Thus, it is not dominated by starbursts, and can be considered a reliable diagnostic for gas processes in the Local Universe.

2.2 Fundamental scaling relations

In Fig. 5, MAGMA galaxies are plotted in the $M_{\rm star}$–SFR plane, forming the SFMS; the lower panels of Fig. 5 show various forms of the correlation between SFR and $M_{\rm gas}$, a global SK law, exploring different gas phases (atomic, molecular and total, i.e., M_{\rm HI}$$+$$M_{\rm H2}) and CO luminosity, $L_{\rm CO}^{\prime}$. Also shown in Fig. 5 are the loci of the SDSS10 galaxies reported in previous figures. The lowest-level contour encloses 90% of the sample, illustrating the extension by MAGMA to lower $M_{\rm star}$ and SFR. We have included in Fig. 5 also the parameters from the LVL sample as measured by Hunt et al. (2016a); most metallicities are from the direct Te method (Marble et al. 2010; Berg et al. 2012). This sample is the best approximation available for the number and types of galaxies present in the nearby Universe.

Figure 6 shows the mass-metallicity ($M_{\rm star}$–$Z$) relation of the MAGMA sample; although with some scatter, galaxies lie along the MZR over $>4$ dex in $M_{\rm star}$ and a factor of $\sim 50$ in $Z$ (1.7 dex in 12$+$log(O/H)). As in Fig. 5, we have also included SDSS10 and LVL; the MAGMA sample is consistent with both, implying that there are no significant selection effects from our gas detection requirements. Also shown as a blue curve in Fig. 6 is the direct-Te calibration for the SDSS obtained by Andrews & Martini (2013); the MAGMA PP04N2$+$direct Te metallicities are in good agreement with this calibration illustrating that PP04N2 is a good approximation to Te methods (see also Hunt et al. 2016a).

Fig. 6 illustrates that the flattening of the MZR that frequently emerges at high $M_{\rm star}$ is present in the MAGMA and SDSS10 samples. This curvature is clearly seen in the contours of SDSS10 where 90% of the SDSS10 galaxies are enclosed in the lowest contour. Again, the overall extension of MAGMA to lower $M_{\rm star}$ and O/H is evident. In what follows, we focus on applying linear scaling relations to this curved MZR. A single linear relation is not altogether appropriate, given the flattening of the MZR at high $M_{\rm star}$. Indeed, as shown by the linear trend in the upper panel of Fig. 6, it can only roughly approximate the overall MAGMA distribution. In what follows, we investigate the best approach to approximate non-linear trends in the data.

A combination of the scaling relations described above produce the result reported in the lower panel of Fig. 6, namely a tight correlation between the specific SFR (sSFR = SFR/$M_{\rm star}$) and metallicity over $\sim 2$ dex in sSFR and $\ga 1.5$ dex in 12$+$log(O/H). Fig. 6 demonstrates that galaxies more actively forming stars (i.e., with a high sSFR) are less enriched (and also more gas-rich; see Mannucci et al. 2010; Cresci et al. 2012; Hunt et al. 2016b; Cresci et al. 2018, for a discussion).

In the next section we focus on the MZR (Fig. 6), exploring its secondary and tertiary dependencies on SFR and $M_{\rm H2}$, $M_{\rm HI}$, and $M_{\rm gas}$. Since $Z$ is the only intensive 131313Intensive properties are physical properties of a system that do not depend on the system size or the amount of material in the system (e.g., the metallicity of a galaxy does not depend on its size). They differ from extensive properties, that are additive for subsystems. For instance, the total $M_{\rm star}$, $M_{\rm gas}$ and SFR in a galaxy are the sums of the parts composing the galaxy; in other words these quantities depend on the size of a galaxy. quantity among the available integrated physical properties of our galaxies, the MZR is, among the others described above, the most sensitive to internal physical mechanisms.

3 Mutual correlations: a PCA analysis

Principal Component Analysis (PCA) is a parameter transformation technique that diagonalizes the covariance matrix of a set of variables. Consequently, a PCA gives the linear combinations of observables, the eigenvectors, that define the orientations whose projections constitute a hyper-plane; these eigenvectors minimize the covariance and are, by definition, mutually orthogonal. In other words, a PCA performs a coordinate transformation that identifies the optimum projection of a dataset and the parameters that are most responsible for the variance in the sample. The most common use of PCA is to search for possible dimensionality reduction of the parameter space needed to describe a dataset. For instance, a PCA approach has shown that galaxies lie roughly on a 2D surface in the 3D parameter space defined by $M_{\rm star}$, SFR and $Z$ (e.g., Hunt et al. 2012, 2016a) or $M_{\rm star}$, $M_{\rm H2}$ and $Z$ (e.g., Bothwell et al. 2016a).

We use the MAGMA sample to expand upon and re-examine previous trends found with PCAs of $M_{\rm star}$, SFR, metallicity and gas. In addition to the “classical” algorithm for PCA (an unweighted matrix diagonalization), we also apply two additional PCA methods which give the uncertainties associated with the best-fit parameters: the “bootstrap PCA” (BSPCA) first introduced by Efron (1979, 1982) and the “probablistic PCA” (PPCA) described by Roweis (1998). BSPCA is a specific example of more generic techniques that resample the original data set with replacement, to construct “independent and identically distributed” observations. PPCA is an expectation-maximization (EM) algorithm which also accommodates missing information. For the PPCA, we randomly remove 5% of the individual entries for each galaxy; in practice, this means that we omit the SFR for 5% of the galaxies, $M_{\rm star}$ for a different 5%, and so on. For both methods, we generate several realizations of 100–1000 repetitions, and calculate the means and standard deviations of the resulting PC coefficients. For all statistical analysis, we rely on the R statistical package141414R is a free software environment for statistical computing and graphics (https://www.r-project.org/)..

To estimate uncertainties, other groups have used Monte Carlo methods with resampling, injecting Gaussian noise into the nominal measurement values (e.g., Bothwell et al. 2016a, b). We performed several detailed tests using this technique and found that it introduces systematics in the results, depending on the amplitude of the noise and the $M_{\rm star}$ and SFR distributions of the samples; hence we prefer resampling techniques in order to avoid potential unreliability of the results. This point will be discussed further in Sect. 4.3 and Appendix C.

Thus, we perform three kinds of PCAs on: (1) a 4D parameter space defined by $M_{\rm star}$, SFR, $Z$, and a gas quantity (either $M_{\rm gas}$, $M_{\rm H2}$, $M_{\rm HI}$, or $L_{\rm CO}^{\prime}$); and (2) a 3D space defined by $M_{\rm star}$, SFR, and either metallicity $Z$ or a gas quantity (as for 4D). We then assess whether two, three, or four parameters are statistically necessary to describe the variance of these quantities in the MAGMA sample; this is decided by comparing the change in scatter produced by adding an additional variable. $L_{\rm CO}^{\prime}$ has been included as a separate gas quantity in order to separate the effects of true CO luminosity from the effects of a metallicity-dependent $\alpha_{\rm CO}$; this point will be discussed further below. We have also performed a five-dimensional PCA on $M_{\rm star}$, SFR, $Z$, $M_{\rm HI}$, and $L_{\rm CO}^{\prime}$ (or $M_{\rm H2}$), but results do not differ significantly from the 4D case, so we do not discuss it here.

3.1 4D PCA

Results for the 4D PCA are given in Table 2; the coefficients of the PC with the least variance (by definition PC4) are reported, together with the fraction of variance contained in PC4, PC3, and the sum of PC1$+$PC2. There are two separate rows for the PPCA and the BSPCA; these are the different methods used to infer uncertainties, and demonstrate that the coefficients of all the methods agree to within the uncertainties. We find that the proportion of variance contained in only the first two eigenvectors, PC1$+$PC2, is generally large, $\ga$ 94%. Because most of the variance of the sample is contained within the first two eigenvectors, the dimensionality of the parameter space of the MAGMA galaxies is two-fold. They are distributed on a 2D plane in the 4D space; the remaining $\la$ 6% of the variance (shared between PC3 and PC4) produces a scatter perpendicular to this plane.

PC4, the eigenvector with the least variance ($\sim$1%), is always dominated by metallicity, $Z$ (see Table 2). Because very little of the sample variance is contained in PC4, it can be set to zero and inverted to give a useful prediction for the dominant term, namely metallicity $Z$ (see Hunt et al. 2012, 2016a, for a discussion). The estimate for the metallicity obtained by setting PC4 equal to zero is formally accurate to the 1–2% level, i.e., the variance associated with PC4; however, a more robust assessment of the accuracy is obtained by fitting the residuals to a Gaussian as described below.

Interestingly, the contribution of $M_{\rm HI}$ to PC4 is consistent with zero within the uncertainties, and the coefficients for $M_{\rm H2}$ and $M_{\rm gas}$ are small, determined to be non-zero only at a 2$\sigma$ level or less. The PC4 coefficient for all gas components is smaller than that for SFR. The only gas PC4 coefficient strongly different from 0 is $L_{\rm CO}^{\prime}$, determined at $\sim 5-7\sigma$, and comparable in magnitude to the SFR coefficient. This result implies that the 2D plane does not depend significantly on gas properties (except for possibly CO luminosity $L_{\rm CO}^{\prime}$). The expression for 12$+$log(O/H) obtained by inverting the expression based on $M_{\rm H2}$ is:

[TABLE]

where $h_{2}$, $m$, $o$, and $s$ are defined as the centered variables, i.e., log($M_{\rm H2}$), log($M_{\rm star}$), 12$+$log(O/H), and log(SFR) minus their respective means in the MAGMA sample as given in Table 3. The accuracy of this expression is $\sim$0.12 dex, assessed by fitting a Gaussian to the residuals of this fit compared to the observations of 12$+$log(O/H). Eqn. (LABEL:eqn:pca4dh2), in which the uncertainties from Table 2 have been propagated, tells us that:

• –

the gas-phase metallicity of galaxies in the MAGMA sample is primarily dependent on $M_{\rm star}$ (a confirmation of the existence of the well-known MZR);

• –

there is a strong secondary dependence on the SFR, whose contribution in determining $Z$ is $\sim$40% as strong as the dependence on $M_{\rm star}$;

• –

the tertiary dependence on $M_{\rm H2}$ is negligible, virtually zero, given that the accuracy with which the coefficient is determined is $\la 2\sigma$.

We have explored the behavior of the other gas quantities, and as suggested by Table 2, it is similar to the behavior of H2; with the possible exception of $L_{\rm CO}^{\prime}$, the gas content is not influencing metallicity $Z$. However, in the case of $L_{\rm CO}^{\prime}$, the $M_{\rm star}$ coefficient is significantly smaller than that for the other gas quantities and the $L_{\rm CO}^{\prime}$ coefficient has the same sign. It seems that, in some sense, the CO content (not necessarily the H2 content which depends also on metallicity as we have inferred it, see Sect. 2.1) is linked to $Z$, and can partially substitute the influence of $M_{\rm star}$. If we express 12$+$log(O/H) in terms of $L_{\rm CO}^{\prime}$, as we have done in Eqn. (LABEL:eqn:pca4dh2), we obtain:

[TABLE]

where $o$, $m$, $s$ are as in Eqn. (LABEL:eqn:pca4dh2), and $\ell$ is the centered variable log($L_{\rm CO}^{\prime}$). This expression is accurate to $\sim$0.11 dex assessed, as above for $M_{\rm H2}$, by fitting a Gaussian to the residuals of the fit. In reality, this fit should not be interpreted rigorously, since the gas content, rather than CO luminosity, is the quantity of interest. The relation between $L_{\rm CO}^{\prime}$  and the molecular gas content is almost entirely governed by metallicity (e.g., Hunt et al. 2015; Accurso et al. 2017); thus the strong dependence by $L_{\rm CO}^{\prime}$ in Eqn. (LABEL:eqn:pca4dlco) is a reflection of the strong metallicity dependence of the conversion factor $\alpha_{\rm CO}$. We will explore this notion more in detail in a future paper.

3.2 3D PCA

Section 3.1 showed that the 4D parameter space can be approximated by a planar surface, with $\ga$ 94% of the variance contained in the first two eigenvectors, PC1$+$PC2. Here we examine the 3D parameter space (in log space) by retaining $M_{\rm star}$ and SFR as the two main observables, and considering 12$+$log(O/H) as one of the variables together with the four gas quantities described above: $M_{\rm HI}$, $M_{\rm H2}$, $M_{\rm gas}$, and $L_{\rm CO}^{\prime}$. The aim of this exercise is to assess whether any of the gas parameters can be described only by $M_{\rm star}$ and SFR, and to investigate the implication of our 4D PCAs that metallicity 12$+$log(O/H) can be adequately described by $M_{\rm star}$ and SFR alone.

Using the same methodology as for the 4D PCA (“classic” PCA without uncertainties, PPCA and BSPCA with uncertainties), we have performed 3D PCAs on the MAGMA sample, and obtain the results reported in Table 4. Like the 4D PCA, the 3D-PCA component with the least variance is dominated by the metallicity, 12$+$log(O/H) (see column 4 in Table 4). Inverting the PC3 dominated by O/H as before for the 4D PCA, we find:

[TABLE]

The coefficients multiplying log($M_{\rm star}$) and log(SFR) in Eqn. (4) are the same to within the uncertainties as those from the 4D PCA given in Eqn. (LABEL:eqn:pca4dh2). This expression describes 12$+$log(O/H) for the MAGMA sample with an accuracy of $\sim 0.12$ dex, again obtained by fitting the residuals to a Gaussian. The scatter of this expression is comparable to the scatter obtained from the 4D PCA, leading to the conclusion that only $M_{\rm star}$ and SFR are necessary to describe metallicity. The MAGMA coefficient for log($M_{\rm star}$) of $0.37$ is the same as that found by Hunt et al. (2016a).

There are two considerations here: one is that a PCA, by definition, must pass through the multi-variable centroid of the dataset. That is why here, in contrast to Hunt et al. (2012) and Hunt et al. (2016a), we have defined the PCA results in terms of centered variables. This is important when applying a PCA determined with one sample to another sample; if the means of the two samples are significantly different, then the PCA will not pass through the barycenter of the data for the second sample, and will apparently not be a good fit. *Thus a PCA must be applied using the centered variables associated with a particular data set. * The second consideration is that despite the similarity in $M_{\rm star}$ coefficients, the curvature of the MZR present in MAGMA (and SDSS10) is somewhat more pronounced than for the sample analyzed by Hunt et al. (2016a). Fig. 7 shows the coefficient of $M_{\rm star}$ for the subsample with $M_{\rm star}$ $\leq$ $M_{\rm break}$ and $M_{\rm star}$ $>$ $M_{\rm break}$, where $M_{\rm break}$ is the value of $M_{\rm star}$ where one PCA ends and the other one starts. The slopes for $M_{\rm star}$ are systematically smaller for increasing $M_{\rm break}$, because of the flattening of the MZR. In the following section, we explore a remedy for this using an approach more appropriate for data showing non-linear relationships.

3.2.1 3D PCA, a non-linear approach

Several methods have been developed to assess mutual dependencies and dimensionality in a dataset that shows non-linear behavior. In particular, curvature in a dataset can be first approximated by a piecewise linear approach (e.g., Hastie & Stuetzle 1989; Strange & Zwiggelaar 2015; Xianxi et al. 2017). In the case of the curved MZR and its relation with SFR (e.g., Mannucci et al. 2010; Cresci et al. 2018), this is a fairly good approximation as we show below. The fit to the MZR given by Andrews & Martini (2013), shown in Fig. 6, consists of a mainly linear portion toward low $M_{\rm star}$, connected smoothly to a roughly flat regime at high $M_{\rm star}$ (see also Curti et al. 2019). Thus, as a simplified solution to the problem of MZR curvature, we approximate its behavior with two linear segments, and perform a PCA separately on each. Such a procedure is a specific example of a more complex piecewise linear approach, and we postpone a more detailed analysis to a future paper.

The only “free parameter” in the piecewise linear PCA exercise is the break mass, $M_{\rm break}$, namely the value of $M_{\rm star}$ where we establish the transition from one PCA to the other. We have investigated $M_{\rm break}$ between $3\times 10^{9}$ M⊙ and $3\times 10^{10}$ M⊙ (see Fig. 7) and find a “sweet spot” around $M_{\rm break}$ = $2\times 10^{10}$ M⊙, where the overall variance is minimized. The best-fit piecewise 3D-PCA for MAGMA is as follows:

[TABLE]

and the averages of the parameters in the two bins are given in Table 3.

Figure 8 shows the piecewise 3D-PCA results [hereafter “3DPCA2(OH)”] where we compare the predictions of 12$+$log(O/H) from Eqn. (5) and the means given above to the observed values (vertical axis). The parent samples of the individual MAGMA galaxies are given in the legend in the middle panel. The standard deviation of the PC3 component is slightly smaller (0.11 dex vs. $\sim$0.12 dex) than that in the metallicity-dominated 4D PCA with, however, gas content taken into account. The piecewise PC3 standard deviation of 0.11 dex is also slightly smaller than the continuous 3D PCA result without gas, $\sim$0.12 dex. The Gaussian fit to the 3DPCA2(OH) residuals is shown in the right panel of Fig. 8. We expect that the degree to which the piecewise is better than the single PCA depends on the number of galaxies more massive than $\sim 2\times 10^{10}$ M⊙, i.e., the amplitude of the curvature in the MZR.

Also reported in Fig. 8 is the SDSS10 sample, to which we have applied the 3DPCA2(OH) determined from MAGMA; the grey contours enclose 90% of the sample. The mean (median) SDSS10 residuals are $0.000$ ($0.002$)  dex with a standard deviation of 0.08 dex over 78579 galaxies. Thus the MAGMA 3DPCA2(OH) applied to SDSS10 represents the metallicities in that sample with accuracy comparable to the scatter found for the new formulation of the FMR for SDSS by Curti et al. (2019), and with low systematics given the zero mean. For LVL, the scatter is slightly worse: mean (median) LVL residuals are $0.000$ ($-0.02$) dex, with a standard deviation of $\sim$0.2 dex over 135 galaxies. Nevertheless, the small mean (median) residuals indicate that the LVL metallicities are also fairly well approximated by the MAGMA 3DPCA2(OH) even for the low masses in LVL, with a median log($M_{\rm star}$/M⊙) of 8.1, and 25% of the galaxies less massive than log($M_{\rm star}$/M⊙) = 7.3.

Ultimately, comparison of the 4D and 3D PCAs shows that there is no need to include gas content, either $M_{\rm gas}$, $M_{\rm HI}$, or $M_{\rm H2}$, in the description of $Z$; it is statistically irrelevant, since a similar scatter is obtained without any gas coefficient. **This is a clear confirmation that metallicity in field galaxies in the Local Universe can be determined to $\la$0.1 dex accuracy using only $M_{\rm star}$ and SFR. ** However, it is not a statement that metallicity is independent of gas content; on the contrary, in a companion paper, we describe how gas content shapes the MZR through star-formation-driven outflows. As we shall see below, the point is that gas content, like metallicity, can be described through $M_{\rm star}$ and SFR dependencies.

4 Comparison with previous work

Our results are in stark contrast with those of Bothwell et al. (2016b) who, as mentioned above, in a 4D PCA found that H2 mass had a stronger link with metallicity than SFR. Bothwell et al. (2016a) found a similar result based on a 3D PCA, namely that gas content drives the relation between $M_{\rm star}$ and metallicity, and that any tertiary dependence on SFR is merely a consequence of the Schmidt-Kennicutt relation between gas mass and SFR. In a similar vein, Brown et al. (2018) through stacking and Bothwell et al. (2013) found that Hi mass is strongly tied to $Z$, more than to SFR, similarly to later results for $M_{\rm H2}$. We conclude, instead, that metallicity is more tightly linked with stellar mass and SFR than with either $M_{\rm gas}$, $M_{\rm HI}$, or even $M_{\rm H2}$. There are several possible reasons for this disagreement, and we explore them here, with additional details furnished in Appendix B.

4.1 Metallicity calibration and CO luminosity-to-molecular gas mass conversion

We first examine how our results change if we use the same metallicity calibration as Bothwell et al. (2016a, b). This is potentially an important consideration because the KD02 O/H calibration used by Bothwell et al. tends to give metallicities that are too high (e.g., Kewley & Ellison 2008), relative to direct-Te estimates; as shown in Fig. 6, the PP04N2 is a better approximation of these (Andrews & Martini 2013; Hunt et al. 2016a; Curti et al. 2017). Together with using the KD02 calibration, we have also assessed the effect of applying the $\alpha_{\rm CO}$ conversion factor used by Bothwell et al.. The exponential metallicity dependence proposed by Wolfire et al. (2010); Bolatto et al. (2013) depends more steeply on metallicity than the power-law dependence we have used above, as formulated by Hunt et al. (2015). Thus it is possible that the effects of metallicity are enhanced for the molecular gas mass $M_{\rm H2}$ with this approach.

We have thus applied these calibrations to the MAGMA sample, and performed a 4D PCA, as in Bothwell et al. (2016b). The results of this exercise are reported in the lower portion of Table 2. With the KD02 calibration and the exponential metallicity dependence of $\alpha_{\rm CO}$, we find that the PCA4 coefficients are slightly altered: the $M_{\rm star}$ and SFR coefficients are larger in amplitude and O/H coefficient is smaller. The H2 term is even smaller than with our original formulation, and zero to within the uncertainties. In agreement with our original formulation, we would have concluded that H2 has negligible impact relative to SFR. Thus, the different approaches for $\alpha_{\rm CO}$ and the metallicity calibration are probably not the cause of the disagreement.

4.2 Differences in sample sizes and properties

Here we examine whether the larger MAGMA sample, its significant low-mass representation, and different SFR relations can influence PCA results. Our MAGMA sample of 392 galaxies is nominally twice as large as the sample studied by Bothwell et al. (2016a, b). However, if we consider only the CO detections in their low-$z$ sample (141 galaxies), judging from Table 2 in Bothwell et al. (2016a), our sample is almost three times larger. Moreover, MAGMA contains a much larger fraction of low-mass galaxies, as it includes the HeViCs dwarf galaxies (Grossi et al. 2015), the DGS (Cormier et al. 2014), the BCDs not yet published by Hunt et al., and DDO 53, Sextans A, Sextans B and WLM, the extremely metal-poor galaxies studied by Shi et al. (2015), Shi et al. (2016) and Elmegreen et al. (2013). The MAGMA mean log($M_{\rm star}$/M⊙) = 9.7 is $\sim$ 3 times lower than the mean log($M_{\rm star}$/M⊙) = 10.2 of the 158 (including high-$z$) detections in the Bothwell et al. (2016a) sample; while 24% (94) of the MAGMA galaxies have M_{\rm star}$$\leq 10^{9.3} M⊙, this is true for only 7% (11) of the Bothwell et al. (2016a) galaxies, and for $\sim$11% (18) of those in Bothwell et al. (2016b).

Nevertheless, the most important difference between the MAGMA sample and the Bothwell et al. sample(s) is the inclusion of galaxies at high redshift in the latter. As shown in Appendix B, the addition of these galaxies significantly increases the amplitude of the $M_{\rm H2}$ term in the 4D PCA, and reduces that of SFR. When the $z\sim 2$ galaxies are not included, the results of a 4D PCA on the Bothwell et al. sample are ambiguous, because the metal content is found to increase with increasing SFR, similarly to the increase with $M_{\rm star}$. However the statistical significance of this result is low, and the sample is ill conditioned because of the behavior of SFR with $M_{\rm star}$ in the sample.

4.3 Methodology comparison and parameter uncertainties

In Appendix C, we assess the consequences of introducing Gaussian noise on an observing sample that is to be subject to a PCA. After constructing several sets of mock samples based on well-defined input scaling relations, we conclude that the accuracy with which the original relations can be retrieved depends on the amount of noise injected. It is fairly common to calculate uncertainties on fitted parameters by injecting noise in a sample and repeating the exercise several times (e.g., Bothwell et al. 2016a, b). However, our results show that this process skews the data because of the broader range in the parameter space, and relative importance of outliers in a PCA.

Another important consideration is the importance of the $M_{\rm star}$ distribution of a sample like the one considered here. At a given level of noise injection $\sigma$, we found that the broader the range of $M_{\rm star}$, the more consistent with the original “true” scaling relations will be the results.

In some sense, as we show in Appendix C, observing samples such as MAGMA already contain noise, and adding more will skew results, compromising reliability. Ultimately, these are the reasons we chose to apply probabilistic PCA and boot-strap PCA with sample replacement, rather than perturb the parameters of the sample by injecting noise.

5 Summary and conclusions

With the aim of investigating the role of gas on the mass-metallicity relation, we have compiled a new ‘MAGMA’ sample of 392 galaxies covering unprecedented ranges in parameter space, spanning more than 5 orders of magnitude in $M_{\rm star}$, SFR, and $M_{\rm gas}$, and almost 2 orders of magnitude in metallicity. Basic galaxy parameters, $M_{\rm star}$ and SFR, have been recalculated using available data from IRAC, WISE, and GALEX archives, and all O/H values have been converted to a common metallicity calibration, PP04N2. All stellar masses and SFRs rely on a common Chabrier (2003) IMF, and the combined sample has been carefully checked for potential systematics among the sub-samples.

Applying 4D and piecewise 3D PCAs to MAGMA confirms previous results that O/H can be accurately ($\lesssim 0.1$ dex) described only using $M_{\rm star}$ and SFR. However, our findings contradict earlier versions of PCA dimension reduction on smaller samples, as we find that the O/H depends on SFR more strongly than on either Hi or H2. Thus, even though a PCA shows mathematically that only a 2D plane is necessary to describe metallicity $Z$ or $M_{\rm gas}$ (or $M_{\rm H2}$), the dependence of $Z$ on gas content is not well constrained with a PCA.

In Sect. 3.1, a 4D PCA showed that the four parameters $M_{\rm star}$, SFR, 12$+$log(O/H), and $M_{\rm H2}$ or $M_{\rm gas}$ are related through a 2D planar relation, with metallicity as the main primary dependent variable. This implies that O/H depends primarily on $M_{\rm star}$ and SFR, but also that $M_{\rm gas}$ must depend primarily on these two variables because of the physical connection between gas content and metallicity.

The observational scaling relations here for O/H are applicable to isolated (field) galaxies in the Local Universe, over a wide range of stellar masses, SFR, and metallicities 12$+$log(O/H) $\ga$ 7.6. They can be used as a local benchmark for cosmological simulations and to calibrate evolutionary trends with redshift. Future papers will consider relations among gas content, star-formation, and metal-loading efficiencies, as well as detailed comparison with evolutionary models.

Acknowledgements

We are grateful to the anonymous referee whose critical comments improved the paper. The authors would also like to thank A. Baker, F. Belfiore, S. Bisogni, C. Cicone, M. Dessauges-Zavadsky, Y. Izotov, S. Kannappan, R. Maiolino, P. Oesch and D. Schaerer for helpful discussions. We acknowledge funding from the INAF PRIN-SKA 2017 program 1.05.01.88.04. We thank Yong Shi for passing us their results in digital form. We have benefited from the public available programming language Python, including the numpy, matplotlib and scipy packages. This research made extensive use of ASTROPY, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), and glueviz, a Python library for multidimensional data exploration (Beaumont et al. 2015). This research has also made use of data from the HRS project; HRS is a Herschel Key Programme utilizing Guaranteed Time from the SPIRE instrument team, ESAC scientists and a mission scientist. The HRS data was accessed through the Herschel Database in Marseille (HeDaM - http://hedam.lam.fr) operated by CeSAM and hosted by the Laboratoire d’Astrophysique de Marseille.

Appendix A Overall comparison of galaxy parameters in MAGMA

Figure 9 shows the median differences of each individual parent samples in MAGMA relative to the sample as a whole as a function of discrete $M_{\rm star}$ bin; the error bars are the standard deviations of the estimate. The grey regions give the standard deviations of the entire MAGMA sample, and the horizontal dashed lines the median difference, zero by definition. If there were any systematic differences as a function of parent sample and/or $M_{\rm star}$, they should stand out in Fig. 9. However, it is clear from the figure that there are no systematic differences in any of the parameters shown. The only exception could be the 12$+$log(O/H) (bottom) panel, where the BCD, DGS, Virgo SFDs have O/H lower than would be expected; there is a similar corresponding excess, although within the spreads, in SFR. This is a real, physical difference, among the samples, due to the anti-correlation between SFR and O/H, the FMR (see, e.g., Mannucci et al. 2010). These dwarf samples at low $M_{\rm star}$ are slightly starburst biased because of our requirement for CO detections. Thus, with this exception, we conclude that there are no apparent systematic differences among the parent samples in MAGMA.

Appendix B Details of comparison with previous work: a cautionary tale

There are two salient differences between the MAGMA sample and the sample from Bothwell et al. (2016b). One is the significant low-mass coverage of the MAGMA sample, and the other is the SFRs. Both these are illustrated in Fig. 10 where the MS for the two samples are plotted; the median points show only those $M_{\rm star}$ bins with $\geq\,3$ data points. The Bothwell et al. (2016b) sample is the same as that from Bothwell et al. (2016a), but with 8 BCDs from Hunt et al. (2015); Bothwell et al. (2016b) altered the $M_{\rm star}$ values, and to best approximate their sample, for the BCDs from Hunt et al. (2015) we have arbitrarily lowered the $M_{\rm star}$ values by a factor of three151515Bothwell et al. (2016b) has used a different technique to estimate $M_{\rm star}$, and find roughly a factor of 3 times lower values, but there is no tabulation of their modified values.. Even with the addition of the BCDs from Hunt et al. (2015), the mass distribution is significantly more extended for MAGMA. The low-mass bin in the Bothwell et al. sample of Fig. 10 are the galaxies from Hunt et al. (2015) that are not present in the sample used by Bothwell et al. (2016a), and the highest-mass bin is dominated by the high-$z$ sub-millimeter galaxies (SMGs) and MS galaxies at $z\sim 2$. MAGMA shows slight downward inflection associated with quiescent or passive high-mass galaxies (see also Fig. 5), while the Bothwell et al. (2016b) sample shows a steep upturn, because of the inclusion of high-$z$ massive galaxies with consequently higher SFR.

There are also differences in the SFRs: while MAGMA adopts COLDGASS galaxies with parameters from Saintonge et al. (2017), Bothwell et al. (2016a, b) use the COLDGASS parameters reported by Saintonge et al. (2011a, b). However, a subsequent study (Huang & Kauffmann 2014) showed that the SFRs in those papers are overestimated by $\sim 0.2$ dex because of aperture effects from the CO single-dish IRAM beam. Roughly half (115 galaxies) of the Bothwell et al. sample is from COLDGASS which means that these values are also discrepant with respect to MAGMA. This can be seen in Fig. 10 where, at a given $M_{\rm star}$, the Bothwell et al. sample tends to have higher SFRs than MAGMA.

Possibly the most important difference in this context relative to MAGMA is the inclusion by Bothwell et al. of the 17 high-$z$ galaxies (9 main-sequence galaxies and 8 SMGs). These galaxies are at $z\sim 2$, and have significantly higher SFRs than the local galaxies of similar stellar mass, because of the increasing normalization of the star-formation main sequence with redshift (e.g., Speagle et al. 2014). The minimum Log($M_{\rm star}$/M⊙) value in the high-$z$ sample is 9.78, and the mean SFR $\sim 274$ M⊙ yr*-1*; for the low-$z$ sample over the same mass range, the mean SFR $\sim 2.4$ M⊙ yr*-1*, more than 100 times smaller. A similar difference applies to the ratios of $M_{\rm H2}$ in the two samples, where the mean $M_{\rm H2}$ in the high-$z$ sample is $\sim 85$ times higher than in the low-$z$ galaxies over the same range in $M_{\rm star}$. To assess the impact of these galaxies on the results by Bothwell et al. (2016a, b), we have performed 4D PCAs on the Bothwell et al. (2016b) sample, with and without the 17 high-$z$ galaxies, using only the galaxies with CO detections in their Table 2 (and Hunt et al. 2015). Results are reported in Table 5.

The upper part of Table 5 shows the 4D PCA for H2 for the Bothwell et al. (2016b) sample including the 17 high-$z$ MS galaxies and SMGs and the 8 BCDs from Hunt et al. (2015). The results are in fairly good agreement161616This sample is not quite the same as that analyzed by Bothwell et al. (2016b) because the 8 galaxies from Hunt et al. (2015) are included with an arbitrary factor of 3 lower $M_{\rm star}$, since Bothwell et al. (2016b) recalculated their $M_{\rm star}$ values but did not tabulate the new $M_{\rm star}$ values. Moreover the numbers of galaxies do not apparently match; here we are analyzing only the 158$+$8 CO detections given in Table 2 of Bothwell et al. (2016a) with 8 BCDs (with $M_{\rm star}$/3 from Hunt et al. 2015). with those of Bothwell et al. (2016b). The dependence of 12$+$log(O/H) on $M_{\rm H2}$ is larger than that on SFR; our probabilistic PCA estimates of the uncertainties show, however, that the coefficient for the $M_{\rm H2}$ dependence is determined with $<3\sigma$ significance.

The lower part of Table 5 shows instead the 4D PCA result for the low-$z$ Bothwell et al. (2016a) sample, without the 17 high-$z$ galaxies, but with the 8 BCDs from Hunt et al. (2015). The $M_{\rm star}$ dependence is significantly reduced, the SFR dependence is increased, and the $M_{\rm H2}$ dependence is consistent with 0.0 to within the errors. Interestingly, the SFR and $M_{\rm star}$ coefficients have the same sign, implying that increasing both SFR and $M_{\rm star}$ will increase $Z$; this is contrary to the “normal” 3DPCA1–OH behavior in which at a given $M_{\rm star}$, increasing SFR tends to reduce $Z$. The comparison in Table 5 of the two versions of the Bothwell et al. sample, including or omitting the high-$z$ galaxies, shows that the PCA significantly changes, and that the most probable driver of the lack of metallicity SFR dependence relative to $M_{\rm H2}$ is caused by the inclusion of $z\sim 2$ galaxies which, at a given $M_{\rm star}$, have a significantly higher molecular gas content and SFR than galaxies in the Local Universe.

Appendix C Assessment of Monte Carlo error injection on PCA fits

To explore the effect of injecting Gaussian noise on a dataset subject to a PCA, we generated several “mock” samples of the 3DPCA1–OH for galaxies at $z\sim 0$. To do this, we first distributed numbers of galaxies in $M_{\rm star}$ bins with redshifts ranging from $z\,=\,0.0$ to $z\,=\,0.02$, according to either a constant $M_{\rm star}$ distribution or one that resembles the GSMF given by Baldry et al. (2012). Within each mass and (small) redshift bin, we selected randomly $M_{\rm star}$ in order to more or less reproduce the assumed distribution. Then we derived a MS of star formation by imposing Speagle et al. (2014) at $z\,=\,0$, or equivalently adopting the relation given by Hunt et al. (2019) for the KINGFISH sample. We added a small (realistic) fraction of starbursts using the formulation of Sargent et al. (2012); this approach separates galaxy populations according to sSFR, and approximates the distribution with two Gaussians. Sargent et al. (2012) further assume that the starburst fraction is independent of mass and redshift, which may or may not be correct (see e.g., Bisigello et al. 2018). Finally, we related 12$+$log(O/H) to $M_{\rm star}$ and SFR via the 3DPCA1–OH (FPZ) reported by Hunt et al. (2016a). This means that the basic input parameter is $M_{\rm star}$, which defines SFR through the MS with the addition of a small fraction of starbursts; then 12$+$log(O/H) is calculated based on $M_{\rm star}$ and SFR. We adopted a total mock sample size for both $M_{\rm star}$ distributions of $\sim$12 000 simulated galaxies.

This construction of the mock samples may seem arbitrary, but in truth the details are not important; we only want to compare what we get out with what we put in, in the case of varying levels of noise injection. To this end, we took this initial “noiseless” sample and introduced varying degrees of Gaussian uncertainty, $\sigma$, to SFR and 12$+$log(O/H); for simplicity, we used the same value of $\sigma$ for both SFR and 12$+$log(O/H). We repeated this procedure several times, and applied a PCA to each of the noise-injected samples. The results are shown in Fig. 11, where the 3D PCA coefficients are plotted against the injected noise level $\sigma$. The left-hand part of the curves for $\sigma\,=\,0$ are the input values of the 3DPCA1–OH by which the mock sample was generated. There are seven separate curves in Fig. 11 for each of the three PCs, corresponding to seven different realizations of the noise injection for the mock samples; the closeness of the curves evident in the figure indicates that the statistical results are quite stable. The discrepant curves for the constant $M_{\rm star}$ case correspond to different lower-mass limits.

The idea here is to simulate an observed sample such as MAGMA, and assess the accuracy of the resulting PCA, compared to the input “true” values. The implicit assumption is that MAGMA, or similar samples, are governed by an underlying 3DPCA1–OH or 2D plane, but suffer from uncertainty in the measurement of the observables. We cannot know whether or not this is true; we can only ascertain how far the observed data set could differ from the underlying relation if it were present.

Also shown in Fig. 11 are six data points corresponding to a PCA on the $\sigma\,=\,0.15$ mock sample, but to which an additional perturbation has been applied. We have chosen $\sigma\,=\,0.15$ dex as a starting point, in order to best approximate the behavior of the MAGMA sample which shows a dispersion in 12$+$log(O/H) of roughly this amplitude around the best-fit 3DPCA1–OH. While the original mock samples are intended to reproduce observed samples, this additional injection of Gaussian noise is aimed at simulating a Monte Carlo perturbation of an observed sample. Thus, we injected Gaussian noise of various amplitudes on the variables of our mock sample: Log($M_{\rm star}$) was varied by 0.3 dex; SFR was varied by 30% [corresponding to $\sim 0.2$ dex on Log(SFR)]; 12$+$log(O/H) was varied by 0.1 dex. This is a similar technique to that described by Bothwell et al. (2016a, b), and in principle helps to establish uncertainties in the final PCA results. However, Fig. 11 shows clearly that the injection of additional noise on the mock sample carries the PCA results even further from the input true 3DPCA1–OH. The amplitude of the noise injected $\sigma$ already masks the 3DPCA1–OH that was the basis for the mock samples, but the additional perturbation exacerbates even more the capacity of the PCA to diagnose the underlying relation.

Fig. 12 illustrates in another way the process of this subsequent perturbation on the mock sample; the orthogonal projections of the PCs are shown in three different panels. The color maps show the original mock sample with $\sigma\,=\,0.15$ injected Gaussian noise, while the contours show the same data as the individual points (with the GSMF from Baldry et al. 2012), but here with the PCs calculated according to the loadings of the original mock sample. This is done to highlight the change in orientation relative to the original sample, thus illustrating that the injection of additional noise alters the orientation of the components. This can be seen in particular in the upper right panel which compares PC3, dominated by O/H, to PC2, which is dominated by $M_{\rm star}$, and to a lesser extent SFR; the contours are oriented at a different angle relative to the underlying color map.

We conclude that:

• •

The injection of Gaussian noise in a noiseless sample changes the PCA characteristics, because of the resulting change of the orientation of the derived PCs (see the curves in Fig. 11). Further injection of Gaussian noise moves the PCA even further from the input relation, as shown by the individual points in Fig. 11. Even though the introduction of noise does not change the means of the parameters, it skews the orientation because PCAs consider the entire distribution of data, including outliers.

• •

The distribution of $M_{\rm star}$ in a sample also impacts the results of a PCA (see varying curves in Fig. 11). This is because a PCA calculates the orthogonal distance from an orientation, and requires a broad distribution in parameters in order to better take into account eventual outliers.

Ultimately, because of the mathematical nature of the PCA, the addition of noise to a sample for which a PCA is to be performed is deleterious for the reliability of the results. For this reason, we have preferred here to use resampling techniques, rather than alter the noise characteristics of the sample.