On the origin of the mass-metallicity gradient relation in the local Universe

Almost all the analytic models that reproduce the observed MZR do not have spatial information of the distribution of metallicities in a galaxy – these are typically developed to study global metallicities in galaxies. Although the primary focus of our work is to explain metallicity gradients by making use of the spatial information of metallicity, our model also reproduces the MZR as a proof of concept.

To produce an MZR from the model, we need an estimate of the mean metallicity in galaxies as a function of $M_{\star}$. For this purpose, we use the SFR-weighted mean metallicity given by Sharda et al. (2021b, equation 46)

where $\dot{s}_{\star}(x)=1/x^{2}$ is the radial distribution of star formation per unit area (Krumholz et al., 2018). We use the SFR-weighted $\overline{\mathcal{Z}}$, because it can be directly compared against available MZRs since they are inherently sensitive to the SFR as the nebular metallicities are measured in H ii regions around young stars (Zahid et al., 2014). Additionally, semi-analytic models and simulations too use SFR-weighted metallicities to construct MZRs (e.g., Tissera et al. 2019; Torrey et al. 2019; Forbes et al. 2019; Yates et al. 2021).

In order to derive results in terms of $M_{\star}$, we treat the rotational velocity, $v_{\phi}$, as the primary quantity that we vary. For each $v_{\phi}$, we can estimate the corresponding halo mass $M_{\rm{h}}$ and halo accretion rate $\dot{M}_{\rm h}$ at $z=0$ (Sharda et al., 2021b, equations 34–35). We convert the halo mass to $M_{\star}$ following the $M_{\mathrm{h}}-M_{\star}$ relation from Moster et al. (2013) for the local Universe. Following Sharda et al. (2021b), we keep the yield reduction factor, $\phi_{y}$, as a free parameter and vary it between 0.1 and 1, though we note that, based on both theory and observations, $\phi_{y}$ is expected to be close to unity in massive galaxies. For all other parameters, in particular the velocity dispersion $\sigma_{g}$, we use the fiducial values listed in Sharda et al. (2021b, Tables 1 and 2). Specifically, we use local dwarf values for galaxies with $M_{\star}\leq 10^{9}\,\rm{M_{\odot}}$, and local spiral values for $M_{\star}\geq 10^{10.5}\,\rm{M_{\odot}}$. For intermediate stellar masses, we linearly interpolate in $\log_{10}M_{\star}$ between these two limits for all parameters. For example, the velocity dispersions we adopt for spirals and dwarfs are $10\,\rm{km\,s^{-1}}$ and $7\,\rm{km\,s^{-1}}$ respectively, so we adopt $\sigma_{g}=\left(2\log_{10}M_{\star}/\mathrm{M_{\odot}}-11\right)\,\mathrm{km\,s^{-1}}$ for intermediate-mass galaxies with $10^{9}\,\mathrm{M_{\odot}}<M_{\star}<10^{10.5}\,\mathrm{M_{\odot}}$. We have verified that the resulting MZR and MZGR are not particularly sensitive to the choice of the $M_{\star}$ boundaries invoked to classify dwarfs and spirals; we also discuss this further in Section 4. We set $\mathcal{Z}_{r0}$ to its equilibrium value (Sharda et al., 2021b), and set the circumgalactic medium metallicity to $\mathcal{Z}_{\rm{CGM}}=0.2$ for all galaxies²²2This is slightly lower than the median $\mathcal{Z}_{\rm{CGM}}=0.3$ found by Prochaska et al. (2017) for $z\sim 0.2$ galaxies (see also, Wotta et al. 2016, where the authors find a bimodal distribution of $\mathcal{Z}_{\rm{CGM}}$); however, these surveys do not cover the entire range in galaxy masses we are interested in, and we expect $\mathcal{Z}_{\rm{CGM}}$ to be lower in low mass galaxies. In any case, this difference does not have a significant effect on the MZGR., which sets $c_{1}$. The MZR (as well as the MZGR discussed below) is insensitive to $\mathcal{Z}_{r_{0}}$ and only weakly sensitive to $\mathcal{Z}_{\rm{CGM}}$ as compared to $\phi_{y}$, so we do not vary $\mathcal{Z}_{\rm GCM}$ separately. Finally, we follow van der Wel et al. (2014) to estimate $r_{\rm{e}}$ as a function of $M_{\star}$, and set $x_{\mathrm{min}}=0.5\,r_{\mathrm{e}}$ and $x_{\mathrm{max}}=3\,r_{\mathrm{e}}$ as the range of radii $x$ over which our model solution applies. This range of radii roughly mimics that over which metallicities are measured.

Figure 1 shows the resulting MZR from our model, color-coded by the ratio $\mathcal{P}/\mathcal{A}$ that describes the relative strength of advection to cosmic accretion. We remind the reader that both $\mathcal{P}$ and $\mathcal{A}$ (as well as $\mathcal{S}$) are normalised by diffusion in the model. The vertical spread in the model MZR is a result of varying $\phi_{y}$. We also overplot the parameter space of observed MZRs from several other works based on the direct $T_{\rm{e}}$ method (Pettini & Pagel, 2004; Andrews & Martini, 2013; Curti et al., 2017, 2020) and photoionization modeling (Kewley & Dopita, 2002; Tremonti et al., 2004; Mannucci et al., 2010), all of which we adopt from Maiolino & Mannucci (2019, Figure 15). We see that the model is able to reproduce the MZR of the local Universe albeit with a large spread due to $\phi_{y}$. There are several factors behind quantitative differences between the model MZR and MZRs in the literature. From the perspective of the model, these differences are attributed to the choice of the metal yield $y$, excluding the galaxy nucleus while finding the mean metallicities, and the absolute size of the galaxy disc. From the perspective of the MZRs we compare the model with, these differences are due to calibration and observational uncertainties, as well as limited coverage of the galaxy discs.

In order to match with the measured MZRs, the model prefers higher $\phi_{y}$ for massive galaxies and lower $\phi_{y}$ for low-mass galaxies. This implies that metals are well-mixed in the ISM in massive galaxies before they are ejected through outflows, whereas in dwarf galaxies, some fraction of metals are ejected directly before they can mix in the ISM; in other words, the best match between the model MZR and the literature MZRs predicts that dwarf galaxies have more metal-enriched winds than massive galaxies. This finding is not new and has been theorized in several works (e.g., Larson, 1974; Dekel & Silk, 1986; Dalcanton, 2007; Finlator & Davé, 2008; Lilly et al., 2013; Dayal et al., 2013; Forbes et al., 2019), simulations (Creasey et al., 2015; Ma et al., 2016; Christensen et al., 2018; Emerick et al., 2018; Emerick et al., 2019), and also has some observational evidence (Martin et al., 2002; Chisholm et al., 2018).

To further treat the question of how $\phi_{y}$ scales with $M_{\star}$ quantitatively, we also plot two models for this scaling. We obtain the first of these from available observations that directly constrain the ratio of wind metallicity to ISM metallicity (Chisholm et al., 2018), and the second simply by forcing the model to reproduce the observed MZR provided by Curti et al. (2020). Appendix A describes how we obtain these scalings (and the associated uncertainties) in detail. While the shape of the first scaling is consistent with observed MZRs, the second is almost identical to the direct $T_{e}$ based MZRs by construction; we include the second scaling nonetheless because there is no guarantee that the scaling we have enforced to produce the MZR will also yield the correct MZGR, a question we explore below.

Figure 1 shows that the MZR bends roughly where the ratio $\mathcal{P}/\mathcal{A}$ passes through unity. We can understand this behaviour as follows: the total metallicity is set by a competition between metal production (the term $\mathcal{S}$) and dilution by metal-poor gas, which can be supplied either by direct cosmological accretion onto the disc ($\mathcal{A}$) or advection of gas from the weakly-star-forming outskirts to the more rapidly-star-forming centre ($\mathcal{P}$). Each of these terms varies differently with rotation curve velocity $v_{\phi}$, which in turn correlates with stellar mass; as shown in Sharda et al. (2021b), $\mathcal{P}$ is independent of $v_{\phi}$,³³3Recall that each of these terms is expressed as the relative importance of a particular process compared to metal diffusion; thus, $\mathcal{P}\propto v^{0}_{\phi}$ does not mean that advection is equally rapid in all galaxies independent of stellar mass, just that the ratio of advection to diffusion does not explicitly depend on stellar mass. while $\mathcal{S}\propto v^{2}_{\phi}$ and $\mathcal{A}\propto v^{3.3}_{\phi}$. In the low-mass regime, corresponding to small $v_{\phi}$, we have $\mathcal{P}>\mathcal{A}$, implying that the metallicities are primarily set by the balance between source and advection. Since $\mathcal{P}\propto v^{0}_{\phi}$ and $\mathcal{S}\propto v^{2}_{\phi}$, as we go to smaller $M_{\star}$ and $v_{\phi}$, the equilibrium metallicity drops because of lower $v_{\phi}$ and lower $\phi_{y}$ as compared to massive galaxies. On the contrary, in the high-mass regime $\mathcal{A}>\mathcal{P}$, implying that the metallicities are set by the balance between $\mathcal{A}$ and $\mathcal{S}$. Since $\mathcal{A}\propto v^{3.3}_{\phi}$, which is stronger than the dependence of $\mathcal{S}$ on $v_{\phi}$, the metallicity, which is proportional to $\mathcal{S}/\mathcal{A}$, ceases to rise with $M_{\star}$, and instead reaches a maximum and starts to decrease. However, the decrease is rather mild, because shortly after passing the value of $v_{\phi}$ where we move into the $\mathcal{A}>\mathcal{P}$ regime, galaxies become so massive that they cease to be star-forming altogether. Thus, among star-forming galaxies, the trend of $\overline{\mathcal{Z}}$ with $M_{\star}$ is simply that $\overline{\mathcal{Z}}$ ceases to increase and reaches a plateau. For less massive galaxies the dominant source of metal-poor gas is advection rather than accretion. However, this only holds as long as advection is non-zero; for low mass galaxies where there is no advection (i.e., there is no turbulence due to gravity), it falls upon cosmic accretion to balance metal production. Since cosmic accretion is much weaker in low mass galaxies, it can take a long time for this balance to approach a steady-state, which can push the gradients out of equilibrium (Sharda et al., 2021b, Section 5.1).

The existence of a local gas phase MZR has been known since early analysis of data from the Sloan Digital Sky Survey (SDSS, Tremonti et al., 2004), although the absolute normalisation of the MZR remains an unsolved issue due to systematic calibration uncertainties (Kewley & Ellison, 2008; Pilyugin & Grebel, 2016; Brown et al., 2016; Curti et al., 2017; Barrera-Ballesteros et al., 2017; Teimoorinia et al., 2021). Despite these uncertainties, however, it is clear both that a relationship exists, and that it has a characteristic mass scale of $\sim 10^{10.5}$ M_⊙ at which the curvature of the relation changes (Blanc et al., 2019). Not surprisingly, there have been numerous attempts to explain these relations theoretically, and it is interesting to put our model in the context of these works. However, we caution that what follows is only a partial discussion of the (vast) literature on this topic, and refer readers to the comprehensive review by Maiolino & Mannucci (2019, Section 5.1).

The basic result from theoretical models to date is that galaxies tend to approach equilibrium between inflows, accretion, star formation and outflows, which naturally gives rise to the observed MZR (Finlator & Davé, 2008; Davé et al., 2012; Lilly et al., 2013; Dayal et al., 2013; Forbes et al., 2014). Our results are broadly consistent with this picture. However, there are some subtle differences among published models, and between existing models and ours. One important point of distinction is the extent to which outflows are metal-enriched relative to the ISM (i.e., $\phi_{y}<1$ in the language our model), and whether this enrichment varies as a function of galaxy mass or other properties (as is the case for our two possible scalings). As already discussed, many authors simply assume that outflows are not metal-enriched (i.e., the outflow metallicity is the same as the ISM metallicity, $\phi_{y}=1$ in our notation; e.g., Finlator & Davé 2008; Davé et al. 2012; Schaye et al. 2015; Hirschmann et al. 2016; Davé et al. 2017; Collacchioni et al. 2018; De Lucia et al. 2020), and produce an MZR based on this assumption. Others explicitly contemplate values of $\phi_{y}<1$ (e.g., Dalcanton, 2007; Spitoni et al., 2010; Peeples & Shankar, 2011; Lu et al., 2015; Forbes et al., 2014; Forbes et al., 2019; Yates et al., 2021; Kudritzki et al., 2021). Our conclusion that reproducing the full shape of the MZR requires $\phi_{y}<1$, particularly in low-mass galaxies, is consistent with the findings of the latter group of investigators. However, many of these authors do not study the relative importance of metal-enriched outflows for dwarfs versus spirals, which we find to be important.

It is also debated whether the MZR really has a curvature at intermediate stellar masses, and if it does, whether it simply flattens out or starts to bend. While some simulations do find curvature in the MZR around $10^{10}-10^{10.5}\,\rm{M_{\odot}}$ (e.g., Davé et al., 2017; Torrey et al., 2019), others do not (e.g., Torrey et al., 2014; De Rossi et al., 2015; Ma et al., 2016). Our model is consistent with the former, especially if we look at the empirical scalings of $\phi_{y}$ with $M_{\star}$. Moreover, recent results also show that the curvature is physical and persists in the data even after observational uncertainties are accounted for (Blanc et al., 2019). However, the cause behind the curvature is not completely understood, and factors like Active Galactic Nuclei (AGN) feedback (De Rossi et al., 2017), gas recycling (Brook et al., 2014), effective gas fraction (Torrey et al., 2019), chemical saturation in the ISM of massive galaxies (Zahid et al., 2013), and a transition in galaxy regimes together with metal-enriched outflows as we show in this work can all play a role.

In addition to the models above, to which our results are directly comparable, a number of authors have studied the dependence of the MZR on factors not included in our work, like downsizing, time-dependent outflows, variations in star formation efficiencies and IMF, presence of satellites, environmental effects, etc. (e.g., Köppen et al., 2007; Cooper et al., 2008; Maiolino et al., 2008; Calura et al., 2009; Spitoni et al., 2010; Bouché et al., 2010; Hughes et al., 2013; Peng & Maiolino, 2014; Genel, 2016; Wu et al., 2017; Bahé et al., 2017; Lian et al., 2018a, b). However, unlike the current work, most models only study the MZR and not the MZGR, thus it is difficult to reconcile whether their conclusions hold or are self-consistent with spatially-resolved galaxy properties.

We use the same metallicity distributions described in Section 3 to compute metallicity gradients. To be consistent with the procedure most commonly used in analysing observations, we obtain the gradient by performing a linear fit to $\log_{10}\,\mathcal{Z}$ from $0.5-2.5\,r_{\rm{e}}$ (e.g. Sánchez et al., 2012, 2014; Sánchez-Menguiano et al., 2016; Poetrodjojo et al., 2018).⁴⁴4To be consistent with observations, we only utilize metallicities till $2.5\,r_{\mathrm{e}}$ to measure the gradients, as opposed to $3\,r_{\mathrm{e}}$ that we use to measure $\overline{\mathcal{Z}}$. Following the discussion on inverted gradients in Sharda et al. (2021b, Section 5.2.3) and the uncertainty around them being in equilibrium, we restrict the model to produce only flat or negative gradients for the purposes of studying the MZGR. Figure 2 shows the MZGR from our model, again color-coded by the ratio of advection to accretion ($\mathcal{P}/\mathcal{A}$). The top and the bottom panels show the metallicity gradients in $\rm{dex\,kpc^{-1}}$ and $\mathrm{dex}\,r^{-1}_{\rm{e}}$ units, respectively. The spread, as for the MZR, is a result of $\phi_{y}$. The transition from the advection-dominated to the accretion-dominated regime, as in the MZR, is also visible in the MZGR. When the gradients are measured in $\mathrm{dex\,kpc^{-1}}$, this transition corresponds to the slight curvature in the MZGR that appears around $M_{\star}\sim 10^{10}-10^{10.5}\,\rm{M_{\odot}}$ (top panel in Figure 2). When they are measured in $\mathrm{dex}\,r^{-1}_{\mathrm{e}}$, it corresponds to the somewhat sharper curvature around the same stellar mass (bottom panel in Figure 2). This finding is strong evidence for the links between the MZR and the MZGR, and also reveals that it is the same underlying physical mechanism that controls the shape of both.

While the stellar mass of the accretion-advection transition influences the location at which our model curves bend, it is not the only factor that does so. The precise location of the bend is also sensitive to parameters like $\mathcal{Z}_{\rm{CGM}}$ and $\phi_{y}$, and both of the MZGR bend and the mass where $\mathcal{P}/\mathcal{A}=1$ depend weakly on the limits in $M_{\star}$ we select for smoothly interpolating between the dwarf and spiral regimes: for example, if we lower the threshold for spirals from $10^{10.5}\,\rm{M_{\odot}}$ to $10^{10}\,\rm{M_{\odot}}$, both shift to lower stellar mass. Similarly, if we increase the threshold for dwarfs from $10^{9}\,\rm{M_{\odot}}$ to $10^{9.5}\,\rm{M_{\odot}}$, both shift to higher stellar mass. However, irrespective of the interpolation limits in $M_{\star}$, both the curvature of the MZGR and the transition from $\mathcal{P}>\mathcal{A}$ to $\mathcal{P}<\mathcal{A}$ are always present. The existence of these features is a robust prediction of the model independent of uncertain parameter choices.

The physical origin for the behaviour of the MZGR is also the same as for the MZR: gradients are at their steepest when both of the processes for smoothing them – accretion, $\mathcal{A}$, and inward advection of gas, $\mathcal{P}$, are at their weakest compared to metal production, $\mathcal{S}$. Diffusion also helps smooth gradients, but is always subdominant compared to either accretion or advection, as evidenced by the fact that we never have $\mathcal{P}<1$ and $\mathcal{A}<1$ simultaneously. The point where advection and accretion are weakest is roughly where galaxies are transitioning from being advection-dominated, $\mathcal{P}>\mathcal{A}$, to accretion-dominated, $\mathcal{P}<\mathcal{A}$. We emphasise that, while the exact stellar mass at which this transition occurs can be somewhat sensitive to choices of model parameters (for example, the Toomre $Q$ of galactic discs), its existence is not; the bends in the coloured bands in Figure 2 that describe our model always occur irrespective of our parameter choices. Additionally, note that the minimum of the model MZGR is not always coincident with $\mathcal{P}/\mathcal{A}=1$; the position of the minimum is dependent on the model parameters, in particular, $\phi_{y}$.

In Figure 2 we also plot MZGRs from the MaNGA (Belfiore et al., 2017), CALIFA (Sánchez et al., 2014; Sánchez-Menguiano et al., 2016) and SAMI (Poetrodjojo et al., 2018, 2021) surveys, homogenized and corrected for spatial resolution by Acharyya et al. (2021). We adopt the $\rm{dex\,kpc^{-1}}$ values from Acharyya et al. (2021), and convert to $\mathrm{dex}\,r^{-1}_{\rm{e}}$ following the $r_{\mathrm{e}}$–$M_{\star}$ scaling relations from van der Wel et al. (2014) to be consistent with our assumptions elsewhere⁵⁵5The qualitative trend of the MZGR remains the same for the $\mathrm{dex}\,r^{-1}_{\rm{e}}$ gradients reported by Acharyya et al. (2021) as compared to the ones shown in the bottom panel of Figure 2 using the scaling relation between $r_{\mathrm{e}}$–$M_{\star}$, with a change in the overall normalisation of the metallicity. We have also verified that the $r_{\rm{e}}$ we find from van der Wel et al. (2014) is in very good agreement with that measured in, for example, the SAMI sample we use.. We also overplot results from MaNGA based on three different metallicity calibrations by Mingozzi et al. (2020): Pettini & Pagel (2004, PP04), Maiolino et al. (2008, M08), and Blanc et al. (2015, IZI).

The first thing to notice is that the qualitative trend found in the data is in good agreement with that predicted by our model: gradients are steepest at $M_{\star}\sim 10^{10}-10^{10.5}$ M_⊙, and flatten at both lower and higher masses. However, the location of the curvature in the data and the model differ by as much as $0.5-1\,\rm{dex}$ in stellar mass. This is not surprising given the uncertainties in the parameters that affect the curvature, as discussed above (e.g., interpolation limits in $M_{\star}$, our constant adopted value of $\mathcal{Z}_{\rm{CGM}}$, and the scaling of $\phi_{y}$ with $M_{\star}$). Moreover, it is important to recall that the data themselves are not fully secure, due to uncertainties caused by the choice of metallicity diagnostic; Poetrodjojo et al. (2021, their Figure 11) show that the exact mass at which the MZGR bends depends on which diagnostic is used to determine the metallicity, and that these variations are reduced but still persist even after the diagnostics are homogenised. Thus, it is presently difficult to accurately determine the location of the curvature, especially given its mildness. Nonetheless, the presence of a bend seems to be robust in the data, as it is in our model.

Second, we see that similar to the MZR, this comparison of the model to the observed MZGR reveals that low-mass galaxies prefer low $\phi_{y}$. However, the spread due to $\phi_{y}$ in the MZGR at the high-mass end is quite narrow; thus, gradients in massive galaxies are not particularly sensitive to $\phi_{y}$, although the data suggests higher $\phi_{y}$ for the MZGR in massive galaxies (note the inverted arrows for $\phi_{y}$ on Figure 2 as compared to Figure 1). Our findings on $\phi_{y}$ being ineffective at setting gradients in massive galaxies is consistent with earlier works (e.g., Fu et al., 2013). However, our proposed explanation for the flattening of gradients in massive galaxies based on the advection-to-accretion transition differs from these studies that attributed the observed flattening to saturation of ISM metallicities (Phillipps & Edmunds, 1991; Mollá et al., 2017), radially-varying star formation efficiency (Belfiore et al., 2019), or past mergers (Rupke et al., 2010; Perez et al., 2011; Fu et al., 2013).

In Figure 2, we also plot model predictions using the two scalings of $\phi_{y}$ with $M_{\star}$ that we described in Section 3. These scalings are able to reproduce the high mass end of the MZGR, and yield a qualitative trend similar to that seen in the data, but quantitatively the predicted gradients from the scalings are steeper than that observed at the low mass end. In retrospect, this is not entirely unexpected given the uncertainties in the two approaches, and the fact that these scalings are sensitive to the absolute metallicity (see Appendix A). Judging from Figure 2, we slightly prefer scaling 2, since it is closer to the observations at intermediate stellar masses; we revisit the comparison between the two scalings in Section 5. Nevertheless, the fact that both the MZR and the MZGR suggest a qualitatively similar scaling between $\phi_{y}$ and $M_{\star}$ is an encouraging sign of consistency. However, it is difficult to derive quantitative similarities given the uncertainties in these empirical scalings.

Only a handful of models exist in the literature that focus on gas phase metallicity gradients rather than global metallicities (Mott et al., 2013; Jones et al., 2013; Ho et al., 2015; Carton et al., 2015; Kudritzki et al., 2015; Pezzulli & Fraternali, 2016; Schönrich & McMillan, 2017; Kang et al., 2021), and even fewer that actually study the local MZGR or its equivalent (Lian et al., 2018b; Lian et al., 2019; Belfiore et al., 2019). Of these, the models by Lian et al. (2018b) and Belfiore et al. (2019) are closest in spirit to ours.⁶⁶6Lian et al. (2019) focus only on low-mass satellites, so our results are not easily comparable. Quantitative comparison between our results and those of Lian et al. is challenging, because they do not quote measurements in $\rm{dex/kpc}$ or equivalent. Examining their plots, it seems that they also find slightly steeper gradients for intermediate mass galaxies, consistent with our findings. Similarly, Belfiore et al. find that observed gradients in local dwarfs and spirals are best reproduced by a model where the star formation timescale at each radius is proportional to the local orbital period. For massive galaxies, this scaling is quite similar to that in the Krumholz et al. (2018) galaxy model that is embedded in our metallicity model, and thus at first glance is also consistent with our findings. However, there remain substantial differences between our model and those of Lian et al. and Belfiore et al.. Neither of these studies include the effects of radial inflow or metal diffusion. Neither adopt our approach of systematically varying the highly-uncertain yield reduction factor $\phi_{y}$: Lian et al. adopt a parameterised, time-dependent functional form that they tune in order to match stellar and gas metallicity gradient data, while Belfiore et al. assume that the ISM and outflow metallicities are equal ($\phi_{y}=1$ in our terminology), contrary to our findings and inconsistent with the available observational evidence (Martin et al., 2002; Schmidt et al., 2016; Chisholm et al., 2018; Telford et al., 2019; Kreckel et al., 2020). Finally, both sets of authors explicitly fit their model parameters to the data, whereas we do not except while introducing the second scaling in $\phi_{y}$. Thus, it is unclear to what extent the agreement between the models is simply a matter of their being enough adjustable parameters to make them behave similarly.

In addition to analytic models, semi-analytic models like L-Galaxies 2020 have also investigated the local MZGR, finding somewhat flatter gradients for massive galaxies as compared to low mass galaxies (Yates et al., 2021). The authors attribute their findings to inside-out star formation that increases the gas phase metallicity in the inner disc in massive galaxies. In the outer disc in these galaxies, Yates et al. either find metal-rich accretion from the CGM that enhances the metallicity (their ‘modified’ model), or metal-poor accretion that dilutes the metallicity at every radius (their ‘default’ model). The combined effect is to produce flatter metallicity profiles in massive galaxies in each case. They further conclude that flattening of the metallicity profiles in massive local galaxies is expected regardless of the mass-loading factors of outflows. Thus, their explanations for the trends seen in the local MZGR are consistent with the findings of our model. It is worth noting that while working with an earlier version of L-Galaxies, Fu et al. (2013) found relative metal enrichment of outflows to be more important than advection in driving gas phase metallicities. These authors also find a trend in the MZGR consistent with Yates et al. and ours.

It is also helpful to compare our results to simulations that have studied the local MZGR. For example, both Tissera et al. (2016) and Ma et al. (2017) find slightly flatter gradients for massive local galaxies in their simulations, consistent with our model and available observations. The EAGLE simulations (Schaye et al., 2015) find that metallicity gradients in their simulated galaxies are systematically shallower at $z=0$ than those observed in the local Universe due to high star formation efficiency at all radii (Tissera et al., 2019, Figure 11). As a result, the MZGR predicted from their simulations does not show any clear trends with the stellar mass. On the other hand, the local MZGR produced by the IllustrisTNG50 simulations (Pillepich et al., 2019; Nelson et al., 2019) is in very good quantitative agreement with that produced by our model, both in terms of the mean gradient and the scatter in gradients at a given stellar mass (Hemler et al., 2021, Figure 8). These authors suspect that gradients flatten in massive local galaxies due to AGN feedback and increasing galaxy size. While the latter of the two is consistent with the findings of Sharda et al. (2021b), the primary driver of flatter gradients in massive galaxies in our model is due to the increasing role of metallicity dilution by cosmic accretion.