The interplay between feedback, accretion, transport and winds in setting gas-phase metal distribution in galaxies

The recent decade has seen an exponential growth in spatially-resolved metallicity measurements in the interstellar medium (ISM) of galaxies. To first order, these measurements are characterised by the slope of the radial metallicity profile, known as the metallicity gradient. In this work, we model the relative role of star formation feedback, gas transport, cosmic gas accretion, and galactic winds in driving radial metallicity profiles and setting the mass-metallicity gradient relation (MZGR). We include a comprehensive treatment of these processes by including them as sources that supply mass, metals, and energy to marginally unstable galactic discs in pressure and energy balance. We show that both feedback and accretion that can drive turbulence and enhance metal-mixing via diffusion are crucial to reproduce the observed MZGR in local galaxies. Metal transport also contributes to setting metallicity profiles, but it is sensitive to the strength of radial gas flows in galaxies. While the mass loading of galactic winds is important to reproduce the mass metallicity relation (MZR), we find that metal mass loading is more important to reproducing the MZGR. Specifically, our model predicts preferential metal enrichment of galactic winds in low-mass galaxies. This conclusion is robust against our adopted scaling of the wind mass-loading factor, uncertainties in measured wind metallicities, and systematics due to metallicity calibrations. Overall, we find that at $z \sim 0$, galactic winds and metal transport are more important in setting metallicity gradients in low-mass galaxies whereas star formation feedback and gas accretion dominate setting metallicity gradients in massive galaxies.


INTRODUCTION
Metals act as natural tracers of galaxy evolution, and provide some of the most important constraints on the physical properties of baryons in galaxies (Kobayashi et al. 2020).This has resulted in the development of several galactic chemical evolution models over the last three decades to reproduce, fit, or explain trends in the metal ★ sharda@strw.leidenuniv.nl(PS) † omry.ginzburg@mail.huji.ac.il (OG) ‡ mark.krumholz@anu.edu.au(MRK) content of galaxies.Most of the literature thus far has focused on understanding the physics of galaxy-integrated (or global) metallicities, and fundamental trends such as the mass-metallicity relation (MZR, Tremonti et al. 2004;Zahid et al. 2014;Sanders et al. 2021).However, thanks to integral field unit (IFU) spectroscopy, the last decade has seen immense advancements in mapping the spatiallyresolved metal content in galaxies.It is now clear that, at least at  = 0, galaxies typically exhibit a gradient in metallicity whereby their inner parts are more metal-rich than their outskirts.Such a negative metallicity gradient is a natural consequence of inside-out galaxy formation (e.g., Shaver et al. 1983).This simplified picture is complicated by various factors internal and external to galaxies, as is evident from the vast diversity of metallicity gradients discovered in numerous galaxies across cosmic time (see reviews by Kewley et al. 2019;Maiolino & Mannucci 2019;Sánchez 2020).High resolution data from IFU spectroscopy has also revealed the existence of resolved versions of the MZR: the mass-metallicity gradient relation (MZGR), and the resolved massmetallicity relation (rMZR).The former describes the evolution of metallicity gradients with stellar mass (e.g., Sánchez et al. 2014;Belfiore et al. 2017;Sánchez-Menguiano et al. 2018;Poetrodjojo et al. 2018;Mingozzi et al. 2020;Sharda et al. 2021b;Franchetto et al. 2021) whereas the latter describes the evolution of local metallicity with the local stellar surface density (e.g., Rosales-Ortega et al. 2012;Sánchez et al. 2013;Barrera-Ballesteros et al. 2016;Trayford & Schaye 2019).
IFU observations have started to reveal the detailed 2D structure of metal distribution in galaxies, finding evidence not just for radial, but also azimuthal variations (e.g., Ho et al. 2017Ho et al. , 2018;;Kreckel et al. 2019Kreckel et al. , 2020;;Li et al. 2021Li et al. , 2023)).However, progress in theoretical work on modeling these 1 st and 2 nd order metallicity variations for a wide range of galaxies has remained rather limited (Kudritzki et al. 2015;Ho et al. 2015;Lian et al. 2018;Krumholz & Ting 2018;Belfiore et al. 2019;Sharda et al. 2021a;Metha et al. 2021).A key challenge in modeling resolved metal distributions is that they are impacted by both local and global processes in galaxies (Baker et al. 2023;see, however, Boardman et al. 2022).Additionally, while the global metallicity is only sensitive to the total metal content of galaxies, the spatially-resolved metallicity structure is also affected by transport processes within the ISM, which redistribute metals from the time they are produced.Thus, metal dynamics are as important as metal production and ejection in setting spatially-resolved metal distributions.We sub-divide and introduce these two categories below (see also, Figure 1).

Gas and metal flows in/out of galaxies
To understand the physics of spatially-resolved metal distribution and metallicity gradients, it is helpful to first explore the life cycle of baryons in galaxies.Gas regulator models paint a simple yet powerful picture of this life cycle: gas accretes onto the galaxy from the cosmic web in the form of cold or hot streams, or through mergers and interactions (e.g., Fardal et al. 2001;Dekel & Birnboim 2006;Bournaud & Elmegreen 2009;Dekel et al. 2009a;Rupke et al. 2010).The accreted gas is then transported across the galactic disc via angular momentum exchange due to viscous torques or encounters with giant molecular clouds (GMC) and bars or spiral arms (e.g., Pfenniger & Norman 1990;Block et al. 2002;Bournaud et al. 2005;Dekel et al. 2009b;Krumholz & Burkert 2010;Hopkins & Quataert 2011).The accreted gas also acts as a fuel for star formation, and dictates subsequent galactic life cycles (e.g., Kereš et al. 2005;Somerville & Davé 2015).Star formation leads to metal production dictated by the mass distribution of massive stars following a stellar initial mass function (IMF).Some of these metals are locked in low mass stars that live for a long time, others are returned to the ISM on relatively shorter timescales (Tinsley 1980), or are depleted onto dust (e.g., Konstantopoulou et al. 2022).Star formation and active galactic nuclei (AGN) feedback drives winds that eject mass and metals out of the galaxy (e.g., Veilleux et al. 2005;Muratov et al. 2015;Pillepich et al. 2018a;Mitchell et al. 2020a;Veilleux et al. 2020;Pandya et al. 2021).A fraction of the ejected mass makes its way back onto the disc in the form of fountains (e.g., Shapiro & Field 1976;Fraternali & Binney 2008;Grand et al. 2019;Li & Tonnesen 2020), whereas the rest enriches the circumgalactic medium (CGM) with metals (Tumlinson et al. 2017).Sufficiently strong winds can also deposit metals beyond the virial radii of galaxies into the intergalactic medium (IGM; Origlia et al. 2004;Ranalli et al. 2008;Tumlinson et al. 2011).Thus, the life cycle of metals is quite dynamic, and it is now known that all these different physical processes alter the metal content of galaxies (e.g., Dalcanton 2007;Finlator & Davé 2008;Lilly et al. 2013;Grand et al. 2019;Sharda et al. 2021a;Wang & Lilly 2022b;Greener et al. 2022).

Metal dynamics within galaxies
We define metal dynamics as processes that alter the spatial distribution of metals within galaxies without altering the total metal mass.Broadly speaking, we can further classify metal dynamics into categories: (1.) transport of metals with the bulk flow of the gas, and (2.) redistribution of metals.The former results in advection of metals, and is primarily driven by radial gas flows generated due to gravitational instabilities (Dekel et al. 2009b;Krumholz & Burkert 2010;Forbes et al. 2012), magnetic stresses (Wang & Lilly 2022a), stellar bars (Kubryk et al. 2013(Kubryk et al. , 2015)), or mismatch in the specific angular momentum of accreting gas and the disc (Mayor & Vigroux 1981;Pezzulli & Fraternali 2016).The latter is caused by turbulence in galactic discs that mix metals produced in different parts of the galaxy (e.g., due to diffusion or thermal instabilities, Yang & Krumholz 2012;Petit et al. 2015;Beniamini & Hotokezaka 2020).Except for a handful of galactic chemical evolution models, most models to date still lack a treatment for metal advection.The few that do include this process find it to be important for setting resolved metallicity trends in galaxies (Sharda et al. 2021a;Chen et al. 2023b).However, the treatment of radial gas flows is often decoupled from star formation and gas accretion, which creates inconsistency in the treatment of metals and the gas, or allows for tunable parameters that are then matched to an observable to constrain metal advection.This is further complicated by the fact that direct evidence for radial gas flows is scarce (Schmidt et al. 2016;Di Teodoro & Peek 2021;Sharda et al. in prep.), even though it has been long suspected to play an important role in galaxy evolution (Tinsley 1980;Lacey & Fall 1985).
The diffusion of metals due to turbulence warrants a detailed discussion since it is a relatively new component of galactic chemical evolution models (Sharda et al. 2021a), and is usually implemented as sub-grid physics in hydrodynamical simulations (Escala et al. 2018).Galactic discs typically exhibit supersonic turbulence, as is now known from measurements of the multi-phase gas velocity dispersion across cosmic time (Ianjamasimanana et al. 2012;Moiseev et al. 2015;Wisnioski et al. 2015;Varidel et al. 2016Varidel et al. , 2020;;Simons et al. 2017Simons et al. , 2019;;Wisnioski et al. 2019;Förster Schreiber et al. 2018;Übler et al. 2019;Girard et al. 2021).Simulations have shown that, in the absence of strong anisotropies (Hansen et al. 2011;Beattie & Federrath 2020) or strong magnetic fields oriented perpendicular to the disc (Kim & Basu 2013), supersonic turbulence decays very rapidly compared to the dynamical time of the disc (Mac Low et al. 1998;Stone et al. 1998;Lemaster & Stone 2009;Downes & O'Sullivan 2009, 2011;Burkhart et al. 2015;Kim & Ostriker 2015b).Therefore, turbulence needs to be continuously driven in galactic discs in order to explain the observed high gas velocity dispersions that cannot be a result of thermal gas motions.cite Sarrato-Alós et al. (2023) The two most common sources of turbulent energy injection into the interstellar medium (ISM) in galactic discs are star formation feedback due to supernovae (Ostriker & Shetty 2011;Faucher-Giguère et al. 2013) and radial gas flows.Krumholz et al. (2018) (hereafter, KBFC18) present a unified model of galaxies that includes both these processes, and show that star formation feedback alone can maintain a floor velocity dispersion of 6 − 10 km s −1 , explaining why the gas velocity dispersion tends to settle down to this value.KBFC18 also find that the high gas velocity dispersions seen in starburst galaxies are primarily driven by gravity due to radial mass transport.This model has been used to reproduce the observed star formation rate (SFR) -gas velocity dispersion trend in several surveys (Johnson et al. 2018;Yu et al. 2019;Varidel et al. 2020;Yu et al. 2021;Girard et al. 2021; see also, Ejdetjärn et al. 2022).A number of other authors have also published models to explain the role of turbulence in galaxy evolution but without the dual role of feedback and transport while simultaneously considering energy and momentum balance (KBFC18, section 1.2 and references therein).
However, KBFC18 only consider gas accretion from outside the galaxy as a source that supplies mass to the galactic disc, but did not consider the possibility that it might also drive turbulence in the disc (Ginzburg et al. 2022).Such driving is possible because the accreting gas can convert some of its kinetic energy into turbulence, as has been demonstrated in several simulations (Klessen & Hennebelle 2010;Elmegreen & Burkert 2010;Genel et al. 2012;Gabor & Bournaud 2014;Heigl et al. 2020;Forbes et al. 2023).Some works (e.g., Hopkins et al. 2013) argue that this effect is negligible; however, they do not consider all regimes (both in terms of halo mass and redshift) where accretion can be dominant.Some simulations are inconclusive about the importance of accretion in driving turbulence in discs (Orr et al. 2020), while others find it to be quite important, at least for certain halo masses and redshifts (Gabor & Bournaud 2013;Forbes et al. 2023;Jiménez et al. 2023;Goldman & Fleck 2023).Accretion is less efficient at driving turbulence if it is more spherical rather than in the form of thin filaments, or if the accreting streams are coherent and co-rotating with the disc, which might be the case for thin-disc star-forming galaxies (Danovich et al. 2015;Hafen et al. 2022), but this hypothesis is complicated by the multi-phase nature of accreting streams (Cornuault et al. 2018).Regardless of the predictions of these simulations, it is clearly the case that if gas accretion drives turbulence in the disc, several assumptions about the origin of gas velocity dispersions (and other quantities that are influenced by the gas velocity dispersion, like metallicity) need to be revisited (Hopkins et al. 2013).

This work and its motivation
Only a handful of theoretical works have focused on the spatiallyresolved metal distribution and MZGR (Ho et al. 2015;Kudritzki et al. 2015;Lian et al. 2018;Belfiore et al. 2019;Yates et al. 2021), and none of these investigate metal dynamics due to both radial gas flows and turbulent diffusion.Ho et al. (2015), Kudritzki et al. (2015), and Belfiore et al. (2019) further assume that galactic winds are not metal enriched as compared to the ISM, which has been well known to be a severe limitation of chemical evolution models (e.g., Dalcanton 2007;Peeples & Shankar 2011;Chisholm et al. 2018;Cameron et al. 2021).Several studies have discussed the role of winds in the context of chemical evolution models and shown that they are crucial to explaining trends in the MZR (e.g., Finlator & Davé 2008;Feldmann 2015;Chisholm et al. 2018;Choi et al. 2020;Tortora et al. 2022), but the corresponding role of winds in setting the spatially-resolved metal distribution remains largely unexplored.
Cosmological simulations have also investigated the MZGR, with varying degrees of success in reproducing the observed trends (Tissera et al. 2016(Tissera et al. , 2019;;Ma et al. 2017;Collacchioni et al. 2020;Hemler et al. 2021;Porter et al. 2022).While these simulations provide an excellent overview of the impacts of large-scale physical processes (including the role of CGM/IGM) and environment on metallicity profiles, they rely on subgrid recipes for star formation and metal diffusion.This is due to limitations imposed by finite resolution and the adopted temperature floor that does not capture the cold ISM.It is also difficult to control for important variables that influence metallicity profiles in these simulations, which leads to uncertainty around the relative contribution of different physical processes in setting the MZGR.
In a series of previous works, Sharda et al. (2021a) introduced a first principles model for the physics of gas-phase metallicity and metallicity gradients in galaxies.The authors built upon the unified disc model of KBFC18 by including metal dynamics to show how various galaxy properties govern the gas-phase metallicity profiles in galaxies (Sharda et al. 2021a,b,c).However, Sharda et al.only incorporated gas accretion as a process that adds mass to the disc, and did not self-consistently model gas turbulence and its relationship to accretion.In this work, following the updates to the Krumholz et al. (2018) model by Ginzburg et al. (2022), we self-consistently model the role of gas accretion by also including it as a source of turbulence, properly capturing the relation between metal dilution by accretion of metal-poor gas and driving of turbulence by the kinetic energy of incoming material.Such a comprehensive treatment of gas accretion has not been included in existing chemical and metallicity evolution models, which renders our understanding of the importance of gas accretion in setting ISM metallicities incomplete.As we will show in this work, gas velocity dispersion plays a central role in setting spatially-resolved gas-phase metallicities, and metal diffusion due to turbulence becomes dominant over other processes in some galaxies.
We arrange the rest of the paper as follows: Section 2 presents the model for the evolution of gas-phase metallicity, Section 3 presents resulting metallicity gradients, Section 4 presents the local MZGR we obtain from the fiducial model, and Section 5 discusses the role of galactic winds in setting the MZGR.Finally, we present our conclusions in Section 6.For the purpose of this paper, we use  ⊙ = 0.0134, corresponding to 12 + log 10 O/H = 8.69 (Asplund et al. 2021), Hubble time  H(0) = 13.8Gyr (Planck Collaboration et al. 2018), and follow the flat ΛCDM cosmology: Ω m = 0.27, Ω Λ = 0.73, ℎ = 0.71, and  8 = 0.81.

MODEL FOR SPATIALLY-RESOLVED GAS-PHASE METALLICITY
In this work, we piece together models for galactic discs and gasphase metallicities from Forbes et al. (2014a), KBFC18, Sharda et al. (2021a), and Ginzburg et al. (2022) to present a semi-analytic model for gas-phase metallicity gradients that self-consistently includes the role of accretion, feedback, and transport.We list all the mathematical symbols representing different physical parameters for the galactic disc in Table 1.Similarly, we provide parameters for the metallicity model in Table 2.For the remainder of this work, we classify galaxies according to their stellar mass,  ★ , as follows: low-mass galaxies -log 10  ★ /M ⊙ ≤ 9, and massive galaxieslog 10  ★ /M ⊙ ≥ 10.5.
Any chemical evolution model should also be able to reproduce (or, be consistent with) gas-phase galaxy scaling relations since the prescription of metals in these models is inherently tied to that of the gas.The galactic disc model that we incorporate from KBFC18 and Ginzburg et al. (2022) has been verified against numerous scaling relations, including the (resolved and unresolved) Kennicutt-Schmidt relation, and the trend between star formation rate and galaxy kinematics (Yu et al. 2019;Varidel et al. 2020;Yu et al. 2021;Girard et al. 2021;Parlanti et al. 2023;Rowland et al. in, prep.).Similarly, modeling metallicity gradients has an additional constraint as compared to modeling integrated/unresolved metallicities -metallicity gradient models should also be able to simultaneously reproduce the scaling relations of integrated metallicity with galaxy properties.We have shown in previous works how our model also reproduces the MZR (Sharda et al. 2021b).Since the updates we present in this work provide a better description of metal dynamics within galaxies but does not change the global metal content in the disc, we only discuss and apply the new model to study metallicity gradients.
We base our model on the premise that galactic discs remain marginally stable due to gas transport driven by gravitational instabilities (Krumholz & Burkert 2010;Forbes et al. 2012Forbes et al. , 2014a)).In other words, we assume that galactic discs are able to self regulate to a particular value of the Toomre  parameter (Toomre 1964).Following Krumholz et al. (2018, equation 7), we set the Toomre  parameter of the gas where  is the epicyclic frequency given by  = √︁ 2(1 + )Ω; here,  =  ln  / ln  is the rotation curve index of the galaxy,   is the rotational velocity at the outer edge of the disc, and Ω is the angular velocity. g is the gas velocity dispersion, and Σ g is the gas surface density.We assume a flat rotation curve for massive galaxies ( = 0) and a rising rotation curve for low-mass galaxies ( = 0.5).Following Sharda et al. (2021b), we create a linear ramp to interpolate between these two for setting  for intermediate-mass galaxies (9 ≤ log 10  ★ /M ⊙ ≤ 10.5).
Following Romeo & Falstad (2013), we also express the total Toomre  parameter of the gas and the stars as where  g,Q is the effective gas fraction in the disc (KBFC18, equation 9).For discs to be marginally stable, we require that  ≥  min , where  min is the minimum Toomre  parameter needed for stability (see, however, Forbes 2023 for a more sophisticated approach).As a fiducial value, we set  min = 1 as appropriate for relatively quiescent discs. 1 We will discuss the role of  min in detail in Section 2.1.We begin by describing the evolution of the halo mass and gas mass (Section 2.1), from which we infer the gas velocity dispersion (Section 2.2) based on the Toomre criterion.We then incorporate these parameters into the metallicity model, solve for the metallicity at each radius, and check if the radial metallicity distribution is in equilibrium.If it is, we find radial metallicity profiles (Section 2.3) and calculate the metallicity gradient.This way, we can self consistently find the gas velocity dispersion and metallicity at each epoch without violating the total mass, energy and metal budget of the disc.We provide a schematic that captures the most important ingredients of our model in Figure 1.

Evolution of the halo mass and gas mass
To describe the evolution of the halo mass,  h , we first estimate the halo accretion rate,  h , following Neistein & Dekel (2008) and Bouché et al. ( 2010) where  = 0.628,  = 0.14 and the self-similar time variable  is such that We integrate equation 3 to find the halo mass,  h as a function of redshift.Note that we use the stellar mass -halo mass relation from Moster et al. (2013,    = 76.17 where  is the halo concentration parameter that encodes the merger history of dark matter haloes.We follow Zhao et al. (2009, fig. 16) to vary  as a function of  h , ignoring the ∼ 0.1 dex scatter in  at fixed  h (e.g., Jing 2000;Wechsler et al. 2002;Diemer & Joyce 2019), since changes by an order of magnitude in  only impact   by ∼ 50 per cent.We also evolve the size of the galactic disc following (Ginzburg et al. 2022) where we have added the leading factor of two to ensure we cover most of the star-forming part of the disc.
The rate of change of gas mass in the disc is controlled by processes that supply mass to the disc (e.g., accretion) and that remove mass from the disc (e.g., star formation and outflows).Thus where the terms on the right hand side represent, respectively, the accretion rate of gas onto the galactic disc (  g,acc ), radial transport of gas through the disc (  trans , also referred to as radial gas flows), star formation rate (  SF ), and mass outflow rate given by the product of the mass-loading factor  w and  SF .In writing equation 7, we have assumed that the gas is transported down the potential well of the galaxy and ends up in a bulge that forms a distinct component of the galaxy.We, however, do not track the evolution of the bulge in the model; the bulge is only intended to act as a boundary condition for the galactic disc rather than a separate reservoir.We have neglected the recycling of metals in the form of galactic fountains (Spitoni et al. 2013;Christensen et al. 2016;Tollet et al. 2019;Grand et al. 2019;Veilleux et al. 2020).Assuming it is the gas expelled by star formation feedback that gets recycled through a fountain, including galactic fountains would add a source term of the form  GF  SF in equation 7, where  GF describes mass-loading of material returned to the galaxy (e.g., Lapi et al. 2020).However, it is unclear if this is sufficient to describe galactic fountains because the mass and dynamics of fountains are not just set by star formation-driven winds, but also by material already present in the CGM (e.g., Pandya et al. 2022).Further, there is no consensus on the timescales on which fountains re-supply mass and metals to the galaxy (e.g., Spitoni et al. 2009;Anglés-Alcázar et al. 2017;Grand et al. 2019;Mitchell et al. 2020b).

Gas accretion
To define the gas accretion rate, we re-write it as where  B = 0.17 is the universal baryonic fraction (White & Fabian 1995;Burkert et al. 2010 In our present model, we take the SFR to be (Krumholz et al. 2018, equation 60) where  sf is the fraction of gas in the star-forming molecular phase,  a is a parameter of order unity that accounts for the integral of the star formation rate per unit area for galaxies with varying  (Krumholz et al. 2018, section 3.1.2), ff is the star formation efficiency per freefall time,  g,P denotes the ratio of the mid-plane pressure due to self gravity of the gas to the total mid-plane pressure,  mp denotes the ratio of the total to the turbulent pressure at the disc mid-plane,  orb,out is the orbital timescale at the outer edge of the disc, and  SF,max ≈ 2 Gyr is the maximum gas depletion timescale.Thus, terms in parentheses represent star formation feedback in the Toomre and the GMC regimes, respectively.In practice, we set  g,Q =  g,P ,  mp = 1.4,and  a = 2. Following numerous results (Krumholz et al. 2012;Federrath & Klessen 2012;Salim et al. 2015;Sharda et al. 2018Sharda et al. , 2019;;Hu et al. 2022), we keep  ff = 0.015, ignoring possible intrinsic scatter (Hu et al. 2021(Hu et al. , 2022)) f 2 B 9 / g C W g Z S l < / l a t e x i t > ṀSF < l a t e x i t s h a 1 _ b a s e 6 4 = " I 5 q w s i S g a o 2 f 0 i t 6 c z H l x 3 p 2 P e e u K U 8 w c o T 9 w P n 8 A X e K T Z g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " I 5 q w s i S g a o 2 f 0 i t 6 c z H l x 3 p 2 P e e u K U 8 w c o T 9 w P n 8 A X e K T Z g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " I 5 q w s i S g a o 2 f 0 i t 6 c z H l x 3 p 2 P e e u K U 8 w c o T 9 w P n 8 A X e K T Z g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " I 5 q w s i S g a o 2 f 0 i t 6 c z H l x 3 p 2 P e e u K U 8 w c o T 9 w P n 8 A X e K T Z g = = < / l a t e x i t > ra di al flow s gas L S A g 4 R l e 4 c 0 z 3 o v 3 7 n 3 M W 0 t e M X M I f + B 9 / g D 3 M J E F < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " variable gas reservoir: Ṁg < l a t e x i t s h a 1 _ b a s e 6 4 = " y a b N A H 1 2 e P p r X g 5 V w o n r e W + 2 z 6 0 r m 6 A E 9 o W f r 3 n q 0 X q z X V W v J K m Z q 6 A e s t 0 8 W P Z W s < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / r m 6 A E 9 o W f r 3 n q 0 X q z X V W v J K m Z q 6 A e s t 0 8 W P Z W s < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / r m 6 A E 9 o W f r 3 n q 0 X q z X V W v J K m Z q 6 A e s t 0 8 W P Z W s < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / / y e t 4 7 J j l 5 3 L S q l u z + P I o z 2 0 j w 6 R g 6 q o j s 5 R A z U R R X f o A T 2 h Z + v e e r R e r N d Z a 8 6 a z + y i H 7 D e P g G o w 5 r r < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " G 3 9 m J z e p J o 5 h c Y K F n 7 / y e t 4 7 J j l 5 3 L S q l u z + P I o z 2 0 j w 6 R g 6 q o j s / y e t 4 7 J j l 5 3 L S q l u z + P I o z 2 0 j w 6 R g 6 q o j s / y e t 4 7 J j l 5 3 L S q l u z + P I o z 2 0 j w 6 R g 6 q o j s We consider marginally unstable galactic discs (quantified by the Toomre  parameter) in pressure and energy balance.The evolution of gas mass,  g , is described by equation 7. It depends on the gas accretion rate (  g,acc ), star formation rate (  SF ), radial gas flows (or, gas transport,  trans ), and galactic winds ( w  SF , where  w is the mass-loading factor).Gas turbulence (quantified by the gas velocity dispersion,  g ) is driven by energy injected into the ISM via accretion, star formation feedback, and gas transport.Metals are produced as star formation leads to supernovae (SNe).Some fraction of the metals produced is returned almost instantaneously to the ISM (  R,inst ) whereas the rest is locked in long-lived stars.Metals are also ejected via galactic winds with metallicity  w , which can be greater than or equal to .  y denotes the preferential metal enrichment of winds.Once produced, metals are also advected with radial gas flows, and diffused into the ISM due to turbulence and thermal instabilities.© AAS.Reproduced with permission.
et al. 2022) or redshift (Tacconi et al. 2020).We express  sf as the ratio of the molecular to the molecular and atomic gas mass, and find it following a compilation of observations by Tacconi et al. (2020, table 3) and Saintonge & Catinella (2022, figure 4).We also consider a theoretical model that can specify radial variations in  sf (Krumholz 2013) in Appendix A.

Gas transport
As gravitational instabilities and non-axisymmetric torques transport angular momentum outwards from the galaxy centre, mass is transported inwards to conserve angular momentum (a process similar to that observed in accretion discs around protostars and black holes - Shakura & Sunyaev 1973;Hartmann et al. 1998).The effective transport rate can be estimated either from considerations of energy conservation (Krumholz & Burkert 2010) or, in the case of gas-rich discs that produce massive clumps, from analysis of clump kinematics (Dekel et al. 2009b).Since the two approaches differ only by a small factor in their predictions (and not at all for galaxies with flat rotation curves), we adopt the transport rate by KBFC18, who adopt the energy conservation approach.Their formulation is based on the premise that star formation feedback can also drive turbulence in the disc in addition to gas transport.Ginzburg et al. (2022) extend the KBFC18 formulation to include gas accretion as another source of turbulence in the disc.We combine the two approaches set in (KBFC18, equation 49) and Ginzburg et al. (2022, equation 39) to define the gas transport rate as where  trans > 0 corresponds to inward mass flow,  = 1.5 defines the dissipation of turbulence across the scale height of the disc in one crossing time,  Q is given by and  nt parameterises the contribution of thermal motions to the total gas velocity dispersion, where  th is the thermal gas velocity dispersion.Following Krumholz et al. (2018, section 2.4.3),we estimate  th =  sf  th,mol + (1 −  sf  th,WNM ).Here,  th,mol ≈ 0.2 km s −1 is the thermal velocity dispersion of molecular gas assuming a gas temperature ≈ 10 K, and  th,WNM ≈ 5.4 km s −1 is the thermal velocity dispersion of the warm neutral medium (WNM, Wolfire et al. 2003).
The last term in equation 10,  ( g ), warrants a detailed discussion.It represents the fractional contribution of star formation feedback and gas accretion as compared to transport in driving turbulence in the disc.Below, we analyze different cases where we only include star formation feedback or accretion as drivers of turbulence in the disc, to study their impact on the gas-phase metallicity gradients.
(i) Turbulence driven by feedback and transport.First, we only consider the case where star formation feedback and gas transport due to radial flows drive  g , and accretion only adds mass to the disc but does not drive turbulence in the disc.This is equivalent to setting where  SF is the gas velocity dispersion that can be maintained by star formation feedback.We define  SF as (Krumholz et al. 2018, equation 39) where we maximise over terms that represent star formation feedback in the Toomre and the GMC regimes, respectively.⟨/⟩ ★ ≈ 3000 km s −1 represents the average momentum per unit stellar mass formed that is injected into the ISM by non-clustered core-collapse supernovae (Cioffi et al. 1988;Thornton et al. 1998 (ii) Turbulence driven by accretion and transport.Next, we consider the case where only accretion and transport drive  g .Here, where  acc is the gas velocity dispersion due to turbulence induced by gas accretion where  a is the fraction of kinetic energy of the accreted gas that drives turbulent motions in the disc. 2 Klessen & Hennebelle (2010) argue that  a is set by the density contrast between the accreting material and the material in the disc.The density of the accreting material depends on the clumpiness of the accreting streams (Mandelker et al. 2020).Taking inspiration from models that find the clumpiness of accreting streams to strongly vary with redshift (Mandelker et al. 2018), Ginzburg et al. (2022) propose  a = 0.2(1 + ).Further, Forbes et al. ( 2023) also find  a ≈ 0.1 − 0.2 on average across the disc at  ≈ 0 for a  ★ = 5×10 9 M ⊙ galaxy they study in the Illustris TNG50 suite of simulations (Nelson et al. 2019).However, there are no direct observational measurements of  a to date.Keeping these facts in mind, we fix  a = 0.2(1 + ), and show how our results remain unchanged for a higher  a in Appendix B.
(iii) Turbulence driven by feedback, accretion and transport.In this case, star formation feedback, gas accretion and gas transport all drive turbulence in the galaxy.The equivalent expression for  ( g ) then becomes (Ginzburg et al. 2022) It is worth noting that   equilibrates to a higher value when both accretion and feedback drive turbulence as compared to the cases above.In the analysis that follows, we use different versions of  ( g ) to understand the impacts of various drivers of turbulence on gas-phase metallicity gradients.

Winds
The role of galactic winds for the evolution of  g is described by  w (see equation 7).Several works have focused on developing models for  w (Creasey et al. 2013(Creasey et al. , 2015;;Forbes et al. 2014b;Hayward & Hopkins 2017;Torrey et al. 2019;Tacchella et al. 2020), as well as measuring it in simulations (Muratov et al. 2015;Christensen et al. 2016;Pillepich et al. 2018b;Kim et al. 2020;Pandya et al. 2021).These works find  w that spans more than three orders of magnitude (Mitchell et al. 2020a, fig. 13), which is a substantially larger range than that estimated from direct observations of galactic winds.This scatter is a complex combination of varying details of (star formation and/or AGN) feedback models, accurately resolving the multi-phase nature of galactic winds, the location at which  w is measured (in isolated versus cosmological simulations), treatment of the CGM, and burstiness of star formation.Observations estimate  w ≈ 0 − 30 (Bouché et al. 2012;Newman et al. 2012;Kacprzak et al. 2014;Schroetter et al. 2015Schroetter et al. , 2019;;Chisholm et al. 2017;Davies et al. 2019;Förster Schreiber et al. 2019;McQuinn et al. 2019), but are plagued by systematics originating from the location as well as the thermal phase in which it is measured.
In our fiducial model, we leave  w undefined because no models exist that self-consistently treat the role of outflows in driving turbulence in galactic discs together with the other sources we describe above. 3In practice, this means that we do not consider the wind term in equation 7 in the fiducial model, but we explore the full possible range of preferential metal ejection via winds below in Section 2.3.To rectify this inconsistency, we will later consider three different scalings of  w derived from theoretical models and simulations to discuss the effects of a non-zero  w in Section 5. Note that we do not consider AGN feedback in our model.In massive haloes, AGN feedback can boost  w (Mitchell et al. 2020a), so it is likely that some of the scalings we consider based on studies excluding AGN feedback underestimate  w at the high mass end.

Gas velocity dispersion
With all the terms in equation 7 defined, we can now integrate the equation to obtain the evolution of the gas mass as a function of redshift.The only input parameter that we need is the halo mass at some high redshift, which we use to evolve the galaxy down to  = 0. 4 The exact choice of redshift is not important since the solution quickly converges to its steady-state value within a few orbital times (Ginzburg et al. 2022, fig. B1).
In cases where  ( g ) ≥ 0, we use the Toomre  criterion (equation 1) to obtain  g ().However, we also encounter cases where  ( g ) < 0; physically, this occurs when mass transport is not needed to drive the required level of gas velocity dispersion, and the disc self-regulates the Toomre  parameter.In such cases,  ≥  min , and we find  g by solving for  ( g ) = 0. 5 This is a major improvement over Sharda et al. (2021a) because we selfconsistently derive  g based on the overall mass and energy budget in galactic discs instead of setting it to ad-hoc values.Although only a handful of theoretical studies have focused on possible correlations between  g and gas-phase metallicity (Ma et al. 2017;Krumholz & Ting 2018;Sharda et al. 2021c;Hemler et al. 2021), it is expected that the ISM metal distribution is sensitive to  g (Queyrel et al. 2012;Gillman et al. 2021), so it is crucial to self-consistently find  g .Solving for  g in this manner means that we do not take any radial variations in  g into account; however, these are expected to be minor, since both modeling (KBFC18) and observations (e.g., 3 Leaving  w undefined is not the same as setting it to 0 since the latter would constrain  y to unity as per equation 26, not permitting variations in  y .Since our aim in this section is to explore the full range of  y , we leave  w undefined. 4We use SCIPY solve_ivp to numerically integrate equation 7 from high- to  = 0. Given the stiffness in the differential equation that arises from the different timescales of accretion and star formation, we employ the backwards differentiation formula (BDF, Shampine & Reichelt 1997) method to search for stable solutions at each . 5 The equation we solve is a sixth degree polynomial in  g , as compared to the cubic equation in Ginzburg et al. (2022)    The three curves corresponds to models where turbulence is driven by either feedback due to star formation, gas accretion, gas transport, or a combination of the three.Wilson et al. 2011;Mogotsi et al. 2016) show that  g varies with radius by at most a factor of two in local galaxies.We do not make a distinction between  g in the star-forming molecular phase versus the ionized phase, though we caution that recent observations and simulations have shown the former can be systematically lower (Girard et al. 2021;Ejdetjärn et al. 2022;Rathjen et al. 2022;Lelli et al. 2023).
Figure 2 shows the resulting  g for the three cases where gas turbulence is driven by feedback + transport (blue), accretion + transport (orange), and feedback + accretion + transport (green).We find that, in energy balance, accretion-driven turbulence (with efficiency  a set to 0.2) plays a minor role in setting  g in lowmass galaxies, but it is the dominant contributor to turbulence in massive galaxies.The high  acc at high  ★ is partially caused by the overestimation of  g,acc in massive haloes as we mention in Section 2.1.Nevertheless, the value of  g at the high-mass end is in good agreement with that observed in local galaxies (Epinat et al. 2008;Moiseev et al. 2015;Varidel et al. 2020).
Conversely, feedback plays a major role in driving turbulence in low-mass galaxies but becomes sub-dominant compared to accretion in massive galaxies.This is seemingly contrary to the findings of both KBFC18 and Ginzburg et al. (2022) who find that mass transport is the dominant driver of turbulence in massive local galaxies.However, KBFC18 did not consider accretion as a possible source of turbulence that would reduce the amount of transport needed to maintain  g , and Ginzburg et al. (2022) did not consider the GMC regime of star formation, which compared to their assumption that all galaxies are in the Toomre regime yields a higher  sf and hence reduces the transport rate  trans required to maintain energy equilibrium (cf.equation 7).Mass transport is thus very sensitive to model parameters that govern  sf and  w at  = 0.As we will discuss later, this has important consequences for the metallicity distribution and gradient in local galaxies.

Other potential drivers of turbulence
We pause here to briefly mention other potential sources of turbulence in the disc that we have not taken into account in this work.
Feedback from supermassive black holes is another potential Constant of integration in the equation 29 metallicity equation source of turbulent energy injection in galaxies (e.g., Churazov et al. 2004;Zhuravleva et al. 2014;Li et al. 2020).For instance, both AGN-driven winds (Wagner et al. 2013) and radio jets (Mukherjee et al. 2018;Zovaro et al. 2019;Nesvadba et al. 2021) can enhance turbulence throughout the galactic disc on kpc scales.We also ignore turbulence injected by magnetorotational instabilities (MRI) which could be important in the outskirts of galactic discs (Sellwood & Balbus 1999;Piontek & Ostriker 2005), since the current evidence remains rather inconclusive (Tamburro et al. 2009;Utomo et al. 2019;Bacchini et al. 2020).
Another class of potential turbulence drivers are stellar components of the galaxy that can transfer energy to the gas, for example, spiral arms and bars.While these can drive turbulence in the gas (Kim & Ostriker 2007;Haywood et al. 2016;Khoperskov et al. 2018;Orr et al. 2023), simulations spanning a variety of redshifts and halo masses find that they remain sub-dominant compared to feedback, transport, and accretion (Bournaud et al. 2007;Bournaud & Elmegreen 2009;Ceverino et al. 2010;Goldbaum et al. 2015Goldbaum et al. , 2016)).Nevertheless, the impact of these drivers of turbulence on the gas-phase metallicity remains largely unexplored (e.g., Zurita et al. 2021;Chen et al. 2023a).

Evolution of gas-phase metallicity
The metallicity model of Sharda et al. (2021a) is a standalone model that uses a galaxy evolution model as an input to predict the metallicity profile in a wide variety of galaxies.The model includes several key physical processes -gas accretion, radial gas flows, metal advection and diffusion, star formation, metal consumption in long-lived stars, and outflows.In its primitive form, the evolution of spatially-resolved gas-phase metallicity in this model is given by (Sharda et al. 2021a where  and  are the dimensionless galactocentric distance and time variables, respectively. = / 0 , where  0 = 1 kpc is the reference radius that we assume to be the inner edge of the disc.Similarly,  =  Ω 0 , where Ω 0 is the angular velocity at  0 .Z is the gas-phase metallicity  normalised to Z ⊙ .The parameters  g (),  (),  ★ (), and  ★ () respectively represent the functional forms of spatial variations in the gas surface density Σ g , epicyclic frequency , star formation rate surface density Σ ★ , and gas accretion rate surface density Σ g,acc .
We pause here to mention a key difference between our approach and that of Sharda et al. (2021a).These authors find  g () = 1/,  () =  by assuming that   is the same throughout the galactic disc.However, this is only valid as long as  ∼ 0. We further generalize their results by re-writing   =  ,r (/) − where  ,r is the rotational velocity at a distance  from the galaxy center, and  is the disc size we define in equation 6.The advantage of re-writing the rotational velocity in this way is that we retain the meaning of   as that defined at the outer edge of the disc.With this generalization, we obtain It is expected that  ★ declines with  (Chiappini et al. 2001;Fu et al. 2009;Courty et al. 2010;Forbes et al. 2014b;Pezzulli & Fraternali 2016;Mollá et al. 2016;Stevens et al. 2017).We retain  ★ () = 1/ 2 from Sharda et al. (2021a).While the motivation for such a radial variation in Σ acc is that it can reproduce the present-day stellar surface density of the Milky Way (Colavitti et al. 2008), there is no reason why all galaxies should follow this profile.Nevertheless, (Sharda et al. 2021a, appendix A) show that a different functional form does not qualitatively change the resulting metallicity gradients in local spirals, and does not matter at all for local dwarfs, and we therefore adopt a single functional form for simplicity.
Since we include both the Toomre and the GMC regime of star formation, we update the definition of  ★ such that it is a continuous function of where  b is the critical location in the disc where Σ ★ in the GMC regime equals that in the Toomre regime (Krumholz et al. 2018, equation 32) Note that the star formation rate profile is less steep in the GMC regime, which will lead to more metal production in the outer regions of the galaxy as compared to the Toomre regime.As we will discuss later in Section 4, this feature is partially responsible for giving rise to steady-state inverted metallicity gradients in the model.
The four dimensionless ratios that connect the metallicity model to the evolution of gas and stars are: (1).T − the ratio of the orbital and metal diffusion timescales, (2).P− the ratio of metal advection to metal diffusion, well known as the Péclet Number in fluid dynamics (Patankar 1980), (3).S− the ratio of metal production (star formation) to metal diffusion, and (4).A− the ratio of gas accretion to metal diffusion.As in Section 2.1, we have not taken into account the effect of galactic fountains in returning metals to the disc.As in equation 19 and equation 20, we also generalize the expressions for T , P, S and A from Sharda et al. ( 2021a) where  max = / 0 is the disc size normalized by  0 , and  y is a fractional quantity that describes the preferential metal enrichment of galactic winds. 6Most galactic chemical evolution models to date assume that the metallicity of galactic winds is the same as that of the ISM.However, this is not necessarily the case for all galaxies. y accounts for the fact that the wind metallicity can be higher than the ISM metallicity in a galaxy if the ejecta from supernovae are imperfectly mixed with the ISM (Dekel & Silk 1986;Mac Low & Ferrara 1999;Recchi et al. 2008;Peeples & Shankar 2011;Fu et al. 2013;Emerick et al. 2018Emerick et al. , 2019;;Taylor et al. 2020;Sharda et al. 2021a;Yates et al. 2021;Andersson et al. 2022), as has now been observed in galactic winds of low mass galaxies (Chisholm et al. 2018;Lopez et al. 2020;Cameron et al. 2021;Xu et al. 2022;Tortora et al. 2022).Following Sharda et al. (2021b, equation A1), we write where  = 0.028 is the yield of newly formed Type II metals (Forbes et al. 2019) following the instantaneous recycling approximation (Tinsley 1980) and Z w is the metallicity of galactic winds normalized to Z ⊙ .It is given by (Forbes et al. 2019, equation 41) where  R,inst = 0.77 is the fraction of metals locked in low mass stars (Tinsley 1980), and  w is a fractional parameter bounded between 0 and 1 that describes the preferential enrichment of galactic winds if some of the newly produced metals are ejected with winds before they mix with the ISM.Thus, equation 27 implies that common assumption of wind metallicity being equal to the ISM metallicity is likely only a lower bound in reality.Further, since a fraction  R,inst of metals are locked in stars, the minimum metal mass ejected is 1 −  R,inst .We see from equation 26 that  y ∼ 0 corresponds to all newly produced metals being entrained in the winds, whereas if  y ∼ 1, newly produced metals are perfectly mixed with the ISM before they are ejected such that the wind metallicity equals the ISM metallicity.In the absence of any scaling relations of  y with galaxy properties, we explore  y in the range 0.1 − 1, and show how the resulting metallicity gradients are quite sensitive to the choice of  y .In principle, we can specify  y as a spatially-varying parameter.However, we treat it as an integrated quantity for this work since variations in  y with  remain unexplored.As we will discuss below, P, S and A form the cornerstone of our metallicity analysis and act as useful diagnostics to understand what sets metallicity gradients in galaxies.
Note that we only consider elements produced by core collapse supernovae in this work, since the observable we study is the oxygen abundance gradient in the ISM.Therefore, we cannot comment on the evolution of  y for elements produced by other nucleosynthetic sources (e.g., AGB stars or Type Ia supernovae).In fact, simulations of low mass galaxies that explicitly resolve feedback find metal enrichment of galactic winds to be dependent on the nucleosynthetic source (Emerick et al. 2018).Similarly, metal mixing in the ISM is also sensitive to the nucleosynthetic source (Krumholz & Ting 2018).

Metal equilibration timescale
Because we are interested in searching for steady-state solutions to equation 18, we first estimate the time it takes for a given metal distribution to reach equilibrium,  eqbm .Following Sharda et al. (2021a, equation 19),  eqbm is given by The physical interpretation of the above equation is that we compare the time taken by gas accretion, star formation, outflows, metal advection and metal diffusion to induce changes in metallicity with the rate at which metallicity evolves with time as given by the leading term in equation 18.Under equilibrium, Z/ → 0. This is a common practice in chemical evolution models as it allows us to estimate if metallicity reaches equilibrium simultaneously with other galaxy properties (Davé et al. 2012;Lilly et al. 2013;Feldmann 2015;Krumholz & Ting 2018).
If the metal distribution does not reach equilibrium within a reasonable time, we consider the resulting gradients to be in nonequilibrium, in which case we cannot apply our model to study gradients.A necessary condition is that the metal equilibration timescale be less than the Hubble time.However, this is not sufficient since physical processes relevant for the model may change on timescales much shorter than the Hubble time.Out of all the physical processes we include, the timescale for star formation, quantified by the molecular gas depletion timescale,  dep , is usually observed to be the shortest (Genzel et al. 2015;Scoville et al. 2017;Utomo et al. 2017;Tacconi et al. 2018;Liu et al. 2019).In nearby galaxies, observations find  dep ≈ 2 Gyr (Leroy et al. 2008;Bigiel et al. 2011;Saintonge et al. 2011;Sun et al. 2023).Keeping this in mind, we consider the metallicity to be in non-equilibrium if  eqbm ≫  dep .To ensure that local fluctuations in metallicity are not responsible for setting the metallicity gradient (that would otherwise generate a random metallicity gradient set by a stochastic metal field), we investigate if  eqbm is longer than the timescale for such fluctuations to become steady.Typically, the latter is of the order of  fluc ≈ 300 Myr in local galaxies (Krumholz & Ting 2018, fig.7), and we check if  eqbm ≳  fluc .Given the self-consistent evolution of gas parameters with metals, the model we present here allows for inverted gradients to form in equilibrium, in contrast to the Sharda et al. (2021a) model.et al. (2021a) show that most galaxies tend to achieve a steady-state metallicity gradient within a timescale much shorter than the Hubble time, and comparable or shorter than  dep .When this condition is satisfied, we can search for equilibrium solutions for equation 18 by setting Z/ → 0. Since the functional form of the star formation term,  ★ (), is different in the Toomre and the GMC regimes of star formation, the resulting solution for Z() is also different.Z() in the Toomre regime is given by functions that are not illuminating, but they can be directly obtained by the applying the above constraints.We set Z CGM = 0.05 for low-mass galaxies, 0.2 for massive galaxies, and create a linear ramp in log 10  ★ between these two for intermediate mass galaxies.Such a setup roughly reproduces the metallicity of metal-poor inflows onto the disc seen in simulations (Muratov et al. 2017).As more measurements of Z CGM become available (e.g., Prochaska et al. 2017;Kacprzak et al. 2019), it will be possible to refine the prescription we use.

Sharda
This completes the formulation of the model.In Table 3, we provide a summary of key differences between our approach and earlier works we make use of.

EQUILIBRIUM METALLICITY GRADIENTS
We now present resulting ISM metallicity profiles and gradients from the model for two classes of galaxies.Table 4 lists the values of model parameters for the two representative galaxies.The top panel of Figure 3 presents the resulting family of metallicity profiles from the model for a fiducial local massive galaxy with  ★ = 10 10.8 M ⊙ at  = 0 with  y = 1.0 (solid curves) and 0.1 (dotted curves).We include all three sources of turbulence (feedback, accretion, and transport) to create these profiles.The family of curves results from the permissible range of values of the constant  1 in the metallicity equation.However, not all the gradients are in equilibrium.The colorbar in Figure 3 informs on the time taken by the metal distribution to reach equilibrium compared to  dep ; the redder the color, the larger are the deviations from equilibrium.As expected, we see that the metallicity at the inner edge of the disc naturally approaches the value Z  0 set by the terms that dominate at small radii.The inner regions ( < 5) of all the profiles are in the Toomre regime of star formation whereas the outer regions are in the GMC regime, as expected from the KBFC18 model.The bottom panel of Figure 3 presents the family of metallicity profiles for a fiducial local low-mass galaxy with  ★ = 10 9 M ⊙ .In this case, the Toomre regime exists for  < 3.5, beyond which star formation occurs in GMCs.For the lowest mass galaxies we study, the entire star-forming disc is in the GMC regime at  = 0.For the massive galaxy, the radial metallicity profiles for high values of  y are consistent with the average metallicity profile observed in nearby galaxies of similar mass (Belfiore et al. 2017).On the other hand, model profiles with low values of  y in the low-mass galaxy better match the observed profiles for similar mass (see also, Sharda et al. 2021a, figures 2 and 5).
To find a 1D gradient, we fit the logarithmic metallicity profiles with a linear function in  from 1 to  max (see Table 4), the slope of which is the metallicity gradient.This is an oversimplification because the profiles are typically more complex given the non-linearity of Z(), so a gradient returned by the fit is often a poor representation of the underlying profile.However, we use this approach to mimic observational measurements where the 2D maps are reduced to 1D by stacking metallicity profiles as a function of galactocentric distance, and then fit with a linear function in log 10 Z.For the massive galaxy, we find gradients to be in the range −0.18 dex  −1 e to 0.08 dex  −1 e , where we take  e ≈ /2 in the model; adopting  e from van der Wel et al. (2014) (as in Sharda et al. 2021a) gives similar results.Similarly, we obtain gradients in the range −0.23 dex  −1 e to 0.13 dex  −1 e for the low-mass galaxy.Typically, profiles given by the maximum possible value of  1 (especially for massive galaxies) are out of equilibrium, and these profiles lead to inverted gradients.While most inverted metallicity gradients are out of equilibrium in the model, we do find some steady-state profiles where the metallicity in the outer regions is larger than that in the inner regions.This differs from our earlier models in Sharda et al. (2021a), which did not include the GMC regime of star formation and did not find any inverted gradient equilibria.
Given this new finding, it is worth investigating how these inverted gradients form and exist in equilibrium.Let us consider a simpler case of a massive galaxy with no advection, so  = P = 0.In such a case, at large radii, accretion and diffusion compete with star formation to set the metallicity.However, star formation in the GMC regime only declines with radius as 1/, while accretion and diffusion go as 1/ 2 .So, even though fewer metals are being produced at large radii, the rate at which they are diluted or transported is even smaller.This is why metallicity in the outskirts can be higher, leading to inverted gradients in the model in equilibrium.However, there are two important corollaries.First, this result is sensitive to our assumed radial profile of accretion,  ★ = 1/ 2 (Colavitti et al. 2008).If we assume that accretion goes as 1/ instead, for example, the resulting gradient will likely be close to zero.Second, if  sf declines strongly with radius (e.g., Krumholz et al. 2009;Krumholz 2013;Kubryk et al. 2015), the source term S would decrease more rapidly than 1/, and we would again have no or less inversion.We also do not invoke radial variations in  y or  w below (e.g., Johnson et al. 2021).Exploring radial variations in these parameters requires direct measurements of both the molecular and atomic gas surface densities, as well as metal enrichment of galactic winds at different radii.This is an area where a multi-wavelength observations, combining JWST IFU instruments, VLT/MUSE and ALMA (e.g., Leroy et al. 2021;Emsellem et al. 2022) have the potential to make a big impact on our understanding of inverted metallicity gradients.
Note that our model is based on a disc flow model, but (major) mergers can significantly impact the disc, and likely reset the metallicity gradients.We showed in Sharda et al. (2021a) that the metallicity distribution falls out of equilibrium during a major merger (a result corroborated by observations -Pérez-Díaz et al. 2023), and thus our equilibrium model cannot be applied in such cases.It is for this reason that we do not study merging galaxies (e.g., ULIRGs) in this work.The impact of minor mergers is small at resetting the gas-phase metallicity gradient, especially at  = 0, which is the focus of this work.There is an associated question of the impact of past major mergers on the present day metallicity gradient in the disc, answering which requires using a combination of hydrodynamic simulations and SAMs like ours (Sharda et al., in prep.), and is beyond the scope of this work.

MASS-METALLICITY GRADIENT RELATION (MZGR) WITH THE FIDUCIAL MODEL
We show the observed MZGR (corrected for spatial resolution) from three different surveys -MaNGA (Mapping Nearby Galaxies at Apache Point Observatory, Bundy et al. 2015), CALIFA (Calar Alto Legacy Integral Field Area, Sánchez et al. 2012), and SAMI (Sydney-AAO Multi-object Integral-field spectrograph, Bryant et al. 2015) in Figure 4.Each marker denotes the average metallicity gradient (normalized by  e ) in different bins of stellar mass.It is immediately clear that even between the three surveys, average gradients differ by as much as 0.1 dex  −1 e .In fact, the galaxy-togalaxy scatter at fixed  ★ is as high as 0.3 dex  −1 e .The scatter is larger at the low and high mass ends than at intermediate masses (Poetrodjojo et al. 2021, fig. 10).While some of this scatter can be physical (especially at the low-mass end), it is also caused by systematic differences between gas-phase metallicity calibrations (Kewley & Ellison 2008), limited spatial resolution and S/N ratio (Mast et al. 2014;Acharyya et al. 2020), the sensitivity of metallicity diagnostics to physical properties of H ii regions such as ionisation parameter (e.g., Pettini & Pagel 2004;Kobulnicky & Kewley 2004;Pilyugin & Grebel 2016;Curti et al. 2017;Mingozzi et al. 2020), and to poorly-understood contributions from diffuse ionized gas (DIG, Sanders et al. 2017;Zhang et al. 2017;Poetrodjojo et al. 2019;Vale Asari et al. 2019).An important effect of this discrepancy is that it is not clear if the MZGR shows an inflection at 9.5 < log 10  ★ /M ⊙ < 10.5 (Belfiore et al. 2017;Mingozzi et al. 2020); some studies favor such a break whereas others find a constant, characteristic gradient for all galaxies (Sánchez et al. 2014;Sánchez-Menguiano et al. 2016, 2018;Poetrodjojo et al. 2018).
Using the second data release of SAMI data, Poetrodjojo et al. (2021) reduce the discrepancy due to different emission line ratios and metallicity calibrations to 0.02 dex  −1 e by using a carefullychosen combination of various diagnostics that include more than two emission line ratios (e.g., Pilyugin & Grebel 2016).These authors find that there is indeed evidence for a break in the MZGR between 9.5 < log 10  ★ /M ⊙ < 10.5, however, the feature is rather weak as compared to the curvature observed in the MZR.Franchetto et al. (2021) reach the same conclusion based on their analysis of the MZGR in the MaNGA and GASP (GAs Stripping Phenomena in galaxies with MUSE, Poggianti et al. 2017) surveys.Deciphering the origin of the MZGR and determining whether there is an inflection at intermediate masses is important to understanding how various physical processes we describe in Section 2 impact metal distribution in galaxies.
To construct an MZGR from our model, we specify a range of halo masses at high- such that the resulting halo mass at  = 0 span 10.9 < log 10  h /M ⊙ < 14.7 .As we saw in Section 3, the model produces a family of metallicity profiles.Below, we find the MZGR for different cases corresponding to the different drivers of turbulence (feedback, accretion, and transport).

Gas turbulence driven by feedback and transport
In this case, the Péclet number P that describes the ratio of metal advection to diffusion is set only by star formation feedback and gas transport ( g corresponding to the blue curve in Figure 2).This is the case studied in Sharda et al. (2021a) since the KBFC18 model we used only considers feedback and transport as sources of turbulence.However, there are several differences between this and earlier work, as we point out in Section 2.
The top panel of Figure 4 plots the resulting MZGR from the model in this case.We color-code the model MZGR by the ratio (P + A)/S; the reason for this coding is that the processes parameterised by the numerator (advection and accretion) tend to flatten the gradient, while the one parameterised by the denominator (star formation) tends to steepen it.Thus this ratio is effectively the ratio of flattening to steepening processes.We show the resulting values of P, A, and S for this case in the left panel of Figure 5.The scatter in S at fixed  ★ is caused by  y (see equation 24).In line with our expectations, we see that steep negative gradients develop when S > P + A; this process is largely controlled by the parameter  y , since at fixed  ★ increasing  y leads to larger S (as indicated by colour in Figure 4) and to steeper, negative gradients (corresponding to moving downward in Figure 4).The exception to this trend is the least massive galaxies with  ★ < 10 9.1 M ⊙ -these galaxies lie in the GMC regime of star formation, as against galaxies  The colorbar indicates the ratio of the strength of processes that create flat metal distributions (P and A) to those that create steep metallicity gradients (S) in equilibrium, as defined in Section 2. The top panel plots the model MZGR for the case where only gas accretion and transport drive turbulence in the galaxy.The middle panel corresponds to the case where only star formation feedback and gas transport drive turbulence.The bottom panel corresponds to the case where all three factors are included as drivers of gas turbulence.The scatter in the model at fixed  ★ corresponds to the results of varying   , the parameter that describes the preferential ejection of metals via galactic winds, from 0.1 − 1, as indicated.Overplotted (orange markers) are average metallicity gradients in bins of  ★ obtained from IFU surveys -MaNGA (Belfiore et al. 2017;Mingozzi et al. 2020), CALIFA (Sánchez et al. 2014;Sánchez-Menguiano et al. 2016), and SAMI (Poetrodjojo et al. 2018), as well as mean of the three (grey markers).P, A, and S) that set the metal distribution at  = 0 in the fiducial model presented in Section 4, as a function of stellar mass, for the case where turbulence is driven by star formation feedback and gas transport.The value of S is directly proportional to the parameter  y that describes the preferential enrichment of galactic winds (see equation 24), so for this parameter we plot a vertical range corresponding to varying  y from 0.1 − 1, as indicated by the colour bar.Metal transport due to radial gas flows (P) is only active for the lowest and most massive galaxies in this case.The dashed line denotes a value of unity, with values of a parameter less than unity indicating that turbulent diffusion is a more important process than the one described by that parameter.Middle panel: Same as the left panel but now turbulence is driven by gas accretion and transport.Metal advection is for most galaxies in this case.Right panel: Same as the other two panels but turbulence is now driven by all three mechanisms -feedback, accretion, and transport.Metal advection is active only in the lowest mass galaxies.These diagnostic plots aid in understanding the nature of metal distributions and metallicity gradients in each of the corresponding models presented in Figure 4.
with  ★ ≥ 10 9.1 M ⊙ where the inner regions are in the Toomre regime and the outer regions are in the GMC regime.Since the star formation timescale is capped by  SF,max in the GMC regime, the resulting S < P + A even with  y = 1, as we can also read off from the left panel of Figure 5.We also find that A > 1 for all  ★ in this case, implying that metal diffusion is weak as compared to gas accretion.However, metal diffusion dominates over metal advection and star formation in some low-mass galaxies that have P < 1 and S < 1.Since diffusion tends to homogenize the metal distribution by moving metals from high to low concentrations, strong diffusion leads to flatter gradients.The inflection we see in the MZGR at  ★ ≈ 10 9.2 M ⊙ occurs due to mass transport shutting off.
We further see that there is a large scatter in the gradient at fixed stellar mass (especially in low-mass galaxies) that narrows at the high-mass end.This scatter occurs because the gradients are quite sensitive to the choice of  y -low  y lowers S, driving flat or even inverted gradients, whereas high  y drives steep negative gradients.At the low-mass end, the data prefers low values of  y , implying a significant boost in preferential metal ejection in winds of low-mass galaxies.In the absence of accretion as a source of turbulence,  g equilibrates to a lower value that is closer or equal to  SF .This in turn causes P to go to zero, thus shutting off mass transport and driving flatter (and even inverted) metallicity gradients at the highmass end.This is complimented by the dramatic increase in S and an even larger increase in A because S ∝  2 g and A ∝  3 g .We find that this model fails to reproduce the observed MZGR, especially at the high-mass end.
At first look, this result seems contradictory to the MZGR presented in Sharda et al. (2021b), because these authors could also reproduce the high-mass end of the observed MZGR with just feedback + transport.However, their ad-hoc choice of  g was larger than that we derive in this work (blue curve in Figure 2).Larger  g led to P > 0 and smaller A in massive galaxies in Sharda et al. (2021b), yielding flatter gradients in good agreement with those observed.As we will discuss below, our self-consistent solution for  g which includes accretion as a source of turbulence does reproduce the MZGR at the high-mass end because of larger  g .This highlights the importance of accretion and turbulence in setting ISM metallicity distributions and gradients in massive galaxies.

Gas turbulence driven by accretion and transport
Next, we discuss the case where gas turbulence is driven by a combination of accretion and transport (orange curve in Figure 2).We show the predictions from the model in the middle panel of Figure 4, and the corresponding dimensionless ratios in the middle panel of Figure 5. Overall, this model produces flatter metallicity gradients as compared to the feedback + transport model in Section 4.2 above.It can also reproduce the diversity of metallicity gradients in lowmass galaxies, although the range of  y that best matches the data is narrower as compared in the feedback + transport case.Including accretion as a source of turbulence brings the model gradients closer to the data, but agreement is still quite poor for massive galaxies.This is because P > 1 in most galaxies where transport is active, so a combination of strong metal advection and relatively higher gas accretion drives the flat and inverted gradients in this model.Thus, we learn that removing feedback and accretion as sources of turbulence results in model gradients that are too flat or inverted at the high-mass end as compared to the observed MZGR.

Gas turbulence driven by accretion, feedback and transport
We have seen that if feedback and accretion are not included as sources of turbulence in the disc, the resulting model MZGRs do not reproduce the observed MZGR at the high-mass end, even with the uncertainty in  y .Now, we consider the more general case where  g is driven by feedback, accretion, and transport acting together (green curve in Figure 2).The bottom panel of Figure 4 plots the model MZGR with the data.We find that the agreement between the model and the data is much better at all  ★ as compared to the cases above where we exclude either feedback or accretion as sources of turbulence.Higher  g caused by turbulence due to all three sources (under energy balance) leads to smaller S and even smaller A that results in steeper gradients, particularly at the high mass end.We can also discern the same from the right panel of Figure 5.As compared to the other two models, we see that A is smaller due to larger  g .The inflection in the MZGR at  ★ ≈ 10 9.2 M ⊙ that we noticed in the feedback+transport model is also present in this case, due to P → 0 as mass transport shuts off.In equilibrium, the strength of mass transport is very sensitive to  sf at  = 0, and thus to variations in the star-forming gas fraction  sf to which  sf is proportional.We explore alternate models of  sf in Appendix A and show that the model MZGR at the high mass end (where transport may or may not be shut off) remains qualitatively similar.
We also see that low-mass galaxies strongly prefer a low  y , implying significant metal enrichment of galactic winds.Below, we will see that this conclusion holds irrespective of the mass-loading of galactic winds in low-mass galaxies.We remind the reader that we adopt  a = 0.2 for local galaxies, and show in Appendix B how the results change in the unlikely case where  a is higher.As we show in Appendix D using MaNGA data, these interpretations are robust to scatter caused by different metallicity diagnostics used to measure metallicity gradients in the data (see also, Poetrodjojo et al. 2021;Sharda et al. 2021a).

ROLE OF GALACTIC WINDS
As we mention in Section 2.1.4,the fiducial model we have presented so far is limited, and to some extent, inconsistent, because we include mass-loading of galactic winds ( w ) in the evolution of metallicity (cf.equation 27) but not in the evolution of the gas mass.There are several different scalings of  w available in the literature, which can broadly be grouped into two categories: one where  w is related to the large-scale properties of the galaxy (e.g.,  h and virial velocity, Murray et al. 2005;Muratov et al. 2015;Pillepich et al. 2018a;Mitchell et al. 2020a;Pandya et al. 2021), and another where it is described by properties internal to the disc (e.g., Σ g , Σ SFR and  g , Dekel & Krumholz 2013;Hayward & Hopkins 2017;Kim et al. 2020;Pandya et al. 2021;Bassini et al. 2023).To rectify the inconsistency between metals and the gas, we discuss three cases below where we adopt  w from the two groups of works as above.Our choice of adopted models is such that we can cover a broad range in  w at fixed  ★ .While this approach of adopting  w from other works still renders the treatment of outflows decoupled from the rest of the model, the advantage of using these recipes is that they can easily be incorporated into semi-analytical models like ours.

Effects of different mass loading parameterisations
(i) Mass-loading factor from Hayward & Hopkins (2017).Using an analytical model that describes the dual role of star formation feedback in regulating star formation as well as launching winds from the galactic disc, Hayward & Hopkins (2017) show that the gas fraction is a critical factor that determines  w in star-forming galaxies.These authors find that gas-rich galaxies are typically  Mitchell et al. 2020a).The dashed grey line demarcates 1 −  R,inst , the fraction of newly formed metals not locked in low mass stars (Tinsley 1980).
(ii) Mass-loading factor from FIRE-2 simulations.We use the best-fit scaling of the multi-phase ISM  w (measured in a radial shell between 0.1 − 0.2  vir ) as a function of  ★ from FIRE-2 zoom-in simulations (Hopkins et al. 2018) This scaling is based on an analysis of multi-phase winds in central galaxies of 13 haloes spanning a wide range of  h (Pandya et al. 2021, equation 15).A key improvement over similar studies is that Pandya et al. define a more physically-meaningful criteria to define outflowing gas particles that also encaptures the slow-moving hot wind.However, the set of simulations used for this analysis does not include AGN feedback and turbulent metal diffusion.
(iii) Mass-loading factor from EAGLE simulations: Lastly, we adopt the relation between  w measured in EAGLE galaxies (Schaye et al. 2015) as a function of  ★ from Mitchell et al. (2020a, fig. 15).Similar to Pandya et al. (2021), this  w corresponds to mass  leaving the ISM of the galaxy, and not the halo, and is measured for central galaxies.Crucially, the simulations used for this work by Mitchell et al. include AGN feedback.However, the criteria used by EAGLE to select outflowing gas particles is different than FIRE-2 (see section 2.6 of Mitchell et al. 2020a).
These scalings of  w also constrain the range over which  y can vary, as against the fiducial model where we swept across  y = 0−1.This is simply because a fixed value of  w limits the range of  w in equation 27 that in turn limits  y .However, as we can expect from equation 27, this constraint is only used in practice in galaxies where  w < 1 −  R,inst .
Figure 6 unifies these scaling relations to represent them as a function of  ★ , for illustration.We find that, at the low mass end,  w from all the three works are qualitatively consistent with each other.While different simulations qualitatively agree on the behaviour of  w with  ★ for low-mass galaxies, they do not converge on the partition of mass and metal loading in the different   (Hayward & Hopkins 2017).Colorbar denotes the permitted range of  y that describes the preferential enrichment of galactic winds (see equation 26).Middle panel: Same as the left panel but  w is scaled with  ★ following FIRE-2 cosmological zoom-in simulations without AGN feedback (Pandya et al. 2021).Right panel: Same as the other two panels but  w is scaled with  ★ following EAGLE cosmological simulations that include AGN feedback (Mitchell et al. 2020a).
phases of galactic winds as this is sensitive to the distance from the midplane where  w is measured.There also exists tension between  w measured in simulations and observations (e.g., Concas et al. 2022;Marasco et al. 2023).Simulations of isolated galaxies that can afford very high resolution tend to measure mass fluxes few kpc from the midplane (e.g., Kim et al. 2020;Vĳayan et al. 2020;Wibking & Krumholz 2022), whereas cosmological simulations tend to measure much farther away from the disc (see the discussion in Pandya et al. 2021).However, there is a large scatter in  w for simulated massive galaxies.The key reason for this difference is that EAGLE includes AGN feedback that boosts mass-loading in massive galaxies at  = 0.In the absence of AGN feedback, it is difficult to drive strong winds in massive galaxies due to their deep potential wells.While our model does not include AGN feedback and is closer in spirit to Hayward & Hopkins (2017), it is insightful to test how AGN feedback-dominated  w impacts the MZGR.

Results
We know from equation 7 that introducing  w ≠ 0 activates the outflow sink term that takes mass out of the disc and thus impacts  g .Since  g ∝  g , and the dimensionless ratios that govern the metallicity profile are all sensitive to  g , we expect wind mass loading to non-linearly modulate the metallicity profiles and resulting metallicity gradients.Figure 7, Figure 8, and Figure 9 plot the resulting model MZGR where we scale  w following Hayward & Hopkins (2017), Pandya et al. (2021), andMitchell et al. (2020a), respectively.We only plot the model for the case where turbulence is driven by a combination of feedback, accretion, and transport, as we have seen in Section 4 that neglecting either feedback-or accretiondriven turbulence does not reproduce the observed MZGR even in the fiducial case with  w undefined.
We first discuss results for the scaling of  w from Hayward & Hopkins (2017).We notice from Figure 7 that, as compared to the fiducial model (bottom panel of Figure 4), the range of metallicity gradients at the high mass end of the MZGR remains unaffected by this scaling.This is because even at fixed  y , the model produces a family of gradients due to constraints on  1 (see Section 2.3.3).However, the range of gradients at the low mass end is reduced by more than 0.1 dex r −1 e , due to a high  w > 1 for low-mass galaxies.While this scaling of  w can capture the observed MZGR, the mechanism giving rise to the agreement between the observations and the model is significantly different as compared to the fiducial model, as we see from the differences in the ratio (P +A)/S denoted by the colorbar.The transition in color from green to blue as a function of increasing  ★ implies a transition in the dominant mechanism that sets metallicity gradients in low-mass versus massive galaxies, and is responsible for driving the model MZGR.Additionally, the best match between the model and the data reveals a large scatter in the preferred value of  y for galaxies with  ★ < 10 9.4 M ⊙ : low-mass galaxies in the SAMI survey prefer  y ∼ 0.1, whereas MaNGA and CALIFA dwarfs prefer  y ∼ 0.4.The average measurement for the three surveys prefers  y ∼ 0.25.Massive galaxies ( ★ ≳ 10 10 M ⊙ ) prefer  y ∼ 0.8 − 1.0.
To further diagnose the cause for these differences, we plot the dimensionless ratios as a function of  ★ for this scaling in the left panel of Figure 10.At the high mass end,  y is restricted to a very narrow range close to unity because  w in the Hayward & Hopkins model is less than 1 −  R,inst , which limits the allowed range of Z w (cf.equation 27).Physically, we can understand this result very simply: for these galaxies the mass flux carried by the winds is so small that, even if this mass were to consist of pure SN ejecta, winds would carry away at most a small fraction of the metals being produced.Thus most metals must remain in the galaxy, corresponding to  y close to unity.Since the dimensionless ratio S ∝  y , this leads to a decrease in the quantity (P + A)/S, which in turn drives the colorbar to blue at the high mass end in Figure 7.We also see from the left panel of Figure 10 that P > 0 only for galaxies with  ★ < 10 10 M ⊙ , implying that mass transport shuts off in massive galaxies because accretion and feedback can drive sufficient levels of turbulence in these galaxies.Lastly, metal diffusion dominates over star formation for a wider range of  ★ in this case as compared to the fiducial model.
Next, we examine the MZGR using the scaling of  w from FIRE-2 zoom-in simulations.Figure 8 shows the results.The scatter in metallicity gradients at the low mass end is larger than both the fiducial model and the model above where  w scales according to the Hayward & Hopkins model.This correlates with the non-zero P as we show in the middle panel of Figure 10, and occurs because  w measured in FIRE-2 dwarf galaxies is smaller by almost an order of magnitude as compared to Hayward & Hopkins.On the contrary, the scatter at the high mass end remains identical to the other models, although the relative contribution of the dimensionless ratio S is larger than the fiducial model but smaller than the model with  w from Hayward & Hopkins.This occurs because the larger  w in FIRE-2, as compared to Hayward & Hopkins, leads to a wider range of  y , as we read off from the middle panel of Figure 10.Nevertheless,  w from FIRE-2 simulations also results in an MZGR that can reproduce the observed data, but the preferred value of  y for low-mass galaxies is lower as compared to the case above with  w from Hayward & Hopkins.Specifically, models with  y ∼ 0.1 − 0.2 best reproduce metallicity gradients in low-mass galaxies in all the three IFU surveys in this case.The preferred values of  y for massive galaxies remain similar to the case above.
Finally, we examine the case where we scale  w following the EAGLE simulations (Mitchell et al. 2020a).We plot the resulting model MZGR in Figure 9, and the corresponding diagnostic plot for the dimensionless ratios in the right panel of Figure 10.For 10 9 M ⊙ <  ★ < 10 10.7 M ⊙ , the MZGR from EAGLE is similar to that we obtain from FIRE-2 results above.The difference in the range of gradients at  ★ < 10 9 M ⊙ between the two is due to a non-zero P in FIRE-2.Further, the best match between the observed and the EAGLE MZGR also leads to  y ∼ 0.1 − 0.2 for low-mass galaxies.This is not surprising given  w from both these simulations for low-and intermediate-mass galaxies is within a factor of few (Figure 6).While the range of metallicity gradients predicted for massive galaxies is also the same, the underlying mechanism is qualitatively different.As compared to Figure 8, we see that (P + A) > S for some gradients at large  ★ .The reason behind this is the fact that EAGLE includes AGN feedback that boosts  w above 1 −  R,inst for massive galaxies.As a consequence, EAGLE allows  y to be as low as possible, which in turn decreases S. We confirm this from the range of  y and S we observe from the right panel of Figure 10.
We can summarise the overall conclusions from using different scalings of  w as follows: (1.) there is a transition in the dominant process setting metallicity gradients in local galaxies as a function of  ★ ; gradients are set by P (advection) and A (accretion) in lowmass galaxies, and S (production) begins to play an important role in setting the gradients in intermediate-mass and massive galaxies, (2.) discrepancies in  w that exist for massive galaxies due to AGN feedback lead to different mechanisms setting the high-mass end of the MZGR, (3.) mass-loading is less critical for the MZGR as compared to metal-loading, and (4.) low-mass galaxies seem to prefer a low  y irrespective of our adopted scaling of  w , implying winds in low-mass galaxies are metal-enriched compared to their ISMs.

Predicted scaling of metal enrichment of galactic winds
Based on this analysis, we can now predict a scaling of  y that describes the preferential metal enrichment of galactic winds using  w from the three studies we discuss above in Section 5.2.Such an exercise is similar in spirit to Peeples & Shankar (2011) who predict the scaling of the ratio Z w /Z based on the best match between their chemical evolution model and the MZR from Sloan Digital Sky Survey (SDSS, Tremonti et al. 2004).For this purpose, we use the observed MZGR from the three IFU surveys as well as their overall mean, and calculate the range of  y needed to reproduce the average metallicity gradient observed per  ★ in the model with  Boxplot showing the predicted scaling of  y that describes the preferential enrichment of galactic winds, as a function of  ★ for local galaxies. y ≈ 1 typically corresponds to wind metallicities being similar to ISM metallicities (equation 26). y ≈ 0 typically corresponds to winds being very metal enriched as compared to the ISM.The predicted  y is based on the best match between the models with different mass-loading factors ( w , see Section 5) and the observed MZGR from the three local IFU surveys (SAMI, MaNGA, CALIFA) as well as their overall mean.Notches denote confidence interval around the median value that is marked in solid black.Errors denote the 5 th and 95 th percentile range.The x-coordinate for each boxplot is the mean  ★ across the IFU surveys.The width of the boxplot is the 1 deviation from the mean  ★ .Pink markers denote upper limits on  y measured by Sharda et al. (2021b, appendix A) based on the observations of galactic wind mass-loading and metal-loading in a sample of five local galaxies by Chisholm et al. (2017Chisholm et al. ( , 2018)).Orange marker denotes our estimate of  y for M31 based on the results of Telford et al. (2019Telford et al. ( ). et al. (2020a)).Specifically, we start with a uniform distribution of  y from 0 to 1; for each  y , we check if the model in question can reproduce the gradients observed in a given survey.If so, then it ends up in our final  y distribution.It is difficult to put strong constraints on the range of  y (especially for massive galaxies) given the family of gradients produced at fixed  y .Nevertheless, this exercise demonstrates the power of using metallicity gradients as a tracer of metal enrichment of galactic winds, and its subsequent impact on the metallicity of the CGM and the IGM.
Figure 11 shows the predicted scaling of  y as a function of  ★ from our analysis.The median value of  y is denoted by the solid black line within the notched boxplots.The errorbars denote the 5 th and the 95 th percentile range.We take the mean and standard deviation of  ★ in different bins in the three IFU surveys to specify the − coordinates and widths for each boxplot, respectively.The key conclusion we draw from Figure 11 is that  y is low and wellconstrained for galaxies with  ★ < 10 9.4 M ⊙ regardless of the IFU survey data, systematics due to metallicity calibrations, adopted scaling of  w , and boundary conditions on the metallicity solution imposed by  1 .The median value of  y increases with  ★ : it is ≈ 0.5 for intermediate-mass galaxies and ≈ 0.8 for massive galaxies.However, the scatter is high, with  y = 0.3 − 1.0 providing a match between the model and the observed MZGRs.
We also overplot  y for a sample of five local galaxies for which mass-and metal-loading were measured by Chisholm et al. (2017Chisholm et al. ( , 2018) ) from Si iv and Si ii absorption features in the UV using HST/COS spectra.The corresponding  y values were obtained by Sharda et al. (2021b, see their appendix A). 8 The predictions for  y are in good agreement with all the galaxies in the Chisholm et al. sample, however the least massive galaxy (Mrk 1486) prefers a slightly higher  y as compared to the predictions.Cameron et al. (2021) estimate  y ≈ 0.8 − 0.95 from measurements of ionized gas outflows for Mrk 1486.Since Chisholm et al. and Cameron et al. only analyze the ionized phase of winds in these galaxies, it is likely that the  y estimated for these galaxies is an upper limit, as the hot, X-ray emitting phase can also carry a substantial amount of metals (Lopez et al. 2020(Lopez et al. , 2023)).In addition to these data, we also expect  y ∼ 1 based on the results of Telford et al. (2019) for the local massive galaxy M31 with  ★ ≈ 10 11 M ⊙ .Our predicted scaling of  y is in qualitative agreement with the scaling of Z w /Z extracted by Peeples & Shankar (2011) from the MZR.It is evident that observational measurements of  y in diverse galaxies will be key to solving the puzzle of the role of galactic winds in driving integrated as well as spatially-resolved gas-phase metal distribution in galaxies.

CONCLUSIONS AND FUTURE OUTLOOK
In this work, we create a model of spatially-resolved gas-phase metallicities in galaxies by combining models of turbulent galactic discs in pressure and energy balance (KBFC18, Ginzburg et al. 2022) with a model for ISM metallicities (Sharda et al. 2021a).The evolution of gas-phase metallicity profiles in our model is dictated by both the large-scale properties of galaxies and properties local to the ISM.We include a comprehensive treatment of metal dynamics (advection with radial gas flows and diffusion due to turbulence powered by star formation feedback, gas transport, and cosmic accretion), and allow for preferential enrichment of galactic winds whereby wind metallicities can be higher than the ISM metallicity (see the schematic in Figure 1).Crucially, the key parameters in our model are constrained by both local and global equilibrium, which ensures that metals are treated self-consistently with the gas and stars in galactic discs.Previous works have shown that our model also reproduces the observed relationship between  SF −  g , Σ SFR −Σ g ,  g −  SF , and  ★ −Z.Table 3 highlights the key differences in our work as compared to previous models, which are important to selfconsistently model spatially-resolved metallicities and metallicity gradients.
We compare the results of our model with observed metallicity gradients in local galaxies.We show that when the gas mass and velocity dispersion are calculated self-consistently, only models where turbulence is driven by a combination of feedback, transport, and accretion can reproduce the local mass-metallicity gradient relation (MZGR; see Figure 4).Turbulence driven only by feedback or accretion produces metallicity gradients that are flatter than that observed in massive galaxies.Metal transport, if active (due to radial gas flows), plays an important role in setting metallicity gradients in low-mass galaxies.
The prescriptions of wind mass-loading we adopt from theoretical works (Hayward & Hopkins 2017;Pandya et al. 2021;Mitchell et al. 2020a) are naturally decoupled from the rest of the model.Regardless of this decoupling, we find strong evidence for the preferential metal enrichment of winds in low-mass galaxies, in line with earlier works based on modeling the MZR (Dalcanton 2007;Peeples & Shankar 2011;Sharda et al. 2021b;Tortora et al. 2022).We also find metal loading to be more critical than mass loading in driving the MZGR.The key impact of metal-loaded winds in low-mass galaxies is to reduce the overall level of metal content in the ISM, which is necessary to reproduce the observed flatness in metallicity gradients in these galaxies without violating other gasphase scaling relations.However, the extent of enrichment of winds in intermediate-mass and massive galaxies remains unclear because metal-loaded winds are less important (compared to feedback and accretion) for driving metallicity gradients in these galaxies (Figure 11).
The current sample of direct measurements of wind metallicities and outflow rates is very limited and often only probes a single phase of galactic winds.Given the importance of winds in explaining both the MZR and the MZGR, it is crucial to obtain direct measurements of mass, metal, and energy loading in diverse galaxies.JWST/NIRSpec and VLT/MAVIS are expected to deliver resolved metallicity measurements for a large number of galaxies at high redshifts (Rigaut et al. 2021;Rigby et al. 2023).On the theoretical front, it is desirable to self-consistently model the CGM and the ISM so that the decoupling that currently exists between the two in chemical evolution models such as ours can be removed (Carr et al. 2022;Pandya et al. 2022).Equally important is understanding multiphase metal mixing at the CGM -ISM interface (Kim et al. 2020;Schneider et al. 2020;Vĳayan et al. 2023) that occurs at pc scales not currently resolved in cosmological simulations (Gentry et al. 2019).Such a unification of ISM and CGM models can also provide a self-consistent formalism for galactic fountains that is missing from models like ours.Particularly interesting is the discovery of inverted gradients in galaxies across redshift (e.g., Cresci et al. 2010;Curti et al. 2020;Wang et al. 2019Wang et al. , 2022)).Explaining the origin of inverted gradients within the context of gas regulator models remains a critical task for future theoretical work.several non-linear terms to model spatial variations in  sf , which is why we refrain from adopting it as our standard choice.
As a workaround, we now use the Krumholz (2013) model to find  sf where we estimate the input metallicity from the MZR.We adopt the value of  sf at  = /2 that roughly represents the radiallyaveraged  sf across the disc.Figure A1 shows the resulting MZGR for the fiducial model we plot in the bottom panel of Figure 4. We also plot the corresponding dimensionless ratios that govern the metallicity in Figure A2.We find that transport is now active even in intermediate-mass galaxies because the average  sf is lower than the one we use in the main text, which lowers  sf .However, the gradients produced by the model for  ★ ≈ 10 10.6 M ⊙ are slightly shallower than that observed, because a lower  sf decreases S (cf.equation 24).

APPENDIX B: IMPACT OF UNCERTAINTIES IN THE EFFICIENCY OF ACCRETION-INDUCED TURBULENCE
The efficiency with which accreting gas can convert its kinetic into driving turbulence in the disc is described by the parameter  a in our work.Throughout the main text, we have assumed  a = 0.2 for all galaxies.However, as Ginzburg et al. (2022) point out,  a is essentially a free parameter as no constraints or measurements are available for it (except from the analysis of one low-mass galaxy in the IllustrisTNG simulations by Forbes et al. 2023 that gives  a = 0.2).It is difficult to constrain  a because it depends on the angular momentum of the accreting gas as well as its clumpiness (Mandelker et al. 2018(Mandelker et al. , 2020)).In this section, we explore higher values of  a = 0.6 and  a = 1.0 to understand its impact on the MZGR, noting that  a = 1 in particular is an extreme value since it implies all the kinetic energy of accreting gas goes into driving turbulent motions.We only show results for the fiducial model where  g is driven by feedback, accretion, and transport.
The top and bottom panels of Figure B1 plots the MZGRs for the case with  a = 0.6 and 1, respectively.The key difference between these plots and the ones we present in the main text is that the scatter due to  y at large  ★ is higher.We can understand this trend with the help of  g -higher  a leads to more turbulence due to gas accretion, and thus, higher  g .Since S ∝  2 g but A ∝  3 g , A decreases much more as compared to S, and a stronger S can drive steeper, negative metallicity gradients even at high  ★ .However, it is unlikely that  a is high in local massive galaxies, especially if a significant fraction of accreting gas is co-rotating with the disc (e.g., Danovich et al. 2015;Trapp et al. 2022).

APPENDIX C: EQUATIONS FOR CENTRAL METALLICITY
Following from Section 2.3.2,we provide here the resulting equations for the central metallicity Z  0 .For the case where Z  0 is set by the competition between source and accretion (e.g., massive galaxies), If Z  0 is set by the competition between advection and diffusion (e.g., low-mass galaxies), then, in the Toomre regime of star forma- (C4)

APPENDIX D: SYSTEMATIC DIFFERENCES IN MZGR DUE TO METALLICITY DIAGNOSTICS
The gas-phase metallicity is typically measured using a combination of emission lines commonly found in the spectra of H ii regions.Various diagnostics of ISM metallicity exist that provide a way to convert emission line ratios to metallicity on the absolute scale (12+log 10 O/H; see the reviews by Kewley et al. 2019 andMaiolino &Mannucci 2019).These diagnostic have systematic differences that are partially responsible for the observed scatter in metallicity gradients (e.g., Poetrodjojo et al. 2021).In this section, we study the impacts of these differences on our interpretation of the MZGR.
To do so, we plot the MZGR obtained from the MaNGA survey (Mingozzi et al. 2020) but for three different calibrators: PP04, based on the N ii to H ratio, (Pettini & Pagel 2004), M08, based on O iii and O ii (Maiolino et al. 2008), and IZI (Blanc et al. 2015).Briefly, the modified version of IZI that Mingozzi et al. Figure D1 shows a comparison of the fiducial model MZGR where  g is driven by feedback, accretion, and transport with the MaNGA MZGR derived from these metallicity calibrations.The overall trend in the MZGR remains the same for all the three calibrations, as already noted by (Mingozzi et al. 2020, fig. 12).In fact, the scatter within different surveys (MaNGA, SAMI, CALIFA) is somewhat larger than that within the same survey but different calibrators.The comparison between the model and the observed MZGR also remains unchanged.Thus, we find that systematic differences arising from different metallicity calibrations has no appreciable impact on the comparison between the model and the data.
This paper has been typeset from a T E X/L A T E X file prepared by the authors.
, variations with radius (Leroy et al. 2018; Fisher long-lived stars: 1 −  !,#t e x i t s h a 1 _ b a s e 6 4 = " w b p 2 0 W 2 n r g w u G c e 2 f u P U H C m d K O 8 2 2 t r K 6 t b 2 x W t q r b O 7 t 7 + / b B Y U f F q S S 0 T W I e y 1 6 A F e V M 0 L Z m m t N e I i m O A k 6 7 w e R m 6 n c f q V Q s F g 8 6 S 6 g X 4 Z F g I S N Y G 8 m 3 a 4 N h r N G d n w 9 k l I / O M S F F 4 d t 1 p + H M g J a J W 5 I 6 l G j 5 9 p d 5 h a w b p 2 0 W 2 n r g w u G c e 2 f u P U H C m d K O 8 2 2 t r K 6 t b 2 x W t q r b O 7 t 7 + / b B Y U f F q S S 0 T W I e y 1 6 A F e V M 0 L Z m m t N e I i m O A k 6 7 w e R m 6 n c f q V Q s F g 8 6 S 6 g X 4 Z F g I S N Y G 8 m 3 a 4 N h r N G d n w 9 k l I / O M S F F 4 d t 1 p + H M g J a J W 5 I 6 l G j 5 9 p d 5 h a w b p 2 0 W 2 n r g w u G c e 2 f u P U H C m d K O 8 2 2 t r K 6 t b 2 x W t q r b O 7 t 7 + / b B Y U f F q S S 0 T W I e y 1 6 A F e V M 0 L Z m m t N e I i m O A k 6 7 w e R m 6 n c f q V Q s F g 8 6 S 6 g X 4 Z F g I S N Y G 8 m 3 a 4 N h r N G d n w 9 k l I / O M S F F 4 d t 1 p + H M g J a J W 5 I 6 l G j 5 9 p d 5 h a w b p 2 0 W 2 n r g w u G c e 2 f u P U H C m d K O 8 2 2 t r K 6 t b 2 x W t q r b O 7 t 7 + / b B Y U f F q S S 0 T W I e y 1 6 A F e V M 0 L Z m m t N e I i m O A k 6 7 w e R m 6 n c f q V Q s F g 8 6 S 6 g X 4 Z F g I S N Y G 8 m 3 a 4 N h r N G d n w 9 k l I / O M S F F 4 d t 1 p + H M g J a J W 5 I 6 l G j 5 9 p d 5 h a R e r N d Z a 8 6 a z + y i H 7 D e P g G o w 5 r r < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " G 3 9 m J z e p J o 5 h c Y K F n 7 P C x d 6 n 2 R e r N d Z a 8 6 a z + y i H 7 D e P g G o w 5 r r < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " G 3 9 m J z e p J o 5 h c Y K F n 7 P C x d 6 n 2

Figure 1 .
Figure 1.Bathtub model of chemical evolution, adapted from Lilly et al. (2013, figure2), and modified to reflect the spatially-resolved metallicity model presented in this work.The galaxy is modeled as a two component reservoir: stellar (red) and gas-phase (blue).The gas reservoir is made up of H, He and metals (denoted by the ISM metallicity ).We consider marginally unstable galactic discs (quantified by the Toomre  parameter) in pressure and energy balance.The evolution of gas mass,  g , is described by equation 7. It depends on the gas accretion rate (  g,acc ), star formation rate (  SF ), radial gas flows (or, gas transport,  trans ), and galactic winds ( w  SF , where  w is the mass-loading factor).Gas turbulence (quantified by the gas velocity dispersion,  g ) is driven by energy injected into the ISM via accretion, star formation feedback, and gas transport.Metals are produced as star formation leads to supernovae (SNe).Some fraction of the metals produced is returned almost instantaneously to the ISM (  R,inst ) whereas the rest is locked in long-lived stars.Metals are also ejected via galactic winds with metallicity  w , which can be greater than or equal to .  y denotes the preferential metal enrichment of winds.Once produced, metals are also advected with radial gas flows, and diffused into the ISM due to turbulence and thermal instabilities.© AAS.Reproduced with permission.

Figure 2 .
Figure 2. Gas velocity dispersion,  g , as a function of stellar mass, at  = 0.The three curves corresponds to models where turbulence is driven by either feedback due to star formation, gas accretion, gas transport, or a combination of the three.

Figure 3 .
Figure 3. Top panel: ISM metallicity profiles (Z = / ⊙ ) as a function of the dimensionless radius ( = / 0 with  0 = 1 kpc) for the fiducial case of a massive galaxy ( ★ = 10 10.8 M ⊙ ) with  y = 1.0 (solid) and 0.1(dotted), respectively. y describes the preferential metal enrichment of galactic winds.The family of curves for each  y arise from constraints on the constant of integration  1 in the solution for metallicity.Colorbar represents the difference between the metal equilibration timescale,  eqbm , and the gas depletion timescale,  dep , with increasingly red color representing larger deviations from equilibrium.Bottom panel: Same as the top panel but for a low-mass galaxy with  ★ = 10 9 M ⊙ .

Figure 4 .
Figure 4.The local mass-metallicity gradient relation (MZGR) as predicted from the fiducial model.It plots metallicity gradients (∇ log 10 Z) as a function of the stellar mass ( ★ ).The colorbar indicates the ratio of the strength of processes that create flat metal distributions (P and A) to those that create steep metallicity gradients (S) in equilibrium, as defined in Section 2. The top panel plots the model MZGR for the case where only gas accretion and transport drive turbulence in the galaxy.The middle panel corresponds to the case where only star formation feedback and gas transport drive turbulence.The bottom panel corresponds to the case where all three factors are included as drivers of gas turbulence.The scatter in the model at fixed  ★ corresponds to the results of varying   , the parameter that describes the preferential ejection of metals via galactic winds, from 0.1 − 1, as indicated.Overplotted (orange markers) are average metallicity gradients in bins of  ★ obtained from IFU surveys -MaNGA(Belfiore et al. 2017;Mingozzi et al. 2020), CALIFA(Sánchez et al. 2014;Sánchez-Menguiano et al. 2016), and SAMI(Poetrodjojo et al. 2018), as well as mean of the three (grey markers).

Figure 5 .
Figure5.Left panel: The three dimensionless ratios (P, A, and S) that set the metal distribution at  = 0 in the fiducial model presented in Section 4, as a function of stellar mass, for the case where turbulence is driven by star formation feedback and gas transport.The value of S is directly proportional to the parameter  y that describes the preferential enrichment of galactic winds (see equation24), so for this parameter we plot a vertical range corresponding to varying  y from 0.1 − 1, as indicated by the colour bar.Metal transport due to radial gas flows (P) is only active for the lowest and most massive galaxies in this case.The dashed line denotes a value of unity, with values of a parameter less than unity indicating that turbulent diffusion is a more important process than the one described by that parameter.Middle panel: Same as the left panel but now turbulence is driven by gas accretion and transport.Metal advection is for most galaxies in this case.Right panel: Same as the other two panels but turbulence is now driven by all three mechanisms -feedback, accretion, and transport.Metal advection is active only in the lowest mass galaxies.These diagnostic plots aid in understanding the nature of metal distributions and metallicity gradients in each of the corresponding models presented in Figure4.

Figure 7 .
Figure 7. Same as the bottom panel of Figure 4 but with the mass-loading factor ( w ) estimated from the analytical model of (Hayward & Hopkins 2017).

Figure 8 .
Figure 8. Same as the bottom panel of Figure 4 but with the mass-loading factor ( w ) estimated from the FIRE-2 cosmological zoom-in simulations without AGN feedback (Pandya et al. 2021).

Figure 9 .
Figure 9. Same as the bottom panel of Figure 4 but with the mass-loading factor ( w ) estimated from the EAGLE cosmological simulations including AGN feedback (Mitchell et al. 2020a).

Figure 10 .
Figure10.Left panel: Strength of the dimensionless ratios in the metallicity model (P, A, and S) as a function of  ★ for the case where turbulence is driven by feedback, accretion, and transport, and the mass-loading factor,  w , scales with the gas fraction and  ★ , following an analytical model of star formation-driven winds(Hayward & Hopkins 2017).Colorbar denotes the permitted range of  y that describes the preferential enrichment of galactic winds (see equation 26).Middle panel: Same as the left panel but  w is scaled with  ★ following FIRE-2 cosmological zoom-in simulations without AGN feedback(Pandya et al. 2021).Right panel: Same as the other two panels but  w is scaled with  ★ following EAGLE cosmological simulations that include AGN feedback(Mitchell et al. 2020a).
Figure11.Boxplot showing the predicted scaling of  y that describes the preferential enrichment of galactic winds, as a function of  ★ for local galaxies. y ≈ 1 typically corresponds to wind metallicities being similar to ISM metallicities (equation 26). y ≈ 0 typically corresponds to winds being very metal enriched as compared to the ISM.The predicted  y is based on the best match between the models with different mass-loading factors ( w , see Section 5) and the observed MZGR from the three local IFU surveys (SAMI, MaNGA, CALIFA) as well as their overall mean.Notches denote confidence interval around the median value that is marked in solid black.Errors denote the 5 th and 95 th percentile range.The x-coordinate for each boxplot is the mean  ★ across the IFU surveys.The width of the boxplot is the 1 deviation from the mean  ★ .Pink markers denote upper limits on  y measured bySharda et al. (2021b, appendix A)  based on the observations of galactic wind mass-loading and metal-loading in a sample of five local galaxies byChisholm et al. (2017Chisholm et al. ( , 2018)).Orange marker denotes our estimate of  y for M31 based on the results ofTelford et al. (2019).

Figure A1 .
Figure A1.Same as the fiducial model in the bottom panel of Figure4but for the fraction of star-forming gas,  sf , estimated using the approach laid out inKrumholz (2013).

Figure A2 .
Figure A2.Same as Figure 5 but for the model presented in Figure A1.

Figure B1 .
Figure B1.Same as the fiducial model in the bottom panel of Figure 4 but with the  a = 0.6 (top panel) and  a = 1.0 (bottom panel), implying 60 per cent and 100 per cent of the kinetic energy of accreting gas goes into driving turbulent motions in the disc, respectively.The mass-loading factor ( w ) is estimated from the EAGLE cosmological simulations (Mitchell et al. 2020a).
create uses Bayesian inference to predict the joint posterior probability distribution function (PDF) of the metallicity, ionization parameter, and extinction along the line of sight by comparing the observed emission line ratios with the Dopita et al. (2013) photoionization models.
table 1) to estimate  ★ .Once we determine Navarro et al. (1997) al. (2013)to find the rotational velocity at the outer edge of the disc,   , assuming aNavarro et al. (1997)profile for the dark matter density as(Krumholz et al. 2018, equations  69−71)

Table 1 .
List of parameters in the galaxy evolution model.
(Karpov et al. 2020)a;Gentry et al. 2017Gentry et al. , 2019;;Kim et al. 2017) 2017).Although simulations of clustered supernovae that lead to the formation of superbubbles (e.g.,Keller et al. 2014) suggest a value higher by a factor of few(Kim & Ostriker 2015a;Gentry et al. 2017Gentry et al. , 2019;;Kim et al. 2017), discrepancies up to an order of magnitude exist between different simulations (see the discussion inHu et al. 2023).These discrepancies stem from a mixture of numerical (treatment of shocks and contact discontinuities in Lagrangian versus Eulerian codes) and physical (impact of magnetic fields) issues.Supernova clustering is also sensitive to the local ISM metallicity due to metallicity-dependent cooling and associated expansion of superbubbles(Karpov et al. 2020).It is also not yet clear what frac- Orr et al. 2022)entum drives turbulence in the ISM as compared to driving winds (e.g.,Orr et al. 2022).Given these caveats, we continue to use momentum injection from non-clustered supernovae in our model, but emphasize that further exploration of supernova clustering (especially at non-Solar metallicities) is highly desirable to more accurately model the role of star formation feedback for ISM metal distribution.

Table 2 .
List of parameters in the ISM metallicity model.

Table 4 .
List of parameter values for metallicity gradients in two representative galaxies as shown in Figure3.