key: cord-0320695-2xv82sr9
authors: Sharda, Piyush; Krumholz, Mark R.; Wisnioski, Emily; Forbes, John C.; Federrath, Christoph; Acharyya, Ayan
title: The physics of gas phase metallicity gradients in galaxies
date: 2021-02-02
journal: nan
DOI: 10.1093/mnras/stab252
sha: 46ffa04516ab422490faba6a62298707837ec55a
doc_id: 320695
cord_uid: 2xv82sr9

We present a new model for the evolution of gas phase metallicity gradients in galaxies from first principles. We show that metallicity gradients depend on four ratios that collectively describe the metal equilibration timescale, production, transport, consumption, and loss. Our model finds that most galaxy metallicity gradients are in equilibrium at all redshifts. When normalized by metal diffusion, metallicity gradients are governed by the competition between radial advection, metal production, and accretion of metal-poor gas from the cosmic web. The model naturally explains the varying gradients measured in local spirals, local dwarfs, and high-redshift star-forming galaxies. We use the model to study the cosmic evolution of gradients across redshift, showing that the gradient in Milky Way-like galaxies has steepened over time, in good agreement with both observations and simulations. We also predict the evolution of metallicity gradients with redshift in galaxy samples constructed using both matched stellar masses and matched abundances. Our model shows that massive galaxies transition from the advection-dominated to the accretion-dominated regime from high to low redshifts, which mirrors the transition from gravity-driven to star formation feedback-driven turbulence. Lastly, we show that gradients in local ultraluminous infrared galaxies (major mergers) and inverted gradients seen both in the local and high-redshift galaxies may not be in equilibrium. In subsequent papers in this series, we show that the model also explains the observed relationship between galaxy mass and metallicity gradients, and between metallicity gradients and galaxy kinematics.

Metals act as tracers of the formation and assembly history of galaxies. Tracking their evolution is crucial to understanding the various pathways a galaxy takes while it forms (Freeman & Bland-Hawthorn 2002) . Metals are produced in galaxies through supernovae (Hillebrandt & Niemeyer 2000; Woosley et al. 2002) , asymptotic giant branch (AGB) stars (van Winckel 2003; Herwig 2005) , neutron star mergers (Thielemann et al. 2017) , etc., and are consumed by low mass stars and retained in stellar remnants (Kobayashi et al. 2006; Sukhbold et al. 2016) . Apart from this in-situ metal production and consumption (Pagel & Patchett 1975) , metals can also be lost through outflows in the form of galactic winds (Heckman et al. 1990; Veilleux et al. 2005; Rupke 2018; Chisholm et al. 2018 ), ★ piyush.sharda@anu.edu.au (PS) † mark.krumholz@anu.edu.au (MRK) ‡ emily.wisnioski@anu.edu.au (EW) or transported into the galaxy from the circumgalactic (CGM) and the intergalactic medium (IGM, Prochaska et al. 2017; Tumlinson et al. 2017) , or during interactions with other galaxies, like flybys or mergers (e.g., Torrey et al. 2012; Grossi et al. 2020 ). All of these processes can be classified into four main categories: metal production (through star formation and supernovae), metal consumption (through stellar remnants and low mass stars), metal transport (through advection, diffusion, and accretion), and metal loss (galactic winds and outflows). The distribution of metals within galaxies places important constraints on galaxy formation (Sánchez Almeida et al. 2014 ). One of the strongest pieces of evidence for the inside-out galaxy formation scenario is the existence of negative metallicity gradients (in the radial direction) in both the gas and stars in most galaxies. The presence of such negative radial gradients is easy to understand: in the inside-out scenario, the centre, i.e., the nucleus of the galaxy forms first, and the disc subsequently forms and evolves in time. The nucleus undergoes greater astration, leading to the pres-ence of more metals in the centre as compared to the disc, thus establishing a negative gradient. Such gradients were first observed and quantified through nebular emission lines in H regions in the interstellar medium (ISM) by Aller (1942) , Searle (1971) and Shaver et al. (1983) . The decrease in metallicity is approximately exponential with galactocentric radius (Wyse & Silk 1989; Zaritsky 1992) , yielding a linear gradient in logarithmic space, with units of dex kpc −1 . Since these early works, metallicity gradients have been measured for thousands of galaxies, both in stars and gas (see recent reviews by Kewley et al. 2019a , Maiolino & Mannucci 2019 and Sánchez 2020 . The stellar metallicity gradients are typically characterised by the abundance of iron, and are written in the form log 10 (Fe/H)/ , whereas the gas phase metallicity gradients are characterised by the abundance of oxygen, and written as log 10 (O/H)/ , where is the galactocentric radius. Hereafter, we will only discuss the gas phase metallicity gradients in galaxies.

In the local Universe, samples of metallicity gradient measurements have been dramatically expanded by three major surveys: CALIFA (Calar Alto Legacy Integral Field Area, Sánchez et al. 2012) , MaNGA (Mapping nearby Galaxies at Apache Point Observatory, Bundy et al. 2015) , and SAMI (Sydney-AAO Multi-object Integral-field spectrograph, Bryant et al. 2015) . These surveys show that local galaxies contain predominantly negative metallicity gradients, with typical values ranging between 0 and −0.1 dex kpc −1 . Measurements in high-( 3) galaxies are more challenging, and the samples are correspondingly smaller, but rapidly expanding (Queyrel et al. 2012; Swinbank et al. 2012; Jones et al. 2013; Stott et al. 2014; Troncoso et al. 2014; Leethochawalit et al. 2016; Wuyts et al. 2016; Förster Schreiber et al. 2018; Wang et al. 2019 Wang et al. , 2020 Curti et al. 2020; Simons et al. 2020; Gillman et al. 2021) .

Theoretical efforts to understand these observations are still in their infancy, with most effort thus far dedicated to understanding galaxies' global metallicities (e.g., Erb 2008; Peeples & Shankar 2011; Davé et al. 2012; Lilly et al. 2013; Dayal et al. 2013; Hunt et al. 2016; Wang & Lilly 2020; Furlanetto 2021 ) rather than their metallicity gradients. Early work on metallicity gradients was tuned to reproduce the present-day Milky Way metallicity gradient (e.g. Chiappini et al. 1997 Chiappini et al. , 2001 Prantzos & Boissier 2000) , and thus offers relatively little insight into how metallicity gradients have evolved over cosmic time and in galaxies with differing histories. More recent work has attempted to address the broader sample of galaxies, using physical models that include a range of processes: cosmological accretion, mass-loaded galactic winds, in situ metal production by stars, and radial gas flows (e.g., Mott et al. 2013; Jones et al. 2013; Ho et al. 2015; Carton et al. 2015; Kudritzki et al. 2015; Pezzulli & Fraternali 2016; Belfiore et al. 2019; Kang et al. 2021 ). However, these works generally suffer from a problem of being under-constrained: the models generally involve multiple free functions (e.g., the radial inflow velocity or mass loading factor as a function of radius and time) that are not constrained by any type of independent physical model, and fits of these free functions to the data are often non-unique, leaving in doubt which physical processes are important in which types of galaxies.

Moreover, not all models include all possible physical processes, making them difficult to compare. For example, some of the models above assume that galactic wind metallicities are equal to ISM metallicities, contrary to observational evidence (e.g., Martin et al. 2002; Strickland & Heckman 2009; Chisholm et al. 2018 ). Many models do not include metal transport processes like advection (flow of metals that are carried by a bulk flow of gas) or turbulent diffusion (metal flow that occurs due to turbulent mixing of gas with non-uniform metal concentration; e.g., de Avillez & Mac Low 2002; Greif et al. 2009; Yang & Krumholz 2012; Aumer et al. 2013; Kubryk et al. 2013; Petit et al. 2015; Armillotta et al. 2018; Rennehan et al. 2019) , despite modelling showing that these effects play an important role in setting metallicity gradients (Forbes et al. 2012 ). These problems have been partly alleviated by recent radiallyresolved semi-analytic models (Kauffmann et al. 2013; Fu et al. 2013; Forbes et al. 2014a Forbes et al. , 2019 Henriques et al. 2020; Yates et al. 2020 ) and cosmological simulations with enough resolution to capture disc radial structure (Pilkington et al. 2012; Gibson et al. 2013; Ma et al. 2017; Tissera et al. 2019; Hemler et al. 2020; Bellardini et al. 2021) , which do at least attempt to model the dynamics of gas in galaxies self-consistently. The general result of these models (as summarised in Figure 8 of Curti et al. 2020) , is that there is only a mild evolution in metallicity gradients between 0 ≤ ≤ 4, with a slight steepening toward the present day. However, the physical origin of these results, and insights from the numerical results in general, have yet to be distilled into analytic models that we can use to understand the overall trends in the data. Thus, to this date, we lack a model that can explain the occurrence of metallicity gradients in a diverse range of galaxies from first principles. This leaves many interesting questions around gas phase metallicity gradients unanswered.

Motivated by this, we present a new theory of gas phase metallicity gradient evolution in galaxies from first principles. As with all theories of metallicity gradients, ours requires a galaxy evolution model that describes the gas in galactic discs as an input. For the purposes of developing the theory, we use the unified galactic disc model of , which has been shown to reproduce a large number of observations of gas kinematics relevant to metallicity gradients, including the radial velocities of gas in local galaxies (Schmidt et al. 2016) , the correlation between galaxies' gas velocity dispersions and star formation rates (SFRs; e.g., Johnson et al. 2018; Yu et al. 2019; Varidel et al. 2020) , and the evolution of velocity dispersion with redshift (e.g., Übler et al. 2019) . However, our metallicity model is a standalone model into which we can incorporate any galaxy evolution model. In this paper, we present the basic formalism and results of the model, and use them to explain the evolution of metallicity gradients with redshift; in two follow-up papers we first use the model developed here to explain the dependence of metallicity gradients on galaxy mass that is observed in the local Universe (Sharda et al. 2021a) , and then use the model to predict the existence of a correlation between galaxy kinematics and metallicity gradients, which we validate against observations (Sharda et al. 2021b ).

We arrange the rest of the paper as follows: Section 2 describes the theory of metal evolution, Section 3 describes the equilibrium metallicity gradients generated by the theory for different types of galaxies both in the local and the high-Universe, Section 4 combines the local and high-predictions of the model to describe the cosmic evolution of metallicity gradients, and Section 5 discusses the limitations of the model, including special cases where the metallicity gradients may not be in equilibrium and thus the model may or may not apply. 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. 2009 ), Hubble time H(0) = 13.8 Gyr (Planck Collaboration et al. 2018) , and follow the flat ΛCDM cosmology: Ω m = 0.27, Ω Λ = 0.73, ℎ = 0.71, and 8 = 0.81 (Springel & Hernquist 2003) .

For convenience, we collect all of the symbols we define in this section in Table 1 and Table 2 .

Let us start by defining to be the volume density of metals at some point in space; this is related to the metallicity, , and gas density, , by =

.

(1)

The density of the metals can change due to transport -via advection with the gas or diffusion through the gas -and due to sources and sinks (e.g., production of new metals by stars or consumption during star formation). The conservation equation for metal mass is then

Here, v is the gas velocity, j is the flux density of metals as a result of diffusion, and represents the source and sink terms. The central assumption of diffusion is that the diffusive flux is proportional to minus the gradient of the quantity being diffused (e.g., Yang & Krumholz 2012; Krumholz & Ting 2018) . The slight subtlety here is that what should diffuse is not the density of metals, but the concentration of metals, i.e., the flux only depends on the gradient of . We can therefore write down the diffusive flux as

where is the diffusion coefficient (with dimensions of mass/length 2 ). Inserting this into the continuity equation, we now have

We can now specialise to the case of a disc. Firstly, we assume that the disc is thin, so we can write in terms of the surface density as Σ = ∫ . We choose our coordinate system so that the disc lies in the plane. Integrating all quantities in the z direction, the equation of mass conservation becomes

where Σ is the metal surface density, ∇ contains only the derivatives in the plane, and = ∫ . Assuming cylindrical symmetry, this reduces to,

where represents the radial component of the velocity. It is helpful to rewrite the velocity in terms of the inward mass flux across the circle at radius , which is

where we have adopted a sign convention whereby > 0 corresponds to inward mass flow 1 . This gives

Similarly, since star formation is the process that is responsible 1 This is the opposite of the sign convention used in Forbes et al. (2012 Forbes et al. ( , 2014b ), but consistent with the one used in .

for the source term, it is convenient to parameterize in terms of the star formation rate. We adopt the instantaneous recycling approximation (Tinsley 1980) , whereby some fraction, R,inst , of the mass incorporated into stars is assumed to be left in long-lived remnants (compact objects and low-mass stars), and the remainder of the mass is returned instantaneously to the ISM through Type II supernovae, enriched by newly formed metals with a yield . Under this approximation, we have

where Σ ★ is the star formation rate surface density. The last term in equation 9 represents loss of metals into a galactic wind; here is the mass loading factor of the wind (i.e., the wind mass flux is Σ ★ ) and is the metallicity of the wind. Following Forbes et al. (2019, equation 41) , we further parameterize the wind metallicity as

where the 1 − R,inst limit specifies the minimum mass that can be ejected if some metals are ejected directly after production. The parameter , which is bounded in the range 0 ≤ ≤ 1, specifies the fraction of metals produced that are directly ejected from the galaxy before they are mixed into the ISM. So, = 0 corresponds to a situation when the metallicity of the wind equals the metallicity of the ISM, whereas = 1 corresponds to the regime when all the metals produced in the galaxy get ejected in winds. Forbes et al. (2014a,b) introduced to relax the assumption that metals fully mix with the ISM before winds are launched, so that = . A number of authors have shown that this assumption leads to severe difficulties in explaining observations, particularly in lowmass systems (Pilyugin 1993; Marconi et al. 1994; Mac Low & Ferrara 1999; Recchi et al. 2001 Recchi et al. , 2008 Martin et al. 2002; Robles-Valdez et al. 2017) .

We can further simplify by writing down the continuity equation for the total gas surface density Σ , which is equation 8 with fixed to unity and = 0, with an additional term for cosmic accretion, 2

where Σ cos is the cosmic accretion rate surface density onto the disc (Oppenheimer et al. 2010; Benson 2010 ). If we now use this to evaluate

2 Note that equation 11 is identical to equation 1 of Forbes et al. (2019) except that Forbes et al. adopt instantaneous recycling only for Type II supernovae, and not for metals returned on longer timescales (e.g., Type Ia or AGB winds). While this approach is feasible in simulations and semianalytic models, it renders analytic models of the type we present here intractable. However, this does not make a significant difference for our work because the most common gas phase metallicity tracer, O, comes almost solely from Type II supernovae. One area where our approach might cause concern is at high redshift, where the gradients are often measured through the [N ]/H emission line ratio, because most of the N comes from AGB stars and is released over Gyr or longer timescales (Herwig 2005) . Table 1 . List of fiducial parameters in the model that are common to all galaxies. All of these parameters are adopted from Table 1 in , except for and R,inst , which we adopt from Forbes et al. (2019) . We refer to , which is bounded in the range 0 ≤ ≤ 1, as the yield reduction factor. Note that R,inst only appears in , implying that metals locked in stars are unimportant for the radial profile of metallicity as long as > 1 − R,inst .

From left to right, we can interpret the terms in equation 12 as follows: the first is the rate of change in the metallicity at fixed gas surface density; the second represents the change due to advection of metals through the disc; the third represents the change due to diffusion of metals; finally, the terms on the right hand side are: (1.) the change in metallicity due to metal production in stars, with an effective yield that is reduced relative to the true yield by the factor , and (2.) the change in metallicity in the disc due to cosmic accretion of metal-poor gas.

The term represents the factor by which the effective metal yield is reduced because some fraction of metals directly escape the galaxy before they mix with the ISM. Higher values of imply that metals are well-mixed into the ISM, whereas lower values imply that the yield is significantly reduced by preferential ejection of unmixed metals. However, does not equate to the mass loading factor : galaxies with heavily mass-loaded winds (high ) may still have close to unity if metals mix efficiently before the winds are launched; conversely, galaxies with weakly mass-loaded winds (low ) may still have small if those winds preferentially carry away metals. We discuss the possible range of values for in more detail in Section 2.2 (see also, Sharda et al. 2021a) .

At this point it is helpful to non-dimensionalise the system. We choose a fiducial radius 0 , which we will later take to be the inner edge of the disc where the bulge begins to dominate; for now, however, we simply take 0 as a specified constant. We measure position in the disc with the dimensionless variable = / 0 and time with = Ω 0 , where Ω 0 is the angular velocity at 0 . We further write out the profiles of gas surface density, diffusion coefficient, star formation surface density and cosmic accretion rate surface density as Σ = Σ 0 ( ), = 0 ( ), Σ ★ = Σ ★0 ★ ( ), and Σ cos = Σ cos0 ★ ( ) respectively. Here, the terms subscripted by 0 are the values evaluated at = 0 , and ( ), ( ), ★ ( ), and ★ ( ) are dimensionless functions that are constrained to have a value of unity at = 1. Note that, in principle, we could introduce a similar scaling function for ; we do not do so because both observations (Schmidt et al. 2016 ) and theoretical models suggest that, in steady state, is close to constant with radius within a galactic disc. We express the metallicity as Z = / .

Using these definitions, we can rewrite equation 12 as a form of the Euler-Cauchy equation (Arfken 1985; Kreyszig et al. 2011) ,

In the above equation, we have suppressed the -dependence of , , ★ and ★ for compactness, and we have defined,

The four quantities T , P, S and A have straightforward physical interpretations: T is the ratio of the orbital and diffusion timescales, P is the Péclet number of the system, which describes the relative importance of advection and diffusion in fluid dynamics (e.g., Patankar 1980) , S measures the relative importance of metal production (the numerator) and diffusion (the denominator), and A measures the relative importance of cosmic accretion and diffusion. T dictates the time it takes for a given metallicity distribution to reach equilibrium in a galaxy, whereas the other three quantities govern the type and strength of the gradients that form in equilibrium. We will only look for the steady-state or equilibrium solutions to equation 14, so we drop the Z/ term. This approach is reasonable because, as we will show below, the equilibration timescale for metals is less than the Hubble time, H(z) , for most galaxies. In our model, the time it takes for the metallicity gradient to approach an equilibrium state, eqbm , is based on the time it takes for the metal surface density to adjust to changes in metallicity triggered by each of the terms in equation 14,

If eqbm > H(z) , the metallicity gradient in the galaxy cannot attain equilibrium within a reasonable time, and the model we present below does not apply. While this is a necessary condition for metal equilibrium, it may not be sufficient. This is because if input parameters to the metallicity model (e.g., accretion rate, surface density, etc.) change on timescales much shorter than H(z) , the equilibrium of metals will depend on that timescale. For a steady-state model like ours, it is safe to assume this is not the case, since the input galactic disc model in the next Section we use is an equilibrium model. We discuss this condition in more detail in Section 3 and Section 5.2 where we also compare eqbm with the molecular gas depletion time that dictates the star formation timescale. The accretion of material from the CGM can also impact metallicity in the galactic disc (Wise et al. 2012; Forbes et al. 2012; Taylor & Kobayashi 2015; Tumlinson et al. 2017; Schaefer et al. 2019) . While this is an important consideration, in the absence of which 'closed-box' galaxy models overestimate metallicity gradients (e.g., Dalcanton 2007; Zahid et al. 2013; Kudritzki et al. 2015) , it typically adds a floor metallicity at the outer edge of the galactic disc, and is of concern for simulations where the entire (star-forming as well as passive) disc up to tens of kpc is considered. CGM metallicity can also be important for long term 0.1 − 1 H(z) wind recycling (Henriques et al. 2013; Pandya et al. 2020 ), which we do not take into account in this model. As we show later in Section 2.3, we make use of this effect only as a boundary condition on the metallicity at the outer edge of the disc, and do not include it directly in the evolution equation.

This completes the basic formulation of the theory of metallicity gradients in galaxies. To further solve for the equilibrium metallicity, we now need a model of the galactic disc. We use the unified galactic disc model of for this purpose. However, we remind the reader that the metallicity evolution described by equation 14 can be used with other galactic disc models as well.

We use the unified galactic disc model of to further solve for metallicity. This model self-consistently incorporates all of the ingredients that we require as inputs: profiles of Σ , , and Σ ★ , and the relationship between them. We refer the reader to for full details of the model, and here, we simply extract the portions that are relevant for this work.

Firstly, note that the angular velocity at 0 is simply,

where is the rotational velocity of gas in the galactic disc. We can solve for the gas surface density Σ by requiring that the Toomre parameter for stars and gas is close to 1; formally, following Forbes et al. (2014a) , we take = min , where min ≈ 1 − 2 is the minimum parameter below which gravitational instability prevents discs from falling (e.g., Martin & Kennicutt 2001; Martin et al. 2002; Genzel et al. 2010; Meurer et al. 2013; Romeo & Falstad 2013; Inoue et al. 2016; Stott et al. 2016; Romeo & Mogotsi 2017) . This can be re-written as (Krumholz et al. 2018, equation 8 

where is the Toomre parameter for the gas alone, and , is the effective gas fraction in the disc (Krumholz et al. 2018, equation 9) , which, based on the estimates of Σ (McKee et al. 2015) and gas velocity dispersion (Kalberla & Kerp 2009 ) is ≈ 0.5 in the Solar neighbourhood. Writing down the Toomre equation (Toomre 1964) , this becomes,

Here, is the epicyclic frequency given by = √︁ 2( + 1)Ω = √︁ 2( + 1) / , where is the index of the rotation curve given by = ln / ln . Following and results from time-dependent numerical solutions for energy equilibrium in galactic discs (Forbes et al. 2014a) , we can assume that in the steady-state, and are independent of radius. Thus, we obtain

This solution provides a 1/ dependence for Σ that is somewhat at odds with observations that find an exponential dependence of Σ (Bigiel & Blitz 2012) . However, these observations trace the entire disc (using CO as well as H ) and the Σ profiles show a large scatter in the inner disc, which is the focus of our work. Given these findings, we cannot conclude that a 1/ profile of Σ is unrealistic, and therefore continue to use it for our work. The quantities Σ 0 and ( ) that we defined in Section 2.1 are thus given by

We can express the diffusion coefficient due to turbulent diffusion as ≈ ℎ /3, where ℎ represents the gas scale height (Karlsson et al. 2013; Krumholz & Ting 2018) given by (Krumholz et al. 2018, equations 24 and 27) ,

where Σ ★ and ★ is the stellar surface density and velocity dispersion, respectively, and − 1 is the ratio of gas to stellar Toomre parameters. This gives

Hence, the factors 0 and ( ) that we defined in Section 2.1 are given by

Thus, the product 0 Σ 0 ∝ 3 / describes an effective metal flow rate in the disc due to diffusion.

To derive Σ ★ , we can use equations 31 and 32 of ,

where ff is the star formation efficiency per free-fall time (Krumholz & McKee 2005; Federrath 2013; Sharda et al. 2018 Sharda et al. , 2019 , sf is the fraction of gas in the cold, molecular phase that is not supported by thermal pressure, and thus forms stars (Krumholz et al. 2008 (Krumholz et al. , 2009 Krumholz 2013) , , is the fraction of the mid-plane pressure due to self-gravity of the gas only, and not stars or dark matter , and mp is the ratio of the total to the turbulent pressure at the mid-plane. Following equation 30, we can derive Σ ★0 and ★ ( ) as,

Next, we consider the cosmic accretion of gas onto the disc. The functional form of ★ ( ) is not provided in the model. Within the framework of inside-out galaxy formation, Σ cos decreases with radius, as has been noted in several works (Chiappini et al. 1997 (Chiappini et al. , 2001 Fu et al. 2009; Courty et al. 2010; Forbes et al. 2014b; Pezzulli & Fraternali 2016; Mollá et al. 2016) . In particular, we find from Colavitti et al. (2008, see their Figure 2 ) that ★ ≈ 1/ 2 is necessary to reproduce the present day total surface mass density along the disc in the Milky Way. Additionally, a 1/ 2 accretion profile is also identical to ★ , implying a direct correlation between star formation and accretion, as has been noticed in simulations (Davé et al. 2011) . Such a profile also means that more accretion is expected in more massive parts of the disc due to higher gravitational potential (Prantzos & Boissier 2000) . Keeping these results in mind, we set ★ ( ) = 1/ 2 . However, we show in Appendix A that changing the functional form of ★ ( ) has only modest effects on the qualitative results. Following Forbes et al. (2014a) , we define

where B ≈ 0.17 is the universal baryonic fraction (White & Fabian 1995; Burkert et al. 2010; Planck Collaboration et al. 2016) , and in is the baryonic accretion efficiency given by Forbes et al. (2014a, equation 22) , which is based on cosmological simulations performed by Faucher-Giguère et al. (2011) . h is the dark matter accretion rate (Neistein & Dekel 2008; Bouché et al. 2010; Dekel et al. 2013) given by Krumholz et al. (2018, equation 

where the halo mass, h , can be written in terms of by assuming a Navarro et al. (1997) density profile for the halo as , equations 69 to 71),

where is the halo concentration parameter (Mo et al. 2010 , section 7.5). It is now known that scales inversely with halo mass (Macciò et al. 2007; Zhao et al. 2009; Dutton & Macciò 2014) . For the purposes of this work, we simply adopt = 10, 15 and 13 for local spirals, local dwarfs and high-galaxies, respectively, rather than adopting more complex empirical relations (e.g., Forbes et al. 2019) . Finally, note that the numerator in equation 33 is simply the baryonic accretion rate, ext . The inflow rate required to maintain a steady state is given by the balance between radial transport, turbulent dissipation and star formation feedback (Krumholz et al. 2018, equation 49) = 4(1 + )

Here sf is the gas velocity dispersion that can be maintained by star formation feedback alone, is the scaling factor for the rate of turbulent dissipation (Krumholz & Burkert 2010) , and nt is the fraction of gas velocity dispersion that is turbulent as opposed to thermal. While a cosmological equilibrium dictates that ext (and also ★ ext , with the former being the star formation rate), it is unclear if these conditions in fact hold for observed galaxies at high redshift. We discuss this in detail in Appendix B, showing that these uncertainties do not affect our qualitative results on metallicity gradients.

Finally, we revisit the yield reduction factor that we introduced in equation 13. Both the mass loading factor and the direct metal ejection fraction that are incorporated into are largely unknown (Creasey et al. 2013 (Creasey et al. , 2015 Christensen et al. 2018) . A number of authors have proposed models for (e.g., Creasey et al. 2013; Forbes et al. 2014b; Torrey et al. 2019; Tacchella et al. 2020) , and it is believed to scale inversely with halo mass. However, there are no robust observational constraints, with current estimates ranging from 0 to 30 (Bouché et al. 2012; Newman et al. 2012; Kacprzak et al. 2014; Schroetter et al. 2015 Schroetter et al. , 2019 Chisholm et al. 2017; Davies et al. 2019; Förster Schreiber et al. 2019; McQuinn et al. 2019) . is even less constrained by observations and theory, although observations and simulations suggest non-zero values in dwarf galaxies (e.g., Chisholm et al. 2018; Emerick et al. 2018 Emerick et al. , 2019 . For this reason we leave as a free parameter in the model and present solutions for metallicity evolution for a range of values. As we show in a companion paper (Sharda et al. 2021a) , galaxies tend to prefer a particular value of based on their stellar mass, ★ .

We list fiducial values of all the parameters used in the model in Table 1 and Table 2 . Plugging in these parameters in equations 15−18, we get,

where we have explicitly retained the dependence of the radial profile of cosmic accretion rate surface density in A. Note that none of these ratios depend on 0 . Some of these parameters are dependent on other parameters: e.g., h can be expressed as a function of as is clear from equation 34 and equation 35.

Now, we can combine the metallicity evolution model from Section 2.1 and the galactic disc model from Section 2.2 to obtain an analytic solution to equation 14 in steady-state ( Z/ = 0). The solution is

where 1 is a constant of integration and Z 0 ≡ Z( = 0 ). We remind the reader that Z = / and = / 0 as we define in Section 2.1. In writing the above analytic solution, we have assumed that the metallicity at the inner edge of the disc (to which we shall hereafter refer as the central metallicity), Z 0 , is known. We show below (Section 3) that this approach is reasonable, because the solutions naturally tend toward a particular value of Z 0 . Thus, in practice, 1 is the only unknown parameter in the solution. We also show later in Section 3 that 1 can be expressed as a function of the metallicity gradient at 0 .

We now turn to constraining 1 . Firstly, note that Z > 0 for all . In practice, we ask that Z > Z min for some fiducial Z min ≈ 0.01. For 1, this gives

where max is the outer radius of the disc at which we apply this condition 3 . Secondly, the total metal flux into the disc across the outer boundary cannot exceed that supplied by advection of gas with metallicity Z CGM into the disc, since otherwise this would imply the presence of a metal reservoir external to the disc that is supplying metals to it, which is only true in special circumstances, e.g., during or after a merger Hani et al. 2018) , or due to long term wind recycling through strong galactic fountains 3 The inequality is such that applying this condition at max ensures that it is also satisfied everywhere else in the disc. (Grand et al. 2019) . Mathematically, this condition can be written as

For 1, this translates to,

Thus, we find that 1 is bounded within a range dictated by the two conditions above. Given a value of 1 , we can also calculate the Σ -weighted and Σ ★ -weighted mean metallicity in the model,

Finding Z is helpful because we can use it to produce a massmetallicity relation (MZR) that can serve as a sanity check for the model. We show in a companion paper that our model can indeed reproduce the MZR (Sharda et al. 2021a ).

We apply our model to four different classes of galaxies: local spirals, local ultra-luminous infrared galaxies (ULIRGs), local dwarfs, and high-galaxies. The fiducial dimensional parameters we adopt for each of these galaxy types are listed in Table 1 and Table 2 . We remind the reader that the metallicity evolution model can only be applied to those galaxies where the metallicity gradient can reach equilibrium. This condition is approximately satisfied if eqbm < H(z) , where H(z) is the Hubble time at redshift . We also compare eqbm with the molecular gas depletion timescale dep,H 2 , since we expect that dep,H 2 controls the metal production timescale (hence, S) and can potentially impact metallicity gradients. Thus, the metallicity gradients may also not be in equilibrium if eqbm dep,H 2 . An exception to this is for local ULIRGs, where we compare eqbm with merge , the merger timescale. This is because the dynamics of the galaxy (as dictated by its rotation curve and orbital time) are dictated by merge for local ULIRGs.

Before checking whether equilibrium is satisfied for each individual galaxy class, it is helpful to put our work in context. Considering galaxies' total metallicity (as opposed to metallicity gradient), Forbes et al. (2014b, see their Figure 15 ) predict that galaxies with halo masses h ≥ 10 10.5 M (corresponding to ★ ≥ 10 9 M - Moster et al. 2010, their Figure 4 ) reach equilibrium by ≈ 2.5. Feldmann (2015) use a linear stability analysis to show that the metal equilibration time is at most of order the gas depletion time dep , which is small compared to H(z) for all massive main sequence galaxies. Similar arguments have been made by Davé et al. (2011 Davé et al. ( , 2012 and Lilly et al. (2013) where the authors find that the metallicity attains equilibrium on very short timescales as compared to dep , and is thus in equilibrium both in the local and the high-Universe. In contrast, Krumholz & Ting (2018) study metallicity fluctuations, and find that these attain equilibrium on an even shorter timescale, ∼ 300 Myr. Our naive expectation is that equilibration times for metallicity gradients should be intermediate between those for total metallicity and those for local metallicity fluctuations, and thus should generally be in equilibrium. We show later in Section 5 that, while these expectations are in general satisfied, some galaxy classes, namely, local dwarfs with no radial inflow, local ULIRGs, and galaxies with inverted gradients, can be out of equilibrium. Thus, our model cannot be applied to these galaxies.

For the rest of the galaxies where the equilibrium model can be applied, we use the fiducial parameters that we list in Table 2 , and solve the resulting differential equation to obtain Z( ), for different yield reduction factors. We list the resulting values of T , P, S and A for different galaxies in Table 3 . To mimic the process followed in observations and simulations (e.g., Carton et al. 2018; Collacchioni et al. 2020 ) as well as existing models (e.g., Fu et al. 2009 ), we linearly fit the resulting metallicity profiles using least squares with equal weighting in logarithmic space log 10 Z ( ) = log 10 Z 0 + ∇ log 10 Z ( ) ,

between = 1 and max , thereby excluding the innermost galactic disc where the rotation curve index is not constant, and where factors such as stellar bars can affect the central metallicity (Florido et al. 2012; Zurita et al. 2021 ). While it is clear from equation 41 that the functional form of Z is such that log 10 Z may not be a linear function of in certain cases, we will continue to use the linear fit as above in order to compare with observations. We show in Appendix C how the gradients change if we vary min or max . For each class of galaxy that we discuss in the subsections below, we plot a range of gradients that results from the constraints on the constant of integration 1 (see Section 2.3), as well as the weighted mean metallicities, Z Σ and Z Σ ★ .

For local spirals, we select the outer boundary of the star-forming disc to be 15 kpc, thus max = 15, reminding the reader that = / 0 where 0 = 1 kpc. We first study the metallicity equilibration time ( eqbm ) to see if metallicity gradients in these galaxies have attained equilibrium, so that the model can be applied to them. Figure 1 shows the value of eqbm we find from equation 19 as a function of for local spirals for different values of the yield reduction factor, ; the bands shown correspond to solutions covering all allowed values of the integration constant 1 . It is clear from Figure 1 that eqbm < H(0) for all possible and 1 , so we conclude that the gradients in local spirals are in equilibrium. Additionally, eqbm ∼ dep,H 2 for local spirals (1−3 Gyr, e.g., Wong & Blitz 2002; Bigiel et al. 2008; Saintonge et al. 2012; Leroy et al. 2013; Huang & Kauffmann 2014) , implying that the metallicity distribution reaches equilibrium on timescales comparable to the molecular gas depletion timescale. The model also predicts that central regions of local spirals should achieve equilibrium earlier than the outskirts, however, this is somewhat sensitive to the choice of 1 and as we can see from Figure Our equilibrium timescales are also consistent with our naive expectation as stated above: long compared to the timescale for local fluctuations to damp, but shorter than the time required for the total metallicity to reach equilibrium. Figure 2 presents the family of radial metallicity distributions we obtain from the model for local spirals; the different lines correspond to varying choices of the outer boundary condition 1 , from the minimum to the maximum allowed. We report in the text annotations that accompany these curves the range of gas-and SFR-weighted mean metallicities Z Σ and Z Σ ★ , and metallicity gradients ∇(log 10 Z), spanned by the models shown. To aid in the interpretation of these results, in Figure 3 we also show the magnitudes of the various terms in the numerator on the right hand side of equation 19, which represent, respectively, the relative importance of advection, diffusion, metal production (reduced by metal ejection in outflows), and cosmological accretion in determining the metallicity gradient. We use this figure to read off which processes are dominant in different parts of the disc. While the source and the accretion terms fall off in the outermost regions due to the 1/ 2 dependence, the advection and diffusion terms slightly increase with , thereby resulting in a shorter metal equilibration time in the outermost regions as compared to intermediate regions, as we see in Figure 1 . Thus, transport processes in the outer regions play an important role in establishing metal equilibrium in local spirals.

There are several noteworthy features in these plots. First, note how the solution asymptotically reaches a particular value of the central metallicity. We choose to set Z 0 to this value, but we emphasise that the behaviour of the solution does not depend on this choice except very close to = 1: if we choose a different value of Z 0 , the solution is (by construction) forced to this value close to = 1, but returns to the asymptotic limit for 1.1. Indeed, we shall see that this is a generic feature for all of our cases: the limiting central metallicity is set by a balance between two dominant processes, and can be deduced analytically by equating the two dominant terms in equation 41. For the case of local spirals, the two dominant terms throughout the disc are production and accretion, as we can read off from Figure 3 . The balance between these two processes gives

This matches the conclusions of Finlator & Davé (2008) regarding the total metallicity. However, we show below in Section 3.2 that this conclusion holds only for local, massive galaxies, since other processes like metal transport also play a significant role in low mass galaxies as well as at high redshift. Using the above definition of Z 0 , we can now express 1 in a more physically-meaningful way 

Thus, for local spirals, 1 essentially describes the metallicity gradient at 0 . Second, both the central metallicity Z 0 and the mean metallicity Z decrease with decreasing , as expected; we obtain mean metallicities close to Solar, as expected for massive local spirals, for fairly close to unity. Thus our models give reasonable total metallicities for local spirals if we assume that there is relatively little preferential ejection of metals, consistent with the results of recent simulations (Du et al. 2017; Tanner 2020; Taylor et al. 2020) . Note that some semi-analytic models find a high metal ejection fraction for spirals, but self-consistently following the evolution of the CGM subsequently leads to high re-accretion of the ejected metals ). In the language of our model, this essentially implies a high when averaging over the metal recycling timescale for local spirals, consistent with our expectations.

Third, and most importantly for our focus in this paper, the value of has little effect on the metallicity gradient, as is clear from the similar range of gradients produced by the model for different . Our models robustly predict a gradient ∇(log 10 Z) ≈ −0.07 to 0 dex kpc −1 , in very good agreement with the range observed in local spirals (e.g., Zaritsky et al. 1994; Sánchez et al. 2014; Ho et al. 2015; Sánchez-Menguiano et al. 2016; Belfiore et al. 2017; Erroz-Ferrer et al. 2019; Mingozzi et al. 2020) , and within the range provided by existing simpler models of metallicity gradients (Chiappini et al. 2001; Fu et al. 2009 ).

Apart from the mean gradient, we can also study the detailed shape of the metal distribution with the model. For the given input parameters as in Table 2 , the model features a nearly-flat metal distribution in the inner galaxy for all allowed values of 1 . Such flat gradients in the inner regions are commonly observed in local spirals (Moran et al. 2012; Belfiore et al. 2017; Mingozzi et al. 2020) , and have been attributed to metallicity reaching saturation in these regions (Zinchenko et al. 2016; Maiolino & Mannucci 2019) , although the flatness depends on the metallicity calibration used (Yates et al. 2020, Figure 4 ). This is also the case for our models of spirals, since the flat region corresponds to the part of the disc where the metallicity is set by the balance between metal injection and dilution by metal-poor infall (c.f. Figure 3 ). For comparison, we also show in Figure 2 the measured average metallicity profiles in local spirals observed in the MaNGA survey (Belfiore et al. 2017) using two different metallicity calibrations (Pettini & Pagel 2004; Maiolino et al. 2008) , where we have adjusted the overall metallicity normalisation by 0.02 dex so that the model profiles overlap with the data. We see that the profiles produced by the model are in reasonable agreement with that seen in the observations (see also, Sánchez-Menguiano et al. 2018) .

Several works have also noted that local spirals with higher gas fractions (at fixed mass) show steeper metallicity gradients (Carton et al. 2015; De Vis et al. 2019; Pace et al. 2020 ). In the language of the model, a higher gas fraction implies a higher value of , and , . Increasing these parameters leads to an increase in the source term S, which gives rise to steeper metallicity gradients in the model, consistent with the above observations. Moreover, a higher gas fraction (i.e., higher , and , ) also results in a rather steep metallicity profile in the inner disc, thus giving slightly lower metallicities in the inner disc as compared to the fiducial case above, consistent with the standard picture of galaxy chemical evolution (Tinsley 1972 (Tinsley , 1973 ; see also, Pace et al. 2020) .

It is difficult to provide robust predictions for the metal distribution in the outer parts of the galaxy without further constraining 1 . The outer-galaxy metal distribution in the model is also sensitive to parameters like the galaxy size and the CGM metallicity. The result of these uncertainties is that depending on the choice of 1 , the model can produce both nearly-flat and quite steep metal distributions in the outer parts of the galaxy. A steep drop in the metallicity in the outer disc has been observed in several local spirals (Moran et al. 2012 ), but is dependent on the metallicity calibration used (Carton et al. 2015) . In our models, this region corresponds to where cosmological accretion of metal-poor gas onto the disc becomes less important than inward advection of metal-poor gas through the disc -a process whose rate we would expect to be correlated with the available mass supply in the far outer disc, as measured by H . Note that the gradient can also flatten again in the outermost regions in the disc (Werk et al. 2011; Sánchez et al. 2014; Sánchez-Menguiano et al. 2016; Bresolin 2019) ; however, these regions typically have insufficient spatial resolution (Acharyya et al. 2020a ) as well as significant diffused ionised gas emission, both of which can cause the gradients to appear flatter than their true values (Kewley et al. 2019a, Section 6) . Given the uncertainties in the model as well as observations of metallicities in the outer discs in spirals, it is not yet obvious if the metal distribution in the outer disc in the model can be validated against the available observations. Thus, we do not study these regions with our model. This analysis also shows that linear fits to the metallicity profiles is a crude approximation to the true underlying metallicity distribution in local spiral galaxies.

Our model can also be applied to local dwarf galaxies that can be classified as rotation-dominated, e.g., the Large Magellanic Cloud (LMC), for which ∼ 60 km s −1 and ∼ 7 km s −1 (Alves & Nelson 2000) . Such galaxies typically lie at the massive end of dwarfs ( ★,LMC = 2 × 10 9 M , as reported in van der Marel 2006; Skibba et al. 2012) , and possess an equilibrium gas disc to which the unified galaxy evolution model of can be applied. We set the outer disc radius to 6 kpc to find the gradient in Table 1 and Table 2 , for different values of the yield reduction factor, . The analytic solution to the metallicity evolution equation is given by equation 41. The slope of the linear fit to the model gradients between = 1 − 15 (black, dashed lines) gives the metallicity gradient that can be compared against simulations and observations. The blue coloured curves show the acceptable parameter space of the gradients based on the constraints on the constant of integration, 1 , using the boundary conditions criteria described in Section 2.3. The metallicity at the inner edge of the disc (referred to as the central metallicity in the text), Z 0 , is set by the balance between source and accretion for local spirals (see equation 48). Z Σ and Z Σ★ represent the range of massweighted and SFR-weighted mean equilibrium metallicities produced by the solution, respectively (see equation 46). We expect closer to unity for local spirals, implying that metals in these galaxies are well-mixed with the ISM before they are ejected. Finally, in the top panel we overplot the average metallicity profiles observed in local spirals in the MaNGA survey by Belfiore et al. (2017) using the PP04 (Pettini & Pagel 2004 ) and M08 (Maiolino et al. 2008 ) calibrations, adjusting the normalisation to overlap with the model profiles. that collectively build the metallicity gradient in local spirals, for a fixed Z CGM = 0.1 and fixed 1 for different yield reduction factors, . These terms are defined in equation 14. The leading terms that set the gradients in local spirals are metal production and accretion of gas onto the galaxy, whereas advection and diffusion play a subdominant role in local spirals, due to the small velocity dispersion, . Note that the sharp feature in the diffusion term near = 1.3 corresponds to the location where this term passes through zero as it changes sign; the term in fact behaves smoothly everywhere, but this behaviour appears as a sharp feature when plotted on a logarithmic axis. the fiducial model, in line with the estimated gas disc size of local dwarfs ( LMC ∼ 4.3 kpc, Westerlund 1990) . Figure 4 shows the metal equilibration time, eqbm , for local dwarfs based on the parameters we list in Table 1 and Table 2 . It is clear that metallicity gradients are in equilibrium in dwarfs, since eqbm < H(0) as in the case of local spirals (see, however, Section 5.2.1 where we show that this may not be the case under certain circumstances). Contrary to local spirals, local dwarfs show a wide range of dep,H 2 , from a few hundred Myr to several Gyr (e.g., Bolatto et al. 2011; Bothwell et al. 2014; Hunt et al. 2015; Jameson et al. 2016; Schruba et al. 2017) , similar to the scatter we find in eqbm (see also, Section 5.2.1) 4 .

Having established metal equilibrium in local dwarfs, we can now study the gradients produced by the model. Figure 5 shows the resulting metallicity versus radius for different (analogous to Figure 2 ), and Figure 6 shows the relative importance of the various processes (analogous to Figure 3 ).

In the case of local dwarfs, we see that Z 0 is set by the balance between advection and diffusion, giving

4 While it is often quoted that dep,H 2 is smaller by a factor of 2 − 5 in local dwarfs as compared to local spirals, Schruba et al. (2017) point out that this may not necessarily be true. This is because it is difficult to trace the entire molecular gas content in dwarfs, and a significant fraction of the molecular gas can be 'CO-faint' or 'CO-dark' (Bolatto et al. 2011; Jameson et al. 2018 Using the above definition of Z 0 , we can express 1 as

Central metallicities are in the range Z 0 ≈ 0.2 − 0.6 depending on the choice of , in good agreement with that observed in local dwarfs, e.g., in the SMC and the LMC (Russell & Dopita 1992; Westerlund 1997) , and M82 (Origlia et al. 2004) . While Z 0 depends only on S/A in local spirals, it also depends on the choice of 1 for local dwarfs, implying that it is independent of the disc properties in the former case but not in the latter. 5 Similarly, mean metallicities range from Z ∼ 0.1 − 0.5 as varies from ≈ 0.1 − 1; both observations (Martin et al. 2002; Strickland & Heckman 2009; Chisholm et al. 2018 ) and numerical simulations (Emerick et al. 2018 (Emerick et al. , 2019 suggest that dwarfs suffer considerable direct metal loss, so considerably smaller than unity seems likely. As opposed to spirals, our models predict that gradients are not necessarily flat in the inner regions of dwarfs, which is also consistent with observations (Belfiore et al. 2017; Mingozzi et al. 2020) . The reason for this difference is due to different physical processes dominating in the two types of galaxies: accretion versus metal production in spirals, and advection versus production in dwarfs. Consequently, we predict linear gradients for local dwarfs that are steeper than the ones for local spirals at fixed and 1 . For the smaller values of expected in local dwarfs, we expect gradients in the range ∼ −0.01 to −0.15 dex kpc −1 , implying a larger scatter in the gradients measured in local dwarfs as compared to that in local spirals, consistent with observations (Mingozzi et al. 2020 , Figure 12 ). The metallicity profiles produced by the model for smaller values of are also in agreement with that observed in the MaNGA survey (Belfiore et al. 2017 ), as we show in Figure 5 , where we have adjusted the overall metallicity normalization by 0.15 dex to facilitate a comparison of the data and the model profiles. Further, the larger range of gradients in low mass local galaxies as compared to massive galaxies allowed within the framework of our model is also relevant and necessary for reproducing the observed steepening of gradients with decreasing galaxy mass (Bresolin 2019, Figure 10 ).

Although this is not illustrated in Figure 5 , we also find that the magnitude of the gradient is quite sensitive to both the "floor" velocity dispersion supplied by star formation, sf , and the Toomre parameter, since these two jointly set the strength of advection and in this case, sf ∼ . Thus, we expect that gradients for local dwarfs will show more scatter than those for local spirals. It is interesting to note that there is a similarly large scatter in simulations of dwarf galaxies, with some groups (e.g., Tissera et al. 2016) finding steeper gradients for dwarfs as compared to spirals whereas others (e.g. Ma et al. 2017 ) finding the opposite. This difference between the simulations has been attributed to the strength of feedback, which, in the language of our model, corresponds to variations in sf and ; thus the sensitivity of our model is at least qualitatively consistent with the strong dependence of feedback strength observed in simulations.

in both the diffusion term and the metallicity profile. For the purposes of plotting, we have chosen a single value of 1 , which in turn forces all models to converge to a single Z 0 . While we could correct this by choosing different values of 1 for different models so that they remain smooth, since the sharp feature does not affect the metallicity gradient that is our main focus in this paper, we choose for reasons of simplicity to retain the fixed 1 . Figure 1 , but for local dwarfs. Here, eqbm < H(0) , implying that the metallicity gradients in local dwarfs are also in equilibrium, even in the case of low (see the text for a discussion on dep,H 2 for local dwarfs). The corresponding metallicity gradients are plotted in Figure 5 .

Massive galaxies at high-are primarily rotation-dominated with underlying disc-like structures (Weiner et al. 2006; Förster Schreiber et al. 2009 Wisnioski et al. 2011 Wisnioski et al. , 2015 Wisnioski et al. , 2019 Wuyts et al. 2011; Di Teodoro et al. 2016; Simons et al. 2017; Übler et al. 2019 ). Thus, we can apply the model to these galaxies. For high-galaxies, we set the outer disc radius to 10 kpc to find the gradient in the fiducial model, acknowledging that galaxies at higher redshifts are smaller than that in the local Universe (e.g., Queyrel et al. 2012; van der Wel et al. 2014 ). Hereafter, we work with = 2 as a fiducial redshift. Figure 7 shows the metal equilibration time for highgalaxies. It is clear that eqbm < H(z) , so that the equilibrium solution can be applied to these galaxies. Following Tacconi et al. (2018 Tacconi et al. ( , 2020 , if we assume that a main sequence high-galaxy follows dep,H 2 ∝ (1 + ) −0.6 , it implies that dep,H 2 ∼ 0.5 − 1.5 Gyr for high-galaxies, which is comparable with eqbm as above. Figure 8 shows the equilibrium metallicity distributions we obtain for a fiducial high-galaxy with parameters listed in Table 1 and Table 2 , and Figure 9 shows our usual diagnostic diagram comparing the importance of different processes. Examining this diagram near = 1, it is clear that, as is the case for local dwarfs, the central metallicity Z 0 is set by the balance between advection and diffusion, which gives

It varies between Z 0 = 0.3-0.7 depending on the value of , in good agreement with observed metallicities in high-galaxies in the mass range we consider (Erb et al. 2006; Yabe et al. 2012) , with 1 same as that in equation 51. While the absolute metallicity depends on , the metallicity gradients for the most part do not -we find ∇(log 10 Z) ≈ −0.15 to −0.05 dex kpc −1 , with order-of-magnitude variations in only altering these values by a few hundredths. The gradients we find for high-galaxies are steeper than for local spirals, and the distributions are steeper at small radii than at larger radii, the opposite of our finding for local spirals. Figure 9 shows why this is the case: gradients over most of the radial extent of high-galaxies are set by the balance between source and advection, whereas accretion, which dilutes the gradients in local spirals, is sub-dominant. The fundamental reason for this change is due to Figure 2 , but for local dwarfs. Here, Z 0 is set by the balance between advection and diffusion, whereas metallicities in the disc are set by the balance between advection and source. The sharp rise and fall in the profile at = 1 is an artefact of the choice of the constant of integration 1 used to calculate Z 0 (see equation 51). The gradients are particularly sensitive to the strength of advection for local dwarfs since turbulence due to star formation feedback is comparable to that due to gravity, sf ∼ . When they are exactly equal, advection vanishes, and the gradients may not be in equilibrium (see Section 5.2.1). In the last panel we also plot (purple lines) the average metallicity profiles observed in local dwarfs in the MaNGA survey; see Figure 2 .

the vastly higher velocity dispersions of high-galaxies, which increase the importance of the advection term (P ∝ (1 − sf / )) while suppressing the accretion term (A ∝ −3 ); this effect is partly diluted by the higher accretion rates found at high-(equation 34), but the net change at high redshift is nonetheless toward a smaller role for accretion onto discs and a larger role for transport through them. We discuss this further in detail in Section 4. 

A significant advantage of our model compared to previous analytic efforts, is that it makes meaningful predictions for how galaxy metallicity gradients have evolved across cosmic time. This is the case because we do not have the freedom to adjust parameters such as radial inflow rates and profiles of star formation rate to match any given observed galaxy. Instead, these parameters are either prescribed directly from our galaxy evolution model or depend on parameters that are directly observable (e.g., galaxy velocity dispersions). The basic inputs to our model are the halo mass h and the gas velocity dispersion as a function of . We consider three different ways of selecting galaxies that yield different tracks of ℎ ( ) (see below for details), while we take the evolution of ( ) from the observed correlation obtained by Wisnioski 

6 As opposed to Wisnioski et al. (2015) , we have explicitly retained the dependence of ( ) on . where gas is the molecular gas fraction of the galaxy (Genzel et al. 2011; Genel et al. 2012; Tacconi et al. 2013; Genzel et al. 2015; Faucher-Giguère 2018) . This scaling is subject to considerable observational uncertainty, the implications of which we explore in Appendix B. We follow Wisnioski et al. (2015) to find gas as a function of ★ and from Tacconi et al. (2013) and Whitaker et al. (2014) , as it is now known that gas decreases with cosmic time and stellar mass Saintonge et al. 2011; Geach et al. 2011; Davé et al. 2012; Tacconi et al. 2013 . Same as Figure 3 , but for high-galaxies. Here, the metallicities in the disc are set by the balance between source and advection, due to efficient radial transport of the gas. produce from the model in this section are in equilibrium across the redshifts we use, since eqbm < H(z) .

We first study how the gradient in a Milky Way-like galaxy has evolved over time using our model. We only need one parameter to begin with: at = 0. We set this to 220 km s −1 (Bland-Hawthorn & Gerhard 2016). Then, we use equation 35 to calculate h ( = 0) for a fixed = 15. Using h ( = 0) as boundary condition, we integrate equation 34 to find h ( ), keeping in mind that this equation represents an average evolution of h ( ) that may not necessarily apply to the Milky Way. Then, we utilize h ( ) to find ( ), by changing the concentration parameter ( ) as an empirical third-order polynomial fit, following Zhao et al. (2009) . This ensures that as we change , we self-consistently find h and . We adopt a simple linear variation for the outer edge of the star-forming disc, max , as a function of such that it is 15 at = 0 and 10 at = 2. Similarly, we vary sf between 0.5 and 1 across redshift, keeping in mind that sf cannot be more than 1 at any redshift. For simplicity, we fix the other parameters as follows: = 0, , = , = 0.5, sf = 7 km s −1 , min = 1.5 and Z CGM = 0.1.

We show the resulting evolution of the gradient in Figure 10 . The model predicts a steepening of the gradient in Milky Way-like galaxies over time, with the exception of a very recent flattening, between ≈ 0.15 and 0. We can understand these trends in terms of the dimensionless parameters S, P, and A that describe the relative importance of in situ metal production, radial advection, and cosmological accretion with diffusion, respectively. The source term S will always make the gradients steeper because of the steep radial profile of Σ ★ , and it is either P or A that balances S to give rise to flatter gradients. The steepest gradients at ≈ 0.15 correspond to when both P and A are at their weakest compared to S. We can understand the trends on either side of this maximum in turn.

First, let us focus on the recent epoch, 0.15. During this period, cosmological accretion (A) is more important than radial transport (P), and accretion and metal production depend on the galaxy rotational velocity as A ∝ 3.3 and S ∝ 2 , respectively. Thus, as the galaxy grows in mass, dilution by accretion gets stronger compared to metal production, leading to the recent flattening in the model. However, this can change if the metal production is underestimated, e.g., due to ignoring the contribution from longterm wind recycling (Leitner & Kravtsov 2011) .

During this epoch advection is more important than accretion, P > A. The ratio of the two effects, P/A, is large at high redshift, and decreases systematically towards the present day, ultimately reaching P/A ≈ 1 at ≈ 0.15. This transition is ultimately driven by the systematic decrease in galaxy velocity dispersions with redshift, as already discussed in the context of our high-galaxy models (Section 3.3): higher velocity dispersions are strongly correlated with higher rates of radial inflow through a galaxy, so that for a Milky Way progenitor at 1, radial inflow transports metal-poor gas into galaxy centres ∼ 10× faster than cosmological accretion (P/A ≈ 10) -despite the fact that the absolute accretion rate is higher at 1 than it is today. Similarly, the ratio of radial inflow to metal production, P/S, scales with velocity dispersion as 2 (for sf ), so radial inflow also becomes more important relative to metal production as we go to higher redshift and higher velocity dispersion. This explains the flatness of gradients at high redshift 7 . This transition from radial advection being dominant to being unimportant is mirrored in the transition from gravity-driven to star formation feedback-driven turbulence from high-to low- , as we noted earlier in Section 3.3.

Lastly, we find that diffusion is sub-dominant compared to both advection and accretion at all cosmological epochs, because P and A are never both less than unity at the same time. Thus, while diffusion can have some effects on the metallicity distributions, particularly towards galaxy centres (cf. Figure 6) , as well as on metal equilibrium timescales (cf. Figure 1) , it is generally unimportant for setting galaxy metallicity gradients.

There is extensive data on the history of the Galaxy's metallicity gradient, as summarised by Mollá et al. (2019, see their Table 1) , and on the history of the gradients in a number of other nearby galaxies. The general outcome of these studies is that gradients measured in H regions (which trace the current-day metal distribution) are steeper than those measured in planetary nebulae or open clusters (which trace older populations) (Stanghellini & Haywood 2010; Stanghellini et al. , 2014 Sanders et al. 2012; Stasińska et al. 2013; Magrini et al. 2016 ). This implies a steepening of the gradient with time in Milky Way-like galaxies, however, this should be treated with caution because measured metallicity gradients in the Galaxy are subject to large errors arising from uncertainties in estimating the ages of the planetary nebulae (Maciel et al. 2010; Cavichia et al. 2011) , and due to radial migration that could result in a movement of the planetary nebulae away from their origin (Minchev et al. 2013) 

To allow a quantitative comparison of these observations with 7 Note that this is a qualitatively different outcome than our comparison of local spirals and high-galaxies in Section 3.1 and Section 3.3, where high-galaxies were found to have steeper gradients. The difference can be understood by recalling that in Section 3.1 and Section 3.3 we were comparing galaxies with comparable rotation curve speeds , whereas here we are following a single growing galaxy, so is much smaller at high-than at = 0. This reduces S at high-. 8 Some earlier work reported the opposite trend, whereby the metallicity gradient in the Galaxy was initially steep and has flattened over time (Maciel et al. 2003; Mollá & Díaz 2005) , while other work found little or no evolution in the gradient over time (Maciel & Costa 2013) . This is a difficult our model, we show measurements of the metallicity gradient for the Milky Way as a function of lookback time from as yellow circles in Figure 10 . The data for the Milky Way (as well as other local spirals, see Stanghellini et al. 2014) are in qualitative agreement with the predictions from our model. However, we also note that for our model to agree quantitatively with the measurements, we would need to be lower at high redshift, and increase towards unity today. Such a change in is plausible and is consistent with our expectation that should be close to unity in more massive galaxies like the present-day Milky Way, and smaller than unity in less massive galaxies with shallower potential wells, such as the Milky Way's high-progenitors. However, the exact form of this evolution is not independently predicted by our model.

On Figure 10 , we also overplot results from Feedback In Realistic Environments (FIRE) simulations (Hopkins et al. 2014 (Hopkins et al. , 2018 ) of a Milky Way-like galaxy (m12i) discussed in Ma et al. (2017) . This simulation finds that metallicity gradients are unstable until ∼ 1, after which they steepen and stabilise to an equilibrium value. This transition is primarily due to the formation of a robust galactic disc that cannot be disrupted again due to internal or external feedback. While the quantitative trends slightly differ at some redshifts between our model and the simulation, which is not unexpected given that the exact implementation of the feedback and measurements of the gradients are different, there is a very good qualitative match. This match also implies that Milky Way-like galaxies would have had lower in the past as compared to the present day, as outflows were more common and stronger in the past due to higher SFR and could have ejected a larger fraction of metals not mixed with the ISM (Muratov et al. 2015; Ma et al. 2017) ; such a scenario has received support from recent high-resolution simulations that spatially resolve multi-phase galactic outflows, and find that the metal enrichment factor in both the cold (< 2 × 10 4 K) and hot (> 5 × 10 5 K) outflows increases with the SFR surface density (Kim et al. 2020) . We can also compare our results with those of Gibson et al. (2013) , where the authors study two identical simulation suites with either weak or enhanced stellar feedback, called MUGS and MaGICC, respectively (Stinson et al. 2010) . The authors find that gas phase metallicity gradients are steep at high redshift in MUGS, whereas they are flat in MaGICC, clearly revealing the close correlation between feedback and metallicity gradients in galaxies. One of their simulated galaxies, MaGICC g1536, resembles the Milky Way in terms of its stellar mass, so we also compare our model results to that simulation in Figure 10 . Again, we find qualitative similarities between the simulations and the model.

In this section, we study the mass-averaged trends of metallicity gradients across cosmic time. For this purpose, we use a compilation of observations of metallicity gradients in (lensed and un-lensed) galaxies spanning 0 ≤ ≤ 2.5 (Queyrel et al. 2012; Swinbank et al. 2012; Stott et al. 2014; Leethochawalit et al. 2016; Wuyts et al. 2016; Molina et al. 2017; Carton et al. 2018; Förster Schreiber et al. 2018; Wang et al. 2020; Curti et al. 2020) , and we also include results measurement, and the error bars and uncertainties are large (Maciel et al. 2010; Cavichia et al. 2011; Minchev et al. 2013 , while symbol colour shows the ratio of the dimensionless numbers P/A that describe the relative importance of radial transport and cosmological accretion, respectively. The grey curve is taken from FIRE simulations of a Milky Way-like galaxy (Ma et al. 2017) whereas the dashed, black curve is from the MaGICC g1536 simulation by Gibson et al. (2013) . The orange points are from observations of H regions, planetary nebulae and open clusters by Stanghellini & Haywood (2010) , with horizontal errorbars representing the uncertainties in the ages of planetary nebulae and open clusters. The data, simulations and the model all qualitatively show that gradients in Milky Way-like galaxies have steepened over time, with the model predicting a mild flattening between = 0.15 and present-day. In the model, this evolution is driven by a transition from the advection-dominated regime (P/A > 1) to the accretion-dominated regime (P/A < 1) around ≈ 0.15. Such a transition in metallicity gradients is mirrored in the transition in gravity-driven turbulence at high to star formation feedback-driven turbulence at = 0 from local surveys Sánchez-Menguiano et al. 2016; Belfiore et al. 2017; Mingozzi et al. 2020; Acharyya et al. 2020b) .

Before proceeding, we warn the reader that there are many uncertainties inherent in comparing metallicity gradients across samples and across cosmic time. For example, most studies in the compiled dataset rely on strong line calibrations that use photoionisation models or electron temperature-based empirical relations to measure metallicity gradients, and the variations between different calibrations can be as high as 0.1 dex per effective half-light radius (Moustakas et al. 2010; Poetrodjojo et al. 2019; Poetrodjojo et al. 2021; Mingozzi et al. 2020 ). Further, since many high-metallicity gradient measurements rely on nitrogen whereas low-measurements use a larger set of (optical) emission lines, we also expect some systematic differences in these measurements with redshift (Carton et al. 2018; Kewley et al. 2019a) . Using nitrogen can also lead to systematically flatter gradients due to different scalings of N/O with O/H in galaxy centres and outskirts .

Lastly, it is not yet clear if strong line metallicity calibrations developed for the ISM properties of local galaxies are also applicable at high-, where ISM electron densities, ionisation parameters, N/O ratios, or other conditions may differ from those in local galaxies (e.g., Shirazi et al. 2014; Sanders et al. 2016; Onodera et al. 2016; Kashino et al. 2017; Kaasinen et al. 2017; Kewley et al. 2019b; Davies et al. 2020) . We acknowledge these biases and uncertainties in the measured sample due to different techniques and calibrations or the lack of spatial and/or spectral resolution (Yuan et al. 2013; Mast et al. 2014; Carton et al. 2017; Acharyya et al. 2020a ). We do not attempt to correct for these effects or homogenize the sample because our goal here is simply to get a qualitative interpretation of the data with the help of the model, and not to obtain precise measurements from these data. Future facilities like JWST and ELTs will provide more reliable metallicity measurements, thereby enabling a more robust comparison of the model with the data (Bunker et al. 2020 ).

We bin the data into three bins of ★ : 9 ≤ log 10 ★ /M < . Trends in metallicity gradients as a function of redshift and lookback time. Colored markers represent individual galaxies within the three ★ bins as shown in the legend, with bigger markers representing binned averages of non-positive gradients across different redshift bins, and errorbars representing the scatter in the data within each redshift bin. The averages at = 0 are taken from local surveys Sánchez-Menguiano et al. 2016; Belfiore et al. 2017; Mingozzi et al. 2020) . The high-redshift compilation data is taken from Queyrel et al. (2012) Curti et al. (2020) , and is inhomogeneous, with systematic issues within the different measurements (see Section 4.2). The colored bands represent models at three ★ values, with the spread resulting from different yield reduction factors , as marked by the arrow besides the shaded region. This spread in the model is largest for the low mass galaxies. While the general trend of mild evolution of gradients across redshift holds true, the models uncover the underlying variations due to galaxies transitioning from advection-to accretion-dominated regimes between = 2.5 and 0, as is visible in the binned data averages. Some data points lie outside the range of the plot, and we do not include those for the purposes of studying the average trends of the data with the model. 10, 10 ≤ log 10 ★ /M < 11 and log 10 ★ /M ≥ 11. Figure 11 shows the individual data as well as the binned averages of nonpositive gradients (represented by bigger markers) with errorbars representing the scatter in the data within different redshift bins. We only select galaxies that show non-positive gradients while estimating the average gradient in different mass bins because our model may not apply to galaxies with positive gradients, as we explore in Section 5.2.3. We bin the data in redshift such that we can avoid redshifts where there is no data due to atmospheric absorption; such a bin selection in redshift also ensures that the binned averages reflect the true underlying sample for which the averages are calculated. We have verified that our results are not sensitive to the choice of binning the data. For simplicity, we do not overplot measurements for individual galaxies at = 0.

For the model, we select three representative ★ values corresponding to the mean of the three stellar mass data bins as above. Specifically, we use: log 10 ★ /M = 9.6, 10.4 and 11.1 M for the model. We start the calculation by selecting rotation curve speeds ( ) corresponding to each of these ★ values based on the ★ − h relation at all (Moster et al. 2013) . Given values of ℎ ( ) and ( ) corresponding to each stellar mass ★ at each redshift , we use our model to predict the equilibrium metallicity gradient exactly as in Section 4.1.

We plot the resulting range of metallicity gradients from the model points in Figure 11 . As in other figures, the spread in the model represents different between 0 and 1 (note the arrow besides the shaded regions corresponding to the models). While there is a large scatter within the individual data points, the binned averages are in good agreement with the model. Note that almost one-third of the observed galaxies show inverted gradients, which may not be in metal equilibrium and thus may fall outside the domain of our model, as we explore in detail in Section 5.2.3. For the most massive galaxies, the model predicts a mild steepening of the gradients from = 2.5 to 1, followed by an upturn (due to the transition from advection-to accretion-dominated regime) and flattening from = 1 to 0. The available data, despite the large scatter and inhomogeneties, also seem to follow the same trend. However, the location where this upturn occurs is unknown because of the lack of data in the most massive galaxy bin around = 0.5. Upcoming large surveys like MAGPI (Foster et al. 2020 ) that will observe massive galaxies between ≈ 0.3 − 0.5 will provide crucial data that can be compared against our model in the future to establish whether this upturn is indeed real.

Additionally, we can compare our results with those from the IllustrisTNG50 simulation (Hemler et al. 2020 , Figure 6 ). While our results match theirs at low redshifts, there are certain differences at high redshifts where IlustrisTNG50 fails to reproduce the observed flattening, as already noted by the corresponding authors. We explain in a companion paper (Sharda et al. 2021b ) that this difference could primarily be due to the gas velocity dispersion ( ). At high redshift, IllustrisTNG50 systematically under-predicts galaxy velocity dispersions as compared to, for example, the EAGLE simulations (Pillepich et al. 2019 , Figure 12a) , and the empirical relation we use from Wisnioski et al. (2015) .

There is a large diversity of gradients at all redshifts (Curti et al. 2020 ), particularly at low stellar mass. This observed scatter can be explained in part due to the range of in our model. For example, we notice from Figure 11 that the scatter in the model due to for the most massive galaxies is lower at low than at high . This is consistent with the trend of larger scatter in the gradients of massive galaxies at higher redshift observed in the IllutrisTNG50 simulations (Hemler et al. 2020 , Figure 6 ). On the other hand, the scatter in the model is the largest near the upturn, where galaxies transition from advection-dominated to accretiondominated regime. Between the three models, the scatter due to is the highest for the lowest ★ , thus reflecting the diverse variety of gradients that can form in low-mass galaxies. This prediction of the model is consistent with observations that find strong evidence for increased scatter in the metallicity gradients in low mass galaxies (Carton et al. 2018; Simons et al. 2020 ).

Finally, we also study the evolution of metallicity gradients across an abundance-matched sample of dark matter haloes spanning a range in 9 . Abundance-matching is based on the premise that the number density of halo progenitors should nearly remain constant across within a comoving volume in the Universe Mo & Fukugita 1996; van Dokkum et al. 2010) . It has been used to study a range of properties in local galaxies together with their high-progenitors (e.g., Marchesini et al. 2009; Papovich et al. 2011; Trujillo-Gomez et al. 2011; Leja et al. 2013; Read & Erkal 2019) , which is not possible with other 9 Abundance in the context of Section 4.3 refers to the abundance of galaxies in a given comoving volume in the Universe, and not the metallicity. selection criteria of galaxies (e.g., selecting galaxies with identical stellar mass, as we do in Section 4.2) as such galaxies evolve in time themselves (Conroy & Wechsler 2009) .

Abundance matching involves assigning more massive galaxies to more massive haloes at every ; this means selecting galaxies at each with h ( ) that satisfy

where 0 is the target number density 10 , and ( h , ) is the number of galaxies per unit mass per unit comoving cubic Mpc given by Mo & White (2002, equation 14 ) based on the Sheth & Tormen (1999) modification of the Press & Schechter (1974) formalism for the number density of haloes across . Thus, using the functional form for , we can deduce the required h at each that would correspond to an abundance-matched sample for a given 0 . Following Marchesini et al. (2009) and Papovich et al. (2011) , we study four sets of log 10 0 /Mpc −3 = −3, −3.325 − 3.5 and −4.0, respectively. For each of these 0 , we find ( ) and ★ ( ) using h ( ) from equation 54, and ( ) from equation 53. We fix = 0 for all galaxies since our choice of 0 results in massive galaxies with ★ > 10 10 M for all 0 ≤ ≤ 2.5. For simplicity, we fix , = , = 0.5 and sf = 7 km s −1 , the same as that for local spirals. Given that sf varies between 0.5 and 1 as increases, we use a cubic interpolation to vary it between = 0 and 4. We also fix Z CGM = 0.1. Figure 12 shows the cosmic evolution of gradients for an abundance-matched sample of galaxies, each panel representing a different 0 . Similar to what we have seen in prior sections, the scatter in the model is the largest at the upturn where gradients start flattening. To the best of our knowledge, there are no existing abundance-matched samples of galaxies across redshift that also contain information on metallicity gradients. However, we can construct an abundance-matched sample from the available data. We caution that constructing an abundance-matched sample from existing observations ex post facto is not as accurate as properly constructing the sample to start with. In the absence of the latter, we use our constructed sample to compare against the model to learn about the kinds of metallicity gradients that existed in progenitors of local galaxies. For this purpose, we construct our pseudo-abundance matched sample as follows: for each target value of 0 , we first select a redshift, and use equation 54 to estimate the halo mass h corresponding to the target 0 at that redshift. We then estimate the stellar mass of that galaxy ★ using the stellar mass-halo mass relation of Moster et al. (2013) . To construct our sample set at that redshift, we then take the data collection described in Section 4.2 and select galaxies that have stellar masses within ± 0.05 dex of the ★ from above; this constitutes our pseudo-abundance matched sample for that redshift, from which we then measure the mean and dispersion of metallicity gradient at that redshift bin. We plot these values in Figure 12 , along with model predictions of the metallicity gradient, which we compute from the halo mass and redshift as in previous sections. The data we obtain in this manner have considerable scatter (shown by the errorbars), but the general trends are reasonably well reproduced by the model. However, given the uncertainties in the procedure we are forced to use to construct the observed sample, 10 This approximation of a fixed 0 breaks down if certain galaxies in the abundance-matched sample do not follow the stellar mass rank order, for example, due to an abrupt increment in stellar mass because of mergers, or abrupt decrement due to quenching (Leja et al. 2013 . Trends in metallicity gradients as a function of (and lookback time) for four different abundance-matched galaxy samples given a fixed comoving number density of galaxies, 0 , color-coded by ★ . Abundance matching leads to the selection of more massive galaxies at lower redshifts, and can be used to collectively study gradients in local spirals and their high-progenitors. The orange data points reflect mean gradients for a constructed abundance-matched sample from available observations, which are the same as that reported in Figure 11 , with errorbars representing the scatter within the data. There is considerable scatter in the data, and the sample is not entirely robust given the ex post facto construction. Nonetheless, the model matches the observations reasonably-well.

it is wiser to regard the model points in Figure 12 as a prediction for future abundance matching measurements, rather than a rigorous comparison to existing data.

In this section, we describe the limitations of the model, first focusing on physical processes that we have excluded, and then discussing galaxies to which we cannot always apply our assumption of equilibrium.

Our model omits three possibly-important physical effects: bars, galactic fountains and long term wind recycling. With regard to the first of these, there is some evidence that gas phase metallicity gradients in the presence of bars in local spirals can be systematically shallower than those non-barred galaxies (Vila-Costas & Edmunds 1992; Zaritsky et al. 1994; Zurita et al. 2021 ; see, however, Sánchez-Menguiano et al. 2016 . We have not included metal redistribution due to bar-driven flows, and for this reason we limit our study to gradients in the parts of a galaxy where the rotation curve slope ( ) is a constant, which excludes bar-dominated regions (Bland-Hawthorn & Gerhard 2016; Martinez-Medina et al. 2020) . In fact, even if we wished to include bar-driven mixing, the galaxy formation model that we use as an input in Section 2.2 is itself not applicable in regions where the bar dominates the dynamics of the galaxy, since it does not include the effects of bar-driven torques on gas and SFR surface density profiles (Sun et al. 2018 (Sun et al. , 2020 .

With regard to the second issue: we do not explicitly incorporate metal redistribution via galactic fountains (Bregman 1980) . However, the combination of an enriched Z CGM and low essentially constructs a fountain process in the model that we can exploit. Semi-analytic models where the evolution of the CGM is self-consistently followed find that the CGM plays a larger role in the evolution of galaxy metallicity as it gets enriched due to outflows . We also note that galactic fountains, owing to their short fall-back timescale (∼ 100 − 300 Myr, Anglés-Alcázar et al. 2017) and short fall-back distance from the starting point (∼ 1 kpc, Spitoni et al. 2008 ) have been shown to play an insignificant role in the metallicity evolution of the local spiral M31 ; see, however, simulations by Grand et al. 2019, where fountains are thought to transport metals to the edge of the Local dwarf ( = 0) y = 1.0 y = 0.3 y = 0.1 Figure 13 . Same as Figure 4 , but without radial inflow such that P = 0.

Here, eqbm H(0) , implying that the metallicity gradients in such cases in local dwarfs may or may not be in equilibrium. Thus, our equilibrium model does not necessarily apply.

star-forming disc). Fountains possibly have a significant effect on the far outskirts of the discs, where there are few or no local sources of metals.

There is some evidence from simulations that long term wind recycling can provide metals to the disc as it re-accretes the ejected material. These simulations also show that this recycling is independent of the halo mass (Christensen et al. 2016; Tollet et al. 2019) , and can be the dominant mode of accretion of cold gas at late times. However, this recycling occurs much farther out in the disc than that we consider in our work, thus, its basic features are captured within Z CGM in the model. Additionally, while the above simulations find long term wind recycling timescale to be of the order of a Gyr, results from the EAGLE simulations find it to be comparable to H(z) . Thus, given these findings from simulations and the lack of direct observations, it is currently difficult to determine the importance of wind recycling for metallicity gradients.

Finally, we caution that our model is intended to apply mainly to metals whose production is dominated by type II SNe, and thus where the injection rate closely follows the star formation rate. We have not attempted to model elements produced by type Ia SNe or AGB stars. This is not a substantial problem for our intended application, however, since type II SNe do dominate production of the elements that are most easily observable in the gas phase (Nomoto & Leung 2018) . The one exception to this statement, where some caution is warranted, is nitrogen, to which AGB stars make a substantial contribution (Meynet & Maeder 2002; Herwig 2005) . This matters because many of the strong-line diagnostics used at high redshift rely on the N 6584 line. While observers who rely on these diagnostics usually attempt to derive the O abundance by calibrating out variations in the N/O ratio (Pettini & Pagel 2004; Masters et al. 2016; Belfiore et al. 2017; Schaefer et al. 2020; see also, Vincenzo et al. 2016) , it is nevertheless the case that variations in N abundance may influence the metallicities derived at high-, and that our model does not capture this effect.

There are certain classes of galaxies where we find that the metallicity distribution can be out of equilibrium, i.e., eqbm H(z) Local ULIRG y = 1.0 y = 0.3 y = 0.1 Figure 14 . Same as Figure 1 , but for local ultraluminous infrared galaxies (ULIRGs). Here, eqbm ∼ merge , where the latter is the merger timescale of the order of ∼ 0.3 − 1 Gyr as seen in models (Jiang et al. 2008; Torrey et al. 2012) . Thus, the metallicity gradients may not be in equilibrium throughout the merger process. In such a case, our equilibrium model for metallicity gradients cannot be applied to local ULIRGs, and the observed gradients, if any, are transient and subject to change as the merger progresses, in line with observations (Rupke et al. 2010b; Rich et al. 2012 ).

1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 Figure 3 , but for ULIRGs, which are known to be major mergers. The non-equilibrium metallicity distribution is set by advection of gas due to tidal inflows during a merger. metallicity gradients in such galaxies. Nonetheless, the limitation of the equilibrium model provides interesting constraints on the evolution of such galaxies. We discuss three such cases below.

The balance between metal production (source) and radial transport of metals through the disc (advection, diffusion) sets the metallicity gradients in local dwarfs (cf. Section 3.2). It has also been suggested that turbulence in these galaxies is mainly driven by star formation feedback and not gravity (Moiseev et al. 2015; , which gives rise to sf ∼ , and the low gas velocity dispersions observed in dwarfs (Yu et al. 2019; Varidel et al. 2020 High-z dwarf [inverted gradients] y = 0.05 CGM = 0.5 3 Figure 16 . Metallicity equilibration timescale eqbm as a function of in galaxies with inverted gradients. The first panel represents eqbm in local dwarfs. The second panel on high-discs is identical to the class of highgalaxies we discuss in Section 3.3. The third panel plots eqbm in the case of high-dwarfs that we create by combining the fiducial parameters for local dwarfs and high-galaxies (see Section 5.2.3 for details). The colors correspond to the different ways that can give rise to an inverted gradient in a galaxy: reduction in metal yield due to high preferential metal ejection ( = 0.05), enrichment of the CGM due to fountains or metal-rich flows (Z CGM = 0.5), and excessive cosmic accretion (A → 3A). The scatter in the model is due to 1 . This plot shows that inverted metallicity gradients may or may not be in equilibrium.

Here, we investigate the case where sf = such that there is no radial inflow of gas through the disc (see equation 36) 11 . Figure 13 shows the radial profile of eqbm in this case. It is 11 sf > is not possible in equilibrium in the Krumholz et al. (2018) model. clear that eqbm H(0) and eqbm dep,H 2 , especially at low , however, the exact values are sensitive to the choice of 1 . The reason for long metal equilibration timescales in this case is that, in the absence of advection, only diffusion and accretion are available to balance the source term. However, diffusion is weak due to the low gas dispersion ( 0 Σ 0 ∝ 3 ), and accretion is weak due to the low halo mass ( Σ cos0 ∝ 1.1 h ). Thus, metallicity gradients may not attain equilibrium in the absence of radial gas inflows in local dwarfs, whereas even a small amount of advection is sufficient to restore metallicity equilibrium (cf. Figure 4) . In the case where there is no accretion, one can expect a diverse range of metallicity gradients that are not constrained by the model. Therefore, caution must be exercised while studying metallicity gradients in such dwarfs with an equilibrium model.

At face value, this result might seem consistent with that of Forbes et al. (2014b) , where the authors find that dwarf galaxies do not attain statistical equilibrium within a Hubble time (see their Figure 15 ; see also , Feldmann 2015; Dashyan & Dubois 2020) . However, the equilibrium scenarios considered by Forbes et al. and us are not necessarily the same, and one is not a precondition of the other. Forbes et al. discuss the equilibrium for the total amount of gas or metals in a galaxy, which is a balance between inflow and outflow. The time required to reach this equilibrium is not necessarily the same as the time to equilibrate the distribution of metals within the galactic disc, for a given total metal content. Thus, one equilibrium time can be longer or shorter than the other.

Similarly, comparing eqbm with the metal correlation timescale for local dwarfs from Krumholz & Ting (2018) , which is the time required for diffusion alone to smooth out the metallicity distribution in the azimuthal direction, reveals that azimuthal metal distribution in these galaxies reaches equilibrium substantially more quickly than the radial distribution that we study here. This is consistent with the findings of Petit et al. (2015) , who also find that metal distributions equilibrate much more quickly in the azimuthal than the radial direction.

Local ULIRGs are very dynamically active, and are well-known to be undergoing major mergers or have companions (Lawrence et al. 1989; Melnick & Mirabel 1990; Leech et al. 1994; Clements et al. 1996; Veilleux et al. 1999 ). These galaxies are often characterized by strong starburst and/or AGN-driven outflows (Veilleux et al. 1995 (Veilleux et al. , 2013 Soto et al. 2012; Arribas et al. 2014) . They also have extremely short orbital timescales (of the order of ∼ 5 Myr, . Local ULIRGs are very compact, with discs extending out only to 2-3 kpc (Downes & Solomon 1998; Rujopakarn et al. 2011) . It is quite challenging to extract gas metallicities in these galaxies because the ionised gas emission lines are often dominated by shocks (Monreal-Ibero et al. 2006 and AGN activity (Ellison et al. 2013) , which interfere with traditional photoionisation-based metallicity diagnostics. In addition, high levels of dust obscuration make it difficult to model the emission line spectra (García-Marín et al. 2009; Nagao et al. 2011; Piqueras López et al. 2013; Stierwalt et al. 2014) . For these reasons, there are only a handful of studies that have been able to extract gas metallicities in local ULIRGs (e.g., Kewley et al. 2006; Monreal-Ibero et al. 2007; Arribas et al. 2008; Rupke et al. 2008; Westmoquette et al. 2012; Kilerci Eser et al. 2014; Pereira-Santaella et al. 2017) , and to the best our knowledge the only published studies of the metallicity gradient in ULIRGs are those of Rich et al. (2012, see their Figure 2 ) and Thorp et al. (2019) .

The short orbital timescales of ULIRGs ensure that they return to dynamical equilibrium quickly compared to their merger timescales, which based on simulations are estimated to be merge ∼ 0.3 − 1 Gyr (Jiang et al. 2008; Torrey et al. 2012 ). Thus our dynamical equilibrium model from is applicable to them. We investigate whether the metallicity distribution is also in equilibrium in Figure 14 , which shows eqbm for local ULIRGs. It is clear that eqbm ∼ merge , thus, metallicity may or may not be in equilibrium during the entire process of a merger. Our results corroborate those of Davé et al. (2012) , who argue that merging galaxies should not be in equilibrium because tidal flows will fuel star formation (Barton et al. 2000; Kewley et al. 2006; Reichard et al. 2009; Perez et al. 2011; Ellison et al. 2013; Moreno et al. 2020) , making cosmic accretion irrelevant. We show this quantitatively in Figure 15 , where advection (radial transport of gas due to tidal inflows) is the dominant term that sets the non-equilibrium metallicity distribution, and cosmic accretion is insignificant in comparison. Our results are also in line with those from simulations and observations where metallicity gradients in local ULIRGs are observed to continuously evolve and flatten as the merger progresses (Rich et al. 2012 , Figure 4) , implying that the metallicity distribution is not in a steady-state. This also implies that non-equilibrium gradients in local ULIRGs are transient; assuming the galaxy settles back to being a quiescent disc after the merger, the metallicity gradient will return to the equilibrium value for a spiral galaxy on the ∼ few Gyr equilibrium timescale for local spirals (cf. Figure 1) .

Given a merger rate, we can estimate the fraction of galaxies as a function of redshift that are expected to be out of metal equilibrium as 1 − − , where is the product of the merger rate and the metallicity equilibration timescale. Following Rodriguez-Gomez et al. (2015, Figure 9 ), we see that the observed average merger rate for massive galaxies ( ★ ≥ 10 10 M ) at = 0 is less than 0.06 Gyr −1 ), so we expect less than 20 per cent of massive galaxies to be out of metal equilibrium at redshift zero. Similarly, based on available observational results that find a merger rate of 0.5 Gyr −1 at ≈ 2 (Bluck et al. 2009 (Bluck et al. , 2012 Man et al. 2012) , we expect less than 40 per cent of the most massive galaxies ( ★ ≥ 10 11 M ) to be out of metal equilibrium at redshift two. The larger fraction of galaxies that are expected to be out of metal equilibrium at high redshift could explain the inverted gradients seen in high-observations, a topic we explore in Section 5.2.3.

Recent observations have discovered the presence of inverted (positive) gas phase metallicity gradients in galaxies Belfiore et al. 2017; Mingozzi et al. 2020) , especially at high redshift (Cresci et al. 2010; Queyrel et al. 2012; Stott et al. 2014; Carton et al. 2018; Wang et al. 2019; Curti et al. 2020; Wang et al. 2020; Simons et al. 2020) . Inverted gradients reflect the possibility of galaxies deviating from the classical, inside-out formation picture, at least temporarily. The three leading mechanisms that are believed to give rise to an inverted gradient are: (1.) substantial metal mass loading or merger-induced tidal flows of metal-poor gas that deprives the galaxy centre of metals, especially in dwarfs (Kereš et al. 2005; Kewley et al. 2006; Emerick et al. 2018 Emerick et al. , 2019 Chisholm et al. 2018; Tissera et al. 2019; see, however, Wilson et al. 2019 ), (2.) re-accretion of ejected metals at the outer edge of the disc from the CGM through cold, metal-rich flows or galactic fountains (Birnboim & Dekel 2003; Dekel & Birnboim 2006; Dekel et al. 2009a; Cresci et al. 2010; Crighton et al. 2013; Suresh et al. 2015) , and (3.) cosmic accretion of metal-poor gas at the centre that dilutes the central metallicity (Cresci et al. 2010 ).

Corresponding to these three scenarios, we can produce inverted gradients in our model by coupling a moderate or high value of Z GCM (i.e., addition of metal-rich gas to galaxy outskirts) with small values of or large values of A (corresponding to depressed central metallicity due to heavy metal loss or rapid dilution by metal-poor gas, respectively). However, any inverted gradients that we get from the model are sensitive to our choice of Z CGM , in the sense that we never get an inverted gradient for a sufficiently low value of Z CGM . Nevertheless, regardless of the value of Z CGM that we adopt, the resulting inverted gradients may or may not be in equilibrium. We illustrate this in Figure 16 , where we plot eqbm for local dwarfs, high-discs (identical to highgalaxies we discuss in Section 3.3), and high-dwarfs. We introduce the latter category by combining fiducial parameters for local dwarfs and high-galaxies from Table 2 in the following manner:

= 0.5, sf = 7 km s −1 , = 40 km s −1 , = 80 km s −1 , , = , = 0.9, sf = 0.4, and max = 4 at = 2. The three colors in all the panels in Figure 16 correspond to low = 0.05 (with Z CGM = 0.1), high Z CGM = 0.5 (with = 0.1), and high accretion where we multiply our fiducial values of A by 3 (with = 0.1, Z CGM = 0.1), respectively. The shaded regions correspond to the allowed values of 1 based on the constraints we introduced in Section 2.3.

We see that whether galaxies with inverted gradients are likely to be in equilibrium or not depends largely on what produces the inversion. Galaxies where the gradient inverts due to rapid accretion (high A) have relatively short values of eqbm , and may be in equilibrium as long as the accretion lasts, while those that invert due to an influx of metal-rich gas at their outskirts (high Z GCM ) are almost certainly out of equilibrium; galaxies with extremely efficient metal loss (low ) are intermediate, and may or may not be in equilibrium. Regardless of these details, the fact that many inverted gradients are not in equilibrium also hints at the possibility of them being transient (see also , Schönrich & McMillan 2017) . This is because subsequent star formation in the galaxy centre (due to cold gas accretion or re-accretion of enriched gas from the CGM) will replenish the metal supply on timescales comparable to the star formation timescale, thus leading to the formation of a negative gradient again. Hence, we expect inverted gradients to be erased within a star formation timescale ( 2 Gyr for massive galaxies, Leroy et al. 2008 ) unless they are re-established on a similar timescale. Since the processes that can cause inverted gradients (strong fountains, mergers, sudden accretion events, etc.) tend to wane with redshift, we expect that most massive galaxies will establish negative gradients by = 0, though some dwarfs, which have longer equilibration (and star formation) timescales, might retain their inverted gradients to = 0 or close to it.

In this work, we present a new theoretical model to explain the occurrence and diversity of gas phase metallicity gradients in galaxies. Starting from the conservation of metal mass, we incorporate major physical processes that can impact the distribution of metals in galaxies, namely, metal production, consumption, loss, advection, accretion and diffusion. Our first-principles based model shows that the radial metallicity gradients observed in galaxies are a natural consequence of inside-out galaxy formation. The equilibrium metallicity evolution model we present is a standalone model, but it requires inputs from a galaxy evolution model to set the galaxy properties that control metallicity. This intricate link between gas and metallicity lets us directly predict the evolution of metallicity gradients without ad hoc assumptions about galaxy properties.

The evolution of metallicities in our model depends on four dimensionless ratios: T , P, S, and A. These describe the ratio of the orbital timescale to the diffusion timescale, advection to diffusion, production (and metal ejection) to diffusion, and cosmic accretion to diffusion, respectively. Based on the input galaxy evolution model , we show how these ratios depend on various properties of the gas (cf. equation 37 − equation 40). The resulting second order differential equation of the radial distribution of metallicity has a simple analytic solution given by equation 41 that we use to predict a possible range of metallicity gradients as a function of galaxy properties. We use this capability to predict the metallicity gradients of local spirals, local dwarfs, and high-redshift disc galaxies, and to predict the evolution of metallicity gradients in galaxies with redshift. Below, we list our main results:

(i) The time required for the metal distribution within a galaxy to reach equilibrium is smaller than the Hubble time and comparable to the molecular gas depletion time in local spirals, (most) local dwarfs, and rotation-dominated high-galaxies. Thus, for most galaxies over most of cosmic time, the gas phase metallicity gradient is in equilibrium. Exceptions to this general trend can include merging galaxies, galaxies with inverted metallicity gradients, and some very low-mass local dwarf galaxies.

(ii) Galaxies tend to approach a particular value of central metallicity, dictated by the balance between the two dominant processes that depend on the properties of the galaxy (see below). The central metallicities we predict agree well with observations.

(iii) In local spirals, the two dominant processes shaping the metallicity gradient are metal production (S), which tries to steepen the gradient, and accretion of metal-poor gas (A), which tries to flatten it. On the other hand, metallicity gradients in local dwarfs and high-galaxies are set by the balance between S and advection of metal-poor gas from the outer to the inner parts of galaxies (P).

(iv) One crucial free parameter that emerges from our model is the "yield reduction factor" , defined as the fraction of supernovaproduced metals that mix with the ISM rather than being lost immediately in metal-enhanced galactic winds. While metallicity gradients in local spirals are not tremendously sensitive to , it has a significant effect on the metallicity gradients in local dwarfs and high-galaxies. also impacts the absolute metallicities in all galaxies. Comparison of the model with observations reveals that massive galaxies prefer a high value of , whereas low-mass galaxies prefer a lower value of . Thus, the model predicts that low-mass galaxies undergo more preferential metal ejection, and should have more metal-enriched winds than massive galaxies. Future work should thus focus on constraining from observations.

As a first application of our model, we study the evolution of metallicity gradients with redshift, both within a single galaxy and over samples of galaxies at different redshifts selected to have matching stellar masses or comoving densities. Our model shows that gradients in Milky Way-like galaxies have steepened over time, in qualitative agreement with recent observations; quantitative agreement between the model and the data requires a scaling of such that was low for the Galaxy in the past as compared to today, consistent with that seen in simulations. We also predict the existence of specific signatures for the evolution of metallicity gradient with redshift as a function of stellar mass that can be tested with future surveys. We show that both the Milky Way in particular and disc galaxies in general transition from the advection-dominated (P > A) to the accretion-dominated (P < A) regime from high to low redshifts. This transition mirrors the transition from gravitydriven to star formation-driven turbulence from high to low redshifts . In companion papers, we show that this transition (along with ) is also responsible for driving the shape of the mass-metallicity relation and the mass-metallicity gradient relation (Sharda et al. 2021a) in the local Universe, and we also apply our model to explain the relationship between metallicity gradients and gas kinematics in high redshift galaxies (Sharda et al. 2021b ).

We dedicate this work to the medical services personnel who have been working tirelessly to combat COVID-19 across the world; this work would not have been possible without their sincere efforts in keeping the community safe. We thank the anonymous reviewer for a thorough referee report that helped improve the manuscript. We also thank Meridith Joyce, Stephanie Monty, J. Trevor Mendel, Lisa Kewley, Kenneth Freeman, and Raymond Simons for several useful discussions. We are grateful to Xiangcheng Ma for sharing results of the FIRE simulations, and to Mirko Curti for sharing their data compilation on metallicity gradients in high-galaxies. PS is supported by the Australian Government Research Training Program (RTP) Scholarship. MRK and CF acknowledge funding provided by the Australian Research Council (ARC) through Discovery Projects DP190101258 (MRK) and DP170100603 (CF) and Future Fellowships FT180100375 (MRK) and FT180100495 (CF). MRK is also the recipient of an Alexander von Humboldt award. PS, EW and AA acknowledge the support of the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. JCF is supported by the Flatiron Institute through the Simons Foundation. Analysis was performed using numpy (Oliphant 2006; Harris et al. 2020 ) and scipy (Virtanen et al. 2020); plots were created using Matplotlib (Hunter 2007) and astropy (Astropy Collaboration et al. 2013 . This research has made extensive use of NASA's Astrophysics Data System Bibliographic Services. This research has also made extensive use of Wolfram|Alpha 12 and Mathematica for numerical analyses, the Wright (2006) cosmology calculator, and the image-to-data tool WebPlotDigitizer 13 .

No data were generated for this work. Figure A1 but for high-galaxies. Similar to local dwarfs, changing the functional form of ★ ( ) has no impact on the metallicity distributions in high-galaxies.

Our predictions of metallicity gradient evolution depend on a scaling of with derived from high-galaxy observations (Wisnioski et al. 2015) , and a scaling of h with derived from cosmological simulations (Neistein & Dekel 2008; Faucher-Giguère et al. 2011) . Since these two methods of deriving the scaling are very different, it is important to comment on any possible discrepancies between the two, and if they affect our results. Indeed, there is what appears at first glance to be an inconsistency: as we noted in Section 2.2, cosmological equilibrium demands that the inflow rate through the disc (and, the star formation rate ★ ) be similar to or less than the accretion rate onto the galaxy ext , in order to conserve the total mass. In terms of our model, the above condition translates into P/A ln max for ★ = 1/ 2 . However, in many cases, our adopted scalings of and ℎ with give considerably larger values of P/A at high-. This discrepancy is simply a manifestation of the known problem that galaxies at ∼ 2 have star formation rates ★ > ext (Erb 2008; Behroozi et al. 2013; Scoville et al. 2017 ); high-galaxies obey the same observed scaling between star formation rate and velocity dispersion as local galaxies Varidel et al. 2020) , and since the inflow rate is directly set by , * > ext directly implies > ext . The discrepancy between star formation rates (and velocity dispersions) and expected cosmological accretion rates has several possible explanations, but from the standpoint of our model for metallicity gradients, we can divide these into two main categories. One is that galaxies near the epoch of peak star formation do in fact form stars and move mass inward faster than their mean cosmological accretion rates, either because mass that was ejected at an earlier epoch falls back onto the galactic outskirts (e.g., Christensen et al. 2016; Tollet et al. 2019) , or because of large angular momentum mismatch between the infalling material and the disc that triggers a sufficiently large radial inflow (Mayor & Vigroux 1981; Bilitewski & Schönrich 2012; Pezzulli & Fraternali 2016) , or because galaxies accumulate large gas reservoirs at 2, which then flow into the star-forming portion of the disc due to compaction events (Dekel et al. 2009b; Dekel & Burkert 2014) , interactions (Rupke et al. 2010a) or mergers (Hani et al. 2018) at ∼ 2. In these cases, the model we present is sufficient and we do not need to make any changes, since in such cases galaxies can maintain a large P/A for a considerable time.

The other possibility is that the ★ measured at high-are overestimated (e.g., Leja et al. 2019 Leja et al. , 2020 , thus altering the − * relationship; this is functionally equivalent to overestimating at fixed ★ . Such an overestimate in high-redshift galaxies could plausibly be due to beam smearing, inclination uncertainty, and similar resolution-dependent (and thus redshift-dependent) factors (Davies et al. 2011; Di Teodoro & Fraternali 2015) . This would have a significant impact on P/A because P/A ∝ 3 . To study the effects this can have on our results, we reproduce the expected cosmic evolution of metallicity gradient for the Milky Way in the model (cf. Figure 10 ), but using obtained from equation 53 reduced by a factor of (1 + ), in line with the redshift scaling of the limits of spectral resolution of different instruments (McDermid et al. 2020, Figure 2 .3). Figure B1 presents the resulting metallicity gradients from the model, where the maximum log P/A is only 0.6. The qualitative shape of the model changes slightly for > 0.5. However, remains the primary factor that drives the model gradients closer to the observed gradients across cosmic time. The main Figure B2 . Metallicity gradients as shown in the main text in Section 3, but with an alternate value of the edge of the star-forming disc, max . The profile of the distribution is preserved in each case, with slight variations in the absolute metallicities and metallicity gradients, with diminishing differences for decreasing . Larger galaxies (in each galaxy class) show higher mean metallicity and flatter gradients.

difference from our fiducial model is that steep gradients become possible at high-if is close to unity, because a smaller P/A implies weaker homogenisation of the ISM by inward advection of metal-poor gas.

In the main text, we non-dimensionalize the solution in terms of , where = / 0 and we adopt a fiducial value of 0 = 1 kpc, finding solutions for in the range ( min , max ), with min = 1 and max chosen based on observations of the sizes of the star-forming region for different galaxies. To explore the sensitivity of our results to the choice of our range in , we show in Figure B2 how the equilibrium gradients change in local spirals if we increase max from the fiducial 15 used in the main text to 20, noting that the qualitative trends remain the same across all types of galaxies. Increasing max leads to higher metallicities at each location in the disc, with slightly higher mean metallicities and shallower metallicity gradients. This is because increasing max decreases A since it means that there is a larger disc for the same total cosmic accretion rate (see equation 34 and equation 40). Thus, the ratio S/A that appears in Z (see equation 41) increases, giving higher Z( ). Further, this increment in S/A is reduced if < 1 because S ∝ . Thus, for lower , changing max does not lead to any appreciable change in Z( ). Similarly, if we shift the inner edge of the galactic disc (where the rotation curve flattens out) by decreasing it to min = 0.5 (or increasing it to 2), the solution allows for slightly higher (lower) mean metallicities, and steeper (shallower) gradients.

Thus, this analysis implies that in galaxies where the transition from the star-forming disc to the bulge occurs at smaller galactic radius, the galaxy could potentially build slightly steeper metallicity gradients. Conversely, in galaxies where the star-forming disc is larger, we expect slightly shallower metallicity gradients. However, these variations remain insignificant compared to the scatter introduced due to other parameters, particularly . Thus we do not regard variations in min or max as a substantial uncertainty in the model. This paper has been typeset from a T E X/L A T E X file prepared by the authors.

Mathematical Methods for Physicists, third edn

NIRSPEC Instrument Science Team JAESs Collaboration

American Institute of Physics Conference Series

Advanced Engineering Mathematics, tenth edn

Galaxy Formation and Evolution

Numerical heat transfer and fluid flow

The Magellanic Clouds

The Local Group as an Astrophysical Laboratory

Here, we describe how the solutions change if we pick a different functional form for the radial profile of cosmological accretion, ★ ( ). Note that we must numerically solve for these functional forms, because analytic solutions either do not exist or are so complex in functional form that a numerical integration is preferable. Specifically, we experiment with ★ ( ) = 1/ and ★ ( ) = 1. Figure A1 , Figure A2 , and Figure A3 show metallicity profiles with different ★ ( ) for local spirals, local dwarfs, and highgalaxies with = 1, respectively; we use = 1 because this maximises the dependence on ★ ( ) -smaller values suppress variations. Note that the dimensionless parameters P and S are identical to that used in the Section 3 for the corresponding galaxies, but that A differs due to its dependence on ★ (equation 40). We see that changing the profile of ★ ( ) has no noticeable effect for local dwarfs or high− galaxies. This is because cosmic accretion is not a dominant term in the metallicity model for these galaxies, and the metallicities are instead mainly set by source, advection and diffusion.The profile of ★ ( ) does matter for local spirals; as ★ ( ) flattens, metallicities in the inner regions of the disc reach higher values whereas the outer regions of the disc become more metal poor, thus leading to somewhat steeper gradients. For the most extreme case of constant ★ ( ), gradients are ≈ −0.05 dex kpc −1 steeper than our fiducial model. Thus, if the cosmic accretion profiles in local spirals are flatter than that we use in the main text, we expect slightly steeper gradients from the metallicity model, which is largely due to the SFR that we input from the galaxy evolution model. This is because under the input galaxy evolution model of where the SFR varies as 1/ 2 , flatter accretion profiles will dilute the metallicity in the metal-deficient outer regions by the same amount as that in the metal-rich inner regions, thus giving a larger difference between metallicities in the inner and outer regions in the disc.