3-D gas-phase elemental abundances across the formation histories of Milky Way-mass galaxies in the FIRE simulations: initial conditions for chemical tagging
Matthew A. Bellardini, Andrew Wetzel, Sarah R. Loebman, Claude-André Faucher-Giguère, Xiangcheng Ma, Robert Feldmann
MMNRAS , 1–23 (2021) Preprint 15 February 2021 Compiled using MNRAS L A TEX style file v3.0
Matthew A. Bellardini, (cid:63) Andrew Wetzel, Sarah R. Loebman, , † Claude-Andr´e Faucher-Gigu`ere, Xiangcheng Ma, and Robert Feldmann Department of Physics & Astronomy, University of California, Davis, One Shields Ave, Davis, CA 95616, USA Department of Physics, University of California, Merced, 5200 Lake Road, Merced, CA 95343, USA Department of Physics and Astronomy and CIERA, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, CA 94720, USA Institute for Computational Science, University of Zurich, Zurich CH-8057, Switzerland
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We use FIRE-2 simulations to examine 3-D variations of gas-phase elemental abun-dances of [O/H], [Fe/H], and [N/H] in 11 MW and M31-mass galaxies across theirformation histories at z ≤ . t lookback ≤ . (cid:46) .
008 dex at all z ≤ .
5. We find negative radial gra-dients (metallicity decreases with galactocentric radius) at all times, which steepenover time from ≈ − .
01 dex kpc − at z = 1 ( t lookback = 7 . ≈ − .
03 dex kpc − at z = 0, and which broadly agree with observations of the MW, M31, and nearbyMW/M31-mass galaxies. Azimuthal variations at fixed radius are typically 0 .
14 dex at z = 1, reducing to 0 .
05 dex at z = 0. Thus, over time radial gradients become steeperwhile azimuthal variations become weaker (more homogeneous). As a result, azimuthalvariations were larger than radial variations at z (cid:38) . t lookback (cid:38) . (cid:46) .
05 dex) across aradial range of ∆ R ≈ . z (cid:38) R ≈ . z = 0. We also measurefull distributions of elemental abundances, finding typically negatively skewed normaldistributions at z (cid:38) z = 0. Ourresults on gas abundances inform the initial conditions for stars, including the spatialand temporal scales for applying chemical tagging to understand stellar birth in theMW. Key words: galaxies: abundances – galaxies: formation – galaxies: ISM – ISM: abun-dances – stars: abundances – methods: numerical
Many current and future observational surveys of starsacross the Milky Way (MW) seek to unveil the MW’s forma-tion history in exquisite detail. Current surveys, such as theRAdial Velocity Experiment (RAVE; Steinmetz et al. 2006),the Gaia-ESO survey (Gilmore et al. 2012), the Large AreaMulti-Object Fiber Spectroscopic Telescope (LAMOST; Cuiet al. 2012), GALactic Archaeology with Hermes (GALAH;De Silva et al. 2015), and the Apache Point Galactic Evo- (cid:63)
E-mail: [email protected] † Hubble Fellow lution Experiment (APOGEE; Majewski et al. 2017) havemeasured elemental abundances of hundreds of thousands ofstars. Future surveys, such as the WHT Enhanced Area Ve-locity Explorer (WEAVE; Dalton et al. 2012), the SubaruPrime Focus Spectrograph (PFS; Takada et al. 2014), theSloan Digital Sky Survey V (SDSS-V; Kollmeier et al. 2017),the 4-metre Multi-Object Spectrograph Telescope (4MOST;de Jong et al. 2019), and the MaunaKea Spectroscopic Ex-plorer (MSE; The MSE Science Team et al. 2019) will in-crease the samples of spectroscopically observed stars intothe millions. A key science driver for these surveys is ‘galac-tic archaeology’: to infer the history of the MW using obser- © a r X i v : . [ a s t r o - ph . GA ] F e b M. A. Bellardini et al. vations of the dynamics and elemental abundances of starstoday.Measurements of stellar dynamics can provide detailedinformation on the MW’s properties and formation history,but the fundamental limitation is that stellar orbits canchange over time via mergers, accretion, scattering, andother dynamical perturbations (e.g. Abadi et al. 2003; Brooket al. 2004; Sch¨onrich & Binney 2009; Loebman et al. 2011).However, a star’s atmospheric elemental abundances will notchange in response to these dynamical processes, providing akey orbital-invariant ‘tag’. ‘Chemical tagging’, introduced inFreeman & Bland-Hawthorn (2002), thus provides tremen-dous potential to infer the formation conditions of a star ofarbitrary age.‘Strong’ chemical tagging represents the most fine-grained scenario, to identify stars born in the same starcluster (e.g. Price-Jones et al. 2020). By contrast, ‘weak’chemical tagging seeks to infer the general location and timewhere/when a stellar population formed, for example, to as-sociate populations of stars to certain birth regions of thegalaxy (e.g. Wojno et al. 2016; Anders et al. 2017) or thataccreted into the MW from galaxy mergers (e.g. Ostdieket al. 2020).Both regimes of chemical tagging rely on sufficientlyprecise measurements of stellar abundances and on assump-tions about the elemental homogeneity (to the measured pre-cision) of the gas from which the stars formed. For example,strong chemical tagging of individual star clusters relies onboth the internal homogeneity of the gas cloud from whichthe stars formed, and on how unique the abundance pat-terns were in that cloud across space and time. Observa-tional evidence of open star clusters suggests the first cri-terion is met (De Silva et al. 2007; Ting et al. 2012; Bovy2016) to measurable precision. Regarding the latter crite-rion, observations of the MW and external galaxies showradial and azimuthal variations in abundances across thedisk (e.g. S´anchez-Menguiano et al. 2016; Moll´a et al. 2019b;Wenger et al. 2019; Kreckel et al. 2020), although more workis needed to understand these spatial variations in the con-text of chemical tagging. Weak chemical tagging is subjectto the same assumptions but applied to larger regions ofgas across the disk (or in accreting galaxies). For example,if all gas in the disk was measurably homogeneous in allabundances at a given time, chemical tagging would offerno spatially discriminating power. Conversely, the limit ofextreme clumpiness, in which each star cluster formed witha measurably unique abundance pattern, would in principleenabled detailed chemical tagging, but it significantly wouldcomplicate the modeling.
Thus, a key question for chemical tagging is: what arethe relevant spatial scales of measurable homogeneity of starsforming at a given time, and how does this evolve acrosscosmic time? . Bland-Hawthorn et al. (2010) previously ex-plored this via a toy model, where they show all star clusters (cid:46) M (cid:12) and a large fraction of clusters with mass be-low ∼ M (cid:12) are expected to be internally homogeneous.Progress in chemical tagging requires addressing these ques-tions regarding stellar birth before examining the subsequentdynamical evolution of stars after they form.Many works have examined abundance variations ofstars across the MW, generally finding a negative radial gra-dient in abundances for stars near the plane of the disk, which flattens or turns positive at larger heights (e.g. Chenget al. 2012; Boeche et al. 2013, 2014; Anders et al. 2014; Hay-den et al. 2014; Mikolaitis et al. 2014; Anders et al. 2017).Furthermore, the MW has an observed negative vertical gra-dient in stars (Cheng et al. 2012; Carrell et al. 2012; Boecheet al. 2014; Hayden et al. 2014), although this slope variessignificantly between observations and depends on radius(Hayden et al. 2014). Both the radial and vertical gradientsvary with stellar age (Wang et al. 2019b). Additionally, Lucket al. (2006); Lemasle et al. (2008); Pedicelli et al. (2009)found evidence for azimuthal variations in the abundancesof young stars, which may result from patchy star forma-tion (Davies et al. 2009; Luck & Lambert 2011; Genovaliet al. 2014). Nieva & Przybilla (2012) also explored the ho-mogeneity of B-type stars in the solar neighborhood, within500 pc, of the sun and found scatter on the order of 0 .
05 dexfor [O/H] which they state is comparable to gas-phase abun-dance scatter out to 1 . (cid:46) .
05 dex) (e.g. Cedr´es& Cepa 2002; Zinchenko et al. 2016).Understanding how these variations change across cos-mic time is imperative for chemical-tagging models. Cur-rently no consensus exists, amongst observations (e.g. Curtiet al. 2020, and references therein) on the redshift evolu-tion of radial elemental abundance gradients, in part becauseof angular resolution limitations (Yuan et al. 2013), whichsome works have addressed via adaptive optics and grav-itational lensing (Jones et al. 2010; Swinbank et al. 2012;Jones et al. 2013). Furthermore, different works use differ-ent calibrators to measure abundances, which often disagree(Hemler et al. 2020), further complicating our understandingof spatial variations.Many theoretical works have used simulations to pre-dict the spatial distribution of gas-phase abundances andtheir evolution. As with observational efforts, there is noconsensus for the redshift evolution of abundance gradientsin theory (e.g. Moll´a et al. 2019a, and references therein).Gibson et al. (2013) compared cosmological simulations withMUGS (‘conservative’ feedback) and MaGICC (‘enhanced’feedback) run with the GASOLINE code and found that thestrength of feedback in simulations is critical for the evolu-tion of radial gradients of abundances, such that strongerfeedback leads to flatter gradients at all times, while galaxieswith weaker feedback have gradients that are steep at highredshift and flatten with time. Ma et al. (2017), studying the
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Figure 1.
Face-on image of all gas within ± (cid:46) .
02 dex at all times. The left panel emphasizes azimuthal variations by showing the deviation from the azimuthally averaged[O/H] at each radius , that is, subtracting off the radial gradient, thus highlighting the enhanced enrichment correlated with spiral arms.The right 3 panels show the deviation from the mean [O/H] of all gas at R ≤
15 kpc at each redshift. White regions have highly diffusegas in which we do not report a measured abundance. The radial gradient in [O/H] dominates over azimuthal variations at late times,but at early times the azimuthal variations are the most significant.
FIRE-1 suite of cosmological simulations, found that galax-ies exhibit a diverse range of radial gradients in abundances,and that these gradients can fluctuate rapidly from steep toshallow (in ∼
100 Myr) at high redshift, so measurements ofhigh-redshift gradients may not be indicative of long-termtrends. They found that galaxies tend to quickly build up anegative gradient once stellar feedback is no longer sufficientto drive strong outflows of gas. By contrast, analyzing star-forming galaxies in the TNG-50 cosmological simulation,Hemler et al. (2020) found that radial gradients in galax-ies are steep at high redshift and flatten with time. Severaltheoretical works also have examined azimuthal variations.Spitoni et al. (2019) developed a 2-D model for abundanceevolution that follows radial and azimuthal density varia-tions in a MW-like disk and found that azimuthal residualsare strongest at early times and at large radii. Using theirS2A model, Spitoni et al. (2019) found azimuthal residuals in[O/H] of ≈ . R = 8 kpc at t lookback = 11 Gyr whichevolve to ≈ .
05 dex at present day. Moll´a et al. (2019b) ex-plored azimuthal variations in a MW-like disk for 5 modelsof 2-D abundance evolution and found [O/H] variations thatare typically small (0 . − . ≈ .
12 dex at z = 0 in galaxies with M star = 10 − . M (cid:12) .Observations of azimuthal variations of abundances ingas in nearby MW-mass galaxies find scatter that is compa-rable to observational measurement uncertainty ( ∼ .
05 dex)(Zinchenko et al. 2016; Kreckel et al. 2019). This implies thatgas in galaxies is well mixed azimuthally at z = 0. Conse-quently, works modeling the abundance evolution of galax-ies, which inform chemical tagging, generally assume thatgas is well mixed azimuthally in the disk at all times (e.g.Minchev et al. 2018; Moll´a et al. 2019a; Frankel et al. 2020),such that the key spatial variation is radial. While galaxiesexhibit radial gradients across a range of redshifts (Queyrelet al. 2012; Stott et al. 2014; Wuyts et al. 2016; Cartonet al. 2018; Patr´ıcio et al. 2019; Curti et al. 2020), radialvariations may not always dominate. At early times, in par-ticular, azimuthal variations may be more important. Someabundance-evolution models have begun to explore both az- imuthal and radial variations across time (e.g. Acharovaet al. 2013; Moll´a et al. 2019b; Spitoni et al. 2019). Kawataet al. (2014); Grand et al. (2015), using N -body simulationsof MW-mass galaxies, and Baba et al. (2016), using bary-onic simulations run with the SPH code ASURA-2, foundthat gas exhibits streaming motion along spiral arms, whichcould contribute to 2-D abundance variations in gas.More detailed 2-D abundance-evolution models, whichaccount for density variations within the disk from spiralarms and bars, result in azimuthal variations in gas-phaseabundances. Moll´a et al. (2019b) found that arm / inter-arm abundance variations quickly dilute through interac-tions with spiral structure. Spitoni et al. (2019) also foundthat azimuthal variations dilute with time, but they foundthat the strength of azimuthal variations at z = 0 approxi-mately agree with the observational results of Kreckel et al.(2019). The simulation analysis of Grand et al. (2015); Babaet al. (2016) showed that, just like stars (Lynden-Bell &Kalnajs 1972), gas experiences radial migration as a resultof spiral structure. Grand et al. (2015) found that this sys-tematic streaming along spiral arms leads to metal-rich gasin the inner galaxy moving to larger radii and metal-poorgas in the outer galaxy moving to the inner galaxy, leadingto non-homogeneous abundances at a given radius. However,they found that gas elements quickly exchange abundancesafter migrating, leading to small azimuthal dispersions inabundance. S´anchez-Menguiano et al. (2020) found that az-imuthal variations in abundances are stronger in galaxieswith stronger bars and grand-design spirals, which sup-ports non-axisymmetric structure driving azimuthal inho-mogeneities.In this paper we use FIRE-2 cosmological simulationsof MW/M31-mass galaxies to explore the cosmic evolutionof 3-D abundance patterns of gas, as a first step towardsunderstanding the spatial and temporal scales of applyingchemical tagging in a cosmological context. Our analysisof gas represents our first step, to characterize the initialconditions for star-forming regions. In future work, we willexamine the resultant trends in stars and their dynamicalevolution across time. Here, we seek to quantify the 3-Dspatial scales over which elemental abundances of gas (andthus the formation of stars) are measurably homogeneous. MNRAS , 1–23 (2021)
M. A. Bellardini et al.
In Section. 2 we describe the simulations used for this anal-ysis. In Section. 3 we first explore the radial gradients andcompare them against observations of the MW, M31, andnearby MW/M31-mass galaxies. Next we examine the cos-mic evolution of radial, vertical, and azimuthal variations ingas-phase abundances, in particular, to understand whichdimension dominates the spatial variations at a given time.We also examine implications of gas (in)homogeneity on cur-rent and upcoming observations of the MW. Finally, we ex-amine full distributions of elemental abundances. Section. 4we summarize the main results of the paper and provide adiscussion of their implications.
We use a suite of MW/M31-mass cosmological zoom-insimulations from the Feedback In Realistic Environments(FIRE) project (Hopkins et al. 2018). We ran these simu-lations using the FIRE-2 numerical implementations of fluiddynamics, star formation, and stellar feedback. These simu-lations use the Lagrangian Meshless Finite Mass (MFM) hy-drodynamics method in Gizmo (Hopkins 2015). The FIRE-2model incorporates realistic gas physics through the inclu-sion of metallicity-dependent radiative heating and coolingprocesses such as free-free, photoionization and recombina-tion, Compton, photo-electric and dust collisional, cosmicray, molecular, metal-line, and fine structure processes, ac-counting for 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca,Fe) across a temperature range of 10 − K. The simu-lations also include a spatially uniform, redshift-dependentUV background from Faucher-Giguere et al. (2009). In calcu-lating metallicities throughout this paper, we scale elementalabundances to the solar values in Asplund et al. (2009).Star particles form out of gas that is self-gravitating,Jeans-unstable, cold (
T < K), dense ( n > − ),and molecular (following Krumholz & Gnedin 2011). Eachstar particle inherits the mass and elemental abundancesof its progenitor gas and represents a single stellar popula-tion, assuming a Kroupa (2001) stellar initial mass function.FIRE-2 evolves star particles along standard stellar popula-tion models from e.g. STARBURST99 v7.0 (Leitherer et al.1999), including time-resolved stellar feedback from core-collapse and Ia supernovae, continuous mass loss, radiationpressure, photoionization, and photo-electric heating. FIRE-2 uses rates of core-collapse and Ia supernovae from STAR-BURST99 (Leitherer et al. 1999) and Mannucci et al. (2006),respectively. The nucleosynthetic yields follow Nomoto et al.(2006) for core-collapse and Iwamoto et al. (1999) for Ia su-pernovae. Stellar wind yields, sourced primarily from O, B,and AGB stars, are from the combination of models fromvan den Hoek & Groenewegen (1997); Marigo (2001); Izzardet al. (2004), synthesized in Wiersma et al. (2009).Critical for this work, these FIRE-2 simulations alsoexplicitly model the sub-grid diffusion/mixing of elementalabundances in gas that occurs via unresolved turbulent ed-dies (Su et al. 2017; Escala et al. 2018; Hopkins et al. 2018). FIRE project web site: http://fire.northwestern.edu
In effect, this smooths abundance variations between gas el-ements, assuming that sub-grid mixing is dominated by thelargest unresolved eddies. (Escala et al. 2018) showed thatthis model yields significantly more realistic distributions ofstellar and gas-phase metallicities in galaxies. We explorethe robustness of our results to variations in the strength ofthe mixing/diffusion coefficient in Appendix C.All simulations assume flat ΛCDM cosmologies withparameters broadly consistent with the Planck Collabora-tion et al. (2018): h = 0 . − .
71, Ω Λ = 0 . − . m = 0 . − .
31, Ω b = 0 . − . σ = 0 . − . n s = 0 . − .
97. For each simulation we generated cos-mological zoom-in initial conditions embedded within cos-mological boxes of length 70 . −
172 Mpc at z ≈
99 usingthe code MUSIC (Hahn & Abel 2011). We saved 600 snap-shots from z = 99 to 0, with typical time spacing (cid:46)
25 Myr.We examine 11 MW/M31-mass galaxies from 2 suitesof simulations. We select only galaxies with a stellar masswithin a factor of ≈ ≈ × M (cid:12) (Bland-Hawthorn & Gerhard 2016). 5 of our galaxies are from the Latte suite of isolated individual MW/M31-mass halos, in-troduced in Wetzel et al. (2016). (We exclude m12w, be-cause it has an unusually compact gas disk at z = 0, with R gas90 = 7 . Latte galaxies have halo masses M =1 − × M (cid:12) , for which M refers to the total masswithin the radius containing 200 times the mean matterdensity of the Universe. These simulations have Dark Mat-ter (DM) particle masses of 3 . × M (cid:12) and initial bary-onic particle masses of 7070 M (cid:12) (though because of stellarmass loss, star particles at z = 0 typically have masses of ∼ (cid:12) ). Star and DM particles have fixed gravitationalforce softening lengths of 4 and 40 pc (Plummer equivalent),comoving at z > M = 1 − × M (cid:12) ,with ≈ × better mass resolution than the Latte suite.In general, we find few systematic differences in any ofour results for the isolated galaxies versus the LG-like pairs,the only notable difference being the relative strength ofazimuthal scatter to radial gradient strength at large radiusand high z , so we combine these suites in all of our results. Fig. 1 shows face-on images of the gas disk of one of oursimulations, Romeo, at several redshifts. We color-code gasby its variation in [O/H], to visualize key trends that weexplore in this work. We do not show results for [Fe/H],because they are qualitatively consistent with [O/H]. Theleft panel shows the deviation of the local [O/H] from themean value at each radius for radial bins of width 200 pc at z = 0, that is, we subtract off the overall radial gradient.This highlights the variations along the azimuthal directionat each radius, showing enhancement in [O/H] along spiralstructure (Orr et al. in prep. will present a detailed analysisof metallicity enhancements along spiral arms).The right 3 panels show the difference between the lo-cal [O/H] and the mean [O/H] across all gas at R ≤
15 kpc
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Table 1.
Properties of the stellar and gas disks of our simulated MW/M31-mass galaxies z = 0. The first column lists the name ofthe galaxy: ‘m12’ indicates isolated galaxies from the Latte suite, while the other galaxies are LG analogues from the ELVIS on FIREsuite. R and R is the radius where the cumulative mass of the disk is 25% and 50%, respectively, within a height ± R and Z as the radius and height where thecumulative mass of the stellar/gas disk are 90% of the total mass of stars/gas within a sphere of 20 kpc. M is the total stellar/gas masscontained within both R and Z . The gas fraction, f gas90 , is the ratio of gas mass to total baryonic mass within R star90 and Z star90 .simulation M star90 [10 M (cid:12) ] R star25 [ kpc] R star50 [ kpc] R star90 [ kpc] Z star90 [ kpc] M gas90 [10 M (cid:12) ] R gas25 [ kpc] R gas50 [ kpc] R gas90 [ kpc] Z gas90 [ kpc] f gas90 ( < R star90 , < Z star90 )m12m 10.0 1.9 4.3 11.6 2.3 2.1 6.6 10.3 15.0 1.2 0.13Romulus 8.0 1.2 3.2 12.9 2.4 2.7 9.0 13.1 18.3 2.3 0.16m12b 7.3 1.0 2.2 9.0 1.8 1.7 6.4 9.6 15.0 1.5 0.11m12f 6.9 1.2 2.9 11.8 2.1 2.3 8.7 12.6 17.8 2.4 0.14Thelma 6.3 6.3 3.4 11.2 3.2 2.6 7.5 12.1 17.6 3.1 0.16Romeo 5.9 1.6 3.6 12.4 1.9 1.8 8.0 12.2 18.1 1.5 0.16m12i 5.3 1.1 2.6 9.8 2.3 1.7 7.1 10.2 16.7 1.7 0.15m12c 5.1 1.3 2.9 9.1 2.0 1.5 5.4 8.4 14.6 2.4 0.15Remus 4.0 1.2 2.9 11.0 2.2 1.5 7.4 11.8 18.0 1.8 0.22Juliet 3.3 0.8 1.8 8.1 2.2 1.5 6.9 11.3 18.6 3.1 0.14Louise 2.3 1.2 2.8 11.2 2.2 1.4 8.0 12.6 18.5 2.0 0.34 and within ± z (cid:38) . t lookback (cid:38) . relative abundance varia-tions, including spatial variations, evolution, and the shapesof abundance distributions. Future work (e.g. Gandhi et al.in prep.), will explore the effects of varying stellar models,such as supernova Ia rates, on abundances in FIRE-2 galax-ies. First we examine the radial profiles of [O/H], [Fe/H], and[O/Fe] in gas for all 11 galaxies at z = 0. We time-averageeach galaxy’s profile across ∼
50 Myr by stacking 3 snap-shots to reduce short-time fluctuations. We present all re-sults in physical radii; in Appendix A, we examine thesetrends scaling to various galactic scale radii, finding thatthe host-to-host scatter in our suite is minimized when ex-amining gradients in physical units. These profiles containall gas within a vertical height Z ± .
25 kpc. We calculate the mass-weighted mean of the gas-phase abundance in each bin. Weshow profiles out to R = 15 kpc; our gas disks generally ex-tend beyond this radius, but because our primary motivationis chemical tagging, we examine only regions with significantstar formation (see Table 1).Fig. 2 (top 2 panels) shows that [O/H] and [Fe/H]decrease monotonically with radius. The mean gradient is ≈ − .
03 dex kpc − for both [O/H] and [Fe/H], and the meanchange in abundance from 0 −
15 kpc is ≈ .
51 dex for [O/H]and ≈ .
56 dex for [Fe/H]. These negative gradients in gasreflect the decreasing ratio of stars (the sources of enrich-ment) to gas towards the outer disk, and show that thesegas disks are not radially well mixed at z = 0.Across our 11 galaxies, the host-to-host standard de-viation is ≈ .
09 dex for [O/H] and ≈ .
07 dex for [Fe/H].The legend of Fig. 2 lists the host galaxies in decreasingorder of stellar mass, highlighting that the abundance ata given radius correlates strongly with the galaxy’s mass.Table. 1 shows that stellar mass drops by a factor of ∼ ≈ . ≈ .
24 dex, almost exactly that in Fig. 2. In other words,the scatter across our suite primarily reflects the mass-metallicity relation (see Lequeux et al. 1979; Tremonti et al.2004; Mannucci et al. 2010; Andrews & Martini 2013).In Fig. 2, dashed lines show the LG-like hosts, and whilethey show typically lower abundance at a fixed radius thanthe isolated hosts, this is because they have somewhat lowerstellar mass on average. We find no systematic differencesbetween LG-like and isolated hosts beyond this, despite thefact that the LG-like hosts form their in-situ stars systemati-cally earlier than the isolated hosts (Santistevan et al. 2020).Thus, this difference in formation history does not imprintitself on gas-phase abundances at z ≤ .
5. As a result, wewill combine these samples in all subsequent analyses.Fig. 2 (bottom panel) shows profiles for [O/Fe], whichare nearly flat at all radii. The mean change in [O/Fe] from0 −
15 kpc is ≈ − .
046 dex. [O/Fe] shows the strongest (posi-tive) gradient in the inner ≈ ≈ ≈ .
027 dex for [O/Fe].We measure this for [Mg/Fe] (not shown here) as anothertracer of core-collapse vs Ia supernovae enrichment. The
MNRAS000
MNRAS000 , 1–23 (2021)
M. A. Bellardini et al.
Figure 2.
Radial profiles of gas-phase elemental abundanceswithin the disks of our 11 galaxies at z = 0, listed by decreasingstellar mass. Each simulation includes all gas, averaging across 3snapshots ( ≈
50 Myr) within a disk height of ± −
15 kpc of ∼ .
51 dexfor [O/H] and ∼ .
56 dex for [Fe/H]. The [O/Fe] profile is approx-imately flat at radii (cid:38) mean change in [Mg/Fe] from 0 −
15 kpc is ≈ − .
111 dexand the host-to-host standard deviation is ≈ .
02 dex. Thesedifferences are likely attributable to our stellar wind modelhaving a metallicity dependent yield for O and not for Mg.This leads to more O production at small radii where themetallicity is higher, thus a flatter profile.Fig. 3 shows the radial gradients of [O/H] in our simu-lations at z = 0 and includes observations of the MW, M31, and nearby MW/M31-mass galaxies. We fit the gradients inour simulations using a least-squares fit of the [O/H] abun-dance across 4 −
12 kpc. As Fig. 2 shows, including the innerregion of our disks, where the bulge dominates ( R (cid:46) −
12 kpccovers the inner and outer disk and generally exhibits asingle power-law profile. The solid blue line shows the me-dian ( − .
028 dex kpc − ) across our 11 galaxies, while theshaded regions show the 68th percentile and the full distri-bution. The latter ranges from − .
042 to − .
024 dex kpc − .[Fe/H] gradients show similar results, with the full distribu-tion spanning − .
044 to − .
024 dex kpc − .Fig. 3 shows [O/H] gradients observed in nearbyMW/M31-mass galaxies as box-and-whisker plots, with thebox showing the 68th percentile and the whiskers show-ing the full observed range. We apply a cut on the stel-lar masses of these observed samples to be comparable toour simulations. The Kreckel et al. (2019, K19) sample in-cludes 5 galaxies from the PHANGS-MUSE survey with10 . ≤ log M star /M (cid:12) ≤ .
6, the Zinchenko et al. (2019,Z19) sample includes 7 galaxies from CALIFA DR3 with10 . ≤ log M star /M (cid:12) ≤ .
8, the Belfiore et al. (2017, B17)sample includes 13 galaxies from the MaNGA survey with10 . ≤ log M star /M (cid:12) ≤
11, and the S´anchez-Menguianoet al. (2016, S16) sample includes 20 galaxies from the CAL-IFA survey with 10 . ≤ log M star /M (cid:12) ≤
11. In additionto the mass cut, we select galaxies that have gradients mea-sured across a radial range comparable to the range in ouranalysis (the measured ranges all fall within 2 −
14 kpc ex-cept for Kreckel et al. (2019) which falls within 1 −
11 kpc).While all of these observed samples show almost exclusivelynegative gradients in [O/H], their abundance gradients aretypically flatter than in our simulations. Our simulationsare consistent at the 1- σ level with K19 and B17, and atthe 2- σ level with Z19. However, our sample does not over-lap with S16. Note that the calibrator used for determiningthe abundances varies from survey to survey. Using differentcalibrators can give drastically different abundance measure-ments (Hemler et al. 2020), which could contribute to dis-crepancies between the different surveys, and to differenceswith our simulations. Note that the difference between oursimulations and these observations are comparable to thedifferences between surveys themselves.Fig. 3 also shows observed abundance gradients in M31and the MW from HII regions. The orange points show ob-served gradients in M31 from Zurita & Bresolin (2012, Z12)and Sanders et al. (2012, S12). These gradients are slightlyshallower than in our simulations, though they agree within2- σ . This may be a consequence of M31 gradient measure-ments spanning ≈ −
25 kpc: from our analysis, includingthe outer regions of a gas disk flattens the inferred gradient.The black points show measured gradients of the MW.Moll´a et al. (2019a), shown in grey, is a best-fit measure-ment of the MW abundance gradient based on the combineddata of Rudolph et al. (2006); Balser et al. (2011); Estebanet al. (2013, 2017); Fern´andez-Mart´ın et al. (2017). We showuncertainties for all samples. The red error bars for Balseret al. (2011); Wenger et al. (2019) show the impact of mea-suring the radial gradient along different galactic azimuths.Balser et al. (2011) finds gradients ranging from − .
03 to
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Figure 3.
Radial gradients in gas-phase [O/H] across our 11 galaxies and observed in the MW, M31, and in nearby MW-mass galaxies.The blue horizontal line shows the median across our 11 galaxies, with the dark shaded region showing the 68th percentile and thelight shaded region the full distribution. We also show observations of radial gradients in external galaxies, from Kreckel et al. (2019,K19), Zinchenko et al. (2019, Z19), Belfiore et al. (2017, B17), and S´anchez-Menguiano et al. (2016, S16), via box-and-whisker, wherethe box displays the 68th percentile, the whiskers display the full distribution, and the orange horizontal line is the median. Orangecircles show observed abundance gradients for M31 derived from HII regions by Zurita & Bresolin (2012, Z12) and Sanders et al. (2012,S12). Black circles show observed abundance gradients for the MW derived from HII regions from Rudolph et al. (2006, R06), Balseret al. (2011, B11), Esteban et al. (2013, E13), Esteban et al. (2017, E17), Fern´andez-Mart´ın et al. (2017, F17), Wenger et al. (2019,W19), and Arellano-C´ordova et al. (2020, A20). Moll´a et al. (2019a, M19), shown in grey, is the gradient derived from a compilationof the data from R06, B11, E13, E17, and F17. We show uncertainties for all points. For W19 and B11 the red shows the variation ingradient observed by looking along different azimuths. The dashed line shows the best-fit MW gradient ( − .
046 dex kpc − ), based on thegradients the observations presented here (excluding M19). The median [O/H] gradient across our galaxies is − .
028 dex kpc − and thestandard deviation is 0 .
005 dex kpc − , in agreement with K19 and B17 to within 1 σ , and with Z19 to within 2 σ , but not in agreementwith S16. Our simulations also agree with observations of M31 and some of the MW observations. − .
07 dex kpc − and Wenger et al. (2019) find gradientsranging from − .
035 to − .
075 dex kpc − which highlightsthat measurements of the MW radial gradient are stronglysensitive to azimuthal variations. The different samples in-clude different radial ranges, so they are not exactly com-parable to each other or our analysis. Most measurementsof [O/H] gradients in the MW overlap with our simulations,though our simulations generally have shallower gradients.While not included in Fig. 3, Hernandez et al. (2020)recently measured the radial [O/H] gradient in neutral andionized gas in M83. The gradients were measured out to ≈ . ≈ − .
02 dex kpc − and for ionized gas is ≈− .
03 dex kpc − , in good agreement with our values.As a whole, the radial gradients in our simulations aresomewhat steeper than in external galaxies but somewhatshallower than in the MW. The MW may be an outlier: asBoardman et al. (2020) note, its gradient is typically steeperthan those observed in MW analogs. These differences arelikely the result of a combination of different factors, such as: measuring over different radial ranges or using differ-ent calibrators. For example, B17 also measure the gradientin their MaNGA observations using a different calibratorfor [O/H] (O3N2, not shown here, as opposed to R23, asFig. 3 shows), which results in a median gradient that is ≈ .
008 dex kpc − shallower. Thus, given that S16 used theO3N2 calibrator applied to their CALIFA observations, thismay explain the discrepancy between S16 and B17. We de-fer a more detailed comparison via synthetic observationsof our simulations, tailored to each observation, to futurework. Rather, Fig. 3 provides a broad comparison, highlight-ing that the radial gradients of gas-phase [O/H] within oursimulations lie within the scatter across the MW, M31, andnearby MW-mass galaxies. We next explore the evolution of gas-phase radial gradi-ents of [O/H], [Fe/H], and [N/H] at z ≤ .
5, over the last ∼
10 Gyr, during the primary epoch of disk assembly, to un-derstand the initial conditions for star formation and chemi-cal tagging of stars. In summary, we find that at earlier timesthe gas disk was more radially homogeneous (flatter gradi-ents), so chemical tagging offers less discriminating powerfor radial birth location at earlier times.Similar to Fig. 2, Fig. 4 shows radial profiles of [O/H],
MNRAS000
MNRAS000 , 1–23 (2021)
M. A. Bellardini et al. [Fe/H], and [O/Fe] in gas at different redshifts. The solid lineshows the mean across our 11 galaxies, while the shaded re-gion shows the 1- σ scatter. At all radii, [O/H] and [Fe/H]increase with time, as the gas mass declines while more starsenrich the interstellar medium (ISM). This evolution agreeswith the observed gas-phase galaxy mass-metallicity rela-tion (Tremonti et al. 2004). Ma et al. (2016) explored thisevolution across a wide galaxy mass range in the FIRE-1simulations: they found that as galaxies grow more massive,the mass-loading factor of their winds decreases, and metalsare more easily held in/near the galaxy as opposed to be-ing driven into the halo (see also Muratov et al. 2015, 2017;Angl´es-Alc´azar et al. 2017).Fig. 4 also shows that both [O/H] and [Fe/H] have nega-tive radial gradients at all times. Ma et al. (2017) also foundprimarily negative gradients in the FIRE-1 suite, becausethe high star-formation efficiency in the inner disks of galax-ies with well ordered rotation leads to sustained negative ra-dial gradients. At z = 0, our average change in [O/H] from0 −
15 kpc is ≈ .
51 dex, while this is ≈ .
56 dex for [Fe/H].At z = 1 . t lookback = 9 . −
15 kpc is ≈ .
24 dex for [O/H] and 0 .
28 dexfor [Fe/H]. Furthermore, as expected given the scatter information history, we find larger host-to-host scatter (aver-aged over all radii) at earlier times: 0 .
09 dex for [O/H] and0 .
07 dex for [Fe/H] at z = 0, but at z = 1 . . z (cid:46) t lookback (cid:46) . ≈ α -elements like O. At later times, the (moredelayed) Ia supernovae preferentially enrich the galaxy inFe, driving down [O/Fe]. However, the typical change ingas-phase [O/Fe] at fixed radius from z = 1 . t lookback =9 . ≈ . − .
03 dex. The [O/Fe] radialgradients are positive at all times, because the outer diskis always younger than the inner disk/bulge, though thegradients are weak at larger radii. The [O/Fe] radial pro-files steepen at small radii at later times, at least at z (cid:46) t lookback (cid:46) . .
015 dex at z = 1 . .
037 dexat z = 0. Overall, [O/Fe] does not provide strong discrimi-nation power for chemical tagging at any radii or time thatwe examine.While not shown here, we also measure the evolution Figure 4.
Radial profiles of gas-phase elemental abundances inour simulations at different redshifts. The solid lines show themean and the shaded regions show the 1- σ scatter across our 11galaxies. The top and middle panels show steepening radial pro-files of [O/H] and [Fe/H] with time, as the gas disk becomes morerotation dominated, with less radial turbulence, thus sustainingstronger radial gradients. The bottom panel shows that the in-nermost regions of the gas disk have lower [O/Fe] than the outerdisk, this indicates that the inner disk is more evolved, that is, hashad more Ia supernovae that produce more Fe than α -elements(like O), than the outer disk. of [Mg/Fe], which is more significant than [O/Fe]. In theouter disk ( R = 12 kpc) [Mg/Fe] decreases from ≈ . z = 1 . ≈ .
22 dex at z = 0. In the inner disk ( R = 4 kpc)the evolution is larger, from ≈ .
29 dex to 0 .
18 dex over thesame redshifts. The stronger evolution seen in [Mg/Fe] likelyresults from stellar winds in our simulations producing rel-atively little Mg. The stellar wind model used in the sim-
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Figure 5.
Evolution of radial gradients in gas elemental abun-dances. The lines show the means for [O/H], [Fe/H], and [N/H].The dark shaded region shows the 1- σ scatter and the light shadedregion shows the total scatter in [O/H] across our 11 galaxies.At each redshift we calculate the gradient via a linear fit acrossa redshift-dependent radial range: 4 < R <
12 kpc for z < t lookback < . < R < z ≥
1. The gradientis flattest ( ≈ − .
015 dex kpc − ) at z = 1 . t lookback = 9 . z ≈
1. At lower redshifts, the gas diskbecomes more rotationally supported and is capable of sustainingstronger radial gradients, becoming steepest ( ≈ − .
03 dex kpc − )at z = 0. The gradients of [O/H] and [Fe/H] are almost identi-cal, despite being sourced primarily by core-collapse and Ia su-pernovae, respectively. However, the [N/H] gradient, sourced pri-marily through stellar winds from massive stars, which have ametallicity dependent mass loss rate, is steeper at all times. ulations has a metallicity dependent yield for O. At lowerredshifts, as the gas disk becomes more enriched, O produc-tion also increases, leading to less [O/Fe] evolution.Fig. 5 shows the evolution of the mean radial gradientsof [O/H], [Fe/H], and [N/H]at redshifts 0, 0 .
2, 0 .
5, 0 .
75, 1,1 .
25, and 1 .
5. The shaded regions show the 1- σ scatter. Wealso checked that the evolution of each galaxy qualitativelyagrees with this mean evolution. We measure radial gradi-ents between 4 < R <
12 kpc at z < t lookback < . < R < z ≥ t lookback ≥ . (cid:46) − R star90 , which is ≈
11 kpc at z = 0and (cid:46) z (cid:38) .
5. At z <
1, we measure the gradient at4 −
12 kpc, because the inner disk is dominated by the bulgeregion and has a steeper profile, as Fig. 4 shows. Thus, atlate times we measure at 4 −
12 kpc, which encompasses boththe inner and outer disk and exhibits a nearly linear profile.We tested our analysis measuring the gradient over varying radial ranges and found that, while the normalization variessomewhat, the shape of the profile and the evolution areconsistent, as Fig. 4.As Fig. 5 shows, the strength of the radial gradientgenerally decreases over time. The minimum magnitude ofthe gradient is ≈ .
01 dex kpc − and occurs at z ≈ . t lookback ≈ . z (cid:46) t lookback (cid:46) . − .
03 dex kpc − at z = 0. The gradients prior to z ≈ z = 0 and is largest at z = 0 .
75 ( t lookback = 6 . . ≈ .
015 dex kpc − )than O or Fe at all times. We next examine the vertical profiles (in absolute height)of elemental abundances, for all gas near the solar circle,within a cylindrical radius of 7 < R < σ scatter; we only show the scatter at z = 1 . t lookback = 9 . z = 0 for clarity.Fig. 6 shows that any systematic trends in abundancewith height to 1 kpc is (cid:46) .
01 dex kpc − absolute on averageat all times, and the 1- σ scatter is (cid:46) .
01 dex kpc − at z = 0and (cid:46) .
02 dex kpc − at z = 1 .
5. Thus, the gas disk is wellmixed vertically. In most of our galaxies, the deviations inabundance increase with distance from the midplane, thatis, they shows a systematic gradient with height. Over time,these vertical gradients become increasingly (if weakly) morenegative, which supports the idea of ‘upside-down’ disk for-mation (e.g. Bird et al. 2013; Ma et al. 2017; Bird et al.2020), such that stars formed in a more vertically extendeddisk at higher redshifts, leading to more enrichment at largerheights, at later times stars formed in a thinner disk and gasfarther out of the midplane became relatively less enriched.At z (cid:38)
1, the absolute strength of this vertical gradient isin fact comparable to the radial gradient (Fig. 5), while at z = 0, the vertical gradient is ≈ × weaker than the radialgradient. This is because the timescale for vertical mixing MNRAS000
1, the absolute strength of this vertical gradient isin fact comparable to the radial gradient (Fig. 5), while at z = 0, the vertical gradient is ≈ × weaker than the radialgradient. This is because the timescale for vertical mixing MNRAS000 , 1–23 (2021) M. A. Bellardini et al.
Figure 6.
The vertical change in [O/H] relative to the midplaneat various redshifts. We measure the vertical profile at R = 7 − σ standard deviation across 11 galaxies. [O/H] showsminimal variation with height: the mean deviation at 1 kpc abovethe midplane is ≈ .
008 dex at z = 1 . t lookback = 9 . ≈ − .
006 dex at z = 0. The majority of stars form within a fewhundred parsecs of the galactic midplane at low z and within ∼ . | Z | ≤
200 pc) the variation in [O/H]is even smaller at all times ( | ∆[O / H] | (cid:46) .
001 dex at z = 1 . | ∆[O / H] | (cid:46) .
002 dex at z = 0). At z < t lookback < . is short, given gas turbulence, and that the vertical scale-height of the gas is itself set by the maximum Jeans lengthat that time. Furthermore, with implications for chemicaltagging, the majority of star formation in our simulations islimited to (cid:46)
500 pc of the midplane at z < .
5, and (cid:46) . z < .
5, and Fig. 6 shows that vertical variations inabundance are minimal on those scales. The vertical trendsin [Fe/H] (not shown here) are consistent with [O/H] within ≈ .
01 dex.
We next investigate azimuthal variations of elemental abun-dances in gas, including its evolution. We thus test a com-mon assumption in galactic evolution, that gas is well mixedazimuthally within a given annulus (e.g. Frankel et al. 2018;Frankel et al. 2020).Fig. 7 shows the standard deviation in [O/H] and [O/Fe]along angular bins of varying length at fixed radius. Specif-ically, we compute the standard deviation within a givenangular bin, and Fig. 7 shows the mean standard deviationacross all bins of a given angular size for all 11 simulations.We stack 3 snapshots (∆ t ≈
50 Myr) for each redshift. Weuse an annulus of gas ± . ± .
15 kpc of the selected cylindrical radius, to ensurethat the angular length dominates over the radial length forour smallest angular bins, that is, to ensure that the radial gradient does not induce significant scatter. We show the1- σ scatter for z = 0 and z = 1 ( t lookback = 7 . z = 1 . φ ≤ ◦ , becauseits angular bins contain too few gas particles.At z = 0 and R = 8 kpc (near the solar circle) the typ-ical azimuthal scatter across the gas disk is (cid:46) .
053 dex for[O/H], (cid:46) .
055 dex for [Fe/H] (not shown), and (cid:46) .
01 dexfor [O/Fe]. This value for [O/H] agrees well with MW ob-servations (Wenger et al. 2019) and observations of externalgalaxies (Sakhibov et al. 2018; Kreckel et al. 2019; Kreckelet al. 2020), though we emphasize that we are not measuringazimuthal scatter in the same way: those observations typ-ically measure differences in abundances between arm andinter-arm regions or measure abundance variations betweenHII regions within an aperture of a given size.Our azimuthal scatter decreases with smaller angu-lar bin size, with a minimum of ≈ .
045 dex for [O/H]( ≈ .
046 dex for [Fe/H]) and ≈ .
009 dex for [O/Fe]at thesmallest angular scales. Interestingly, this minimal scatterremains well above 0 dex as ∆ φ goes to 0. We empha-size that our analysis does not zoom-in on giant molecularclouds (GMC) or individual star-forming regions, but ratherwe examine all of the ISM centered on (effectively) randompositions. Thus, our results on small scales do not immedi-ately inform the homogeneity of individual GMCs, especiallygiven their short lifetimes ( (cid:46) σ host-to-host scatter is approximately indepen-dent of bin size and is (cid:46) .
014 dex for [O/H], (cid:46) .
015 dexfor [Fe/H], and (cid:46) .
005 dex for [O/Fe]. Thus, at z = 0 gaswithin all of our galaxies is well mixed, that is, the azimuthalscatter is comparable to typical measurement uncertainties( ∼ .
05 dex) for elemental abundances.Fig. 7 shows that, at all radii and at all angular bin sizes,the azimuthal scatter was more significant at earlier times,that is, gas was less azimuthally mixed than it is now. Thisis likely because higher accretion and star formation ratescombined with burstier star formation leads to more pro-nounced local pockets of enrichment in gas. At R = 8 kpcand at z = 1 . t lookback = 9 . (cid:46) . (cid:46) .
05 dex for [O/Fe]. The scatter does not drop below ≈ .
15 dex for either [O/H] or [Fe/H] at the smallest az-imuthal scales (0 .
035 dex for [O/Fe]). The 1- σ host-to-hostscatter is (cid:46) .
05 dex for [O/H] and [Fe/H] and (cid:46) .
016 dexfor [O/Fe].Additionally, Fig. 7 shows that the difference in scat-ter between large and small angular scales varies with time.This difference is more significant at earlier times: at z (cid:38) t lookback (cid:38) . ≈ .
042 dex at R = 8 kpc.Thus, at early times, galaxy-scale fluctuations are more im-portant in driving azimuthal scatter (as visible in Fig. 1).However, at z ∼
0, the azimuthal scatter across small versuslarge angular scales differs by only ≈ .
008 dex, so small-scale fluctuations drive most of the azimuthal scatter (alsovisible in Fig. 1). These results at low redshift are usefulfrom an observational perspective, because they means thatone can generalize smaller-scale observations of gas-phaseabundances to overall azimuthal trends at fixed radius.
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Figure 7.
Azimuthal scatter in elemental abundances for [O/H] and [O/Fe] in gas, as a function of angular scale, at different redshiftsand different radii. The solid lines show the mean and the shaded regions show the 1- σ scatter across our suite of 11 galaxies: we showscatter only at z = 1 ( t lookback = 7 . z = 0, near the solar circle ( R = 8 kpc), the average azimuthal scatter across the disk is ≈ .
053 dex for [O/H] and [Fe/H] (not shown)and ≈ .
009 dex for [O/Fe]. For all angular bin sizes, the average scatter increases with redshift: at earlier cosmic times the gas diskswere less well mixed within a given annulus. At z = 0 the scatter across the disk is ≈ . ≈ .
05 dex for [O/Fe]. The scatter also increases with angular bin size at all redshifts, although the increase is minimal at late times.At low z , this means that azimuthal variations are dominated by local (and not global) fluctuations in the disk. Finally, the azimuthalscatter increases with radius for individual abundances: gas is azimuthally better mixed in the inner disk, likely a result of shorter orbitaltimes leading to faster mixing. Fig. 7 also shows that the azimuthal scatter dependson radius. The azimuthal scatter increases with increasingradius for a given angular bin size, and in fact, this is truefor both fixed angular and physical bin size. At z = 0 theazimuthal scatter across the entire disk at R = 4 kpc is (cid:46) .
042 dex for [O/H] ( (cid:46) .
046 dex for [Fe/H], not shown),and this increases to (cid:46) .
062 dex for [O/H] and [Fe/H] at R = 12 kpc. At z = 1 . t lookback = 9 . (cid:46) .
015 dex to 0 .
25 dex for [O/H] and[Fe/H]. We do not find radial dependence in [O/Fe], whichhas a maximal scatter (cid:46) .
01 dex ( (cid:46) .
052 dex) at all radiiat z = 0 ( z = 1 . R = 8 kpc at 3 redshifts. As we discussedabove, metallicity-dependent stellar winds from massivestars, rather than supernovae, primarily source N, so thiscompares the azimuthal mixing of elements sourced by these different processes. The azimuthal scatter of [N/H] is largerthan that of [O/H] and [Fe/H] for all bins at each redshift.On the scale of the entire annulus, the scatter in [N/H] is ap-proximately 0 .
015 dex larger at z = 1 ( t lookback = 7 . .
01 dex larger at z = 0. This discrep-ancy is slightly smaller for smaller angular scales at z = 1(approximately 0 .
01 dex difference), but the difference in az-imuthal scatter is independent of scale at z = 0. As previ-ously stated, our stellar-wind rate depends linearly on pro-genitor metallicity, which likely drives these differences forN as compared with O and Fe. One might expect Fe, beingsourced primarily by Ia supernovae, to have less scatter atsmall scales than O, because Ia are caused by preferentiallyolder stars than core-collapse supernovae, which occur closerto stellar birth location. A possible cause of our similaritycomes from our assumed Ia delay-time distribution (Man-nucci et al. 2006), which has a significant component fromprompt Ia, at ages (cid:46)
100 Myr. A Ia delay-time distributionwith a more significant contribution from older stellar ages(Maoz & Graur 2017) may lead to less small-scale scatter(e.g. Gandhi et al., in prep.).
MNRAS000
MNRAS000 , 1–23 (2021) M. A. Bellardini et al.
Figure 8.
Azimuthal scatter, as in Fig. 7, at R = 8 kpc for [O/H](solid), [Fe/H] (dashed), and [N/H] (dash dot). As with the ra-dial gradients, the azimuthal variations of [O/H] and [Fe/H] arealmost identical, despite being sourced primarily by core-collapseand Ia supernovae, respectively, while the variation in [N/H],which is sourced primarily through metallicity-dependent stellarwinds from massive stars, is higher at all redshifts. We now compare the relative importance of radial gradientsversus azimuthal scatter in gas-phase abundances. Fig. 9shows the evolution of the radial variations in [O/H], frommultiplying the radial gradient (as calculated in § −
12 kpc at z < t lookback < . − z ≥
1, where the radial profile is approximately linear (seeFig. 4). Fig. 9 also shows evolution of the azimuthal scatteraround each disk (360 ◦ ) at 3 radii.Fig. 9 shows that at early times, z (cid:38) t lookback (cid:38) . ◦ angular scatter in [O/H] at all radiiis larger than the radial variation, with the angular scat-ter at large radii being approximately 2 × higher. How-ever, the radial variation dominates at all radii at z (cid:46) . t lookback (cid:46) . z = 0 the radial variations in [O/H]is approximately a factor of 4 higher than the azimuthalvariation. Thus, at z ∼
0, one can approximate the gas diskvariations primarily via the radial gradient, but at earliertimes, the azimuthal variations dominate.Fig. 10 shows the redshift (and lookback time) whenthe radial variation begins to dominate over the azimuthalscatter in [O/H] (the intersection of the solid and dashedlines in Fig. 9), for 3 radii. We measure this transition red-shift separately for each simulation: the horizontal lines showthe median, and the boxes and whiskers show the 68th per-centile and the full distribution. For hosts with no transi-tion redshift (2 at R = 4 kpc and 1 at R = 8 kpc), we as-sign an arbitrarily high transition redshift, which extendsabove the figure. The radial variation starts to dominateearlier at smaller radii, because the azimuthal scatter in-creases with radius at all redshifts (Fig. 7). The mediantransition redshift, after which the radial variation dom-inates, is z ≈ . t lookback ≈ . R = 4 kpc, z ≈ . t lookback ≈ . R = 8 kpc, and z ≈ . t lookback ≈ . R = 12 kpc.While we in general find no systematic differences be-tween our LG-like hosts and our isolated hosts, the tran-sition redshift of the LG-like simulations is systematicallyhigher than that of the isolated simulations at large radii.At R = 12 kpc the the median transition of the LG suiteoccurs at z ≈ . t lookback ≈ . z ≈ . t lookback ≈ . equality , the ratio of the 360 ◦ azimuthal scatter (at a given radius) to the radial gradient.This represents the characteristic radial scale over whichthe radial and azimuthal abundance variations are equal.The lines show the median ∆R equality for [O/H] and [Fe/H],measuring the azimuthal variation at 3 radii, and the shadedregion shows the 68th percentile of [O/H] at R = 8 kpc.The median ∆R equality for [O/H] is (cid:46)
19 kpc at z = 1 . t lookback = 9 . (cid:46) . z = 0 at the so-lar circle. This corresponds to the radial range over whichazimuthal scatter dominates the variations in abundance,rather than the radial gradient. Thus, for the purposes ofchemical tagging, this represents a limit for the radial pre-cision that chemical tagging (of a single element) can placeon the formation location of a star without also modelingazimuthal location. At early times, ∆R equality is compara-ble to or greater than the size of the disk, meaning that theazimuthal coordinate determines the abundance of newlyforming stars more than the radial position.∆R equality is largest at z = 1 . t lookback = 9 . equality on ra-dius. This ratio slightly increases as a function of radius,because the azimuthal scatter increases with radius at alltimes (Fig. 7). This means that modeling chemical taggingof stellar birth radius is more challenging at larger radii.In summary, any models for chemical tagging shouldincorporate azimuthal scatter in abundance especially at z (cid:38) . t lookback (cid:38) . We next explore observational implications of our measuredradial gradients, by comparing them against typical mea-
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Figure 9.
The evolution of variations in [O/H] in gas, both radiallyand azimuthally. (We do not show [Fe/H], its trends are consis-tent with [O/H] to (cid:46) .
01 dex.) The solid lines show the meanscatter, across our 11 galaxies, for the full (360 ◦ ) annulus of gasat R = 4 kpc (blue), R = 8 kpc (yellow), and R = 12 kpc (green).The red dashed line shows the mean radial change in [O/H] acrossa radial distance of 8 kpc. The shaded regions show the 1- σ scat-ter. While the radial gradient dominates the spatial variationsat late cosmic times, azimuthal variations were more significantthan the radial gradient at earlier cosmic times ( z (cid:38) .
9, or t lookback (cid:38) . R = 8 kpc). Larger radii transition to ra-dially dominated abundance variations at later times (see alsoFig. 10). Thus, elemental evolution models should not assume az-imuthal homogeneity of abundances at early times; instead, az-imuthal variations are the primary source of spatial informationfor chemical tagging of stars forming at early times. Figure 10.
Following Fig. 9, the redshift below which radial vari-ations in [O/H] dominate over azimuthal variations at 3 radii inour 11 simulations (the intersection of the dashed and solid lines).The horizontal line shows the median, the box shows the 68th per-centile, and the whiskers show the full distribution. At R = 4 kpctwo hosts have no transition redshift and at R = 8 kpc one hostshas no transition redshift, that is, the radial variation is alwaysstronger than azimuthal scatter for the redshift range we observe.This transition redshift is the last time the azimuthal variation isstronger than the radial variation. This transition occurs earlierat smaller radii, where azimuthal variations are smaller (Fig. 7),given the shorter timescale for mixing at smaller radii. Before thesetransition redshifts, any model of chemical tagging should accountfor azimuthal variations as the primary source of spatial variation. surement uncertainties in elemental abundances, to under-stand the radial scales over which gas is effectively homo-geneous in a measurable sense. Thus, in this sub-section weignore azimuthal variations and focus just on radial gradi-ents. While we examine gas-phase abundances, the chemicaltagging of stars (that form out of this gas) ultimately mo-tivates our work, so we examine measurement uncertaintiestypical of MW stellar surveys.Motivated by observational surveys of stellar abun-dances, we select 3 observational measurement uncertaintiesof δ m = 0 .
01, 0 .
05, and 0 . homogeneous ,the ratio of δ m to the radial gradient in abundance, ver-sus redshift. The solid line shows the median and theshaded region shows the 68th percentile accross our 11 hosts.∆R homogeneous represents the radial scale of measurable ho-mogeneity: below this radial length scale, the change inabundance from the radial gradient is less than this mea-surement uncertainty. The larger ∆R homogeneous is, the lessprecisely measurements can pinpoint a star’s birth radius.For a fiducial measurement uncertainty of δ = 0 .
05 dexat z = 1 . t lookback = 9 . homogeneous ≈ . ≈ . z = 0 for [Fe/H]. [O/H] is con-sistent to within ≈ z = 1 . homogeneous is larger and has large scat-ter. The largest scatter, at z = 1 .
5, comes from the ra-dial gradients being flattest: some galaxies have gradients approaching 0 dex kpc − , as Fig. 5 shows. After z = 0 . t lookback = 6 . homogeneous decreases over time.This means that chemical tagging with measured abun-dances can identify the birth radius of more recently formedstars more precisely.Comparing Fig. 11 with Fig. 12 shows that, in termsof limitations on chemical tagging for a star’s birth radius,at z (cid:38) . t lookback (cid:38) . δ = 0 .
05 dex. In the outer disk( R ≥ δ = 0 .
01 dex, azimuthal variations dominateat all times at all radii. This implies that, if a primary mo-tivation for chemical tagging is inferring the birth locationof a star, there is not much benefit in pushing to higherprecision, because azimuthal variations dominate. In fact,Fig. 9 can indicate the maximum precision in elementalabundance that one should aim to measure stars of a givenage for this purpose, given our predicted azimuthal scatter,unless a given chemical tagging approach includes model-ing azimuthal variations. We will explore possible models infuture work.
MNRAS000
MNRAS000 , 1–23 (2021) M. A. Bellardini et al.
Figure 11. ∆R equality is the ratio of the azimuthal scatter tothe radial gradient: it indicates the characteristic radial scale be-low which azimuthal variations dominate over radial variations,for [O/H] (solid) and [Fe/H] (dashed). We compute the azimuthalscatter using the full (360 ◦ ) annulus of gas at R = 4, 8, and 12 kpc.Lines show the median, and shaded regions show the 68th per-centile across our 11 galaxies. At early times, ∆R equality is compa-rable to the size of the gas disk, when azimuthal scatter dominatedover radial variation, but after z ≈ . t lookback = 5 . z = 0 the median ∆R equality near the solar circle is ≈ . ≈ . z = 1 . t lookback = 9 . ≈ . ≈ . equality , the primary source of inhomogeneity of elementalabundances in gas is azimuthal variations, which complicate/limitthe use of elemental abundances to infer a star’s birth location toa radial precision smaller than ∆R equality . Finally, we explore the full distributions of elemental abun-dances in gas that our simulations predict. Again, we em-phasize that our FIRE-2 simulations explicitly model thesub-grid diffusion/mixing of elemental abundances in gasvia unresolved turbulent eddies, which leads to significantlymore realistic distributions (Su et al. 2017; Escala et al. 2018;Hopkins et al. 2018).We measure [O/H] and [Fe/H] distribution at R =4 kpc, R = 8 kpc, and R = 12 kpc at z = 1 ( t lookback =7 . z = 0 . t lookback = 5 . z = 0 for allgalaxies. We fit these with a skew normal distribution, us-ing the LevMarLSQ fitter in Astropy (Astropy Collaborationet al. 2013; Price-Whelan et al. 2018): dFdx = A × exp (cid:18) − . (cid:16) x − µσ (cid:17) (cid:19) × erf (cid:16) α × x − µ √ σ (cid:17) µ is the mean, σ is the standard deviation, and α is theskewness. Fig. 13 shows representative example distributionsof [Fe/H], for good and bad fits to this distribution, for asingle simulation, and we list the percent of galaxies andradii that fall into each category.As the left panel of Fig. 13 shows, a skew normal dis-tribution reasonably fits these distributions in most cases at z ∼
0. However, there are several common failure modes. We
Figure 12.
Assuming a given measurement uncertainty in an el-emental abundance, ∆R homogeneous is the ratio of this measure-ment uncertainty to the radial gradient in that abundance; itis the characteristic radial scale below which the gas disk is ef-fectively radially homogeneous in a measurable sense. We show∆R homogeneous versus redshift, for 3 measurement uncertainties,motivated by observational surveys of stellar abundances. Thelines show the median across our 11 galaxies for [Fe/H] (solid) and[O/H] (dashed), and the shaded regions show the 68th percentilesfor [Fe/H]. At z = 1 . t lookback = 9 . δ = 0 .
05 dex) across significant radialscales, ∆ r (cid:46) (cid:46) . z = 0.For high-precision abundance measurements ( δ = 0 .
01 dex), thegas disk is measurably homogeneous at ∆ r (cid:46) . z = 1 . r (cid:46) . z = 0. This highlights the limitations fromjust measurement uncertainties on chemical tagging to infer astar’s birth radius; it does not include the additional complica-tions from azimuthal variations (Fig. 11), which can be more im-portant. categorize them as: failing to capture the positive or nega-tive tails of the distribution, failing to capture the width ofthe distribution, or the distribution being multi-modal. Theright panel of Fig. 13 shows examples of each of these fail-ures, along with the percentage of fits ([Fe/H] and [O/H]combined) that we identify to fall into each category ateach redshift, stacking all galaxies and radii at that red-shift. In general, the fit to [O/H] and [Fe/H] at a givenredshift and radius falls into the same category. At z = 0,the vast majority ( ≈ z = 1 ( t lookback = 7 . ≈
11% are well fit. Most common is having a neg-ative underfit or both tails underfit, likely driven by morerapid accretion of low-metallicity gas at earlier times.While not perfect, especially at earlier times, a skewnormal fit does provide a simple characterization of thefull distribution. Thus, Fig. 14 shows the fit parameters for[Fe/H] ([O/H], not shown, is consistent with this) at dif-ferent radii and redshifts. The box-and-whisker plots showthe median, 68th percentile, and the full distribution. Thetop row shows the fitted standard deviation, while the bot-
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Figure 13.
Example distributions of [Fe/H] from our 11 galaxies. Each panel shows the elemental distribution for a single simulation ata single radius. For each distribution, we stack 3 consecutive snapshots (∆ t ≈
50 Myr) and measure all gas within ± . R = 4, 8, and 12 kpc) and within a height ± z = 0, but it provides a worse fit athigher redshift, when the (negative) tails are preferentially underfit. In general, the fit to the distribution of [O/H] and [Fe/H] for thesame R and redshift fall into the same category. tom row shows the fitted skewness, and left to right showsincreasing redshift.At all radii, σ decreases over time. At R = 8 kpc, nearthe solar circle, the median standard deviation decreasesfrom ≈ .
12 dex at z = 1 ( t lookback = 7 . ≈ . z = 0 for [Fe/H] ( ≈ .
11 to ≈ .
06 for [O/H]). Thus, con-sistent with the results for azimuthal scatter, gas at a givenradius becomes more homogeneous over time. Also consis-tent with our results for azimuthal scatter, Fig. 14 showsthat at all redshifts, σ increases with radius, that is, metalsare less well mixed at larger radii. At z = 0 the median is σ ≈ .
05 at R = 4 kpc and increases to ≈ .
07 at R = 12 kpcfor [Fe/H] ( ≈ .
04 to ≈ .
06 for [O/H]).Fig. 14 (bottom row) shows that the distributions arepreferentially negatively skewed at earlier times, but theytrend toward Gaussian over time. At R = 8 kpc the medianskewness is α ≈ − .
54 at z = 1 ( t lookback = 7 . α ≈ − .
21 at z = 0 for [Fe/H] ( α ≈ − .
52 to ≈ − . ≈ − .
19 at R = 4 kpc to ≈ − .
58 at z = 1( ≈ − .
15 to ≈ − . z ∼
0, all radii show abundancedistributions consistent with no skewness at the 1- σ level. We use a suite of FIRE-2 cosmological zoom-in simulationsof 11 MW/M31-mass galaxies to explore the 3-D spatial vari-ations and evolution of elemental abundances [O/H], [Fe/H],and [N/H] in gas at z ≤ . t lookback ≤ . z = 1 .
5, the last ∼
10 Gyr mark the primary epoch ofdisk assembly, which is where we are primarily interested inchemically tagging stars. Our main results are: • Vertical gradients : are negligible. Abundances in gas arewell mixed vertically at all times. At R = 8 kpc, the meandeviation in [O/H] at 1 kpc from the galactic midplane is < .
01 dex at all times. The inner ∼
200 pc of the disks, wherethe majority of star formation for z < . | ∆[ O/H ] | (cid:46) .
002 dex) at alltimes. The inner ∼ . z > . | ∆[ O/H ] | (cid:46) .
01 dex at all times.Thus there is minimal vertical information for chemical tag-ging. • Radial gradients : are negative at all times andfor all abundances, with a maximum steepness of ≈− .
03 dex kpc − at z = 0 and a minimum of ≈− .
01 dex kpc − at z (cid:38) t lookback (cid:38) . MNRAS000
01 dex kpc − at z (cid:38) t lookback (cid:38) . MNRAS000 , 1–23 (2021) M. A. Bellardini et al.
Figure 14. Top : standard deviation of the fitted elemental distribution of gas [Fe/H] at 3 radii, for increasing redshift (left to right).For each distribution, we stack 3 consecutive snapshots (∆ t ≈
50 Myr) and measure all gas within ± . ± Bottom : skewness of the fitted elemental distributionsof gas [Fe/H]. At earlier times, the gas disk had stronger negative skewness, but the disk relaxes to near zero skewness at z = 0. Theskewness shows a slight radial dependence at both z = 1 ( t lookback = 7 . z = 0. At z = 1 the distributions at larger radii weremore negatively skewed, whereas at z = 0 the distributions at smaller radii are more negatively skewed. more rotationally supported and are better able to sustain agradient against radial mixing, as noted in analysis of FIRE-1 simulations in Ma et al. (2017). [N/H] has a steeper gra-dient at all times, because their production is dominated bystellar winds, whose mass-loss rates increase with metallic-ity in our simulations, enhancing the discrepancies betweenmetal-rich and metal-poor regions. [O/Fe] shows little vari-ation with redshift, and is approximately flat across the diskindicating it provides limited discriminating power for chem-ically tagging birth radii. Our [O/H] gradients broadly agreewith most observations of nearby MW-mass galaxies, includ-ing M31, at the 1 or 2- σ level, though our gradients aresomewhat steeper on average. By contrast, our gradientsare somewhat shallower than most observations of the MW,though they agree at the 1 or 2- σ level with 9 of 13 MWobservations. • Azimuthal scatter : systematically decreases over timefor all abundances, from ≈ . z = 1 . t lookback =9 . ≈ .
05 dex at z = 0 for [O/H] and [Fe/H] aroundthe entire disk at R = 8 kpc. This evolution is a result ofhigher gas accretion and also star-formation rates at earliertimes, which lead to stronger variations in abundances onsmall scales, especially at larger radii, where orbital/mixingtimescales are longer. The azimuthal scatter in [N/H] islarger (by ≈ .
01 dex at R = 8 kpc) at all times than in[Fe/H] or [O/H], for the same reasons as above. Azimuthalvariations reduce somewhat with smaller azimuthal aper- ture. However, even in angular bins as small as ≈
350 pc,they remain ≈ .
04 dex at z = 0 and ≈ . z = 1( t lookback = 7 . z = 0( ≈ .
05 dex) and observations of nearby galaxies (Sakhibovet al. 2018; Kreckel et al. 2019; Kreckel et al. 2020). • Azimuthal versus radial scatter : At early times, the az-imuthal scatter was larger than the radial variation for allabundances. We quantify the redshifts when the radial varia-tion (across ∆ R = 8 kpc) first dominates over the azimuthalscatter, finding a median of z ≈ . t lookback ≈ . R = 4 kpc and z ≈ . t lookback ≈ . R = 12 kpc.Before this time, stars born at the same radius could havethe same difference in metallicity as stars born ∆ R (cid:38) R equality . At z ∼
0, ∆ R equality is small at ≈ . z (cid:38) t lookback (cid:38) . R equality is larger than the sizeof the disk. These results indicate that azimuthal variationsin abundances provide the dominant information contentfor chemical tagging for stars formed (cid:38) • Measurable homogeneity : We quantified the radial
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations scales across which our gas disks are effectively homoge-neous in a measurable sense, given representative measure-ment uncertainties. For an uncertainty in elemental abun-dance of 0 .
05 dex, our gas disks are measurably homoge-neous across ∆ R ≈ . z = 0 and ∆ R ≈ . z (cid:38) .
75 ( t lookback (cid:38) . R (cid:38) .
05 dex at all times.Thus, for any measurement uncertainty at or below this,using chemical tagging to measure birth radius is limitednot by measurement uncertainty but instead by azimuthalvariations. These results inform the needed precision for ob-servations, given targeted precision for chemical tagging ofstars across age/time. For example, if one only cares aboutmodeling birth radius, there is little-to-no benefit in mea-suring a stellar abundance to better than ≈ .
05 dex. • Elemental abundance distributions : We measured thefull distributions of elemental abundances in radial annuliand fit them with skew normal distributions. The skewnormal distributions fit these distributions reasonably well,but there are failure modes that become more common athigher redshift, most notably underfitting the negative tailsof the distribution and simultaneously underfitting the pos-itive and negative tails. We find typically negatively skewednormal distributions at z (cid:38) t lookback (cid:38) . z = 0. The primary goal of this paper is to understand the homo-geneity of gas as a proxy for the birth conditions of starsacross space and time, as a first step to understanding theefficacy of chemical tagging in a cosmological context. Thereare caveats to our analysis, though. Namely, our analysis isperformed looking at individual elements, with the excep-tion of [O/Fe]. Examining multi-element abundance distri-butions may well offer more discerning power. Our analysisof [O/Fe] suggests that this may be limited. Furthermore,uncertainty in our fiducial diffusion coefficient leads to un-certainty in our small scale azimuthal abundance scatter, asseen in Appendix. C.Additional complications arise when comparing oursimulations to observational data. We present results in thecontext of constraining chemical tagging in the MW, butour simulations are not exact MW-analogs. Also, when com-paring the redshift evolution of our results to observationsof external galaxies, we track the evolutionary history ofindividual galaxies across time, as opposed to measuringproperties of different galaxies at fixed mass across time.For all of our comparisons to observations, we explore allgas whereas observers typically measure abundances in HIIregions specifically. However, Hernandez et al. (2020) com-pared observations of ionized and neutral gas-phase abun-dance gradients in M83, finding gradients for neutral gasto be − .
17 dex kpc − and gradients for ionized gas to be − .
03 dex kpc − . This might imply that our measured gra-dients are much flatter than one would expect, given obser-vations. Hernandez et al. (2020) did the same analysis ex-cluding the nuclear region of M83 and found the neutral gasto be in much better agreement with ionized gas (a gradientof − .
02 dex kpc − . Additionally, we compare to observa- tions which have measurements in broadly similar physicalregions to those we analyze in the simulations, but they arenot exactly the same.We compared against observations of radial gradients inthe MW, M31, and similar-mass galaxies at z = 0, findingbroad agreement. We also connect our evolutionary trendswith high-redshift observations of gas-phase abundances. Inparticular, we find that our MW/M31-mass galaxies all havenegative radial gradients at z ∼ z (cid:38) M ≈ . × M (cid:12) ). This trend agrees wellwith many observations of comparable mass galaxies at thesehigher redshifts (e.g. Queyrel et al. 2012; Stott et al. 2014;Wuyts et al. 2016; Patr´ıcio et al. 2019; Curti et al. 2020).However, some observational works have found strong neg-ative radial gradients at these masses and redshifts (Wuytset al. 2016; Carton et al. 2018; Wang et al. 2019a). Fur-thermore, while less common than negative radial gradi-ents, some observations find some positive radial gradientsat these redshifts as well (Queyrel et al. 2012; Wuyts et al.2016; Carton et al. 2018; Wang et al. 2019a) which we donot find in any of our galaxies. In general, we find that thesteepening of radial gradients with time in our simulations isconsistent with observational results and follows the trendsin other theoretical analyses (e.g. Ma et al. 2017).One of the most important aspects of our analysisis quantifying azimuthal variations in gas abundances andcomparing their strength relative to radial gradients acrosscosmic time. With the advent of integral field spectroscopy,observations have begun characterizing 2-D abundance dis-tributions in nearby galaxies. These works (S´anchez et al.2015; Vogt et al. 2017; Ho et al. 2017, 2018) all find non-trivial azimuthal variations in nearby galaxies, for example,Kreckel et al. (2019) found variations of ≈ .
05 dex at fixedradius, which agrees well with our results. However, someobservations (e.g. Zinchenko et al. 2016) found no evidencefor large-scale azimuthal variation in nearby galaxies. Oneof our key results/predictions is the evolution of azimuthalvariations, which we predict were stronger at higher red-shifts. Observations of gravitationally lensed systems nowallow sub-kpc measurements at high redshift (Jones et al.2013, 2015), making it possible to test this predicted evolu-tion in more detail.Kreckel et al. (2020) examined azimuthal variationsin gas-phase [O/H] across eight nearby galaxies usingPHANGS-MUSE optical integral field spectroscopy. Whileour technique for measuring azimuthal variations are not ex-actly comparable to their methods, we find similar results.In our analysis we focus on scatter in all gas by measuring amean scatter in angular bins of varying size, so we in effectmeasure the azimuthal inhomogeneities of random patchesof gas at a given radius. In contrast to this, Kreckel et al.(2020) measure abundances specifically in HII regions anddetermine scatter by first subtracting off the radial gradientand then centering apertures of various sizes on individualHII regions and measuring the scatter between the HII re-gions contained within the aperture. They find a slight scaledependence associated with the scatter, which we also see at z = 0, with the scatter on scales larger than ≈ ≈ .
05 dex. The small-scale scatter in Kreckel et al. (2020)( ≈ .
02 dex) is slightly smaller than the z = 0 scatter weobserve, but this could be attributed to the discrepancy in MNRAS , 1–23 (2021) M. A. Bellardini et al. our methods. HII regions are likely better mixed in abun-dances than random patches of gas, so our analysis may beartificially inflating the typical azimuthal scatter of the gasfrom which stars are forming. However, centering on HII re-gions is beyond the scope of our analysis, and in future workwe will examine azimuthal variations in newly formed stars,which may be closer to the values in HII regions.Krumholz & Ting (2018) derived the expected correla-tion function of metal distribution in galaxies across spaceand time using a stochastic diffusion model. While we didnot explore the correlation function of metals, we did exam-ine homogeneity as a function of azimuthal scale, which wecan compare broadly with their work. They found that gas-phase abundances produced primarily through core-collapsesupernovae, in MW-like conditions near the solar circle, arecorrelated on scales of ≈ . − . − . ≈ z = 0 fall between the gradientsHemler et al. (2020) measured in the TNG50 simula-tions ( ≈ − .
02 dex kpc − ) and the gradients Gibson et al.(2013) measured in the MaGICC and MUGS simulations( ≈ − .
04 dex kpc − ). In particular, Hemler et al. (2020)found a gradual flattening of the gradients with time, whichcould come from an ‘inside-out’ growth of galaxies whereinstar formation, hence elemental enrichment, proceeds fromthe inner galaxy to the outer galaxy (e.g. Prantzos & Boissier2000; Bird et al. 2013). The flat(ter) radial gradients at ear-lier times in our galaxies result from higher turbulence andoutflows that frequently eject much of the ISM at thosetimes, perturbations such as mergers and rapid gas infallresult in the velocity dispersion of gas particles dominat-ing over their rotational velocity leading to galaxy-scaleradial mixing (Ma et al. 2017). As the disk settles overtime, it becomes more rotationally supported, so strongerradial gradients can develop/persist. Our results qualita-tively agree with those of the EAGLE simulations (Tisseraet al. 2019), though as with TNG50, Tissera et al. (2019)found [O/H] gradients that are slightly shallower than ours, ≈ − .
011 dex kpc − at z = 0.The evolution of our gas-phase abundance gradients dis-agrees with Agertz et al. (2020), who analyzed the VINTER-GATAN simulation of the m12i initial conditions, performedusing the adaptive mesh refinement (AMR) code RAMSES.They found that the gas-phase profile of [Fe/H] becomesshallower over time (their Fig. 7), compared with our steep-ening with time. One possible explanation is the differencein hydrodynamic solvers: we use the mesh-free finite-mass(MFM) quasi-Lagrangian method in Gizmo , coupled withexplicit modeling of sub-grid mixing, while the AMR simu-lation of VINTERGATAN induces significantly more mixingin gas (complete mixing on the scale of an individual cell),which may contribute to the flattening of their gradient overtime. However, as shown in Appendix. C, the qualitativesteepening of the gradient we observe with time is indepen-dent of the strength of our diffusion coefficient.Galactic evolution models often simplify the abundance distributions of gas in galaxies to a 1-D model (e.g. Minchevet al. 2018; Moll´a et al. 2019a; Frankel et al. 2020), withazimuthal scatter assumed from measurements at z = 0.While this is a useful first step in understanding the abun-dance evolution of galaxies and testing chemical tagging, ourresults mean that this simplifying assumption overestimatesthe radial information content in elemental abundances, in-cluding how well chemical tagging can constrain the birthradius of a star. On the one hand, the non-trivial azimuthalscatter that we find, especially at earlier times, complicatesmodeling the abundances of stars at a given radius. On theother hand, this likely makes individual GMCs more ele-mentally distinct at fixed radius, providing greater discrim-inating power, as we will explore in future work. However,we do not explore the homogeneity of individual GMCs inthis work. Recent work has started to pursue 2-D abundanceevolution models (eg Moll´a et al. 2019b) which may addressthis question, works such as Spitoni et al. (2019) find az-imuthal abundance variations on the order of 0 . z = 0 in our simulation suite.Related to our analysis of a transition epoch, after whichradial variations dominate over azimuthal scatter as the disksettles, Yu et al. in prep. examine the transition epoch from‘bursty’ to ‘steady’ star formation and disk settling in thesame simulations. We checked that their measurement ofthis bursty/steady transition agrees moderately well withthe transition epoch that we present here, at least at smallerradius (4 kpc). We find weaker agreement for our transitiontimes at larger radii (8 and 12 kpc). Furthermore, our tran-sition times are consistently earlier ( ∼ R = 4 kpcand ∼ r = 12 kpc) than those in Yu et al. in prep.,with the transition times on average being more similar atlarger radii. Thus, we find a broad correlation between thetransition from bursty to smooth star formation and thetransition from azimuthal to radial abundance variations,but with significant scatter and some time delay.Our simulations show the importance of consideringazimuthal variations in addition to radial variations whenstudying gas-phase elemental abundance distributions. Thisis particularly important in the context of chemical tagging;in order to accurately identify the birth locations of starsusing elemental abundances the initial conditions of starsneed to be well defined. As we showed in Section. 3.5 az-imuthal variations in abundance are greater than or compa-rable to radial variations at earlier times, so chemical tag-ging models that only account for radial variations will failto accurately capture the scatter in abundances at a givenradius. This could lead to incorrectly assigning stars as co-natal, or vice versa. We also fit the elemental distributionsof our galaxies at different radii at different times, findingthat they shift from negative to zero skewness over time.While skew-normal fits are not a perfect fit for the elementaldistributions of our galaxies at all redshifts, they more ac-curately represent the distributions than a Gaussian. Thus,using our measured distributions would be a useful step inbuilding more accurate abundance evolution models includ-ing for chemical tagging.Next generation telescopes are crucial for testing thepredictions of gas-phase abundance homogeneity presentedin this work, particularly the predictions for azimuthal scat-ter at high redshifts. Current measurements of azimuthalscatter in abundances have been restricted to nearby galax- MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations ies. However, with the advent of JWST NIRSPec IFU andnext-generation adaptive-optic IFUs on telescopes like IRISand TMT, spatially resolved measurements of metallicitiesin distant galaxies are feasible, providing tests of our predic-tions for azimuthal scatter and the transition redshifts whenit becomes sub-dominant.This work is the first step of testing the limits of chem-ical tagging in the FIRE simulations. In the future we willexamine the degree to which these results for gas are mir-rored in newly formed stars across time. Combining those re-sults with measurements of the dynamical evolution of starsin our simulations, we more directly will test the efficacy ofchemical tagging of stars in the FIRE simulations. ACKNOWLEDGEMENTS
We thank Francesco Belfiore for providing abundance gra-dient measurements for galaxies from the MaNGA survey.We also thank Tucker Jones, Ryan Sanders, Joss Bland-Hawthorn, Trey Wenger, and Dana Balser for useful dis-cussions that improved this manuscript.MAB and AW received support from NASA throughATP grants 80NSSC18K1097 and 80NSSC20K0513; HSTgrants GO-14734, AR-15057, AR-15809, and GO-15902 fromSTScI; a Scialog Award from the Heising-Simons Founda-tion; and a Hellman Fellowship. Support for SRL was pro-vided by NASA through Hubble Fellowship grant . DATA AVAILABILITY
Full simulation snapshots at z = 0 are available for m12i,m12f, and m12m at ananke.hub.yt. Some of the python codeused to analyze these data includes the publicly availablepackages https://bitbucket.org/awetzel/gizmo analysis andhttps://bitbucket.org/awetzel/utilities. REFERENCES
Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ,597, 21Acharova I., Gibson B. K., Mishurov Y. N., Kovtyukh V., 2013,Astronomy & Astrophysics, 557, A107Agertz O., et al., 2020, arXiv e-prints, p. arXiv:2006.06008Anders F., et al., 2014, A&A, 564, A115Anders F., et al., 2017, A&A, 600, A70Andrews B. H., Martini P., 2013, ApJ, 765, 140 Angl´es-Alc´azar D., Faucher-Gigu`ere C.-A., Kereˇs D., HopkinsP. F., Quataert E., Murray N., 2017, MNRAS, 470, 4698Arellano-C´ordova K. Z., Esteban C., Garc´ıa-Rojas J., M´endez-Delgado J. E., 2020, Monthly Notices of the Royal Astronom-ical SocietyAsplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A,47, 481Astropy Collaboration et al., 2013, A&A, 558, A33Baba J., Morokuma-Matsui K., Miyamoto Y., Egusa F., KunoN., 2016, MNRAS, 460, 2472Balser D. S., Rood R. T., Bania T. M., Anderson L. D., 2011,ApJ, 738, 27Balser D. S., Wenger T. V., Anderson L. D., Bania T. M., 2015,ApJ, 806, 199Belfiore F., et al., 2017, MNRAS, 469, 151Benincasa S. M., et al., 2020, MNRAS, 497, 3993Bird J. C., Kazantzidis S., Weinberg D. H., Guedes J., CallegariS., Mayer L., Madau P., 2013, ApJ, 773, 43Bird J. C., Loebman S. R., Weinberg D. H., Brooks A.,Quinn T. R., Christensen C. R., 2020, arXiv e-prints, p.arXiv:2005.12948Bland-Hawthorn J., Gerhard O., 2016, ARA&A, 54, 529Bland-Hawthorn J., Krumholz M. R., Freeman K., 2010, ApJ,713, 166Boardman N., et al., 2020, MNRAS, 491, 3672Boeche C., et al., 2013, A&A, 559, A59Boeche C., et al., 2014, A&A, 568, A71Bovy J., 2016, ApJ, 817, 49Brook C. B., Kawata D., Gibson B. K., Freeman K. C., 2004,ApJ, 612, 894Carrell K., Chen Y., Zhao G., 2012, AJ, 144, 185Carton D., et al., 2018, Monthly Notices of the Royal Astronom-ical Society, 478, 4293Cedr´es B., Cepa J., 2002, A&A, 391, 809Cheng J. Y., et al., 2012, ApJ, 746, 149Colbrook M. J., Ma X., Hopkins P. F., Squire J., 2017, MNRAS,467, 2421Cui X.-Q., et al., 2012, Research in Astronomy and Astrophysics,12, 1197Curti M., et al., 2020, Monthly Notices of the Royal AstronomicalSociety, 492, 821Dalton G., et al., 2012, in Ground-based and Airborne Instrumen-tation for Astronomy IV. p. 84460P, doi:10.1117/12.925950Davies B., Origlia L., Kudritzki R.-P., Figer D. F., Rich R. M.,Najarro F., Negueruela I., Clark J. S., 2009, ApJ, 696, 2014De Silva G. M., Freeman K. C., Bland-Hawthorn J., Asplund M.,Bessell M. S., 2007, AJ, 133, 694De Silva G. M., et al., 2015, MNRAS, 449, 2604Escala I., et al., 2018, MNRAS, 474, 2194Esteban C., Carigi L., Copetti M. V. F., Garc´ıa-Rojas J., Mesa-Delgado A., Casta˜neda H. O., P´equignot D., 2013, MNRAS,433, 382Esteban C., Fang X., Garc´ıa-Rojas J., Toribio San Cipriano L.,2017, MNRAS, 471, 987Faucher-Giguere C.-A., Lidz A., Zaldarriaga M., Hernquist L.,2009, The Astrophysical Journal, 703, 1416Fern´andez-Mart´ın A., P´erez-Montero E., V´ılchez J. M., MampasoA., 2017, A&A, 597, A84Frankel N., Rix H.-W., Ting Y.-S., Ness M., Hogg D. W., 2018,ApJ, 865, 96Frankel N., Sanders J., Ting Y.-S., Rix H.-W., 2020, The Astro-physical Journal, 896, 15Freeman K., Bland-Hawthorn J., 2002, Annual Review of Astron-omy and Astrophysics, 40, 487Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Lee K.,2014, MNRAS, 438, 2578Garrison-Kimmel S., et al., 2019, MNRAS, 487, 1380Genovali K., et al., 2014, A&A, 566, A37MNRAS000
Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ,597, 21Acharova I., Gibson B. K., Mishurov Y. N., Kovtyukh V., 2013,Astronomy & Astrophysics, 557, A107Agertz O., et al., 2020, arXiv e-prints, p. arXiv:2006.06008Anders F., et al., 2014, A&A, 564, A115Anders F., et al., 2017, A&A, 600, A70Andrews B. H., Martini P., 2013, ApJ, 765, 140 Angl´es-Alc´azar D., Faucher-Gigu`ere C.-A., Kereˇs D., HopkinsP. F., Quataert E., Murray N., 2017, MNRAS, 470, 4698Arellano-C´ordova K. Z., Esteban C., Garc´ıa-Rojas J., M´endez-Delgado J. E., 2020, Monthly Notices of the Royal Astronom-ical SocietyAsplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A,47, 481Astropy Collaboration et al., 2013, A&A, 558, A33Baba J., Morokuma-Matsui K., Miyamoto Y., Egusa F., KunoN., 2016, MNRAS, 460, 2472Balser D. S., Rood R. T., Bania T. M., Anderson L. D., 2011,ApJ, 738, 27Balser D. S., Wenger T. V., Anderson L. D., Bania T. M., 2015,ApJ, 806, 199Belfiore F., et al., 2017, MNRAS, 469, 151Benincasa S. M., et al., 2020, MNRAS, 497, 3993Bird J. C., Kazantzidis S., Weinberg D. H., Guedes J., CallegariS., Mayer L., Madau P., 2013, ApJ, 773, 43Bird J. C., Loebman S. R., Weinberg D. H., Brooks A.,Quinn T. R., Christensen C. R., 2020, arXiv e-prints, p.arXiv:2005.12948Bland-Hawthorn J., Gerhard O., 2016, ARA&A, 54, 529Bland-Hawthorn J., Krumholz M. R., Freeman K., 2010, ApJ,713, 166Boardman N., et al., 2020, MNRAS, 491, 3672Boeche C., et al., 2013, A&A, 559, A59Boeche C., et al., 2014, A&A, 568, A71Bovy J., 2016, ApJ, 817, 49Brook C. B., Kawata D., Gibson B. K., Freeman K. C., 2004,ApJ, 612, 894Carrell K., Chen Y., Zhao G., 2012, AJ, 144, 185Carton D., et al., 2018, Monthly Notices of the Royal Astronom-ical Society, 478, 4293Cedr´es B., Cepa J., 2002, A&A, 391, 809Cheng J. Y., et al., 2012, ApJ, 746, 149Colbrook M. J., Ma X., Hopkins P. F., Squire J., 2017, MNRAS,467, 2421Cui X.-Q., et al., 2012, Research in Astronomy and Astrophysics,12, 1197Curti M., et al., 2020, Monthly Notices of the Royal AstronomicalSociety, 492, 821Dalton G., et al., 2012, in Ground-based and Airborne Instrumen-tation for Astronomy IV. p. 84460P, doi:10.1117/12.925950Davies B., Origlia L., Kudritzki R.-P., Figer D. F., Rich R. M.,Najarro F., Negueruela I., Clark J. S., 2009, ApJ, 696, 2014De Silva G. M., Freeman K. C., Bland-Hawthorn J., Asplund M.,Bessell M. S., 2007, AJ, 133, 694De Silva G. M., et al., 2015, MNRAS, 449, 2604Escala I., et al., 2018, MNRAS, 474, 2194Esteban C., Carigi L., Copetti M. V. F., Garc´ıa-Rojas J., Mesa-Delgado A., Casta˜neda H. O., P´equignot D., 2013, MNRAS,433, 382Esteban C., Fang X., Garc´ıa-Rojas J., Toribio San Cipriano L.,2017, MNRAS, 471, 987Faucher-Giguere C.-A., Lidz A., Zaldarriaga M., Hernquist L.,2009, The Astrophysical Journal, 703, 1416Fern´andez-Mart´ın A., P´erez-Montero E., V´ılchez J. M., MampasoA., 2017, A&A, 597, A84Frankel N., Rix H.-W., Ting Y.-S., Ness M., Hogg D. W., 2018,ApJ, 865, 96Frankel N., Sanders J., Ting Y.-S., Rix H.-W., 2020, The Astro-physical Journal, 896, 15Freeman K., Bland-Hawthorn J., 2002, Annual Review of Astron-omy and Astrophysics, 40, 487Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Lee K.,2014, MNRAS, 438, 2578Garrison-Kimmel S., et al., 2019, MNRAS, 487, 1380Genovali K., et al., 2014, A&A, 566, A37MNRAS000 , 1–23 (2021) M. A. Bellardini et al.
Gibson B. K., Pilkington K., Brook C. B., Stinson G. S., BailinJ., 2013, A&A, 554, 47Gilmore G., et al., 2012, The Messenger, 147, 25Grand R. J. J., Kawata D., Cropper M., 2015, MNRAS, 447, 4018Hahn O., Abel T., 2011, MNRAS, 415, 2101Hayden M. R., et al., 2014, AJ, 147, 116Hemler Z. S., et al., 2020, arXiv e-prints, p. arXiv:2007.10993Hernandez S., et al., 2020, arXiv e-prints, p. arXiv:2012.12933Ho I.-T., et al., 2017, The Astrophysical Journal, 846, 39Ho I.-T., et al., 2018, A&A, 615Hopkins P. F., 2015, MNRAS, 450, 53Hopkins P. F., et al., 2018, MNRAS, 480, 800Iwamoto K., Brachwitz F., Nomoto K., Kishimoto N., Umeda H.,Hix W. R., Thielemann F.-K., 1999, ApJS, 125, 439Izzard R. G., Tout C. A., Karakas A. I., Pols O. R., 2004, MNRAS,350, 407Jones T., Ellis R., Jullo E., Richard J., 2010, ApJ, 725, L176Jones T., Ellis R. S., Richard J., Jullo E., 2013, The AstrophysicalJournal, 765, 48Jones T., et al., 2015, The Astronomical Journal, 149, 107Kawata D., Hunt J. A. S., Grand R. J. J., Pasetto S., CropperM., 2014, MNRAS, 443, 2757Keller B. W., Wadsley J. W., Wang L., Kruijssen J. M. D., 2019,MNRAS, 482, 2244Kollmeier J. A., et al., 2017, arXiv e-prints, p. arXiv:1711.03234Kreckel K., et al., 2019, The Astrophysical Journal, 887, 80Kreckel K., et al., 2020, MNRAS,Kroupa P., 2001, MNRAS, 322, 231Krumholz M. R., Gnedin N. Y., 2011, ApJ, 729, 36Krumholz M. R., Ting Y.-S., 2018, MNRAS, 475, 2236Leitherer C., et al., 1999, ApJS, 123, 3Lemasle B., Fran¸cois P., Piersimoni A., Pedicelli S., Bono G.,Laney C. D., Primas F., Romaniello M., 2008, A&A, 490, 613Lequeux J., Peimbert M., Rayo J. F., Serrano A., Torres-PeimbertS., 1979, A&A, 500, 145Loebman S. R., Roˇskar R., Debattista V. P., Ivezi´c ˇZ., QuinnT. R., Wadsley J., 2011, ApJ, 737, 8Luck R. E., Lambert D. L., 2011, AJ, 142, 136Luck R. E., Kovtyukh V. V., Andrievsky S. M., 2006, AJ, 132,902Lynden-Bell D., Kalnajs A. J., 1972, MNRAS, 157, 1Ma X., Hopkins P. F., Faucher-Gigu`ere C.-A., Zolman N., Mura-tov A. L., Kereˇs D., Quataert E., 2016, MNRAS, 456, 2140Ma X., Hopkins P. F., Feldmann R., Torrey P., Faucher-Gigu`ereC.-A., Kereˇs D., 2017, Monthly Notices of the Royal Astro-nomical Society, 466, 4780Majewski S. R., et al., 2017, AJ, 154, 94Mannucci F., Della Valle M., Panagia N., 2006, MNRAS, 370, 773Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A.,2010, MNRAS, 408, 2115Maoz D., Graur O., 2017, ApJ, 848, 25Marigo P., 2001, A&A, 370, 194Mikolaitis ˇS., et al., 2014, A&A, 572, A33Minchev I., et al., 2018, Monthly Notices of the Royal Astronom-ical Society, 481, 1645Moll´a M., D´ıaz ´A. I., Cavichia O., Gibson B. K., Maciel W. J.,Costa R. D. D., Ascasibar Y., Few C. G., 2019a, MNRAS,482, 3071Moll´a M., et al., 2019b, MNRAS, 490, 665Muratov A. L., Kereˇs D., Faucher-Gigu`ere C.-A., Hopkins P. F.,Quataert E., Murray N., 2015, MNRAS, 454, 2691Muratov A. L., et al., 2017, MNRAS, 468, 4170Nieva M. F., Przybilla N., 2012, A&A, 539, A143Nomoto K., Tominaga N., Umeda H., Kobayashi C., Maeda K.,2006, Nucl. Phys. A, 777, 424Ostdiek B., et al., 2020, A&A, 636, A75Patr´ıcio V., et al., 2019, MNRAS, 489, 224Pedicelli S., et al., 2009, A&A, 504, 81 Pilyugin L., Grebel E., Kniazev A., 2014, The Astronomical Jour-nal, 147, 131Planck Collaboration et al., 2018, arXiv e-prints, p.arXiv:1807.06209Poetrodjojo H., et al., 2018, MNRAS, 479, 5235Prantzos N., Boissier S., 2000, Monthly Notices of the Royal As-tronomical Society, 313, 338Price-Jones N., et al., 2020, Monthly Notices of the Royal Astro-nomical Society, 496, 5101Price-Whelan A. M., et al., 2018, AJ, 156, 123Queyrel J., et al., 2012, A&A, 539, 93Rennehan D., Babul A., Hopkins P. F., Dav´e R., Moa B., 2019,MNRAS, 483, 3810Rudolph A. L., Fich M., Bell G. R., Norsen T., Simpson J. P.,Haas M. R., Erickson E. F., 2006, ApJS, 162, 346Sakhibov F., Zinchenko I. A., Pilyugin L. S., Grebel E. K., JustA., V´ılchez J. M., 2018, MNRAS, 474, 1657S´anchez-Menguiano L., et al., 2016, A&A, 587, A70S´anchez-Menguiano L., S´anchez S. F., P´erez I., Ruiz-Lara T.,Galbany L., Anderson J. P., Kuncarayakti H., 2020, MonthlyNotices of the Royal Astronomical Society, 492, 4149S´anchez S. F., et al., 2015, A&A, 573, 105Sanders N. E., Caldwell N., McDowell J., Harding P., 2012, ApJ,758, 133Santistevan I. B., Wetzel A., El-Badry K., Bland-Hawthorn J.,Boylan-Kolchin M., Bailin J., Faucher-Gigu`ere C.-A., Benin-casa S., 2020, MNRAS, 497, 747Sch¨onrich R., Binney J., 2009, MNRAS, 396, 203Solar M., Tissera P. B., Hernandez-Jimenez J. A., 2020, MNRAS,491, 4894Spitoni E., Cescutti G., Minchev I., Matteucci F., Silva AguirreV., Martig M., Bono G., Chiappini C., 2019, A&A, 628, A38Steinmetz M., et al., 2006, AJ, 132, 1645Stott J. P., et al., 2014, MNRAS, 443, 2695Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Gigu`ere C.-A.,Kereˇs D., Ma X., Robles V. H., 2017, MNRAS, 471, 144Swinbank A. M., Sobral D., Smail I., Geach J. E., Best P. N.,McCarthy I. G., Crain R. A., Theuns T., 2012, MNRAS, 426,935Takada M., et al., 2014, PASJ, 66, R1The MSE Science Team et al., 2019, arXiv e-prints, p.arXiv:1904.04907Ting Y.-S., De Silva G. M., Freeman K. C., Parker S. J., 2012,MNRAS, 427, 882Tissera P. B., Rosas-Guevara Y., Bower R. G., Crain R. A., LagosC. d. P., Schaller M., Schaye J., Theuns T., 2019, MonthlyNotices of the Royal Astronomical Society, 482, 2208Tremonti C. A., et al., 2004, ApJ, 613, 898Vogt F. P. A., P´erez E., Dopita M. A., Verdes-Montenegro L.,Borthakur S., 2017, A&A, 601, 61Wang X., et al., 2019a, arXiv e-prints, p. arXiv:1911.09841Wang C., et al., 2019b, MNRAS, 482, 2189Wenger T. V., Balser D. S., Anderson L. D., Bania T. M., 2019,ApJ, 887, 114Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Gigu`ere C.-A.,Kereˇs D., Quataert E., 2016, The Astrophysical Journal Let-ters, 827, L23Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tor-natore L., 2009, MNRAS, 399, 574Wojno J., et al., 2016, MNRAS, 461, 4246Wuyts E., et al., 2016, The Astrophysical Journal, 827, 74Yuan T. T., Kewley L. J., Rich J., 2013, ApJ, 767, 106Zinchenko I. A., Pilyugin L. S., Grebel E. K., S´anchez S. F.,V´ılchez J. M., 2016, MNRAS, 462, 2715Zinchenko I., Just A., Pilyugin L., Lara-Lopez M., 2019, Astron-omy & Astrophysics, 623, A7Zurita A., Bresolin F., 2012, MNRAS, 427, 1463de Jong R. S., et al., 2019, The Messenger, 175, 3MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Figure A1.
The 1- σ host-to-host scatter of the scaled radial gra-dients for different scale radii at z = 0. We define R , R , and R for the gas and stars in Table 1. R disk is the exponential scalelength of the stellar disk determined via a 2-component fit to thesurface density, and R phys is the physical radial coordinates ofthe disk, i.e. unscaled coordinates. For each scale radius, we mea-sure the gradient of all galaxies across an equal radial range thatcorresponds to 4 −
12 kpc physical for the galaxy with the me-dian scale length. The 1- σ host-to-host scatter is smallest whenmeasuring the gradients in physical space, which motivates ourchoice for our analysis in this paper.van den Hoek L. B., Groenewegen M. A. T., 1997, A&AS, 123,305 APPENDIX A: SCALED RADIAL PROFILES
Fig. A1 compares the host-to-host scatter in radial gradientsof [O/H] in gas in our simulated galaxies when scaling thesegradients to various galaxy scale radii at z = 0. We scale eachgalaxy’s profile using: R , R , and R for the gas and thestars, along with the exponential scale length, R disk , from a2-component (s´ersic plus exponential) fit to the surface den-sity. Table 1 lists the values for each galaxy. We also com-pare these scaled gradients to the gradients in physical radii(as measured in Section. 3.2. We bin each profile equally inscaled radius, defining the bin width such that the galaxywith the median scale length has a binwidth of 0 .
25 kpc phys-ical. We measure the radial gradient of each galaxy across anequal radial range for each scale radius. We define this radialrange such that we measure the galaxy with the median scalelength across a physical range 4 −
12 kpc. This range corre-sponds to: ≈ . − . R star25 , ≈ . − . R star50 , ≈ . − . R star90 , ≈ . − . R gas25 , ≈ . − . R gas50 , ≈ . − . R gas90 , and ≈ . − . R disk .Measuring the gas abundance radial gradient (from4 −
12 kpc) in physical space minimizes the host-to-host scat-ter, to σ ≈ .
005 dex. R star25 has the next smallest 1- σ scatterwith σ ≈ .
009 dex. The gradients are the least self-similarwhen scaled by R gas90 , for which the 1- σ scatter is ≈ .
1. Theself-similarity of the radial profiles in physical space moti-vates our choice in this paper, because there is no compelling reason to scale the profiles of our galaxies. We emphasize,though, that this may be a result of the small mass range ofour suite (halo masses are M = 1 − × M (cid:12) , stel-lar masses are in Table 1) and may not be generalizable togalaxies across a wider mass range. APPENDIX B: ALL GAS VERSUS STAR-FORMING GAS
In this paper, we examine elemental abundances in all gas, asinitial conditions for chemical tagging of stars. We choose tomeasure all gas in part because star-forming gas representsonly a small fraction of all gas elements at a given snapshot,leading to significant Poisson noise. In principal, we couldattempt to identify photo-ionized (HII) regions near youngstar particles to compare with gas-phase measurements vianebular emission lines, but doing this correctly requires gen-erating synthetic observations via ray-tracing, which is be-yond the scope of our analysis. In future work (Bellardini etal., in prep.) we will compare in detail the spatial variationsin abundance of star particles that form out of this gas tothe gas itself. Here, we explore the impact of measuring onlystar-forming gas instead of all gas.Fig. B1 compares measuring [O/H] in star-forming ver-sus all gas at z = 1 ( t lookback = 7 . z = 0. Foreach galaxy, we select gas elements at 4 < R <
12 kpc and | Z | < ≈
200 Myr), because at any single snapshot there arefew star-forming gas elements. For reference, for these samesimulations at z = 0, Benincasa et al. (2020) find typicalGMC lifetimes, and hence lifetimes of star-forming regions,of 5 − ≈ .
04 dex at z = 1 and ≈ .
01 dex at z = 0. The dif-ference in [O/H] is typically (cid:46) .
02 dex for z = 0 and alwaysless than 0 .
03 dex. The discrepancy is larger at higher red-shift, the difference is typically (cid:46) .
04 dex and always lessthan 0 .
13 dex. This is likely because cosmic accretion andstar-formation rates are higher at earlier times, leading toless efficient small-scale mixing of gas. Of course, a simpleoffset in the [O/H] normalization does not alone mean thatspatial variations are different.Fig. B1 (bottom row) shows the difference in the stan-dard deviation of star-forming versus all gas. Again, theblack line shows the mean value. On average, [O/H] for star-forming gas has slightly smaller standard deviation than forall gas. This difference is larger at z = 1 than at z = 0.However, the difference is typically small, (cid:46) .
05 dex. Thissuggests that the azimuthal variations of star-forming gasmay be smaller than that of all gas, especially if the scatteris driven primarily by radial variations in abundance. Thus,chemically tagging the birth radii of stars may may be com-plicated by azimuthal variations for redshifts higher than weshow in Sec. 3.5, which we will explore further in Bellardiniet al. in prep.We also explore the differences in the radial gradientsfor star-forming versus all gas (not shown). At z = 0 the MNRAS000
05 dex. Thissuggests that the azimuthal variations of star-forming gasmay be smaller than that of all gas, especially if the scatteris driven primarily by radial variations in abundance. Thus,chemically tagging the birth radii of stars may may be com-plicated by azimuthal variations for redshifts higher than weshow in Sec. 3.5, which we will explore further in Bellardiniet al. in prep.We also explore the differences in the radial gradientsfor star-forming versus all gas (not shown). At z = 0 the MNRAS000 , 1–23 (2021) M. A. Bellardini et al.
Figure B1.
A comparison of [O/H] measured in all gas versus only star-forming gas, at 4 < R <
12 kpc and | Z | < ≈
200 Myr to boost the number ofstar-forming gas elements), and we compute the galaxy-wide difference between star-forming gas and all gas. Top panels show histogramsof the difference in the mean [O/H], while bottom panels show histograms of the difference in the standard deviation of [O/H]. Left panelsshow z = 0 and right panels show z = 1 ( t lookback = 7 . ≈ .
04 dex at z = 1 and ≈ .
01 dex at z = 0. Furthermore, star-forming gas is slightlybetter mixed (with less scatter), with σ [O / H] ≈ .
05 dex smaller at z = 1 and ≈ .
008 dex smaller at z = 0. radial gradients of star-forming and all gas are consistent towithin ± .
005 dex kpc − , less than the host-to-host scatterin Fig. 5. At z = 1, for the galaxies with sufficient star-forming gas to measure a reliable radial gradient, they agreewith the gradients for all gas to within ± .
002 dex kpc − .We also compare compare the radial profiles of [O/H] fornewly formed stars (in age bins of 200 Myr) with that ofall gas at the same snapshots. For z (cid:46) .
5, these profilesoverlap to within uncertainty. For z (cid:38) .
5, the stellar andgas abundance profiles start to diverge, such that the youngstars tend to have higher [O/H] and flatter gradients thanall gas, which we will explore further in Bellardini et al., inprep.In summary, while analyzing all gas is a reasonable, ifnot perfect, proxy for star-forming gas in our simulations.Given the short lifetimes of star-forming gas clouds (Benin-casa et al. 2020) and the strict conditions for particles to bestar-forming (Hopkins et al. 2018), only a small fraction ofgas elements are star-forming at a given snapshot ( (cid:46)
1% ofgas particles at the redshifts we observe), so analyzing all gasgreatly reduces the statistical uncertainty. In the future wewill study the abundances of star particles that form fromthis gas, to compare with these results in detail.
APPENDIX C: IMPACT OF DIFFUSION COEFFICIENT
Our FIRE-2 simulations model the sub-grid diffu-sion/mixing of metals in gas via unresolved turbulent eddies(Su et al. 2017; Escala et al. 2018; Hopkins et al. 2018): ∂Z i ∂t + ∇ · ( D ∇ Z i ) = 0 (C1)where Z i is the mass fraction of a metal in gas element i ,and D is the diffusion coefficient. While there is some uncer-tainty in the exact value to choose for this coefficient, ourfiducial value is physically motivated based on tests of themetal diffusion implementation in FIRE-2 on idealized, con-verged turbulent box simulations by Colbrook et al. (2017)and other more extensive studies by Rennehan et al. (2019).Here, we compare our key results using our fiducial diffusioncoefficient D in m12i against a re-simulations of m12i withall identical physics/parameters, except one has a diffusioncoefficient that is 10 times higher (that is, faster mixing) andthe other includes no subgrid mixing.Fig. C1 compares the vertical, radial, and azimuthalvariations for m12i. The left panel shows the vertical gra-dient in gas [O/H], similar to Fig. 6. At z = 0, we find nodifferences within 200 pc and at most ∼ .
015 dex difference
MNRAS , 1–23 (2021)
D gas abundances in FIRE MW-mass simulations Figure C1.
Vertical (left), radial (middle), and azimuthal variations (right) in [O/H] between our fiducial simulation of m12i, a versionwith no subgrid metal diffusion, and a re-simulation increasing the diffusion coefficient by 10 times. The vertical profiles show no clearsystematic variations at a level important for our analysis. We normalize the radial profiles to 4 kpc (the approximate edge of thebulge) for clarity in comparison. The radial gradients (measured from 4 −
12 kpc for z = 0, 2 − z = 1) vary by no more than ≈ .
005 dex kpc − between our fiducial simulation and the simulation with 10 times higher metal diffusion, while the simulation with nometal diffusion has a ≈ .
13 dex kpc − steepe gradient at z = 0. In the right panel, we scaled down the azimuthal scatter in the simulationwith no metal diffusion by a factor of 10, for comparison. Thus, neglecting metal diffusion/mixing leads to 10 × higher azimuthal scatter,and moreover, this scatter does not depend much on azimuthal scale. The enhanced metal diffusion re-simulation shows smaller azimuthalscatter at small azimuthal scales, given the enhanced mixing rate on these small scales. However, that simulation shows similar scatterat large azimuthal scales, indicating that disk-wide azimuthal scatter is not sensitive to the detailed choice of diffusion coefficient. at 1 kpc. The differences are stronger at z = 1 for 10 × higherdiffusion and stronger at z = 0 for the simulation with nodiffusion, though again, not at a significant level to changeour interpretations, especially within 200 pc.Fig. C1 (center) shows the radial profile in gas [O/H],normalized to the abundance at R = 4 kpc (given a strongupturn at smaller R ). The radial gradients, measured overour fiducial radial ranges, vary by (cid:46) .
005 dex kpc − be-tween the 10 × diffusion simulation and the fiducial simu-lation. The gradients vary by less than 0 .
014 dex betweenthe fiducial simulation and the one with no metal diffusionat z = 0. The simulation with no subgrid diffusion has asteeper gradient at z = 0, potentially because the metals areless efficient at spreading from a given radius, in the absenceof subgrid diffusion, once the disk has become rotationallydominated and radial turbulence is no longer efficient atmoving the gas particles. We thus conclude that the radialgradients are reasonably robust to choices of the strengthof the diffusion coefficient, however, in the unphysical caseof no subgrid diffusion, the gradient can be (unphysically)steeper.Fig. C1 (right) compares the azimuthal variations ver-sus angular bin width, at R = 8 kpc. The simulation withno subgrid diffusion has 10 × higher azimuthal scatter, sowe scale down its values in Fig. C1 by 10 × for visual com-parison. Using no subgrid diffusion leads to scatter that islargely independent of azimuthal scale at high z , but that in-creases with azimuthal scale at low z . Without subgrid diffu-sion, a small number of gas elements can absorb most of themetals, while a significant number of (neighboring) elementscan remain nearly un-enriched. This is a patently unphysicalscenario, and it yields azimuthal scatter that disagrees withobservations by an order of magnitude. At both redshifts, us-ing a higher diffusion coefficient leads to smaller azimuthalvariations at small scales, because diffusion smooths varia-tions between nearby gas elements on scales approaching theresolution (Escala et al. 2018). However, the azimuthal varia-tions are nearly unchanged on large azimuthal scales. There- fore, our results on small azimuthal scales are likely sensitiveto the exact choice of diffusion coefficient, but the large-scaleazimuthal variations are robust. An important caveat to thiscomparison is that it is only one simulated galaxy, and in-dividual simulations with the same initial conditions andphysics can show non-trivial stochastic variations from ran-dom number generators, floating-point roundoff, and chaoticbehavior Keller et al. (e.g. 2019). Indeed, we find minor fluc-tuations between these simulations, for example, in the ex-act timing of mergers, which can affect all panels in Fig. C1.We consider it likely that the differences in azimuthal vari-ations on small scales are robust, but any other differencesin Fig. C1 are potentially stochastic. This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000