More fundamental than the fundamental metallicity relation: The effect of the stellar metallicity on the gas-phase mass-metallicity and gravitational potential-metallicity relations

,


Introduction
Scaling relations constitute a fundamental tool to improve our understanding on the formation and evolution of galaxies.A primary scaling relation links the gas-phase metallicity (Z g ) and the stellar mass (MZR, Tremonti et al. 2004;Kewley & Ellison 2008;Wu et al. 2016;Barrera-Ballesteros et al. 2017;Sánchez et al. 2019;Yates et al. 2020).Its shape is characterised by a positive correlation between both parameters that flattens up at high masses.Confirmed up to redshift z ∼ 10 thanks to the advent of JWST (presenting a similar shape but with a downward offset to lower Z g in high-z galaxies, e.g.Curti et al. 2023b;Nakajima et al. 2023), the MZR has been proposed to result from the interplay of different processes including metal removal by outflows, dilution by metal-poor gas infall, or enrichment by previous star formation.A simple evolution by the so-called 'downsizing', where most massive galaxies exhaust their gas reservoir faster, would also explain its origin (Maiolino & Mannucci 2019).
Notable attention has been paid to the study of the scatter in the MZR, whose physical origin can unveil the relative importance of key processes in galaxy evolution.Thus, several secondary dependences in the MZR have been reported, such as that with the star formation rate (SFR) or specific SFR, gas fraction or gas mass, galaxy size, and stellar age (e.g.Ellison et al. 2008;Lara-López et al. 2010;Lian et al. 2015;Bothwell et al. 2016a;Sánchez Almeida & Dalla Vecchia 2018;Sánchez Almeida & Sánchez-Menguiano 2019;Curti et al. 2020;Alvarez-Hurtado et al. 2022;Scholte & Saintonge 2023), to name a few.Particularly relevant is the relation between the gas metallicity, stellar mass (M ), and SFR, called the fundamental metallicity relation (FMR; Mannucci et al. 2010), which shows how for low-and intermediate-mass galaxies (log[M /M ] < 10.5) there exists an anti-correlation between the SFR and Z g at a given M .The FMR is ascribed to the accretion of metal-poor cosmic gas fueling star formation (Mannucci et al. 2010;Davé et al. 2012;Sánchez Almeida et al. 2014).However, there is a current debate on whether the SFR is truly needed to describe the relation between M and Z g (e.g.Izotov et al. 2014;de los Reyes et al. 2015;Barrera-Ballesteros et al. 2017;Sánchez et al. 2017;Duarte Puertas et al. 2022).
Recent studies have also revealed a stronger dependence of Z g with the baryonic gravitational potential (Φ baryon ) than with M (ΦZR, D 'Eugenio et al. 2018;Sánchez-Menguiano et al. 2024).In particular, in Sánchez-Menguiano et al. (2024, hereafter Paper I), we applied a random forest (RF) regressor on a sample of ∼3500 nearby galaxies from the Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey and included, in the model, around 140 galaxy properties.Among the analysed parameters, we can find those that have been previously reported to strongly correlate with Z g , such as M , Φ baryon , Σ and, M gas .We reported the ΦZR to be the tightest relation involving Z g , and so the one with higher chances of being the primary one.
In Paper I, we also show that a model including only Φ baryon or M as the input parameters (ΦZR and MZR, respectively) presented a root mean square error (RMSE) of the difference between the predictions (modelled Z g ) and the targets (measured Z g ) much higher than the complete model with all the input parameters.This result provides evidence that other parameters besides the gravitational potential (and stellar mass) play a significant role in shaping the global gas metallicity.Whether this parameter is the SFR or any other has yet to be found.Therefore, in order to investigate which is the most significant secondary dependence in the ΦZR and MZR, we analyse in this study the residuals of the Z g once the dependence with Φ baryon (and, alternatively, M ) has been modelled and subtracted out.Based on the use of the RF regressor and an extensive set of galaxy properties, we determine which of them present the tightest correlation with these metallicity residuals.As argued, this allows us to unveil the relative importance of key processes in galaxy evolution.
This Letter is organised as follows.Section 2 briefly presents the MaNGA data, the subset of galaxies comprising our study sample, and the physical parameters included in the model.A brief overview of the RF algorithm is given in Sect.3. We describe the outcome of the RF in Sect. 4. Finally, the results are discussed in Sect. 5 and Sect.6 compiles and summarises the main conclusions of the work.Supplementary material includes Appendices A and B. Appendix A reproduces the analysis using alternative calibrations to estimate Z g .Appendix B describes the analysis for the ΦZR based on an alternative tracer of the total gravitational potential (including the effect of the dark matter halo).

Sample and analysed data
This study is based on the data collected by the MaNGA project (Bundy et al. 2015), which is part of the fourth generation Sloan Digital Sky Survey (SDSS-IV).A very brief overview of the data is given in Paper I, and further details on the MaNGA mother sample, survey design, observational strategy, and data reduction are provided in the literature (Yan et al. 2016;Wake et al. 2017;Law et al. 2015Law et al. , 2016Law et al. , 2021)).
In order to investigate the existence of any secondary dependences in the ΦZR and MZR, we explored a list of 148 galaxy parameters extracted from the pyPipe3D Value Added Catalogue (VAC, Sánchez et al. 2022), which is publicly accesible through the SDSS-IV VAC website1 .We refer the reader to Paper I for further details on the selected physical properties, which are listed in Appendix A of that article.This list includes parameters such as M , Φ baryon , Σ , M gas , morphological type, colours, SFR, stellar age and metallicity, or dust extinction.We used the oxygen abundance (O/H) measured at one disc effective radius (R e ) as a proxy for Z g .This estimation is based on the O3N2 index and the empirical calibration proposed by Marino et al. (2013).The results described in this Letter (Sect.4) are confirmed applying alternative calibrators (see Appendix A).
The analysed galaxy sample has been drawn from the 10 010 galaxies comprising the MaNGA mother sample.As described in Paper I, for selecting the galaxies, we have adopted two simple criteria: galaxies have to meet the required quality standards of the analysis pipeline (i.e.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 contain values for all of the analysed parameters.The second condition is quite restrictive.It excludes early-type systems and galaxies in general with minimal or negligible levels of ionised gas, making it impossible to derive gas-related attributes such as oxygen abundance or dust extinction.The resultant sample consists of 3430 galaxies, a sufficiently large number to guarantee the statistical robustness of the results.However, we note that the sample is biased towards large star-forming galaxies (mostly Sa to Sc galaxies).

The RF regressor
The use of RF regressors (Breiman 2001) in the study of scaling relations has grown over the last years, offering a very compelling tool to unveil correlations between galaxy properties and identify the most significant dependencies (e.g.Sánchez-Menguiano et al. 2019;Bluck et al. 2020;Moster et al. 2021;Piotrowska et al. 2022;Baker et al. 2023, among others).For that, this algorithm employs an ensemble of decision trees to identify the input features (galaxy properties in this context) that carry the most complete information on the target feature (the galaxy parameter of interest, the residual gas metallicity in our case).Subsequently, it constructs a predictive model by establishing a set of conditions based on the values of the input features, gathering information on the relative importance of each property in the model in the process.
For this study we used a RF regressor to investigate the existence of any secondary dependences in the ΦZR and MZR.For that, we analysed the residuals of Z g once the dependence on Φ baryon (and, alternatively, M ) was modelled and subtracted.The RF algorithm was implemented in the scikit-learn package for Python (Pedregosa et al. 2011).While in Paper I we provide a brief overview of the basic steps involved and of the selected values for the model parameters, we refer the reader to the scikit-learn User Guide documentation for comprehensive details on the complete algorithmic implementation2 .

Results
The relation between Z g and Φ baryon , the so-called ΦZR studied in Paper I, is represented in the left panel of Fig. 1.We can see that Z g increases when increasing Φ baryon , with a flattening at the high end.In order to study secondary dependences on such a relation, we need to model the primary dependence of Z g on Φ baryon .For that, the median values in ten bins of ∼0.1 dex width (salmon squares) were fitted with a spline function (solid salmon line).We see a reduction of 42% in the dispersion with respect to that of the original relation.The parameter quantifying the reduction (denoted as σ red ) was derived by subtracting the standard deviation of the residuals from the fit (difference between the measured Z g and the fit) to the standard deviation of the measured Z g values.This was then multiplied by 100 and divided by the standard deviation of the measured Z g values to obtain a percentage where Z g is the measured global gas metallicity and Z g,fit is the modelled metallicity.
We subsequently ran the RF using the residual Z g (∆Z g ) as the targets and removing Φ baryon from the input parameters.The algorithm generated a model with T95-[Z/H] Re as the most relevant feature, which corresponds to the stellar metallicity at 1 R e measured at time T95 (the look-back time at which the galaxy formed 95% of its mass, typically within the last 1 Gyr).T95-[Z/H] Re presents an importance value of 0.193 and is followed in the model by the Sersic index (importance value of 0.06).The remaining parameters all have importance values below 0.03.This analysis shows that the ΦZR seems to present a secondary dependence on the stellar metallicity.To quantify this dependence and its significance, we fitted the relation between ∆Z g and T95-[Z/H] Re using again a spline function, and derived the decrease in the dispersion of ∆Z g following Eq.( 1).We obtain a value of 14%.To explore further dependences of Z g , we repeated this exercise fitting the residuals from the previous relation and running a RF with the remaining parameters.Figure 2 (black line, bottom x-axis) represents the decrease in the dispersion of Z g as we modelled its dependence with different galaxy parameters revealed by the RF.We can see that after fitting the dependence on Φ baryon (42%) and the stellar metallicity (14%), the fit with the remaining galaxy properties yields a decrease in the dispersion that falls below 5% in all cases.This result suggests that the global gas metallicity can be unambiguously characterised by two parameters: the galactic gravitational potential and the stellar metallicity (in particular, T95-[Z/H] Re ).In the left panel of Fig. 1, the ΦZR is colour-coded by T95-[Z/H] Re .We clearly see how, at a given value of the gravitational potential, Z g increases with increasing T95-[Z/H] Re (the inset shows this correlation for the particular bin 8.8 ≤ Φ baryon ≤ 9.0), highlighting the existence of this 3D relation between these parameters.The rest of galaxy properties do not seem to have a significant effect on Z g .
Alternatively, we repeated this exercise based on the residuals from subtracting the dependence of Z g with M .This will L11, page 3 of 8 In order to compare the analysed relations with the FMR, the right panels of Fig. 1 and Fig. 3 show the dependence of both ΦZR and MZR on the SFR.At a given value of Φ baryon and M , the colour-coding and the insets reflect a decrease in gas metallicity when SFR increases.We thus checked how much of a reduction in the dispersion would result by fitting the SFR as the secondary parameter in both relations.Whereas the stellar metallicity produces a reduction of ∼11−14%, this amount falls below 3% when considering the SFR.Since the shape of the dependence of the MZR with SFR depends on the considered stellar-mass range (Mannucci et al. 2010;Sánchez Almeida & Sánchez-Menguiano 2019), we performed several tests by dividing the sample into different mass ranges and determining the reduction in the dispersion on these subsamples.Table 1 summarises the results.We can see that the dependence with the SFR is higher in the MZR, and quite significant for masses below ∼10 log M .However, the dependence with LW-[Z/H] Re is stronger in all mass ranges.In the case of the ΦZR, the effect of the SFR on the relation is negligible in all ranges of gravitational potential, whereas the dependence with T95-[Z/H] Re is clear.

Discussion
The scatter in the MZR has been observed to correlate with several galaxy properties.The anti-correlation between the gas metallicity and SFR at a fixed stellar mass was reported in various studies, which defended the existence of a thin surface in 3D space defined by M , Z g , and SFR where star-forming galaxies are confined (Mannucci et al. 2010;Lara-López et al. 2010;Curti et al. 2020).Its initially revealed lack of evolution with redshift4 earned it the name fundamental metallicity relation, making it a key ingredient for chemical and galaxy evolution studies.Some works, nonetheless, argue that the secondary dependence on the SFR does not truly improve the relationship between Z g and M , with the reduction in the scatter being insignificant (e.g.Hughes et al. 2013;Sánchez et al. 2019).Recent studies also suggest that, if it exists, the reported correlation between the MZR scatter and the SFR could be a by-product of a more primordial relation with the gas mass and/or fraction (Bothwell et al. 2013(Bothwell et al. , 2016a,b;,b;Chen et al. 2022;Scholte & Saintonge 2023), due to the well-known link between gas content and star formation activity.
Regarding the ΦZR, no previous works have explored how the scatter might depend on different galaxy properties.In this study, we investigate for the first time this matter, as well as try to provide a definite answer as to the reported correlations of the MZR residuals.Using a RF regressor algorithm, we find that the MZR and the ΦZR residuals are best predicted by two estimators of the global stellar metallicity: LW-[Z/H] Re (for the MZR) and T95-[Z/H] Re (for the ΦZR).These parameters provide a reduction in the scatter of the primary relations of around 11% and 14%, respectively.We note that the reduction in the scatter of the MZR provided by T95-[Z/H] Re , LW-[Z/H] Re , or even MW-[Z/H] Re are very similar, with differences below 0.5%.
Table 1.Decrease in the dispersion (in percent) of the ΦZR and MZR when considering a secondary relation with the SFR in comparison with an indicator of the stellar metallicity for different ranges of mass (for MZR) or gravitational potential (for ΦZR).See main text for details.All [8.2,8.4] [8.6,8.8] [9.0,9.2] [9.4,9.6] [9.8,10.0] [10.2,10.4All [8.7,8.9] [9.1,9.3] [9.5,9.7] [9.9,10.1] [10.3,10.5] [10.7,10.9These differences, although not significant, make the RF prioritise one of these parameters as opposed to the others when developing the model.In any case, it is clear that the stellar metallicity, however estimated, is the secondary parameter with the highest effect on MZR and ΦZR.In contrast, the relative importance of parameters such as SFR or M gas 5 in the model is very low, with the reduction in the scatter for the case of the SFR being barely 0.5−2%.When we restricted the analysis to specific stellar mass ranges, the reduction in the scatter was higher for the stellar metallicity in all cases.Nevertheless, we see an improvement in the reduction of the MZR scatter with the SFR up to 20% for log(M ) ∼ 9 M , in agreement with the original studies of the FMR that find a stronger anti-correlation between Z g and SFR for low-mass galaxies (Sánchez Almeida & Sánchez-Menguiano 2019).This is not the case for the ΦZR, where this percentage is always below 5%.We ascribe this difference to a size effect: while the size of the galaxies is accounted for in the gravitational potential (estimated as M /R e ), it is not in the stellar mass.Therefore, at a fixed stellar mass, a more compact galaxy would present a higher gas density, and thus a higher SFR (Camps-Fariña et al. 2023).This effect would be stronger at lower masses because the dynamical range of R e is larger.At a fixed Φ baryon , since the size effect is accounted for, the dependence with the SFR disappears.

ΦZR
Our results hold independently of the adopted Z g indicator.This is shown in Appendix A, where we reproduce the analysis with 15 alternative estimators, revealing that for most cases an indicator of the global stellar metallicity is the galactic property with the highest effect on the MZR and ΦZR residuals.In the remaining cases, we find a stronger dependence on the Z g gradient that we ascribe to the uncertainty in the determination of R e , leading to a wrong estimation of the characteristic Z g .
The secondary dependence of the baryonic ΦZR on the stellar metallicity that was found is also confirmed when exploring the total ΦZR using M /R 0.6 e as a tracer of the total gravitational potential (the exponent α 0.6 accounts for the inclusion of the dark matter component; see Paper I).In Appendix B we look at the total ΦZR residuals, removing the dependence of Z g on M /R 0.6 e instead of Φ baryon , and we obtain very similar results to the latter, with the stellar metallicity being the parameter with the strongest effect on ∆Z g .
Previous studies attribute the dependence of the MZR scatter found on the SFR to the significant effect of recent gas accretion on the present-day chemical distribution of galaxies.Our results suggest otherwise.The strongest role played by the stel- 5 We note that M gas was indirectly estimated using the dust extinction as a tracer via the dust-to-gas relation, whereas the SFR was determined from dust-corrected Hα luminosity.A more direct or alternative measurement of the parameters might yield different results.lar metallicity in the RF models for predicting the MZR and ΦZR residuals provides evidence that the present-day gas metallicity is mostly driven by the overall chemical enrichment history (ChEH), and not the most recent episodes (of which the present SFR is indicative).In this sense, the fact that T95-[Z/H] Re or LW-[Z/H] Re are the estimators of the global stellar metallicity that best predict ∆Z g reinforces this idea.The mean LW-[Z/H] enhances the contribution of the youngest stars compared with the mean MW-[Z/H], which makes it more sensitive to the ChEH because the content of old stars is similar in all galaxies (in the analysed mass range).In any case, as discussed, the differences obtained by considering different estimators of stellar metallicity are small; what is relevant is the role played by stellar metallicity over other galaxy properties on both relations.

Summary and conclusions
We have analysed the existence of secondary dependences in the scaling relations between stellar mass and gas metallicity (MZR) and between gravitational potential and gas metallicity (ΦZR).Using a RF regressor, and based on a sample of 3430 star-forming MaNGA galaxies, we looked for the galactic property that best predicts the MZR and ΦZR residuals.We paid special attention to the role of the SFR in the model, which is claimed to be a primordial one according to the fundamental metallicity relation.
We conclude that the global stellar metallicity is the galaxy property with the strongest effect on the MZR and ΦZR residuals.This is true independently of the different adopted indicators of the gas metallicity (Appendix A) and considering both the baryonic and the total ΦZR (Appendix B).This parameter reduces the scatter in the MZ and ΦZ relations by 10−15%, whereas the SFR only reduces the scatter by 0.5−2%.For particular mass ranges (specially for log(M ) ∼ 9 M ), the reduction obtained with the SFR can improve up to 20% for the MZR, but it is always below the decrease reached with the stellar metallicity.These results suggest that the present-day gas metallicity distribution is mostly affected by the overall chemical enrichment history of the galaxy, rather than recent events (of which the present SFR is representative) driven by gas accretion, as previously claimed in the literature.
Appendix A: The effect of the adopted metallicity calibration The determination of gas metallicity is a complex process fraught with numerous systematics and sources of uncertainty.
Well-documented discrepancies exist when different strongline calibrators are employed (e.g.Kewley & Ellison 2008;López-Sánchez et al. 2012;Kewley et al. 2019).In this appendix we show that our results are not contingent upon the choice of the adopted gas-metallicity indicator.For that, we conducted a parallel analysis based on the oxygen abundance estimates from 15 alternative calibrations included in the pyPipe3D VAC table 6 .All of these calibrators present a monotonic behaviour in the covered range of abundances of this study (8.2 < 12 + log(O/H) < 8.6), and therefore they allowed us to unambiguously measure the gas metallicity for all galaxies in the sample.Table A.1 lists the calibrations used and the corresponding references.Similar to the procedure followed in Sec. 4 for the main abundance calibrator, for each alternative indicator, we ran a RF using the residual Z g (∆Z g ) from the ΦZR as the targets and removing Φ baryon from the input parameters.In the top panel of Figure A.1, we show a histogram with the most relevant parameter in the model according to its importance value in all 15 cases.We can see that in ten of them, an indicator of the global stellar metallicity is again the galactic property with the highest effect on ∆Z g .In the remaining five cases, a stronger dependence with the gas metallicity gradient (α([O/H])) is found.We attribute this latter finding to the uncertainty in the determination of R e , which can result in the Z g measured at this R e not being completely representative of the global gas metallicity (and thus the need to apply a small correction considering the existing negative gradient appears).Alternatively, we replicated this analysis with the residual Z g from the MZR, whose outcome is represented in a similar way in the bottom panel of Figure A.1.In this case, a secondary dependence of Z g with the global luminosityweighted stellar metallicity is also reported for the majority of cases.We consider that the analyses described in this appendix reinforce the role of the stellar metallicity as the strongest secondary parameter shaping the gas metallicity in galaxies.
Appendix B: The ΦZR relation based on M /R 0.6 e as a tracer of the total gravitational potential In Paper I we showed how the inclusion in the RF of a parameter of the type M /R α e with α = 0.6 (different coefficients for α were explored) performs better than Φ baryon = M /R e when predicting Z g .We investigated the effect of the dark matter (DM) halo on the gravitational potential and we argued how a scale of the form M /R 0.6 e for the total Φ matches very well the theoretical relation between the DM fraction and the baryonic surface density predicted in cosmological numerical simulations of galaxy formation (Nestor Shachar et al. 2023).Thus it is important to confirm that the stellar metallicity plays the same role as the secondary parameter in the total ΦZR.
For this test, we included M /R 0.6 e (as a tracer of Φ T ) as an input parameter in the RF and we repeated the procedure followed in Sec. 4. The black line in Fig. B.1 (bottom x-axis) represents the decrease in the dispersion of Z g as we modelled its relation with different galaxy parameters.Similar to Fig. 2, the secondary x-axis (top) shows the case when we forced the L11, page 7 of 8 stellar mass to be the first parameter to model (orange line).After fitting the dependence on Φ T (44%), the RF again revealed the strong effect of the stellar metallicity (in particular, T95-[Z/H] Re ) on predicting ∆Z g , which is the parameter in the model with the highest importance value and the one that reduces the dispersion of the residual metallicities the most (14%).The subsequent modelling of the residuals with the remaining galaxy properties produces a decrease in the dispersion always below 5%.This test confirms the main conclusion reached in our study, that is, the stellar metallicity is the most important secondary parameter in the ΦZR (both baryonic and total) and the MZR.

Fig. 1 .
Fig. 1.Global gas metallicity as a function of the galactic gravitational potential (ΦZR) colour-coded with the stellar metallicity measured at 1 R e and at the look-back time T95 (T95-[Z/H] Re , left panel), and with the SFR (right panel).Salmon squares represent the median Z g values in ten bins.The averaged-binned values were fitted with a spline function (solid salmon lines).The insets show a scatter plot Z g vs. T95-[Z/H] Re (left) and Z g vs. SFR (right) for the galaxies in the bin 8.8 ≤ Φ baryon ≤ 9.0, which is marked with vertical dashed lines in the main figure.The solid brown lines represent an ODR fit to the scatter points.The Pearson correlation coefficient r is also displayed.

Fig. 2 .
Fig.2.Decrease in the dispersion (in percent) of the original or residual global gas metallicity once the subsequent relation with different galaxy properties (x-axis) was modelled (with a spline function on binned values) and removed.This exercise is accumulative, so the modelled relation was subtracted to the residuals from the previous one.The secondary x-axis (top) shows the case when we forced the stellar mass to be the first parameter to model (orange line).

Fig. 3 .
Fig. 3. Global gas metallicity as a function of the stellar mass (MZR) colour-coded with the LW-stellar metallicity measured at 1 R e (LW-[Z/H] Re , left panel), and with the SFR (right panel).Salmon squares represent the median Z g values in ten mass bins.The averaged-binned values were fitted with a spline function (solid salmon lines).The insets show a scatter plot Z g vs. LW-[Z/H] Re (left) and Z g vs. SFR (right) for the galaxies in the mass bin 9.2 ≤ log(M ) ≤ 9.4, which is marked with vertical dashed lines in the main figure.The solid brown lines represent an ODR fit to the scatter points.The Pearson correlation coefficient r is also displayed.

Fig. A. 1 .
Fig. A.1.Histogram of the most relevant parameter in the RF when studying the residuals of the ΦZR (top panel) and the MZR (bottom panel) using 15 different estimators of Z g (see Table A.1 for references).For the top panel, LW (MW) means light-(mass-)weighted and T95 (T99) represents the look-back time at which the galaxy has formed 95% (99%) of its mass.All of these attributes refer to different estimations of the stellar metallicity measured at 1 R e .

Fig. B. 1 .
Fig. B.1.Analogous to Fig. 2 but including a tracer of the total gravitational potential in the model (Φ T ).See the caption of Fig. 2 for details.

Table A .
1. List of alternative calibrations used to derive the oxygen abundances.