The metallicity gradients of star-forming regions store information of the assembly history of galaxies

The variations in metallicity and spatial patterns within star-forming regions of galaxies result from diverse physical processes unfolding throughout their evolutionary history, with a particular emphasis in recent events. Analysing MaNGA and \textsc{eagle} galaxies, we discovered an additional dependence of the mass-metallicity relation (MZR) on metallicity gradients ($\nabla_{{\rm (O/H)}}$). Two regimes emerged for low and high stellar mass galaxies, distinctly separated at approximately ${\rm M_{\star}}>10^{9.75}$. Low-mass galaxies with strong positive $\nabla_{{\rm (O/H)}}$ appear less enriched than the MZR median, while those with strong negative gradients are consistently more enriched in both simulated and observed samples. Interestingly, low-mass galaxies with strong negative $\nabla_{{\rm (O/H)}}$ exhibit high star-forming activity, regardless of stellar surface density or $\nabla_{{\rm (O/H)}}$. In contrast, a discrepancy arises for massive galaxies between MaNGA and \textsc{eagle} datasets. The latter exhibit a notable anticorrelation between specific star formation rate and stellar surface density, independent of $\nabla_{{\rm (O/H)}}$, while MaNGA galaxies show this trend mainly for strong positive $\nabla_{{\rm (O/H)}}$. Further investigation indicates that galaxies with strong negative gradients tend to host smaller central black holes in observed datasets, a trend not replicated in simulations. These findings suggest disparities in metallicity recycling and mixing history between observations and simulations, particularly in massive galaxies with varying metallicity gradients. These distinctions could contribute to a more comprehensive understanding of the underlying physics.


INTRODUCTION
The presence of oxygen radial profiles in the star-forming regions of galaxies is a well-determined observational fact (Searle 1971;Peimbert 1975).The slope of these profile, hereafter metallicity gradients ∇ (O/H) , in massive, spiral galaxies tend to become shallower with increasing mass, while a large variety of ∇ (O/H) are detected for small galaxies (Zaritsky et al. 1994;Ho et al. 2015).More recently, large surveys of galaxies using integral field spectroscopy (IFS) such as CALIFA (Sánchez et al. 2013), MaNGA (Bundy et al. 2015) and SAMI (Poetrodjojo et al. 2018) have al-★ E-mail: fujara@uc.cllowed a more detailed characterization of the metallicity distribution in nearby galaxies.From this data, a weak anti-correlation between ∇ (O/H) and stellar mass for the low-mass end has been detected, so that that slopes became shallower or even positive 1 (Sánchez et al. 2013;Sánchez-Menguiano et al. 2016;Belfiore et al. 2017).Positive gradients have also been observed at high redshift, particularly, in connection with gas-rich galaxies that are frequently undergoing interactions (e.g Troncoso et al. 2014;Jones et al. 2015;Carton et al. 2018;Wang et al. 2019Wang et al. , 2022).
1 For these galaxies, the central gas-phase metallicity is shallower than in the outskirts.These ∇ (O/H) will be also referred as inverted metallicity gradients.
Metallicity gradients are expected to be shaped by the action of different physical processes (Tinsley 1980;Bland-Hawthorn & Freeman 2003), such as star formation, energy feedback as well as gas inflows, mergers, interaction and probably environmental mechanisms such as ram pressure and tidal striping (Bahé et al. 2017).Hence, they encode information on these processes, which can provide clues for understanding the history of evolution of galaxies (Pilkington et al. 2012;Tissera et al. 2016a;Maiolino & Mannucci 2019;Sharda et al. 2023).
Within the standard inside-out disc formation scenario (Fall & Efstathiou 1980), negative metallicity gradients are expected to form as the star formation takes place in high-density gas clouds, starting from the central regions and moving out progressively to the outskirts of the discs.If the gas accretion continues in a steady fashion, this gas will settle onto the disc, fuelling the ongoing star formation activity and setting off supernovae (SNe) feedback.Energy feedback contributes to modify the thermodynamical properties of the interstellar medium (ISM) (e.g.Wiersma et al. 2009b;Scannapieco et al. 2008), to mix chemical elements or even to eject part of them via galactic winds (e.g.Stinson et al. 2006;Hopkins et al. 2013).Previous works have shown that SNe feedback is one of the main processes that regulate the slope of metallicity gradients together with the availability of gas to form stars and the efficiency of transformation of gas into stars (Pilkington et al. 2012;Lilly et al. 2013;Stinson et al. 2013;Gibson et al. 2013;Ma et al. 2017;Mollá et al. 2017;Belfiore et al. 2019;Tissera et al. 2019;Wang & Lilly 2022).
However, galaxies are rarely isolated entities, but have ongoing interactions with their surroundings.These interactions can lead to the capture of gas and stars from neighbouring galaxies as well as the stripping of material from the galaxies themselves, among other possible effects (e.g.Karademir et al. 2019;Pallero et al. 2022;Rodríguez et al. 2022).These mechanisms could affect the metallicity distribution in the star-forming regions, shaping the metallicity gradients and the level of enrichment of galaxies (Franchetto et al. 2021;Sharda et al. 2021;Maier et al. 2022).Previous works have shown that galaxy-galaxy interactions and mergers can generate tidal torques that affect the distribution of gas and stars, and hence, the star formation activity and metallicity distribution across galaxies (Rupke et al. 2010;Di Matteo et al. 2013;Sillero et al. 2017;Moreno et al. 2019).These tidal torques are also expected to trigger bars, which have been shown to be efficient at transporting gas into the central regions as the gas loses angular momentum (Weinberg 1985;Athanassoula 2003).This gas has the potential to boost star formation activity in these regions and hence, to modify the chemical abundances of the ISM as the galaxies draw closer to each other (Katz et al. 1992;Tissera 2000;Lambas et al. 2003;Torrey et al. 2012;Moreno et al. 2019).
These gas inflows induced during galaxy-galaxy interactions can also explained the dilution of the central oxygen abundance as shown by numerous observational and numerical works (e.g.Kewley et al. 2006;Rupke et al. 2010;Kewley et al. 2010;Perez et al. 2011;Tissera et al. 2016a;Bustamante et al. 2018;Moreno et al. 2019;Bustamante et al. 2020).Depending on the characteristics of the interactions and gas richness of the encountering galaxies, strong starbursts could be generated, which on their turn, can induce metal-loaded outflows (e.g Scannapieco et al. 2009;Stinson et al. 2013).This sequence of events is expected to change the chemical distribution of the ISM.In the first stage, the transport of material from outer regions to the centre will dilute the metallicity but as new stars are formed, the -enrichment will increase by the contribution of Type II SNe (SNeII).Depending on the strength of the starbursts and the potential well of the galaxies, metal-loaded galactic winds can potentially diminish the level of enrichment and modify the thermodynamical properties of the ISM by expelling enriched material or by redistributing it via galactic fountain, regulating the subsequent star formation activity.Sillero et al. (2017) showed that if a gas inflow is triggered in a short time-scale, ∼ 0.5 Gyr, then, it could both induce a starbursts and change metallicity slopes.Active Galactic Nuclei (AGN) activity could be also triggered during these events and have been claimed to be very efficient at quenching star formation in massive galaxies (e.g.Rosas-Guevara et al. 2019).However, it is not yet clear enough how AGN feedback can modify the thermodynamical properties of the ISM and its chemical characteristics across the discs (Hopkins et al. 2013;Rosas-Guevara et al. 2015).
Collectively, the aforementioned physical processes govern the star formation activity, the degree of enrichment, and the distribution of chemical elements within galaxies of varying stellar masses across different evolutionary phases.Consequently, these processes can potentially imprint discernible effects on relations such as the mass-metallicity relation, MZR.
Although, the MZR has been extensively studied (Maiolino & Mannucci 2019, and references there in), it remains to be fully understood which physical mechanisms could drive it and what is the origin of its scatter as a function of time (Curti et al. 2023).Ellison et al. (2008) found that part of the scatter could be related to the size and the star formation activity so that galaxies with larger star formation rates, SFR, tend to have lower metallicity (see also Mannucci et al. 2010;Lara-López et al. 2010).A similar correlation of the scatter with the HI mass has been also reported (Bothwell et al. 2013).
Observations have also found that galaxy interactions could shift the location of galaxies on the MZR.Michel-Dansac et al. (2008) studied pairs of galaxies selected from the SDSS-DR4 and found that merging galaxies exhibited an excess of metallicity in massive galaxies, while smaller galaxies showed a dilution of metallicity.Similarly, Omori & Takeuchi (2022) reported a dilution of metallicity in very close pairs compared to systems with members located at larger distances.Both studies attributed the variations in chemical enrichment in pairs to the effects of the interaction itself, as it is actually predicted by numerical simulations (Perez et al. 2006(Perez et al. , 2011;;Torrey et al. 2012;Tissera et al. 2016b;Sillero et al. 2017).Recent observational studies of the resolved MZR have reported that a key property influencing its scatter is the stellar surface density or stellar compactness, Σ ★ (Boardman et al. 2022).Furthermore, these studies indicate that metallicity gradients may not simply be a byproduct of the resolved properties, which implies a direct increase of the chemical abundances for increasing Σ ★ , but could have an intrinsic origin related to the action of other processes such as gas inflows and mergers (Baker et al. 2023;Boardman et al. 2023).
From a numerical point of view, Torrey et al. (2019) found that the scatter in the simulated MZR estimated for galaxies in Illustris TNG50 is correlated with the gas fraction and star formation activity.The eagle simulations have been extensively used to study the MZR and its scatter (e.g.Crain et al. 2015;Lagos et al. 2016;De Rossi et al. 2017;Trayford & Schaye 2019).In particular van Loon et al. (2021) found that the scatter of the MZR in eagle simulations was related to variations of long time-scales, probably reflecting the combined action of inflows, outflows and star formation.These results agree with previous findings from analytical modelling of the relation between SFR, metallicity and gas fraction and the time variation (e.g.Forbes et al. 2014;Wang & Lilly 2021).Recently, Zenocratti et al. (2020Zenocratti et al. ( , 2022) ) studied the MZR of eagle galaxies with stellar mass in the range [10 9 , 10 10 ]M ⊙ .These authors studied the level of enrichment in relation to galaxy morphology, finding that, at given stellar mass, more disc-dominated galaxies tended to have lower level of enrichment than dispersion-dominated systems.They related this trend to the late assembly history of discdominated galaxies and hence, the latter accretion of low-metallicity gas, compared to spheroidal-systems.
Therefore, the simulated and observational results on the MZR and the the mass-metallicity gradients relation, MZGR, have prompted us to further analyse the information contained in the ∇ (O/H) of star-forming regions in relation to their level of enrichment.In this paper, we explore the potential link between the MZR and MZGR to unravel fossil records of their evolution, possible stored by the metallicity gradients, by exploring a possible secondary dependence of the MZR on ∇ (O/H) .
To achieve this goal, we resort to galaxies from the MaNGA survey (Sánchez et al. 2013) and simulated galaxies from the higher resolution run, Recal-L025N0752, of the eagle project (Schaye et al. 2015).We conduct a comparative analysis between galaxies with different ∇ (O/H) , grouping them according to their metallicity gradient and stellar mass.For this purpose, we follow Tissera et al. (2022, hereafter T22) who studied the metallicity gradient evolution and the relation with the merger and gas inflow histories.These authors found that galaxies with recent strong accretion (i.e.either mergers or gas inflows) tend to be more star-forming and have strong negative or positive gradients and that after about 2 Gyr they recovered their weak metallicity gradients.This work is be taken as reference and, thus we adopt their definitions as stated throughout the paper.
This paper is organised as follows.In Section 2 we summarises the main characteristics of selected galaxies from the MaNGA survey and the eagle simulation.In Section 3, we analyse the relations between the metallicity gradients, specific star formation rate, sSFR, and strong accretion history and AGN activity.Section 4 summarises our main findings.

The MaNGA galaxies
We use the galaxy catalogue built with Pipe3D (Sánchez et al. 2016) on the Data Release 16 of the Mapping Nearby Galaxies at Apache Point Observatory (MaNGA, Bundy et al. 2015).The MaNGA survey is one of the projects included in the Sloan Digital Sky Survey IV (SDSS-IV: Blanton et al. 2017).This survey uses integral field units (IFU) and has generated a database of about 10000 spatially resolved nearby galaxies, with a wide variety of morphologies and a relatively uniform mass distribution.Only galaxies with inclinations lower than 75 degrees have been considered for the analysis.From the catalogue, we also obtained the stellar mass, M ★ , the star formation rate, SFR, and the effective radius,  eff2 .
We selected galaxies based on four criteria: (i) redshift  < 0.05, (ii) stellar masses, M ★ ∈ [10 9 , 10 11 ]M ⊙ , (iii) metallicity gradients between −0.2 and 0.2 dex  eff −1 , and (iv) sSFR within the range covered by the eagle galaxies for each of the three defined mass intervals (see Section 2.2).This fourth selection criterion constrains the sSFR to [10 −10.44 , 10 −9.40 ]yr −1 for low-mass, [10 −11.32 , 10 −9.57]yr −1 for intermediate mass and [10 −11.46 , 10 −9.75 ]yr −1 for high mass galaxies.By applying these criteria, we obtained our main sample of 1765 galaxies, hereafter referred to as the MaNGA sample.We selected the metallicity gradients measured in the radial range [0.5, 1.5] eff , using the abundance of oxygen as a tracer.Our reference metallicities were measured using the O3N2 calibrator (Marino et al. 2013), which is a direct method that depends on two strong emission lines; [NII] 6583 and [OIII] 5007.However, we also explored other calibrators finding that most of them produced similar relations, as shown in Appendix B.
We acknowledge the fact that all the metallicity gradients used in this paper were obtained by fitting a single linear regression to the metallicity profiles within the radial range [0.5 − 1.5]  eff .Recent observational results suggest that more than one linear fit might be required to fully describe the metallicity profiles (e.g.Diaz 1989;Sánchez-Menguiano et al. 2016;Barrera-Ballesteros et al. 2022).Numerical simulations also report the existence of breaks from a single power law, which could be reflecting the action of different physical parameters such as accretion, gas inflows, galactic fountains (e.g.Garcia et al. 2023, Tapia in preparation).Because in this study, we specifically concentrate on the regions dominated by the disc components, we anticipate no significant changes because of this, as breaks are typically identified in the inner regions or outskirts (Sánchez-Menguiano et al. 2016;Tapia et al. 2022).
Additionally, we estimated the mass of central black hole, M BH , by using the known correlation between M BH and the stellar velocity dispersion of the central region,   .The stellar velocity dispersion in the central 2.5 arcsec/aperture was used as given by Pipe3D.The correlation reported by Tremaine et al. (2002) was followed to estimate M BH .The authors proposed empirical values for the coefficients of a linear correlation with the form log(M BH ) =  +  log(/ 0 ), which are  = 8.13 ± 0.06 and  = 4.02 ± 0.32, for  0 = 200 km s −1 .Hence, we obtained M BH with a median error of ∼ log(0.12).
In this work we use the Recal-L025N0752 run, which corresponds to a volume of 25 cMpc (co-moving mega-parsecs) size length, and uses 752 3 initial baryonic and dark matter particles.The mass resolution is 2.26 × 10 5 M ⊙ and 1.21 × 10 6 M ⊙ for the initial gas and dark matter particles, respectively (see Schaye et al. 2015, for more details).This corresponds to the highest resolution within the eagle suite.And, in particular, the recalibrated model predicts a MZR slope in better agreement with observations of low-mass galaxies (Schaye et al. 2015;De Rossi et al. 2017).
The simulations were performed with the so-called ANAR-CHY version of p-gadget-3 (Schaye et al. 2015;Crain et al. 2015, see for more details on the subgrid physics).Briefly, it includes a star formation algorithm which forms stars stochastically, following the model of Schaye & Dalla Vecchia (2008), from cold and dense gas with a metallicity dependent star formation density threshold (Schaye et al. 2010).The cooling and photo-heating models are implemented following Wiersma et al. (2009a).The stellar feedback is treated stochastically, adopting the thermal injection scheme described in Dalla Vecchia & Schaye (2012).The chemical model follows eleven chemical elements (Wiersma et al. 2009b).The eagle simulations adopt the initial mass function of Chabrier ( 2003), with minimum and maximum mass cut-offs of 0.1 M ⊙ and 100 M ⊙ , respectively.
The AGN feedback model locates black holes, BHs, with an initial mass of m seed = 1.48 × 10 5 M ⊙ at the centres of haloes with a dynamical mass exceeding ∼ 10 10 M ⊙ (Springel 2005).The growth of BHs occurs through gas accretion and mergers.The accretion rates of the supermassive black holes, SMBH, are computed using the modified Bondi-Hoyle accretion rate formula as explained by Rosas-Guevara et al. (2015); Schaye et al. (2015).The AGN feedback parameters were calibrated to replicate the stellar mass galaxy function, the M ★ − M BH and the scaling relation between stellar mass and the mass of the central black hole for galaxies observed in the Local Universe as described in Crain et al. (2015).
To explore the possible impact of AGN feedback we use the Eddington rates,  Edd 5 , and M BH .In eagle the BH accretion regimes are controlled by the  Edd as described by Rosas-Guevara et al. (2016), where galaxies with  Edd > 10 −2 are assumed to be highly accretive, X-rays luminous sources; galaxies with 10 −4 <  Edd < 10 −2 are expected to be radiatively inefficient (Narayan & Yi 1994;Abramowicz et al. 1995); and galaxies with  Edd < 10 −4 are considered as inactive.Nevertheless, we stress the fact that  Edd show large variability as discussed by Rosas-Guevara et al. (2019).Thus, whilst still showing the trend with  Edd , we focus the main analysis of the AGN impact on the sSFR and metallicity gradients by using the M BH provided by the eagle database.The detailed estimations of M BH can be found in Rosas-Guevara et al. (2016).

eagle galaxy sample
Galaxies were identified by using the Friends-of-Friends (FoF, Davis et al. 1985) and subfind (Springel et al. 2001;Dolag et al. 2009) algorithms.Our analysis is restricted to central galaxies.We adopted the galaxy catalogue built by T22.These authors identified the discs and spheroidal components based on the angular momentum and the binding energy (AM-E method) of the stars and gas as described by Tissera et al. (2012).This decomposition also allows the estimation of the stellar disc-to-total mass, D/T, which will be used as a morphology estimator (e.g.Pedrosa & Tissera 2015;Correa et al. 2017).
For each simulated central galaxy in the eagle sample, the stellar mass (M ★ ), the stellar half-mass radius ( eff ), the sSFR, the gas fraction (  gas = M gas /(M ★ + M gas )) and the stellar surface density (Σ ★ ) were computed.Additionally, we adopted the parameter  25 defined by T22, which measures the time elapsed since the last major increase of stellar mass, by more than 25 per cent, between two available snapshots.Hence, we are agnostic with respect to the processes that could caused it (e.g.gas inflows or mergers).These 5 The Eddington rate can be estimated as where  is the efficiency (fixed to 0.1),  T is the Thompson cross-section and   is the mass of the proton.
parameters have been estimated by using the gas or stellar particles within 1.5 times the optical radius6 .The ∇ (O/H) of the star-forming regions in the discs components are calculated by using the the oxygen abundances of the gas particles that satisfy all simulation requirements for being converted to stars but have not being yet transformed.As explained by T22, galaxies with stellar masses below 10 9  ⊙ , hence, having typically fewer than 1000 gas particles, are not considered in order to mitigate numerical effects on the estimation of the metallicity gradients due to limited numerical resolution.
To enhance the realism of the comparison with observations, the oxygen profiles are weighted by the SFR of the corresponding regions.Then, the logarithm azimuthal-averaged radial profiles are fitted by a linear regression within the radial range ∈ [0.5, 1.5]  eff in those galaxies with more than 100 star-forming gas particles.Galaxies with  eff < 1 kpc (about three times the gravitational softening) and stellar masses lower than 10 9 M ⊙ were excluded from the sample to mitigate the effects of low numerical resolution on the determination of the metallicity profiles.
The global metallicity of the eagle sample is estimated as the global oxygen abundance defined by the star forming regions in each galaxy (de Graaff et al. 2022) within 1.5  eff .De Rossi et al. (2017) reported that galaxies in the Recal-L025N0752 run have a MZR in good agreement with observations.We note that our estimations are all referred to 1.5  eff , instead of adopting a fixed aperture radius (Zenocratti et al. 2022).Also, our sample is restricted to star-forming galaxies, which satisfy the minimum number of starforming gas particles to reliably estimate the metallicity gradients, so galaxies with low star-formation activity will not be included.
Throughout our analysis, we use results from Tissera et al. (2022).To ensure consistency, we adopt the same criteria for grouping galaxies into stellar mass and metallicity gradient intervals.Three mass intervals are defined: [10 9 , 10 9.75 )M ⊙ , [10 9.75 , 10 10.25 ]M ⊙ and (10 10.25 , 10 11 ]M ⊙ .Additionally, galaxies are analysed according to three ∇ (O/H) subsamples: strong negative (< −0.02 dex kpc −1 ), strong positive (> 0.02 dex kpc −1 ) and weak metallicity gradients.This procedure is applied to both simulated and observational data.While we use metallicity gradients in dex kpc −1 to define the samples for direct comparison with the results of Tissera et al. (2022), we also estimated them normalised by  eff and checked that this does not affect the trends.In Table 1, the median properties of galaxies in each mass interval are summarised.
We note that, in eagle we only consider central galaxies (i.e. the most massive galaxy within a dark matter halo).However, for MaNGA we have not distinguished between central and satellite galaxies.These could introduce same differences considering that satellite galaxies seem to be slightly more enriched and less starforming, particularly in high density regions (Pasquali et al. 2012).However, since we are focusing on star-forming galaxies in both observations and simulations, we expect the comparison not to be strongly affected by the possible presence of satellite galaxies in the MaNGA sample.
Table 1.Summary of the main properties in our eagle subsamples at  = 0.0.From the left to the right columns the fraction with respect to the complete sample, the medians of the gas-phase metallicity gradient, global metallicity, effective radius, specific star formation rate, gas fraction and Eddington ratio are displayed.The medias are shown for low-mass galaxies, M ★ < 10 9.75 M ⊙ , intermediate mass galaxies, M ★ = [10 9.75 , 10 10.25 ] M ⊙ , and high-mass galaxies, M ★ > 10 10.25 M ⊙ .The strong negative group (S-Neg.)corresponds to galaxies with ∇ O/H < −0.02 dex kpc −1 , the gradients within the range of (−0.02, 0.02) dex kpc −1 are defined as weak group, and the strong positive group (S-Pos.)corresponds to galaxies with ∇ O/H > −0.02 dex kpc −1 (Tissera et al. 2022).The associated uncertainties correspond to 95 per cent confidence limits calculated by using a bootstrap technique.From the left to the right columns the fraction with respect to the complete sample, the medians of the gas-phase metallicity gradient, the global metallicity, the effective radius (kpc), and the specific star formation rate (yr −1 )are displayed.Additionally,  rot / and D/T have been included for the MaNGA and eagle samples, respectively.The medias are shown for low-mass galaxies, M ★ < 10 9.75 M ⊙ , intermediate mass galaxies, M ★ = [10 9.75 , 10 10.25 ] M ⊙ , and high-mass galaxies, M ★ > 10 10.25 M ⊙ .The strong negative group (S-Neg.)corresponds to galaxies with ∇ O/H < −0.02 dex kpc −1 , the gradients within the range of (−0.02, 0.02) dex kpc −1 are defined as weak group, and the strong positive group (S-Pos.)corresponds to galaxies with ∇ O/H > −0.02 dex kpc −1 (Tissera et al. 2022).The associated uncertainties correspond to 95 per cent confidence limits calculated by using a bootstrap technique.

ANALYSIS AND RESULTS
In this section, we first investigate the MZR for galaxies exhibiting different metallicity gradients.Figure 1 presents the MZR for both MaNGA and eagle samples in bins of ∇ (O/H) , categorised by the metallicity gradient intervals.The solid black line represents the main MZR relation defined by the MaNGA sample.In order to better compare the behaviour of both samples on the MZR, we renormalised the gas-phase metallicity of the eagle galaxies, which are known to be systematically more metal-rich compared to the MaNGA galaxies (Bahé et al. 2017).In Appendix A1 we explained this in more detail (see also Fig. A1).
As can be seen from Fig. 1 galaxies with weak metallicity gradients follow a similar relation to the full sample as they constitute the dominant population.On the other hand, galaxies with strong positive gradients are consistently less enriched at a given stellar mass, while those with strong negative gradients are more enriched for M ★ ≲ 10 10 M ⊙ .For higher masses, the three subsamples tend to converge to a comparable level of enrichment, albeit with those possessing positive gradients located systematically below the full MZR.Moreover, it is worth mentioning that galaxies with strong positive gradients are less prevalent among intermediate and high mass galaxies.We also notice that there is a large dispersion in the eagle samples, which is partly due to the lower number of members in these subsamples, particularly at the high-mass end, compared to the MaNGA subsamples.
The scatter in the observed MZR has been shown to have a secondary dependence on sSFR and  eff , in the sense that galaxies with higher sSFR tend to have lower levels of enrichment and to be more extended, at a given stellar mass (Ellison et al. 2008).As expected the MaNGA sample is consistent with this trend, as shown in the left panels of Fig. 2. To reveal the correlation between the properties, the Locally Weighted Scatterplot Smoothing (LOESS; Cappellari et al. 2013) technique was applied to the MZR as a function of sSFR and  eff .For low-mass galaxies, we also computed a correlation angle, , performing a PCC analysis (see Appendix A2).This angle reflects how dominant a variable is on a plane formed by two other variables.Our eagle sample shows a comparable trend to that defined by MaNGA galaxies, albeit with larger dispersion (see also Zenocratti et al. 2022, who analysed a larger sample of eagle galaxies and found similar results).Although the distributions are noisier due to the smaller number of galaxies in the simulated sample, the trends are globally reproduced as quantified by  in Fig. 2. As mentioned before we note that low-mass galaxies with low sSFR are not present in the eagle sample probably because of our strict selection criteria (see also Appendix A1). Figure 1.The mass-metallicity relation as a function of the metallicity gradients of galaxies in our MaNGA and eagle samples.The median MZR for galaxies with strong negative (red), weak (pink) and strong positive (blue) metallicity gradients for MaNGA (symbols and lines) and the eagle (shaded areas) samples are denoted according to information in the inset.The median MZR for the full MaNGA sample (black, solid line) and the 95% confidence intervals (black, dotted-dashed lines).Error bars and the widths of the shaded areas were estimated by using a bootstrap technique.

Metallicity gradients, the specific star formation and the stellar surface density
In order to relate the secondary dependence of MZR on ∇ (O/H) found in Fig. 1 with the star formation activity and the size, we investigate the relationship between the sSFR and Σ ★ , defined as M ★ / eff 2 , for galaxies with different ∇ (O/H) in each mass interval.Σ ★ provides a metric for the stellar compactness of a galaxy, indicating that, at a specific stellar mass, there is a spectrum of potential sizes.In Table 2 the medians of the parameters analysed are displayed 7 .
The upper panels of Fig. 3 show the relationship between sSFR and Σ ★ for galaxies in the three defined stellar mass intervals and for the three defined gradient subsamples for the MaNGA galaxies (solid lines and symbols).As can be seen in this figure, lowmass galaxies (left panel) exhibit a high level of sSFR regardless of ∇ (O/H) , with no clear trend present between Σ ★ and sSFR, within the corresponding bootstrap errors.
However, intermediate-mass galaxies exhibit different behaviours depending on the metallicity gradient.As depicted in the middle panel of Fig. 3, galaxies with strong positive gradients are systematically less star-forming and more compact.On the other hand, galaxies with weak and strong negative gradients continue to form stars actively.
In the high stellar mass interval (right panel of Fig. 3), galaxies tend to be more quenched as their ∇ (O/H) becomes shallower and 7 In order to increase the number of eagle galaxies, the simulated sample is extended by including those with  ≤ 0.18.Although this redshift range is wider than that of MaNGA, it allows us to compensate for the low volume size and render a more robust analysis.We confirmed that the dependence on ∇ (O/H) is still present although the full MZR is shallower (De Rossi et al. 2017).more positive.However, galaxies with strong negative gradients continue to form stars actively, regardless of their Σ ★ .These galaxies show only a small decrease in their sSFR, of at most ∼ −0.10 dex, as their stellar mass increases from the low to the high end (see Table 2).On the other hand, galaxies with positive ∇ (O/H) exhibit lower sSFR levels for all Σ ★ values and the decrease of sSFR is ∼ −0.40 dex from low to high mass galaxies (Table 2).
The upper, left panel of Fig. 3 reveals that compact, low stellar mass galaxies tend to have higher frequencies of both strong positive and negative ∇ (O/H) .Conversely, galaxies with weak gradients are more frequent in less concentrated galaxies which suggest they tend to have a large frequency of extended and disc-dominated galaxies.The trend with size can be also seen in Table 2, which displays the median sizes for galaxies in each mass interval and ∇ (O/H) subsample.
The upper, middle and right panels of Fig. 3 show that for galaxies with masses higher than 10 9.75 M ⊙ , weak gradients are also more common and evenly distributed at all Σ ★ .The presence of strong types becomes less significant for increasing galaxy mass.Galaxies exhibiting strong positive gradients are rare and tend to be more concentrated within a given mass interval, whereas galaxies with strong negative gradients are more evenly distributed as a function of Σ ★ in high mass galaxies than in low-mass systems.These trends, along with those found for sSFR as a function of Σ ★ , persisted for most of the metallicity calibrators that we tested, as shown in Fig. B2.
In summary, low-mass galaxies exhibit high sSFR irrespective of ∇ (O/H) and do not show any trend with Σ ★ , in agreement with the findings for MaNGA galaxies However, the dependence of sSFR on Σ ★ by eagle galaxies is significantly different with respect to the MaNGA's trends for intermediate and high mass galaxies.Both intermediate and high-mass galaxies in eagle exhibit a systematic decline in sSFR as Σ ★ increases, independently of ∇ (O/H) .
From Fig. 3, it can be seen that the quenching start to manifest clearly for both observational and simulated galaxies within the intermediate mass interval.This agrees with the mass range where AGN feedback is expected to act more efficiently in eagle simulations (e.g.van Loon et al. 2021).The transition mass between the regimes dominated by SNe feedback and that dominated by AGN feedback is expected to be around 10 10 M ⊙ .The lack of dependence of the sSFR-Σ ★ relation on ∇ (O/H) in the eagle galaxies might suggest that the AGN feedback model used could be more efficient at decreasing the star formation activity by heating the gas than mixing or expelling the chemical elements in the star-forming regions in the simulated discs.Hence, this tension could, indeed, contribute to improve the subgrid physics related with star formation and feedback mechanisms.Our findings suggest that the diminishing of sSFR as a function of the stellar compactness can be also linked to the characteristics of the distribution of chemical abundances in the discs, encoded by the metallicity gradients.

The impact of strong accretion history in eagle galaxies
In this section, we examine the impact of strong accretion events on the location a galaxy occupies on the MZR in the context of the results reported in Fig. 1.As mentioned before, accretion events could be gas inflows, minor or major mergers.We adopt the parameter  25 to quantify the occurrence of the latest strong event that caused a significant change of stellar mass as explained in Section 2.2.1.
In order to better understand the relationship between the MZR and the recent strong accretion history of galaxies, we have computed the MZR as a function of sSFR,  eff , D/T and  gas , as shown To reveal the correlation between the properties, the Locally Weighted Scatterplot Smoothing (LOESS; Cappellari et al. 2013) technique was applied in all the panels.For low-mass galaxies (the dash-dotted black line denotes M ★ = 10 9.75 M ⊙ ) a correlation angle  was computed using PCC analysis (Partial Correlation Coefficients, Bluck et al. 2020b).This angle denotes the direction in which the third dependency of each panel grows.eagle galaxies follow the observational trends for both parameters, with a global trend for galaxies with high metallicity to have lower star-forming and have smaller  eff for M ★ ≲ 10 10 M ⊙ , at a given stellar mass.Higher masses tend to have lower sSFR and to be more extended.The errors on  were estimated by applying a bootstrap technique.
Fig. 4.Those galaxies which had a main accretion event in the last 3 Gyr, i.e.  25 < 3 Gyr, are highlighted by black circles.This value of  25 is taken from T22 where an excess of strong gradients was detected from galaxies with strong accretion events within this time.We applied the LOESS algorithm (Cappellari et al. 2013) in all the panels.From this figure, we can see that galaxies with M ★ ≲ 10 9.75 M ⊙ (dashed vertical line) and recent accretion events are mostly high star-forming galaxies.Those with low-metallicity tend to be more extended and more disc-dominated.However, at the lowest mass end, there are galaxies which are mostly small and spheroidal dominated (see also Table2. Massive galaxies with recent accretion have similar sSFR and  eff than the rest of the populations as shown by the first and second panels of Fig. 4.However, there is a weak trend to have higher metallicity at a given stellar mass (black encircled symbols).In fact, some eagle massive galaxies with recent accretion have strong negative ∇ (O/H) as pointed out by T22.This could be explained by the fact that massive systems might be more stable against external forces and, hence, the accreted low-metallicity material can settle onto a disc, contributing to steepen the negative metallicity gradients (Collacchioni et al. 2019).They could also have more prominent bulges, which can also provide stability to the systems.From the third panel of Fig. 4 we can see that more disc-dominated galaxies populate the intermediate mass region and that for low mass galaxies, there is a population of more disc-dominated systems with low metallicity, which are also more extended and have low  25 , as displayed in the previous panel.The fourth panel shows the trend with  gas which show the expected trend for most the low mass galaxies with recent accretion events to have larger gas fractions.Massive galaxies are more gas-poor than low mass galaxies as expected.
Hence, at a given stellar mass, high metallicity galaxies could be more passively evolving and show negative metallicity gradients or could have been part of a major merger event, which contributes with low metallicity gas that reinforces the negative slopes and/or triggers strong star formation activity.The latter is found to be more common in low stellar mass galaxies and the former, in massive ones.Galaxies with positive metallicity gradients are preferentially found at lower metallicity for a given stellar mass (Fig. 1), and this could indicate that they have experienced inward radial inflows or strong SNe outflows.These mechanisms dilute the metallicity in the inner region as studied in previous works (e.g.Rupke et al. 2010;Perez et al. 2011;Di Matteo et al. 2013;Moreno et al. 2019).
Although we find some trends with morphologies, we did not find clear differences in median of the morphological indicators as can be seen from Table 2.This aspect will be delayed to future work.

The impact of AGN feedback in galaxies
In this section we analyse the possible impact of AGN feedback in the relations analysed in Fig. 3.For this purpose we analyse  Edd and  BH .In the case of the MaNGA sample we only have an estimation of  BH .
In Fig. 5, the  Edd as a function of M ★ is displayed for galaxies in the three different ∇ (O/H) intervals.The upper panel also shows the dependence on sSFR while the lowest ones, on f gas .As can be seen there is a clear correlation in each of the ∇ (O/H) intervals, so Individual values (open blue squares) have been added when number of galaxies in the sample is below 5).Lower panel: Median sSFR as a function of Σ ★ for galaxies in the eagle sample with strong negative (orange triangles and lines), weak (purple crosses and lines) and strong positive (blue squares and lines) ∇ (O/H) .All error bars were estimated by using a bootstrap technique.The top small panels display the fraction of galaxies per bin Σ ★ in each of the subsmaples a function of Σ ★ .that the level of activity increases for increasing M ★ .We recall that more massive galaxies are expected to have more massive SMBHs.However, from this figure we can see that there are very few AGNs according to the adopted classification (see Section 2).Galaxies with intermediate and high masses also have lower star formation activity and lower gas fractions than low mass ones as expected.Only 6.5 per cent of the galaxies in the low-mass interval has relevant accretion activity (ADAFs) and none hosts an AGN as expected.
We have also identified galaxies that have had a recent strong mass growth event quantified by  25 < 3 Gyr as in Fig. 4. As The median trends estimated for galaxies with strong negative (orange triangles and dashed lines), weak (purple crosses and dashed lines) and strong positive (blue squares and dashed lines) metallicity gradients are also shown.The grey shaded area represents the region where galaxies are considered as AGN, while the pink shaded area represents the region were galaxies follow the ADAFs (Advection Dominated Accretion Flows) regime.Both panels contain the threshold adopted to define out stellar mass subsamples (dashdotted line).Galaxies with recent massive mergers or accretion are denoted by a black enclosing circle (i.e.t 25 < 3 Gyr).Additionally, the upper panel shows the percentage of galaxies within each accretion regime per mass bin.
noted in Fig. 5, these galaxies are preferentially low mass (5 over 10 galaxies).This supports the claim that low-mass galaxies are mainly regulated by SN feedback, probably triggered by strong accretions which induced starbursts as these systems are also gas-rich and have the higher sSFR.Interestingly, galaxies with different metallicity gradients occupy different regions in Fig. 5. Galaxies with strong negative gradients tend to be located above the median relation for galaxies with weak metallicity gradients, whereas galaxies with strong positive ∇ (O/H) prefer to be clearly below it.Most of low mass galaxies are inactive as expected.At a given stellar mass, there is a systematic increase of the  Edd from galaxies with strong positive to strong negative ∇ (O/H) .Because the star formation activity also increases in the same fashion, we speculate that the process feeding the star formation, i.e. galaxy mergers, could be also doing the same with the central BH.And, at a given  Edd , galaxies with positive gradients tend to be more massive than the rest.The correlations remains when M ★ is replaced by Σ ★ , so that more massive and compact galaxies show higher  Edd and at a given, Σ ★ , galaxies with strong positive gradients have lower  Edd .
The M BH is known to be more stable and is the result of the BH growth along the history of formation of galaxies.Hence, we analyzed it as a function of both M ★ and Σ ★ for intermediate and high mass galaxies since they are potential affected by AGN feedback.Fig. 6 displays the trends for the MaNGA (upper panels) and the eagle (lower panels) samples.We note that the eagle galaxies have more massive M BH at a given M ★ or Σ ★ .However, in this analysis we are only interested in assessing the relative trends between the relations for galaxies with different ∇ (O/H) within each sample.
At a given stellar mass, galaxies in MaNGA have larger M BH if they have strong positive gradients.Indeed, there is a systematic decrease of the black hole mass from galaxies with strong positive, weak to strong negative gradients.As a function of Σ ★ galaxies with strong negative gradients also tend to have smaller M BH at a given stellar mass and Σ ★ .Those with strong positive gradients tend to have more massive BHs, which is consistent with having already experienced AGN activity which fed the BHs and diminish the star formation activity via AGN feedback as shown in Fig. 3.However, as a function of Σ ★ , there is no significant differences between the trends of galaxies with weak and positive ∇ (O/H) .
For the eagle sample, the correlations are also present but there is no clear difference between galaxies with different gradients as can be seen from the lower panels of Fig. 3.There are only indication that galaxies with strong gradients would tend to have larger M BH than those with weak ones at a given stellar mass.This could be due to the impact of mergers and strong gas inflows (see T22) which can feed the bulge in the eagle simulations.
The results reported in this section could explain the behaviour of the eagle galaxies in the sSFR -Σ ★ relation and their apparent independence on metallicity gradients, suggesting that the AGN feedback might be efficient at diminishing the star formation for more compact galaxies, but without strongly affecting the metallicity distribution in the star forming regions across the disc.Additionally, this could indicate the need to revise the modelling of both SN and AGN feedback.

CONCLUSIONS
We analysed galaxies from the MaNGA survey and from the eagle Project to understand the connection between the metallicity gradient, the global metallicity, and the stellar mass as a function of galaxy size, star formation and AGN activity.The ability to access the evolutionary history of the eagle galaxies enabled us to investigate the information stored in the metallicity gradients of galaxies in relation to their recent history of significant accretion.The metallicity gradients were estimated for the star-forming discs within [0.5, 1.5] eff for both observations and simulations.For the observations we performed several tests by using different metallicity calibrators, finding similar global results for most of them.
Our main results are the following: • Our analysis shows that galaxies with different metallicity gradients tend to occupy different regions of the MZR at a given stellar mass.Galaxies with stellar masses lower than ∼ 10 9.75 M ⊙ and strong negative gradients tend to be more chemically enriched than those with strong positive gradients.However, for higher mass galaxies, the level of enrichment converges to similar values regardless of ∇ (O/H) .
• A clear relation between the sSFR and the stellar compactness is identified, which varies with the metallicity gradients for MaNGA galaxies.In observations and simulations, low mass galaxies are found to be actively forming stars irrespective of their compactness and metallicity gradients.However, galaxies with strong positive and negative gradients tend to be more compact, i.e., higher Σ ★ , while those with weak metallicity gradients tend to be less concentrated.
For intermediate and high mass galaxies, we found differences between the trends in MaNGA and eagle.In MaNGA, we observed a decrease in sSFR with increasing Σ ★ for galaxies with positive gradients.For galaxies with weak gradients, the sSFR decreases with increasing stellar mass, while for galaxies with strong negative gradients, the sSFR remains high regardless of Σ ★ .The strong negative gradient galaxies tend to be more disc-dominated, while the more quenched galaxies are both more concentrated and more dispersiondominated.However, simulated galaxies with M ★ ≳ 10 9.75 M ⊙ are systematically less star-forming for increasing stellar concentration regardless of the metallicity gradient.This is at odd with our observational findings.
• Galaxies in eagle with recent strong accretion tend to be more star-forming and have strong negative or positive gradients, and hence.The results of T22 show that eagle galaxies recover their weak metallicity gradients after 1.4 − 2 Gyr, on average.Hence, galaxies participating in this event could be temporally located at low or high metallicity before they recover a weak gradients and set on the main MZR.At this stage, the contribute to the scatter of the MZR and could show signals of outflows at this stage .
• On the one hand, we find that eagle galaxies with strong positive gradients tend to show weaker BH accretion activity, i.e. lower  Edd , lower sSFR and are more compact than those with weak and negative metallicity gradients.This suggests that those simulated galaxies with strong positive ∇ (O/H) might already passed their more active star-forming and probably AGN phases.Galaxies with strong negative gradients are associated with strong  Edd , larger M BH , high sSFR and larger gas fractions.Regarding M BH , galaxies with strong negative and positive gradients show a weak trend to have larger M BH than those with weaker gradients at a given stellar mass.But, the Σ ★ -M BH does not depend on the metallicity gradients.On the other hand, in the MaNGA sample, galaxies with strong negative gradients tend to have lower M BH for a given M ★ and Σ ★ for M ★ > 10 9.75 .Additionally, galaxies with strong positive gradients have systematically more massive M BH than the rest of the galaxies.
Both observations and simulations agrees that, at low masses, the sSFR is high and shows no dependence on Σ ★ or ∇ (O/H) , indicating a different star formation regulation compared to intermediate and high masses.The latter are expected to be more affected by AGN feedback compared to smaller galaxies, which are affected mainly by SN feedback.At these higher masses, in simulations the sSFR decreases systematically with increasing Σ ★ , suggesting a more efficient quenching mechanism in the most compact systems that we could link to the action of AGN feedback.The lack of dependence of these trends on ∇ (O/H) in the simulations suggests that AGN feedback model could be more efficient at quenching the star formation activity in the central regions than ejecting or mixing chemical elements in the star-forming regions across the discs (see Chen et al. 2023, for a discussion of different feedback modelling).
Our results suggest differences between observations and simulations in the metallicity recycling and mixing history of massive galaxies with different metallicity gradients, in relation to their global parameters.This discrepancies could contribute to further understanding the physical processes that take place in galaxy formation.In this section, we report some general considerations regarding our main samples, eagle and MaNGA.As described in Section 3, we re-normalised the metallicity of eagle galaxies to match the enrichment level of galaxies in the MaNGA sample.We have taken as the normalisation point the median metallicity of MaNGA galaxies at M ★ = 10 10 M ⊙ , x(O/H), manga ≃ 8.46.Compared to the measured value in eagle galaxies, x(O/H), eagle ≃ 8.58, the difference Δx (O/H) ≃ 0.126 was subtracted to all galaxies in the eagle sample.In the Figure A1, we plotted the medians trends of the eagle sample before (purple stars) and after (red crosses) the re-normalisation, along with the median trend of MaNGA galaxies (light-blue circles).As can be seen in this figure, the re-normalised eagle sample reproduces the observable MZR with good accuracy.
Having defined the re-normalised sample as our final eagle sample, in Fig. A2 we show the distribution of four key variables for the analysis carried out in this paper; gas-phase metallicity gradient ∇ (O/H) , specific star formation rate sSFR, effective radius  eff , and stellar mass M ★ .Such distributions are deployed in mass ranges, for our eagle (pink striped regions) and MaNGA (blue lines) samples.In general, we can note that the distribution of these properties agrees globally between both samples, with some exceptions.At low masses, galaxies in the eagle sample are more star-forming, tend to be slightly more extended, and have a saturation at the lower end of the mass range, while galaxies in the MaNGA sample have masses closer to the higher end of the range.At intermediate and high masses, a higher level of agreement can be seen in the distribution of both samples, with a slight trend in eagle galaxies having more positive metallicity gradients.Comparison of the trends obtained from the six additional metallicity calibrators as mentioned in Section 2. Orange, purple and blue shaded lines represent the median trends of galaxies with strong negative, weak and strong positive metallicity gradients, respectively.All error bars were estimated by using a bootstrap technique.

Figure 2 .
Figure2.Mass-metallicity relation for the MaNGA (left panel) and eagle (right panel) samples coloured by sSFR (upper panels) and  eff (lower panels).To reveal the correlation between the properties, the Locally Weighted Scatterplot Smoothing (LOESS;Cappellari et al. 2013) technique was applied in all the panels.For low-mass galaxies (the dash-dotted black line denotes M ★ = 10 9.75 M ⊙ ) a correlation angle  was computed using PCC analysis (Partial Correlation Coefficients,Bluck et al. 2020b).This angle denotes the direction in which the third dependency of each panel grows.eagle galaxies follow the observational trends for both parameters, with a global trend for galaxies with high metallicity to have lower star-forming and have smaller  eff for M ★ ≲ 10 10 M ⊙ , at a given stellar mass.Higher masses tend to have lower sSFR and to be more extended.The errors on  were estimated by applying a bootstrap technique.

Figure 3 .
Figure 3. Upper panel: Median sSFR as a function of Σ ★ for galaxies in the MaNGA with strong negative (orange triangles and lines), weak (purple crosses and lines) and strong positive (blue, filled, squares and lines) ∇ (O/H) .Individual values (open blue squares) have been added when number of galaxies in the sample is below 5).Lower panel: Median sSFR as a function of Σ ★ for galaxies in the eagle sample with strong negative (orange triangles and lines), weak (purple crosses and lines) and strong positive (blue squares and lines) ∇ (O/H) .All error bars were estimated by using a bootstrap technique.The top small panels display the fraction of galaxies per bin Σ ★ in each of the subsmaples a function of Σ ★ .

Figure 4 .
Figure 4. Mass-metallicity relation coloured by specific star formation rate, sSFR (upper panel), effective radius,  eff (second panel), disc-to-total stellar mass ratio, D/T (third panel) and gas fraction,  gas (lower panel).To reveal the correlation between the properties, the Locally Weighted Scatterplot Smoothing (LOESS; Cappellari et al. 2013) technique was applied in all the panels.Galaxies with t 25 lower than 3 Gyr have been highlighted with a black circle.The dotted-dashed line denotes M ★ = 10 9.75 M ⊙ .

Figure 5 .
Figure5.Eddington ratios,  Edd , as a function of M ★ , coloured by the sSFR (upper panel) and the gas fraction,  gas (lower panel).The median trends estimated for galaxies with strong negative (orange triangles and dashed lines), weak (purple crosses and dashed lines) and strong positive (blue squares and dashed lines) metallicity gradients are also shown.The grey shaded area represents the region where galaxies are considered as AGN, while the pink shaded area represents the region were galaxies follow the ADAFs (Advection Dominated Accretion Flows) regime.Both panels contain the threshold adopted to define out stellar mass subsamples (dashdotted line).Galaxies with recent massive mergers or accretion are denoted by a black enclosing circle (i.e.t 25 < 3 Gyr).Additionally, the upper panel shows the percentage of galaxies within each accretion regime per mass bin.

Figure 6 .
Figure6.M BH as a function of M ★ and Σ ★ for galaxies in the MaNGA (upper panels) and eagle sample (lower panels) separated according to their metallicity gradients as displayed in the inset labels.

Figure A1 .
Figure A1.MZR for MaNGA (light-blue circles) and our original (purple stars) and normalized (red crosses) eagle samples.Medians are plotted with the raw data on the background.The black vertical line shows the limit M ★ = 10 10 M ⊙ , which is the matching stellar mass between MaNGA and normalised eagle medians.

Figure A2 .
Figure A2.Comparative distribution of four main galactic properties, metallicity gradient ∇ (O/H) , specific star formation rate, sSFR, effective radius,  eff , and stellar mass, M ★ , between our MaNGA (blue, solid lines) and eagle (pink, shaded areas) samples, displayed for each of the three defined stellar mass intervals as indicated by the labels in the figure.

Figure B1 .
Figure B1.Dependence of MZR on metallicity gradients in our MaNGA sample.Comparison of the trends obtained from the six additional metallicity calibrators as mentioned in Section 2. Orange, purple and blue shaded lines represent the median trends of galaxies with strong negative, weak and strong positive metallicity gradients, respectively.All error bars were estimated by using a bootstrap technique.

Table 2 .
Summary of the main properties in our MaNGA and eagle subsamples.