© The Authors 2024

1. Introduction

The gas-phase metallicity (Z_g) of local galaxies has proven to be a key observable for numerical simulations and theoretical models. As a consequence of recurrent star formation processes, the amount of metals in galaxies increases gradually. In this way, every internal and external process that a galaxy experiences leaves imprints on the chemical distribution of its gas. This parameter has been widely confirmed as a powerful tool for improving our knowledge of the formation and evolution of galaxies (Mosconi et al. 2001; Lia et al. 2002; Brooks et al. 2007; Oppenheimer & Davé 2008; Davé et al. 2011; Pilkington et al. 2012; Mollá et al. 2019; Spitoni et al. 2019, among many others).

Many works have reported significant correlations both globally and locally between the gas metallicity and other galaxy properties, such as the stellar mass (M_⋆, Lequeux et al. 1979; Tremonti et al. 2004; Sánchez et al. 2017), the star formation rate (SFR, Mannucci et al. 2010; Sánchez-Menguiano et al. 2019; Sánchez Almeida & Sánchez-Menguiano 2019), the stellar age (Lian et al. 2015; Sánchez-Menguiano et al. 2020), the rotation velocity (Garnett 2002; Pilyugin et al. 2004), the gas mass fraction (Bothwell et al. 2013; Barrera-Ballesteros et al. 2018), and the galaxy size (Ellison et al. 2008; Sánchez Almeida & Dalla Vecchia 2018). Of all galaxy properties, M_⋆ has been shown to present the strongest correlation with Z_g, following the so-called mass–metallicity relation (MZR). As the stellar mass increases, the metallicity increases, until the relation flattens at high masses (e.g., Tremonti et al. 2004; Kewley & Ellison 2008; Sánchez et al. 2013, 2017, 2019; Wu et al. 2016; Barrera-Ballesteros et al. 2017; Curti et al. 2020). The MZR is found to exist up to z ∼ 3.5 with a similar shape, although it evolves with redshift such that high-z galaxies are less enriched than local ones for a given stellar mass, especially at low masses (Maiolino et al. 2008; Mannucci et al. 2009; Yabe et al. 2012; Maier et al. 2014; Zahid et al. 2014; Sanders et al. 2015; Hidalgo 2017). This relation has been proposed to be the result of the interplay of many processes such as metal removal by outflows, dilution by metal-poor gas infall, enrichment by previous episodes of star formation, or simple evolution by the so-called downsizing, where most massive galaxies exhaust their gas reservoir faster (Maiolino & Mannucci 2019).

Despite the relevance confirmed for the MZR in the literature, a recent study by D’Eugenio et al. (2018) suggested that gas metallicity correlates more tightly with the average gravitational potential (estimated as Φ = M_⋆/R_e, where R_e is the effective radius of the galaxy) than with either the stellar mass or the average mass surface density (${\mathrm{\Sigma }}_{\star }\equiv M_{\star }/R_{\text{e}}^2$). Using a volume-limited sample drawn from the publicly available Sloan Digital Sky Survey Data Release 7 (SDSS DR7), their conclusions were based on (i) a lower scatter in the relation, (ii) a higher Spearman rank correlation coefficient, and (iii) a weaker residual trend with R_e. The three relations were explained in terms of metal enrichment (driven by the gas fraction, and therefore Σ_⋆) and metal loss (proportional to the stellar mass). Alternatively, they could emerge as a consequence of an existing local relation between Z_g and the escape velocity (e.g., reported in Barrera-Ballesteros et al. 2018), a local quantity directly related with the global gravitational potential. In this way, the potential–metallicity relation would be the primary relation, and the MZR and the ΣZR would arise from the local relation between Z_g and Φ. This result contradicts another recent study by Baker & Maiolino (2023), who find a tighter relation between Z_g and M_⋆ than between Z_g and the dynamical mass (M_dyn) or Z_g and the gravitational potential (estimated as Φ = M_dyn/R_e). We note, however, that dynamical masses were derived via Jeans anisotropic modelling, which is a complex method that might increase the final scatter and thus undermine its importance.

In general, studies focus on investigating particular relations between the gas metallicity and other galaxy properties, but we lack works that address the issue in a holistic approach, involving a large number of parameters in the analysis. This is key for assessing with accuracy which property predicts Z_g best, that is, which property yields the tightest relation with Z_g, and so more probably is the primary relation. In this regard, the recent study by Alvarez-Hurtado et al. (2022) explored the impact of 29 physical parameters on the gas metallicity for a sample of 299 star-forming galaxies from the Calar Alto Legacy Integral Field Area survey (CALIFA, Sánchez et al. 2012). The results of the analysis indicate that the stellar mass is the physical parameter that traces the observed Z_g better. However, they did not include the gravitational potential as a parameter in the analysis.

This paper describes our work to provide a robust answer to the question. In order to do this, we applied a random forest regressor algorithm on a sample of more than 3000 nearby star-forming galaxies from the Mapping Nearby Galaxies at Apache Point Observatory survey (MaNGA, Bundy et al. 2015). Using this technique, we explore the effect of 148 parameters on the global oxygen abundance of the galaxies as an indicator of the gas metallicity. We include in the model all properties that have previously been reported to strongly correlate with Z_g, such as M_⋆, Φ, Σ_⋆, and M_gas.

The paper is organised as follows. Section 2 briefly presents the MaNGA data, and the description of the procedure we followed in the analysis is given in Sect. 3, including the physical parameters of the model (Sect. 3.1), the criteria we adopted for the sample selection (Sect. 3.2), and an overview of the random forest algorithm (Sect. 3.3). The outcome of the random forest is described in Sect. 4.1, and a subsequent analysis including different combinations of M_⋆ and R_e is described in Sect. 4.2. All the results are discussed in Sects. 5, and Sect. 6 compiles and summarises the main conclusions of the work. Finally, supplementary material is included in Appendices A, B, and C. Appendix A presents the complete list of all the galactic parameters included in the model. Appendix B reproduces the analysis using alternative calibrations to estimate Z_g. Appendix C describes some tests on how the inhomogeneous mass coverage of the sample might affect the results.

2. MaNGA data

We took advantage of the large statistics provided by the MaNGA survey (Bundy et al. 2015), a now complete project based on integral field spectroscopy (IFS) and part of the fourth-generation SDSS (SDSS-IV). With all the data publicly released (DR17, Abdurro’uf et al. 2022), MaNGA has gathered spatially resolved information for 10 010 galaxies up to redshift ∼0.15.

The data were collected using the two BOSS spectrographs mounted on the Sloan 2.5 m telescope at Apache Point Observatory (Gunn et al. 2006). The 17 simultaneously available fibre bundles are distributed in five different hexagonal configurations, whose field of view (FoV) varies from 12.5″ to 32.5″ in diameter (Drory et al. 2015). The covered wavelength range spans from 3600 Å to 10 300 Å, with a nominal resolution of λ/Δλ ∼ 2100 at 6000 Å (Smee et al. 2013). The resulting spectra for each sampled spaxel of 0.5″ × 0.5″ present a final spatial resolution of FWHM ∼ 2.5″, which corresponds to a physical resolution of ∼1.5 kpc at an average redshift of 0.03.

Details of the MaNGA mother sample, the survey design, the observational strategy, and the data reduction are provided in Law et al. (2015;, 2016, 2021), Yan et al. (2016), and Wake et al. (2017).

3. Analysis

3.1. Galaxy parameters

In order to determine the galactic property that best predicts the gas metallicity, we explored the effect of 148 parameters. These parameters were extracted from the pyPipe3D Value Added catalogue (VAC, Sánchez et al. 2022), which is publicly accessible through the SDSS-IV VAC website¹. This VAC contains the dataproducts of pyPipe3D for the full DR17 MaNGA sample (Lacerda et al. 2022, see also Sánchez et al. 2016a,b). This analysis pipeline was developed to characterise the properties of the stellar populations and the ionised gas. The catalogue comprises a single table with the integrated values, those measured at the effective radius, and the radial gradients of different properties (e.g., stellar mass, star-formation, oxygen and nitrogen abundances, or dust attenuation), as well as a FITS file per galaxy with the corresponding spatially resolved quantities.

For this study, we focused on the global galaxy properties provided in the VAC table. This table comprises 304 parameters (and in most cases, their associated error estimates), many of which correspond to different line intensity values or ratios (e.g., the logarithm of the [NII]6583/Hα line ratio in the central 2.5″ aperture). These quantities were discarded from the study because including them would result in obvious dependences of some galaxy properties and the emission lines from which they are derived (e.g., between the star formation rate and the Hα intensity). In addition, as a proxy for the global gas metallicity (Z_g), we used the oxygen abundance (O/H) measured at one disc effective radius (R_e)² and determined from the empirical calibration proposed by Marino et al. (2013, hereafter M13) for the O3N2 index. Other Z_g calibrators were employed and led to results similar to those described in Sect. 4.1 for O3N2 (see Appendix B for a deeper discussion). All these additional measurements of O/H were also discarded from the main analysis for clear reasons. Finally, the ionisation parameter, the nitrogen abundance, and the N/O abundance were also ignored because their determination involves line ratios in common with the O/H derivation. All in all, the final selection of galaxy properties available for the analysis (N_par) comprises 148 parameters. Table 1 compiles the most representative galactic properties, and the complete list is detailed in Appendix A.

3.2. Sample selection

The galaxy sample we analysed was drawn from the 10 010 galaxies comprising the MaNGA mother sample. In order to select the galaxies, we adopted two simple criteria: Galaxies have to meet the required quality standards of the analysis pipeline (i.e., a QCFLAG field equal to zero in the pyPipe3D VAC table; see Sect. 4.5 of Sánchez et al. 2022 for details), and they must provide values for all of the analysed parameters. The second criterion is quite restrictive. For instance, it removes early-type systems and galaxies in general with little or no ionised gas (for which gas properties such as the oxygen abundance or the dust extinction cannot be derived). The final sample contains 3430 galaxies. This number is high enough to ensure the statistical significance of the results. However, we note that the sample is biased towards large star-forming galaxies.

3.3. Random forest regressor

We applied a supervised machine-learning technique known as random forest (RF) regressor (Breiman 2001). Through a combination of decision trees, the algorithm finds the input features (the galaxy properties in our case) that contain more information on the target feature (Z_g) and creates a model to predict it by defining a set of conditions on the values of the input features. RF regressors have been extensively used in astronomy in the past decade with high levels of success (e.g., Carliles et al. 2010; Carrasco Kind & Brunner 2013; Sánchez-Menguiano et al. 2019; Moster et al. 2021, among many others). The algorithm was implemented using the scikit-learn package for Python (Pedregosa et al. 2011). We briefly describe the basic steps we followed to run the algorithm, but we refer to the scikit-learn User Guide documentation for more details of the complete algorithm implementation³.

We first randomly split the sample of 3430 galaxies into two subsets: One set was used to create the model (training sample; 2572 galaxies, corresponding to two-thirds of the sample), and the other set was used to evaluate its performance (test sample; 858 galaxies, comprising the remaining one-third of the sample). For the training sample, the selected set of galaxy properties was provided to the algorithm to train the model (these are the predictors). Z_g was considered the target feature that was to be predicted by the RF model (this is the solution). The algorithm requires not only the predictors, but also the solutions (called labels) to train the model. When the model was trained, we applied it to the test sample, obtaining a set of predictions for Z_g. The comparison of these predictions with the measured values provides an estimate of the model precision.

Before training the model, we selected the values of the parameters in the algorithm (called hyperparameters, HP) that optimise its performance. The HPs include the number of trees in the forest (n), the number of randomly selected features to consider in each split (m), the maximum depth of the trees (max_depth), the minimum number of samples required in each split (min_split), and the minimum number of samples that remain at the end of the different decision tree branches (min_leaf). We set m = N_par = 148 to fully disentangle the relative interdependences among the input parameters, allowing the algorithm to select any of them at each split. If we choose m < N_par, a variable containing significant information on the target will be replaced by any correlated parameter when it is not included among the randomly selected features to consider in a split. This will falsely reduce the importance of this variable in the model. By setting m = N_par, any small difference in the performance of two correlated parameters will be detected by the RF, breaking possible degeneracies in the relations of the parameters. We refer to Bluck et al. (2022, Appendix B.2) for detailed tests of how well RF can perform with inter-correlated data. Nevertheless, this decision must be balanced with a careful choice of the other HPs to avoid overfitting the data. In order to select the remaining HPs, we used the fivefold cross-validation (CV) method (James et al. 2013), which consists of splitting the training sample into five subsets, called folds. We trained the model on a different combination of four folds and evaluated it on the remaining fold. The reported performance measure was then the average of the values computed in the loop. We performed 100 iterations of the entire 5-Fold CV process, each time using different values for the HPs. We selected the HP values that achieved the highest average performance across the five folds. In our case, this resulted in n = 700, max_depth = 30, min_split = 4, and min_leaf = 2. We note, however, that the results we obtained are not very sensitive to the selected values for the HPs.

Finally, using the selected HPs, we trained the model on the full training sample, and then evaluated it on the test sample. Figure 1 represents the predicted Z_g by the RF algorithm against the measured values for the training (blue circles) and test (brown triangles) samples. The test sample follows a trend around the one-to-one relation (dashed line), similar to the training sample (R² score of 0.87 and 0.98, respectively), showing that the algorithm is able to predict the Z_g of galaxies that are not used to train the model with comparable accuracy. This ensures that the model is not affected by overfitting. The median error in the predictions (predictions-targets) is 0.017 dex, which is about the same as the observational errors propagated into Z_g (∼0.011 dex, which does not consider the systematic error in the calibrator).

Fig. 1. Global gas-phase metallicity – described by the oxygen abundance 12 + log(O/H) – predicted by the RF algorithm vs. the measured values. The blue circles correspond to the training set, and the brown triangles represent the test sample. The dashed line indicates the one-to-one relation.

4. Results

4.1. Feature importance from random forest

The model obtained by the RF algorithm allows us to investigate the connection between different galaxy properties (the features) and their overall gas metallicity (the target). Although the model was fed with more than 100 input parameters that characterise the galaxies, only a few of them are typically important for an accurate prediction for the target. The remaining parameters have little to no impact on Z_g.

In order to find the most relevant features, the RF provides a measurement of the relative importance of each feature in predicting the solution. In rough outlines, the so-called feature importance is a measure of how effective the feature is at reducing variance when the variables are split along the decision trees. The higher the value, the more important the feature. Figure 2 shows the relative importance of the first ten input features (see Table 1 for a description of the labels), ranked by decreasing order of the importance value. In order to test the stability of these values, we ran the RF algorithm 50 times, each time with a different randomly selected training sample. The figure shows the average trend of the values for the 50 realisations together with the standard deviations (shaded area). To our initial surprise, the baryonic gravitational potential (Φ_baryon) was the primary factor determining the global gas metallicity, with a mean importance value of ∼0.51, which is well above that of the other values. The second most relevant parameter was the mass-weighted stellar metallicity measured at 1 R_e (MW–[Z/H]_Re), which presents a relative importance of ∼0.08. The remaining properties had an importance value below 0.05 and are therefore of little or no relevance for determining the gas metallicity. This includes the stellar mass, which occupies the third position in the ranking, with a relative importance of ∼0.04.

Fig. 2. Average relative importance of the first ten input features in predicting the global gas metallicity for 50 runs of the RF algorithm. The features are ranked by decreasing order of importance. The shaded area represents the standard deviation of the mean relative importance among the 50 runs. The meaning of the labels is detailed in Table 1.

The importance value of Φ_baryon is much higher than that of the stellar mass, which means that the gravitational potential is significantly more efficient than the mass to decrease the variance when the variables are split along the decision trees. To confirm this, we ran the RF with the gravitational potential as the only parameter in the model. The root mean square error (RMSE) of the difference between the predictions (modelled Z_g) and the targets (measured Z_g) is 0.052 ± 0.001 (50 runs of the algorithm). If instead we use the stellar mass as the only parameter in the model, the RMSE increases up to 0.054 ± 0.001. This test supports the finding that Φ_baryon better describes the gas metallicity than M_⋆. We would also like to note that although the relative importance of the stellar mass is lower than MW–[Z/H]_Re, we must not forget the correlation between Φ_baryon and M_⋆, which might cause the reduced importance of M_⋆ in the model. We also ran the RF with the stellar metallicity as the only parameter in the model. We obtained an RMSE of 0.056 ± 0.001, confirming that although M_⋆ is the third most relevant parameter in the original model, it regulates the gas metallicity better than MW–[Z/H]_Re.

In order to visualise the shape of the dependence of the global gas metallicity on the stellar mass and on the baryonic gravitational potential, we present in Fig. 3Z_g as a function of both Φ_baryon (hereafter called ΦZR, top) and M_⋆ (the so-called MZR, bottom). The shape of the relations is very similar in both cases: Z_g increases when M_⋆ and Φ_baryon increase, although the flattening at the high end is more pronounced in the MZR. Regarding the dispersion, we find ΦZR to be slightly tighter than the MZR, although the difference is small (2%). In any case, these results agree with the conclusions derived from the RF analysis, supporting the finding that Φ_baryon is a slightly better predictor of Z_g than M_⋆.

Fig. 3. Global gas metallicity as a function of the stellar mass (the so-called MZR, bottom) and the galactic gravitational potential (hereafter called ΦZR, top). The brown squares represent the mean Zg values in ten mass and Φ bins, respectively. The averaged binned values were fitted with a spline function (solid red lines).

4.2. Gravitational potential, or the combined effect of stellar mass and galaxy size

The ratio M_⋆/R_e was used in the analysis as a proxy for the galactic gravitational potential at one effective radius. This first-order approximation ignores the effects of the dark matter (DM) halo on the potential in these inner regions, which is assumed to be dominated by the baryonic component (which in turn is dominated by the stellar component). We address the role of DM below.

The correlation between gas metallicity and both stellar mass and galaxy size has been described in numerous studies (e.g., Ellison et al. 2008; Sánchez Almeida & Dalla Vecchia 2018; Sánchez et al. 2017, 2019). Since Φ_baryon is estimated as M_⋆/R_e, the fact that Φ_baryon resulted as the primary factor determining Z_g could simply reflect therefore the combined effect of both stellar mass and galaxy size on Z_g. In this case, the combination of stellar mass and radius that best fits the metallicity could be in the form of $M_{\star }/R_{\mathrm{e}}^{\alpha }$, with the exponent α similar to but not necessarily equal to one. In order to further test this hypothesis, we ran the RF again, and this time included 20 new input parameters in the shape of $M_{\star }/R_{\mathrm{e}}^{\alpha }$, with α = 0.1, …, 2, in steps of 0.1 (α = 1, corresponding to Φ_baryon, and α = 2 to Σ_⋆). If another value of α ≠ 1 fits better than Φ_baryon, it could suggest that the real predictor of the Z_g is the combined effect of stellar mass and galaxy size and not the baryon gravitational potential of the galaxy alone.

The outcome of this test is shown in Table 2, that is, the average importance values (out of 50 realisations) of the new $M_{\star }/R_{\mathrm{e}}^{\alpha }$ features included in the RF, together with their position in the ranking list. Only the first five parameters are displayed. Interestingly, Φ_baryon (α = 1) is replaced by $M_{\star }/R_{\mathrm{e}}^{0.6}$ as the parameter with the strongest effect on Z_g.

This result seems to suggest that the gravitational potential set by the baryons is indeed not the best predictor of Z_g, but a particular combination of mass and size, $M_{\star }/R_{\mathrm{e}}^{0.6}$. However, as we pointed out above, in our approximation, we ignored the effects of the DM halo. α ≃ 0.6 might account for the inclusion of the DM component in the whole gravitational potential of the galaxy. In this way, $M_{\star }/R_{\mathrm{e}}^{0.6}$ would be a better tracer of the real Φ than M_⋆/R_e, and our main conclusion would remain: The total gravitational potential would be the galactic property with the strongest impact on Z_g. This seems indeed to be the case, as we discuss below.

In order to investigate the conjecture, we worked out the expression for the fraction of the DM,

$$\begin{array}{c}\hfill f_{\mathrm{DM}}=M_{\mathrm{DM}}/(M_{\mathrm{baryon}}+M_{\mathrm{DM}}),\end{array}$$(1)

that is needed for the logarithm of the total potential to scale with the parameter used in the RF, that is, the logarithm of $M_{\star }/R_{\mathrm{e}}^{\alpha }$,

$$\begin{array}{c}\hfill log\mathrm{\Phi }=log(M_{\mathrm{T}}/R_{\mathrm{e}})=a+blog(M_{\star }/R_{\mathrm{e}}^{\alpha }),\end{array}$$(2)

where a and b are two scaling constants, and M_T is the total mass, composed of baryons, M_baryon, and DM, M_DM. Assuming a functional dependence of f_DM of the type

$$\begin{array}{c}\hfill f_{\mathrm{DM}}=1-B{(M_{\star }/2\pi R_{\mathrm{e}}^2)}^{\beta },\end{array}$$(3)

with B an arbitrary factor, we find after some algebra that the functional form in Eq. (3) is consistent with Eq. (2) provided

$$\begin{array}{c}\hfill \beta =(1-\alpha )/(2-\alpha ).\end{array}$$(4)

We considered M_baryon ≃ M_⋆, which is a good approximation for most MaNGA galaxies (Lin et al. 2020). The functional form in Eq. (3) was chosen because it seems to be followed by the galaxies produced in cosmological numerical simulations, as we argue below.

Figure 4 uses Eqs. (3) and (4) to show the baryon fraction for Eq. (2) to hold, depending on the value of α (the coloured lines). The figure also includes the theoretical expectation for f_DM (the symbols) as predicted by the Illustris-TNG100 simulations (Lovell et al. 2018) at redshift 0, a prediction that we extracted from Fig. 5 in Nestor Shachar et al. (2023). The predictions from two additional simulations by Oser et al. (2010) and Remus et al. (2017) are also shown, reflecting that the slope of the relation between f_DM and Σ_⋆ given by β might slightly change depending on the set of simulations that is used. Figure 4 illustrates that α = 0.6 (the solid orange line) matches the theoretical expectation from Illustris-TNG100 (grey symbols and shaded area) extremely well, except for extreme surface densities. The match is encouraging because no freedom or fitting is shown in the representation in the figure (except for a shift in the y-axis given by the constant B). In any case, the match for α = 0.6 is much better than 0.5 or 0.7. We interpret this match as evidence that $M_{\star }/R_{\mathrm{e}}^{0.6}$ traces Φ, and therefore, evidence that our RF reveals a true relation between the total gravitational potential of the galaxy and the gas metallicity, which is tighter and more primordial than the MZR or the baryonic ΦZR.

Fig. 4. Dependence of the dark matter fraction on the baryonic surface density for galaxies in the local Universe predicted by Illustris-TNG100 simulations (Nestor Shachar et al. 2023, Fig. 5). This relation is compared with the one obtained from the five $M_{\star }/R_{\mathrm{e}}^{\alpha }$ parameters with the highest relative importance in the RF. $M_{\star }/R_{\mathrm{e}}^{0.6}$ provides a good fit for the theoretical relation. Two additional relations for the dependence of fDM on Σ⋆ from different simulations are also included (Oser et al. 2010; Remus et al. 2017).

5. Discussion

The search for relations between integrated physical properties in galaxies has been widely addressed in the literature. In particular, scaling relations involving the gas-phase oxygen abundance (as a proxy for the gas metalllicity, Z_g) have received much attention because it is a tracer of the overall chemical evolution of galaxies, modulated by the inflow or outflow of gas in different regions (Sánchez Almeida et al. 2014; Maiolino & Mannucci 2019; Sánchez 2020). For this reason, the dependence of Z_g on other properties (or vice versa) can provide important insights into galaxy formation and evolution.

In this study, we analysed the global properties of a sample of more than 3000 nearby galaxies from the SDSS-IV MaNGA survey (Bundy et al. 2015). The aim of this analysis is to evaluate the galactic property that is the best predictor of Z_g, and to do this, we explored the effect of 148 parameters extracted from the pyPipe3D Value Added catalogue (Sánchez et al. 2022). We include the parameters that were previously reported to strongly correlate with Z_g, such as M_⋆, Φ, Σ_⋆, and M_gas. We used the oxygen abundance measured at 1 R_e as a proxy for the global gas metallicity, determined with the O3N2 calibration proposed by Marino et al. (2013).

For this analysis, we applied a RF regressor, which is a supervised machine-learning technique that has been proven successful in the past in searches for the primary relations between galaxy properties (e.g., Sánchez-Menguiano et al. 2019; Baker et al. 2023a,b). Our trained model was able to predict Z_g with an average error (defined as the difference between the predictions and the targets) of about the propagated uncertainties associated with the measured values.

Based on the relative importance of each galaxy property in predicting Z_g, we obtained that the baryonic gravitational potential (defined as Φ_baryon = M_⋆/R_e) is the most relevant parameter of the model, followed by MW–[Z/H]_Re and M_⋆, with an importance value more than five times lower than that of Φ_baryon. This result was supported by the outcome of a subsequent run of the RF with only one parameter in the model. When that parameter was M_⋆ (MW–[Z/H]_Re), the RMSE of the differences between the predictions and the targets was 0.054 (0.056)±0.001, whereas this value was reduced to 0.052 ± 0.001 for Φ_baryon. The difference between these values and that found for the complete model (0.029) suggests that, assuming there is no overfitting, other secondary parameters play a significant role in shaping the gas metallicity. This will be investigated in a future study (Sánchez-Menguiano et al., in prep.). Furthermore, the scatter of the relations between Z_g and Φ_baryon (ΦZR) and between Z_g and M_⋆ (the so-called MZR) shows that the ΦZR is slightly tighter than the MZR (2%). Although the differences are small in both cases, these tests support the finding that Φ_baryon is better suited than M_⋆ to describe the gas metallicity.

These results are representative of the whole population of analysed galaxies, which spreads over a wide range of stellar masses (8.5 < log(M_⋆/M_⊙) < 11.5). Although our mass coverage is not homogeneous (there is an excess of massive galaxies compared to low-mass systems), we show in Appendix C that this bias does not affect the analysis because using a subsample of galaxies with a uniform mass distribution (constant number of objects per mass bin) leads to the same results. Restricting the analysis to a more limited mass range, however, might lead to different conclusions. For instance, because both the MZR and the ΦZR saturate at high masses, a sample that is highly dominated by massive systems and with a poor coverage of low-mass galaxies could reduce the effect of both M_⋆ and Φ_baryon on Z_g (because Z_g remains almost constant, while both parameters increase). It is therefore important to rely on a wide mass range that allows covering a wide dynamical range for Z_g in order to detect true primary and secondary relations.

It is interesting to note that the stellar mass estimated from photometry (within the MaNGA FoV) based on the use of the M/L ratio and the relations provided in Bell & de Jong (2000) is a better predictor of Z_g than the mass derived from the spectroscopic analysis of the stellar populations performed by pyPipe3D. Whereas the first obtains an importance value of 0.043, the later obtains just 0.001. Another estimation of the integrated stellar mass from the NSA catalog provides an intermediate value (0.010). Photometric stellar masses are considered to be less accurate, but are probably more precise than estimates based on the spectroscopic analysis performed by pyPipe3D. For a more detailed comparison of the mass estimations, we refer to Sánchez et al. (2022).

Due to the well-known discrepancies between the oxygen abundances derived using different strong-line calibrators (for a review on the topic, see Kewley et al. 2019), we reproduced the analysis employing 18 alternative O/H estimators. Consistent with the use of O3N2-M13, in 15 out of the 18 cases, we obtained that Φ_baryon is the best predictor of the global gas metallicity. This strengthens the validity of the result, which is independent of the method used to estimate the abundances. In the remaining three, the discrepancies might be due to the small dynamical range of the calibrators (more details can be found in Appendix B).

This is not the first study to find that the correlation of the gas metallicity with the baryonic gravitational potential is stronger than with the stellar mass. Using a volume-limited sample drawn from the publicly available SDSS DR7, D’Eugenio et al. (2018) analysed the relation between Z_g and three galaxy properties: Φ_baryon, M_⋆, and average Σ_⋆. To minimise aperture bias, the authors also defined three aperture-matched subsamples to ensure an homogeneous coverage of the physical size of galaxies from the fixed-aperture SDSS fibres: R_fib/R_e = 0.5, 1.0, and 1.5. For the three subsamples, they find that ΦZR had the lowest scatter and the highest Spearman correlation coefficient. In addition, this relation had the lowest residual trends with galaxy size. Based on these results, the authors concluded that of the three parameters, Φ_baryon is the best predictor of Z_g. According to D’Eugenio et al. (2018), there are two possible explanations for this. The correlation of Z_g and Φ might arise from the combination of metal loss (which depends on the value of the escape velocity, which in turn is related to M_⋆) and metal enrichment (driven by f_gas, and therefore Σ_⋆). Each of these effects induces a correlation between Z_g and $M_{\star }/R_{\mathrm{e}}^n$, with n = 0 and n = 2, respectively, and their combination could produce an overall correlation with an intermediate value of n. However, because the metallicity of low-mass galaxies is highly affected by outflows, Z_g is predicted to correlate best with M_⋆ at the low-mass end of the relation. Similarly, because the escape fraction for high-mass galaxies is small, Z_g should correlate best with Σ_⋆ at the high-mass end of the relation. These predicted trends are not observed in their analysis, and therefore, the authors pointed to another explanation for the ΦMR: a direct physical link between the average depth of the gravitational potential and the average metallicity, which is a consequence of an existing local relation between Z_g and the escape velocity (this local quantity is directly related with the global gravitational potential). We note that this local relation does exist and has already been reported in Barrera-Ballesteros et al. 2018. The first explanation would imply that the ΦMR is the result of the already known MZR and ΣZR, whereas the latter indicates that the potential–metallicity relation is the primordial relation, and the MZR and the ΣZR would arise from the local relation between Z_g and Φ.

The first explanation provided by D’Eugenio et al. (2018) could be further tested by re-applying the RF and include different parameters in the model in the shape of $M_{\star }/R_{\mathrm{e}}^{\alpha }$. We tried 20 combinations with α = 0.1, …, 2, in steps of 0.1. Interestingly, Φ_baryon (α = 1) is replaced by $M_{\star }/R_{\mathrm{e}}^{0.6}$ as the parameter with the strongest effect on Z_g. When the alternative calibrations were applied to determine Z_g, the coefficients α = 0.7 and α = 0.8 in some cases replaced α = 0.6 as those occupying the highest position in the importance ranking. This may suggest that the ΦMR is indeed not the primordial relation, but results from the MZR and ΣZR with a particular combination of the effects of metal loss and metal enrichment.

However, we provide here another explanation to justify the result. We have ignored the effects of the DM halo in the analysis so far. α ≃ 0.6 might account for the inclusion of the DM component in the whole gravitational potential of the galaxy. In Fig. 4 we showed that a scale of the form $M_{\star }/R_{\mathrm{e}}^{0.6}$ for the total gravitational potential matches the theoretical relation between the DM fraction and the baryonic surface density predicted in the Illustris TNG cosmological numerical simulation very well (Nestor Shachar et al. 2023). In this way, $M_{\star }/R_{\mathrm{e}}^{0.6}$ would be a better tracer of the real Φ than M_⋆/R_e, and therefore, the total gravitational potential (and not its baryonic component) would be the galactic property that predicts Z_g best. The substitution of α = 0.6 by α = 0.7 or α = 0.8 for the parameter with the strongest effect on Z_g when using other abundance calibrations can be accounted for by the large scatter in the prediction by Nestor Shachar et al. (2023; shaded area in Fig. 4), and the fact that the slope of the relation might vary slightly depending on the simulations used (Oser et al. 2010; Remus et al. 2017).

A recent study by Baker & Maiolino (2023) also explored the relation between Z_g and (i) M_⋆, (ii) dynamical mass (M_dyn, DMZR), and (iii) gravitational potential (estimated as Φ = M_dyn/R_e) using MaNGA data. Based on three different methods (average dispersion, partial correlation coefficients, and RF), the authors reported that the gas metallicity primarily depends on the stellar mass, with little or no dependence on either dynamical mass or gravitational potential when the dependence on stellar mass is taken into account. However, their M_dyn were derived via Jeans anisotropic modelling, which is a complex method with its own systematic problems that might increase the intrinsic scatter in DMZR and ΦZR with respect to MZR (e.g., Li et al. 2016; El-Badry et al. 2017; Read et al. 2021). We argue that our derivation of Φ_baryon, although less accurate, produces more precise results that allow us to determine the intrinsic scatter in ΦZR, resulting in a tighter relation than the MZR. This is in contrast to the results reported by Baker & Maiolino (2023).

6. Conclusions

We studied the dependence of the gas-phase metallicity (Z_g) on a large set of galaxy properties, including parameters that have been reported in the literature to correlate strongly with Z_g, such as the stellar mass (M_⋆), the baryonic gravitational potential (Φ_baryon), the average stellar mass surface density (Σ_⋆), or the gas mass (M_gas). In order to do this, we applied a RF regressor algorithm on a sample of more than 3000 nearby galaxies from the SDSS-IV MaNGA survey. Our results are listed below.

(i) M_⋆/R_e, as a proxy for Φ_baryon, is found to be the primary factor determining Z_g (Sect. 4.1; Fig. 2). Because the gravitational potential regulates the ability of a galaxy to retain metals, it is not surprising that ΦMR is the primordial scaling relation involving Z_g.

(ii) When we included a parameter of the type $M_{\star }/R_{\mathrm{e}}^{\alpha }$ in the analysis, we find that α = 0.6 best correlates with Z_g (Sect. 4.2; Table 2). We explored the effect of the DM halo on the gravitational potential, and we showed that a scale of the form $M_{\star }/R_{\mathrm{e}}^{0.6}$ for the total Φ matches the theoretical relation between the DM fraction and the baryonic surface density predicted in cosmological numerical simulations of galaxy formation very well (Sect. 4.2; Fig. 4). This suggests that the total gravitational potential would be the galactic property that predicts Z_g best.

(iii) In the analysis, we used the oxygen abundance measured with the O3N2 calibration by Marino et al. (2013) as an indicator of the gas metallicity. The use of alternative estimators still gives ΦMR as the primordial relation (Appendix B). The inclusion of the DM halo might affect the particular coefficient α, which could be better described by theoretical predictions from other simulations (Sect. 5; Fig. 4).

(iv) We find evidence that Φ_baryon alone cannot predict Z_g and that other secondary parameters might play a significant role in shaping the gas metallicity (Sect. 5). This will be investigated in a future study (Sánchez-Menguiano et al., in prep.).

Acknowledgments

We thank the anonymous referee for suggestions that allowed us to improve the paper. L.S.M. acknowledges support from Juan de la Cierva fellowship (IJC2019-041527-I) and from project PID2020-114414GB-100, financed by MCIN/AEI/10.13039/501100011033. This research was partly funded by the Spanish Ministry of Science and Innovation, projects PID2019-107408GB-C43 and PID2022-136598NB-C31(ESTALLIDOS), and by Gobierno de Canarias through EU FEDER funding, project PID2020010050. S.F.S. thanks the PAPIIT-DGAPA AG100622 project. This project makes use of the MaNGA-Pipe3D dataproducts. We thank the IA-UNAM MaNGA team for creating this catalogue, and the Conacyt Project CB-285080 for supporting them.

References

Appendix A: Galaxy parameters

In this appendix, we list all the galactic parameters extracted from the pyPipe3D VAC that are included in the RF analysis. This includes a brief description of their meaning. The information is displayed in Table A.1.

Appendix B: Effect of the adopted metallicity calibration

The derivation of the gas metallicity involves a large number of systematics and sources of uncertainties, with well-known discrepancies between the use of different strong-line calibrators (e.g. Kewley & Ellison 2008; López-Sánchez et al. 2012; Kewley et al. 2019). In this regard, we replicated the analysis employing 18 alternative calibrators for which the derivation of the oxygen abundance is provided in the pyPipe3D VAC table to assess whether the choice of the gas-metallicity indicator affects the results of this study. In Table B.1 we summarise the adopted calibrations and the corresponding references. We note that although some of the calibrators are multi-valued, they present a monotonic behaviour in the covered range of abundances of this study (8.2 < 12 + log(O/H) < 8.6), and they can therefore be used to unambiguously determine the gas metallicity for all galaxies in the sample.

In Figure B.1 we show the relative importance in the RF of the 15 features with the highest importance values as ranked by the use of the main gas metallicity estimator (O3N2-M13, see Sec. 4.1). Each colour corresponds to the choice of an alternative calibration (see Table 2 for a description of the labels). Consistent with the use of the O3N2-M13 indicator, in 15 out of the 18 cases we obtain that the baryonic gravitational potential is the best predictor of the global gas metallicity. The only exceptions are when the O3O2, the R2, and the O3S2 indicators proposed in Curti et al. (2017) are used, for which the oxygen abundance gradient (two cases) or the stellar metallicity (one case) are found to be more relevant. The small dynamical range of these estimators might cause these differences. Figure 29 of Sánchez et al. (2022) shows that while the other calibrators are clearly correlated with the abundances derived from H19 (arbitrarily selected as the fiducial calibrator), O3O2, R2, and O3S2 show a flatter distribution, which could bias the results and therefore dismiss their validity as abundance estimators for this work.

Fig. B.1. Importance of the 15 features of the RF with the highest importance values as shown in Figure 2. The predicted Zg is estimated using different calibrations (coloured lines). See the caption of Figure 2 for more details.

For the 15 abundance indicators for which Φ_baryon was found to be the most relevant parameter in the RF, we again ran the algorithm now including in the model 20 different combinations of stellar mass and galaxy size in the shape of $M_{\star }/R_{\mathrm{e}}^{\alpha }$, with α = 0.1, ..., 2 in steps of 0.1 (see Sec. 4.2 for more details). We show in Figure B.2 the resulting average importances (out of 50 realisations; the error bars represent the standard deviations) of the 6 α values that present the highest values overall. In general, the coefficients α = 0.7 and α = 0.8 occupy the highest position in the ranking. The coefficient α = 0.6, which was the parameter with the strongest effect on Z_g when estimated with O3N2-M13, presents the highest importance value for 3 out of the 15 calibrations, but it is the third parameter in the ranking overall when the results of all abundance indicators are averaged. In Sec. 4.2 we showed that a scale of the form $M_{\star }/R_{\mathrm{e}}^{0.6}$ for the total gravitational potential matches the theoretical relation between the DM fraction and the baryonic surface density predicted in Nestor Shachar et al. (2023) very well. However, as we argued in Sec. 5, the scatter in the relation, as well as the fact that other simulations provide different slopes for the dependence of f_DM with Σ_⋆, mean that the interpretation is still reasonable that Φ, traced by the parameter $M_{\star }/R_{\mathrm{e}}^{\alpha }$, is the best predictor of Z_g. We refer to that section for further discussions of the matter.

Fig. B.2. Importance of the six combinations of $M_{\star }/R_{\mathrm{e}}^{\alpha }$ with the overall highest importance values in the RF for the alternative abundance calibrations (x-axis labels; see Table B.1 for references). We represent the average values (small circles) for 50 runs of the algorithm together with their standard deviations (error bars).

Appendix C: Effect of a non-homogeneous mass coverage of the sample

The stellar mass distribution of the galaxies in the sample is not uniform, but there is an excess of massive systems compared to low-mass objects. This bias is inherited from the MaNGA mother sample. In order to demonstrate that the dominance of massive galaxies in the analysis does not bias the results, we selected a subsample of 1600 galaxies presenting a uniform distribution of stellar masses (i.e. a constant number of objects per mass bin). To do this, we randomly selected 200 galaxies in eight mass ranges of 0.25 log M_⊙ width from 9.0 to 11.0 log M_⊙. In Figure C.1 we represent the stellar mass distribution for the three samples: the MaNGA mother sample (in grey), the analysed sample (in green), and the subsample of uniform mass distribution (in blue). Outside the range 9-11 log M_⊙, we do not have a sufficient number of galaxies to include in the test, but given the low statistics, it is unlikely that the results are driven by these objects, which are barely 6% of the sample.

Fig. C.1. Stellar mass distribution of the galaxies analysed in this work (green) in comparison with the distribution of the subsample selected for this test (blue), characterised by a uniform distribution for the stellar mass. The numbers represented at the top indicate the number of galaxies from the original sample comprising each mass bin of 0.25 log(M⊙). The MaNGA mother sample is also shown in grey.

We ran the RF using the uniform subsample, and we obtained very similar results as for the original sample. Figure C.2 shows the relative importance of the first ten input features (see Table A.1 for a description of the labels), ranked by decreasing order of the importance value. The solid line shows the average trend of the values for 50 realisations, each of which with a different randomly selected training sample, and the shaded area represents the standard deviation of the trend. As in the original analysis, we find Φ_baryon to be the primary factor determining the global gas metallicity, with a mean importance value of ∼0.51. The next parameters in the ranking are M_⋆ and MW-[Z/H]_Re, with relative importances of 0.07 and 0.05, respectively. This test strengthens the main conclusion reached in this study and proves that this is not affected by any bias induced by a non-homogeneous mass coverage of the sample. All our tests show that the gravitational potential is the most important galaxy property setting the gas metallicity via the ΦZR.

Fig. C.2. Average relative importance of the first ten input features in predicting the global gas metallicity for 50 runs of the RF algorithm. The features are ranked by decreasing order of importance. The shaded area represents the standard deviation of the mean relative importance of the 50 runs. The meaning of the labels is detailed in Table A.1. For this test, the uniform sample was used (see Figure C.1).

All Tables

All Figures

	Fig. 1. Global gas-phase metallicity – described by the oxygen abundance 12 + log(O/H) – predicted by the RF algorithm vs. the measured values. The blue circles correspond to the training set, and the brown triangles represent the test sample. The dashed line indicates the one-to-one relation.
In the text

	Fig. 2. Average relative importance of the first ten input features in predicting the global gas metallicity for 50 runs of the RF algorithm. The features are ranked by decreasing order of importance. The shaded area represents the standard deviation of the mean relative importance among the 50 runs. The meaning of the labels is detailed in Table 1.
In the text

	Fig. 3. Global gas metallicity as a function of the stellar mass (the so-called MZR, bottom) and the galactic gravitational potential (hereafter called ΦZR, top). The brown squares represent the mean Zg values in ten mass and Φ bins, respectively. The averaged binned values were fitted with a spline function (solid red lines).
In the text

Fig. 4. Dependence of the dark matter fraction on the baryonic surface density for galaxies in the local Universe predicted by Illustris-TNG100 simulations (Nestor Shachar et al. 2023, Fig. 5). This relation is compared with the one obtained from the five $M_{\star }/R_{\mathrm{e}}^{\alpha }$ parameters with the highest relative importance in the RF. $M_{\star }/R_{\mathrm{e}}^{0.6}$ provides a good fit for the theoretical relation. Two additional relations for the dependence of fDM on Σ⋆ from different simulations are also included (Oser et al. 2010; Remus et al. 2017).

In the text

	Fig. B.1. Importance of the 15 features of the RF with the highest importance values as shown in Figure 2. The predicted Zg is estimated using different calibrations (coloured lines). See the caption of Figure 2 for more details.
In the text

	Fig. B.2. Importance of the six combinations of $M_{\star }/R_{\mathrm{e}}^{\alpha }$ with the overall highest importance values in the RF for the alternative abundance calibrations (x-axis labels; see Table B.1 for references). We represent the average values (small circles) for 50 runs of the algorithm together with their standard deviations (error bars).
In the text

Fig. C.1. Stellar mass distribution of the galaxies analysed in this work (green) in comparison with the distribution of the subsample selected for this test (blue), characterised by a uniform distribution for the stellar mass. The numbers represented at the top indicate the number of galaxies from the original sample comprising each mass bin of 0.25 log(M⊙). The MaNGA mother sample is also shown in grey.

In the text

Fig. C.2. Average relative importance of the first ten input features in predicting the global gas metallicity for 50 runs of the RF algorithm. The features are ranked by decreasing order of importance. The shaded area represents the standard deviation of the mean relative importance of the 50 runs. The meaning of the labels is detailed in Table A.1. For this test, the uniform sample was used (see Figure C.1).

In the text