# Figuring Out Gas & Galaxies In Enzo (FOGGIE) VIII:
Complex and Stochastic Metallicity Gradients at $z>2$

Ayan Acharyya Department of Physics & Astronomy, Johns Hopkins University,
3400 N. Charles Street, Baltimore, MD 21218 Ayan Acharyya<EMAIL_ADDRESS>Molly S. Peeples Space Telescope Science Institute, 3700 San Martin Dr.,
Baltimore, MD 21218  Department of Physics & Astronomy, Johns Hopkins
University, 3400 N. Charles Street, Baltimore, MD 21218  Jason Tumlinson Space
Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218
Department of Physics & Astronomy, Johns Hopkins University, 3400 N. Charles
Street, Baltimore, MD 21218  Brian W. O’Shea Department of Computational
Mathematics, Science, and Engineering, Department of Physics and Astronomy,
National Superconducting Cyclotron Laboratory, Michigan State University
Cassandra Lochhaas Space Telescope Science Institute, 3700 San Martin Dr.,
Baltimore, MD 21218  Center for Astrophysics, Harvard & Smithsonian, 60 Garden
St., Cambridge, MA 02138  NASA Hubble Fellow  Anna C. Wright Department of
Physics & Astronomy, Johns Hopkins University, 3400 N. Charles Street,
Baltimore, MD 21218  Raymond C. Simons Space Telescope Science Institute, 3700
San Martin Dr., Baltimore, MD 21218  Department of Physics, University of
Connecticut, 196A Auditorium Road Unit 3046, Storrs, CT 06269 USA  Ramona
Augustin Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD
21218  Leibniz Institute for Astrophysics Potsdam (AIP), An der Sternwarte 16,
Potsdam 14482, Germany  Britton D. Smith Institute for Astronomy, University
of Edinburgh, Royal Observatory, EH9 3HJ, UK  Eugene Hyeonmin Lee Space
Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218
University of Texas at Austin, Physics, Math, and Astronomy Building, 2515
Speedway, Austin, TX 78712

###### Abstract

Gas-phase metallicity gradients are a crucial element in understanding the
chemical evolution of galaxies. We use the FOGGIE simulations to study the
metallicity gradients ($\nabla Z$) of six Milky Way-like galaxies throughout
their evolution. FOGGIE galaxies generally exhibit steep negative gradients
for most of their history, with only a few short-lived instances reaching
positive slopes that appear to arise mainly from interactions with other
galaxies. FOGGIE concurs with other simulation results but disagrees with the
robust observational finding that flat and positive gradients are common at
$z>1$. By tracking the metallicity gradient at a rapid cadence of simulation
outputs ($\sim 5$–10 Myr), we find that theoretical gradients are highly
stochastic: the FOGGIE galaxies spend $\sim 30-50$% of their time far away
from a smoothed trajectory inferred from analytic models or other, less high-
cadence simulations. This rapid variation makes instantaneous gradients from
observations more difficult to interpret in terms of physical processes.
Because of these geometric and stochastic complications, we explore non-
parametric methods of quantifying the evolving metallicity distribution at
$z>1$. We investigate how efficiently non-parametric measures of the 2-D
metallicity distribution respond to metal production and mixing. Our results
suggest that new methods of quantifying and interpreting gas-phase metallicity
will be needed to relate trends in upcoming high-$z$ JWST observations with
the underlying physics of gas accretion, expulsion, and recycling in early
galaxies.

interstellar medium — galaxy evolution — metallicity — metallicity gradient

##  1 Introduction

The spatial distribution of metals in the gas phase of the interstellar medium
is a crucial diagnostic for the chemical evolution of galaxies.
Conventionally, the distribution of metals is quantified in terms of the slope
of the azimuthally-averaged metallicity profile, the “metallicity gradient”,
defined as

$\nabla Z=\derivative{\log{Z/Z_{\odot}}}{r}.$ (1)

A negative $\nabla Z$ implies that the metallicity decreases radially, whereas
a positive $\nabla Z$ denotes a metal-poor center of the disk relative to the
outskirts. Metallicity gradient studies have helped us make considerable
progress in understanding galaxy evolution at low-redshift. Thanks to
spatially resolved spectroscopy and advanced models, the last $\sim$two
decades has witnessed a plethora of such studies.

Steep negative metallicity gradients have improved our understanding of the
inside-out star-formation scenario whereby stars are formed first in the
centers of galaxies and subsequently at outer radii (Marino et al., 2016;
Belfiore et al., 2017), allowing more time for metal enrichment in galaxy
centers than galaxy outskirts. Mingozzi et al. (2020) attribute the steepening
of the metallicity gradient with increasing stellar mass to increasing star-
formation activity, and the flattening of gradient at the highest stellar
masses to the increasing prevalence of mergers (see also Belfiore et al.,
2017). Thus, negative gradients can teach us about the star-formation history
and mass assembly of galaxies.

Flat gradients, on the other hand, imply efficient, radially outward mixing of
metals via gas flows (e.g., Sánchez-Menguiano et al., 2018; Simons et al.,
2021). Sánchez-Menguiano et al. (2018) suggest that, for galaxies (or regions
of galaxies) with flat metallicity gradients, this could be an indicator that
the metal budget in these regions is dominated by gas flows, rather than in
situ production. Moreover, interactions and merger events can lead to
efficient gas mixing, and consequently, flat or positive gradients (e.g.,
Krabbe et al., 2008; Kewley et al., 2010; Rich et al., 2012; Miralles-
Caballero et al., 2014; Vogt et al., 2015; Muñoz-Elgueta et al., 2018). Thus,
flat or positive gradients help us learn about gas flows and interactions in
the local Universe.

Typically, a mix of negative (e.g. Vila-Costas & Edmunds, 1992; Garnett et
al., 1997; Bresolin et al., 2002; Kennicutt et al., 2003; Bresolin, 2007;
Moustakas et al., 2010; Cecil et al., 2014; Sánchez et al., 2014) and flat or
positive metallicity gradients (e.g. Sánchez-Menguiano et al., 2016; Belfiore
et al., 2017; Molina et al., 2017; Poetrodjojo et al., 2018; Sánchez-Menguiano
et al., 2018) are observed in the local Universe. However, a significant
fraction of the high-redshift ($z>1$) measurements of metallicity gradients
are flat or positive (Wuyts et al., 2016; Wang et al., 2017, 2019; Curti et
al., 2020; Simons et al., 2021). This is in stark contrast with recent
simulations with a wide range of numerical methods and input physics (Gibson
et al., 2013; Ma et al., 2017; Hemler et al., 2021), which generally show
negative gradients at all redshifts.

With its high angular resolution imaging spectroscopy and IFU capabilities,
the James Webb Space Telescope (JWST) is able to map spatially-resolved
metallicity gradients for statistically large samples of galaxies at $z>2$
(Wang et al., 2022; Venturi et al., 2024). For instance, Wang et al. (2022)
measured the first metallicity gradient with the JWST as part of the GLASS-
JWST-ERS survey (Grism Lens Amplified Survey from Space). They reported a
significantly positive metallicity gradient for a 109 M⊙ galaxy at $z=3$. This
single case already intensifies the tension between observed and theoretical
metallicity gradients, particularly at high redshift. Therefore, it is timely
that we employ theoretical models to aid the interpretation of the forthcoming
high-$z$ metallicity gradient observations.

In this work, we use simulations from the FOGGIE project (“Figuring Out Gas &
Galaxies in Enzo”). FOGGIE ran high spatial resolution ($\sim 270$ comoving
pc) cosmological zoom-in simulations of six Milky-Way mass galaxies. In
addition to very high resolution in the CGM and the disk-halo interface,
FOGGIE generates a high cadence of outputs (every $\sim$ 5 Myr), which allows
us to assess the stochastic variation of metallicity gradients over time. As
we show in this paper, these attributes of FOGGIE demonstrate that the
traditional approach of quantifying spatial distribution of metals via a
single, azimuthally-averaged slope is inadequate in capturing the full
distribution of metals in high-$z$ galaxies. We therefore call for a more
informative measure of the spatial distribution of ISM metallicity.

This paper is organised as follows. In Section 2 we describe the FOGGIE
simulations and the methods we apply in our analysis. We present the
corresponding results in Section 3. We discuss the implications of positive
metallicity gradients in the context of the literature in Section 4. In
Section 5 we discuss what current studies, including our work, hint at the
physics governing the metallicity gradient evolution. We note the caveats and
challenges for this topic of study, and put them in the context of the
upcoming JWST observations in Section 6, followed by the summary in Section 7.

##  2 Analysis methods

This section details our simulations (§2.1) and how we select out the galactic
disk (§2.2). Section 2.3 lays out a detailed description of the radial
metallicity profiles, Section 2.4 considers the effects of mapping 3-D
simulations to 2-D projections, and Section 2.5 describes the non-parametric
quantification of the metallicity distribution in detail.

###  2.1 The FOGGIE simulations

Figure 1: Top: Mass-weighted, projected metallicity field of one of the FOGGIE
halos, Hurricane, at $z=0.07$. Middle: Corresponding radial metallicity
profile out to 10 $h^{-1}$ckpc. The blue 2-D histogram denotes the density of
cells at a given value of metallicity and radius, with lighter blue implying a
higher concentration. The orange circles denote mass-weighted mean metallicity
in radial bins of size 1 $h^{-1}$ckpc. The dashed line is a linear fit to the
orange circles, the slope of which we define as the metallicity gradient.
Bottom: Histogram of the mass-weighted metallicity out to 10 $h^{-1}$ckpc is
denoted in orange. The vertical lines denote various percentiles as indicated
in the text boxes.

We use cosmological zoom-in simulations from the Figuring Out Gas & Galaxies
In Enzo (FOGGIE) project111https://foggie.science/index.html (Peeples et al.,
2019; Corlies et al., 2020; Zheng et al., 2020; Simons et al., 2020; Lochhaas
et al., 2021, 2023; Wright et al., 2023). The simulations have been described
in detail in previous FOGGIE papers; however, we summarize them briefly here.

Firstly, a large cosmological box of 100 $h^{-1}$ comoving Mpc was run with
Enzo at relatively low resolution to provide halo catalogs from which to
select for zooms. Six halos were picked such that they are Milky-Way mass at
$z=0$. Then, these six halos were re-simulated with Enzo (Bryan et al., 2014;
last described in Brummel-Smith et al., 2019) within a “zoom in” region.
FOGGIE is distinguished from most other cosmological simulations by a unique
adaptive mesh refinement (AMR) scheme that allows “forced refinement” within a
refine box that follows the zoomed galaxy through the domain. Inside this
moving box, 200 $h^{-1}$ comoving kpc on a side, there is a minimum level of
refinement, corresponding to a maximum cell size, regardless of density or
other physical properties. As first presented in FOGGIE Paper IV (Simons et
al., 2020), the current production FOGGIE simulations impose a minimum level
of refinement $n_{\rm ref}=9$ ($1{,}100$ comoving pc cells) everywhere within
the “refined box”. In order to ensure that the cooling length is resolved, we
also employ a “cooling refinement” prescription that refines cells larger than
the cooling length—determined by the local sound speed and cooling time—up to
a maximum refinement level $n_{\rm ref}=11$ (274 comoving pc cells). Since
both forced and cooling refinement are active in the moving box, all of the
ISM and almost all of the CGM ($\gtrsim 95\%$) are resolved at the highest
level, i.e., 91 pc at $z=2$ (Simons et al., 2020), and the median cell mass in
the CGM and disk-halo interface is $<200$ $\mathrm{M}_{\odot}$(Lochhaas et
al., 2023). The high resolution enables FOGGIE to capture small scale dynamic
processes such as metal mixing and transport, thereby making these simulations
well suited for our study (see also van de Voort et al., 2019; Hummels et al.,
2019). The six halos were selected to have no major mergers after $z=2$. This
serves as the ideal “control” scenario for our study because we can compare
the behavior of the metallicity gradient with ($z\gtrsim 2$) and without
($z\lesssim 2$) major mergers. We refer the reader to Table 1 for a detailed
list of global properties of all the halos.

We use the $z=0.07$ snapshot of the Hurricane halo to demonstrate our methods.
We present our main results from all six halos—Tempest, Maelstrom, Squall,
Blizzard, Hurricane, Cyclone—for the first part of our analysis where we
compare the cosmic evolution of metallicity gradients of FOGGIE galaxies in
the context of current observations. Thereafter, we focus on one halo,
Hurricane, for the rest of this work where we investigate potential factors
that impact the metallicity gradient evolution seen in the FOGGIE simulations.
We choose this halo as it has the most interactions out of all the FOGGIE
halos, and therefore has a diverse range of gradients. While we present
detailed analysis for only Hurricane in the main body of the paper for
simplicity, we show in Appendix A that our main conclusions remain
qualitatively unchanged for the other halos as well.

###  2.2 Selection criteria for the gas disk

Observed metallicity gradients are limited to the gaseous, star-forming
regions of galaxies and their immediate environment. We therefore limit our
analysis to within 10 $h^{-1}$ comoving kpc of the galaxy center. We choose
this extent because it is equivalent to $\approx 2.5\times r_{e}$, where
$r_{e}\sim 4$ kpc is a typical effective radius of the gas disk of star-
forming galaxies at $z\sim 0$. Additionally, we impose a simple redshift-
dependent density criterion such that $\rho_{\rm cut}=2\times 10^{-26}$ g cm-3
if $z>0.5$, $\rho_{\rm cut}=2\times 10^{-27}$ g cm-3 if $z<0.25$, and a linear
interpolation with respect to physical time (not redshift) between $0.25\leq
z\leq 0.5$. Note that this density criterion is based on proper density, and
not comoving density. We select only the material with density above
$\rho_{\rm cut}$, and apply our analysis methods on this dense material. The
choice of this density criterion was used by Lochhaas et al. (2021) to capture
the dense ISM material. We deliberately do not impose any cylindrical or
geometrical criteria to capture the ISM disk, because the “disk” is likely to
have irregular, clumpy morphology, particularly at high redshifts (see Section
2.5) which would not be captured by a well-behaved cylindrical geometry. It is
worth noting that even at low-$z$ FOGGIE galaxies have significant warps
(Lochhaas et al., 2023, R. Simons et al. in prep).

###  2.3 Radial metallicity gradient

Figure 1 shows (a) the projected gas phase metallicity at the top, projected
along the arbitrary $y$ axis, (b) the radial mass-weighted metallicity profile
in the middle, and (c) the mass-weighted histogram at the bottom. All the
panels correspond to the $z=0.07$ snapshot of the Hurricane halo, and show
metallicities within a 10 $h^{-1}$ comoving kpc sphere (corresponds to 13
physical kpc at $z=0.07$), after employing the density criterion described
above. We provide figures similar to Figure 1, but corresponding to the the
other FOGGIE halos and several redshift epochs, in Figures A1-A6 in Appendix
A.

The top panel reveals distinct “inner” and “outer” disks in the metallicity
space, where there is a steep drop in metallicity in the outer disk, compared
to the inner disk. We note that whenever there is such a steep drop in
metallicity between an “inner” and “outer” disk in the FOGGIE simulations,
active star formation is typically confined to the inner disk only.

The color-coding of the middle panel of Figure 1 represents a 2-D histogram of
the gas cells, with lighter blue implying a higher concentration of cells. In
this case there is a notable “break” in the metallicity profile between a
steeper trend in the central few kpc and a flatter, almost constant,
metallicity profile in the outskirts.

The outer profile is coincident with a diffuse, non-star-forming disk in the
outskirts. Although we see this “break” in the metallicity profile
ubiquitously in our simulations, we highlight a particularly prominent example
here to aid our explanation of the fitting routine.

The orange circles depict the mass-weighted average metallicity in spherical
bins of width $1h^{-1}$ comoving kpc (ckpc) centered at the galaxy center, and
a linear fit to them is shown by the dashed orange line. The resulting slope,
along with its uncertainty, is denoted in the orange text in the middle panel
of Figure 1. Note that only the gas that passes the disk density cut (Section
2.2) is included here. The line is fitted to the mass-weighted mean for each
bin, and the natural, physical variation around each mean value is treated as
an “uncertainty” in the inverse-square weighting, i.e., each bin receives a
weight equal to the inverse square of the “uncertainty”. Owing to mass-
weighting, the central bins dominate the fit because that is where most of the
disk mass is concentrated.

Observed gas-phase metallicity profiles are primarily derived from emission
line maps tracing the ionized ISM. The primary drivers of the ionization in
such H II regions are young, massive stars. Therefore, observed profiles
preferentially trace a luminosity-weighted metallicity. The central regions of
Milky Way type galaxies typically not only have more stars but also more
massive bright stars, leading to a more luminous center than the outskirts
(Minniti & Zoccali, 2008). Therefore, the observed metallicity profile has a
higher weight-contribution from the center. In our case, the fit to the mass-
weighted radial bins preferentially traces the central part of the galaxy,
which is similar in effect to a fit to an observed metallicity profile.

We repeat this radial fitting routine for each snapshot of each halo, from
$z=4$ to the latest snapshot available. We start our analysis at $z=4$ because
at higher redshifts the “disk” morphology is not apparent and mergers are too
frequent and rapid to maintain a radial profile.

###  2.4 Evaluation of projection effects

Figure 2: Comparison of metallicity distribution as seen in projection vs the
full 3-D distribution. Top: Mass-weighted, projected metallicity of Hurricane
at $z=0.07$, as seen from three different lines of sight – $x$, $y$, and $z$
axes, with respect to the simulation box. Bottom left: Corresponding radial
profile out to $10h^{-1}$ckpc for the mass-weighted metallicity projected
along the three axes, with each color corresponding to an axis of projection.
The fitting process is outlined in Section 2.3. Bottom right: Same as the left
panel but for the full 3-D metallicity in the simulation, rather than
projected metallicity. The projected and 3-D metallicity distributions look
qualitatively similar.

Observed metallicity distributions are, necessarily, projected quantities
along the line of sight of the observation. Simulations, on the other hand,
have access to the full 3-D distribution of the physical property concerned.
Therefore, it is important to compare the full 3-D distribution of metals in
the FOGGIE simulations to its projected counterpart in order to assess the
potential impact of projection effects on simulated gradients. Our analysis
for the rest of this paper is based on the intrinsic 3-D metallicities in the
FOGGIE galaxies, but here we show that our qualitative results are not
expected to change significantly when accounting for projection effects.

In Figure 2 we demonstrate the effect of projection on the metallicity maps
and radial profiles as seen in projection from three different lines of sight
(along the cardinal axes), as well as on the full 3-D metallicity profile. The
top row of Figure 2 corresponds to the projected 2-D metallicity maps. The
bottom-left panel shows the azimuthally-averaged fits to the radial
metallicity profile, where each color corresponds to a different projection.
The bottom-right panel depicts the radial profile of the full 3-D metallicity
distribution, same as the middle panel of Figure 1. We do not see any
qualitative difference in the radial profiles due to the different angles of
projection. Indeed, the fitted gradients (colored text in Figure 2) obtained
in cases of projections along various lines of sight and that from the 3-D
distribution all agree with each other within 1$\sigma$ uncertainties.
Therefore, we conclude the projection effects do not impact the main science
results of this paper. Henceforth we proceed with the full 3-D metal
distribution for the rest of the paper, disregarding projection effects.

It is worth noting that forward modeling the simulated data to produce
synthetic observables is generally a good approach in comparing observations
with simulations. However, that is outside the scope of our current work for
the reasons described in Section 6. Therefore, all measurements of metallicity
discussed in this work refer to the intrinsic metallicity of the gas cells in
the FOGGIE simulations.

###  2.5 Non-parametric characterization of the metallicity distribution

Fitting a single, smooth gradient to the entire radial metallicity profile of
the disk, while informative, may not necessarily be the best representation of
the full 2-D distribution of metallicity. For instance, “breaks” in radial
profiles and azimuthal variations are washed out in the conventional gradient
approach. Moreover, defining a gradient makes assumptions about a disk
geometry, which may not hold at high-$z$. We therefore explore a non-
parametric approach to quantify the distribution of metallicity. We discuss
the detailed implications of this approach in Section 5.3. We note that a
similar approach of quantifying the metallicity distribution function (MDF)
has already been utilised by the stellar metallicity community to great
scientific success, ranging from the Milky Way (Hayden et al., 2015; Loebman
et al., 2016) to stellar halos (Durrell et al., 2001) and satellites (Kirby et
al., 2011; Chiti et al., 2020).

We characterize the spatial distribution of metallicity within $10h^{-1}$
comoving kpc from the galaxy center by representing it as a histogram and
computing the 25th, 50th and 75th percentiles. The bottom panel of Figure 1
shows the metallicity distribution for Hurricane at $z=0.07$, along with
vertical lines depicting the different percentiles. We compute the inter-
quartile range (ZIQR) as the difference between the 75th (Z75) and 25th (Z25)
percentiles, and use it, along with the median, to trace the evolution of the
metallicity distribution across cosmic time (Section 3.2). All the percentiles
are in $\log{(\mathrm{Z}_{\odot})}$ units, and therefore ZIQR is the
difference between Z75 and Z25 in log space.

In the snapshot depicted in Figure 1, we see a bimodal distribution: a narrow
peak at low metallicity, followed by a broad distribution at higher
metallicity. This is a manifestation of the ‘break’ in the radial metallicity
profile, as described in Section 2.3. We conclude that the low metallicity
peak in the distribution is contributed by an extended, metal-poor cold gas
disk on the outskirts of the galaxy. Such a bimodal (in metallicity) disk is
likely a result of misaligned disks (see Section 5.4), which survive for a
majority of the cosmic time, until very low redshift ($z\lesssim 0.3$).
Therefore, the metallicity distribution exhibits a double-peaked nature for
most snapshots, except at very low redshift when it becomes single-peaked.

Our attempts to fit these multi-peaked distributions with skewed Gaussian
models have proven to be non-trivial, and not sufficiently informative for
tracing the chemical evolution of the galaxies. This further justifies our
view that it is difficult to capture the underlying physics leading to the
full 3-D distribution of metals in the galaxy disk by parameterizing it with
simple 1D models.

Halo | $z$ | Time (Gyr) |  M⋆/M⊙ |  $\log Z_{\rm total}$/Z⊙ |  $\log Z_{\rm median}$/Z⊙ |  $\log Z_{\rm IQR}$/Z⊙ |  $\nabla Z$ (dex/kpc)
---|---|---|---|---|---|---|---
Tempest | $0.0$ | $13.62$ | $10.72$ | $-0.13$ | $-0.18$ | $0.36$ | $-0.064\pm 0.000$
Tempest | $1.0$ | $5.89$ | $10.56$ | $0.17$ | $\phantom{-}0.17$ | $0.12$ | $-0.096\pm 0.000$
Tempest | $2.0$ | $3.32$ | $10.19$ | $-0.10$ | $-0.31$ | $0.41$ | $-0.699\pm 0.001$
Maelstrom | $0.0$ | $13.62$ | $11.02$ | $0.13$ | $\phantom{-}0.12$ | $0.20$ | $-0.032\pm 0.000$
Maelstrom | $1.0$ | $5.89$ | $10.79$ | $-0.05$ | $-0.09$ | $0.47$ | $-0.143\pm 0.000$
Maelstrom | $2.0$ | $3.32$ | $10.48$ | $-0.20$ | $-0.42$ | $0.61$ | $-0.153\pm 0.000$
Squall | $0.0$ | $13.62$ | $11.07$ | $-0.08$ | $-0.38$ | $0.64$ | $-0.110\pm 0.000$
Squall | $1.0$ | $5.89$ | $10.45$ | $0.09$ | $\phantom{-}0.09$ | $0.21$ | $-0.080\pm 0.000$
Squall | $2.0$ | $3.32$ | $9.18$ | $0.07$ | $\phantom{-}0.00$ | $0.18$ | $-0.241\pm 0.000$
Blizzard | $0.0$ | $13.62$ | $11.14$ | $-0.13$ | $-0.20$ | $0.42$ | $-0.068\pm 0.000$
Blizzard | $1.0$ | $5.89$ | $10.91$ | $0.16$ | $\phantom{-}0.13$ | $0.23$ | $-0.059\pm 0.000$
Cyclone | $1.0$ | $5.90$ | $11.12$ | $-0.04$ | $-0.09$ | $0.09$ | $-0.036\pm 0.000$
Cyclone | $2.0$ | $3.32$ | $10.49$ | $0.16$ | $\phantom{-}0.19$ | $0.26$ | $-0.220\pm 0.000$
Hurricane | $0.0$ | $13.62$ | $11.37$ | $-0.20$ | $-0.45$ | $0.54$ | $-0.080\pm 0.000$
Hurricane | $1.0$ | $5.89$ | $11.04$ | $0.10$ | $\phantom{-}0.06$ | $0.28$ | $-0.124\pm 0.000$
Hurricane | $2.0$ | $3.32$ | $10.42$ | $0.08$ | $\phantom{-}0.08$ | $0.65$ | $-0.301\pm 0.001$
Table 1: A list of the measured quantities within 10 $h^{-1}$ ckpc for each
snapshot of each FOGGIE halo – stellar mass, total metallicity, median and
inter-quartile range (IQR) of the metallicity distribution, and the fitted
radial gradient (see Section 2.3). All the quantities are quoted in log10
scale, with the exception of redshift, time (in Gyr) and metallicity gradient
(in dex/kpc). This table is intended to merely function as an abridged version
of the full table which is available online in a machine-readable format at
https://ayanacharyya.github.io/data/master_table_Zpaper_upto10.0ckpchinv_wtby_mass_fitmultiple_wdencut_islog.txt.

In Table 1 we tabulate the fitted and other measured quantities for all FOGGIE
halos. Here we only include a few snapshots, intended to serve as examples.
The full table is available online in a machine readable format at
https://ayanacharyya.github.io/data/master_table_Zpaper_upto10.0ckpchinv_wtby_mass_fitmultiple_wdencut_islog.txt.

Figure 3: Top: Comparison of the evolution of FOGGIE metallicity gradients and
observed gradients from the literature. The colored lines show the evolution
of metallicity gradients (within 10 $h^{-1}$ckpc) of the six FOGGIE galaxies.
The rest of this figure was adapted from Figure 3 of Wang et al. (2022). The
grey circles denote observed gradients from a variety of literature sources –
ground-based integral-field spectroscopy (IFS) with (Swinbank et al., 2012;
Jones et al., 2013b; Leethochawalit et al., 2016; Förster Schreiber et al.,
2018) and without (Bundy et al., 2015) AO-assist, as well as space-based
slitless spectroscopy with HST/WFC3 (Jones et al., 2015; Wang et al., 2020;
Simons et al., 2021). Bottom: Comparison between FOGGIE (salmon) and
metallicity gradient evolution predicted by other simulations—MUGS (Gibson et
al., 2013, teal), FIRE (Ma et al., 2017, light blue) and TNG50 (Hemler et al.,
2021, darker blue). The density of observed gradients is shown by the hexbins,
with darker bins implying a higher density of observations.

##  3 Results

###  3.1 General patterns in radial gradients

Our starting point for results is the traditional radial metallicity gradient
as seen in data and in models (top panel of Figure 3). The observational data
for comparison is shown in the gray circles, as compiled by Wang et al. (2022)
in their paper presenting the first high-$z$ gradient measurement with JWST.
These data are predominantly ground-based IFU measurements or HST slitless
spectroscopy. Notably, the data at $z=1-2.5$ exhibits just as many flat or
positive gradients as negative ones.

The six FOGGIE galaxies have negative gradients for nearly all of their
history, with only one (Hurricane) briefly approaching or touching a flat
gradient (see Section 4). There is also a broad trend for the FOGGIE galaxies
to start with steeper negative gradients ($\sim-0.5$ to $-0.1$ dex/kpc) and
evolve toward shallower negative gradients ($\sim-0.3$ to $-0.1$ dex/kpc). The
bottom panel of Figure 3 depicts the spread in metallicity gradients predicted
by different simulations—MUGS (Gibson et al., 2013), FIRE (Ma et al., 2017)
and Illustris TNG (Hemler et al., 2021) in comparison to that of FOGGIE (see
also Figure 3 of Wang et al., 2022). This basic comparison between observed
and simulated gradients, in both the panels, shows a strong tension; none of
these state-of-the-art simulation suites yields gradients that are positive as
often as they are negative. This systematic offset is robust: the different
observational techniques employed here generally agree with one another about
a distribution scattered around a flat gradient, while the simulations all
occupy the negative range. This inconsistency poses a challenge to our
interpretation of metallicity gradients in terms of the production and
distribution of metals in galaxies.

FOGGIE generates snapshots every $\sim$5.3 Myr, which makes the rapid
stochastic fluctuations of the gradients readily apparent. This contrasts with
the $50-100$ Myr cadence of most previous simulations (e.g., Gibson et al.,
2013; Ma et al., 2017; Hemler et al., 2021), as shown in the bottom panel of
Figure 3. The considerable degree of stochastic fluctuation in the FOGGIE
metallicity gradient evolution demonstrates that metallicity gradients indeed
vary significantly over short timescales, suggesting that trends deduced from
outputs several Gyr apart may not be robust.

The abundance of observed positive gradients at all redshifts is in contrast
with simulated ones. On par with other simulations, FOGGIE too does not
generally reproduce the observed significantly positive metallicity gradients.
Notably, we do measure positive gradients in FOGGIE on a few
occasions—particularly for Hurricane, as shown in the top panel of Figure 3.
However, such instances in FOGGIE are short-lived, and stem from external
processes such as accretion or mergers, and as such are highly subject to
rapid fluctuations (see Section 4 for detailed discussion).

We do not see long-lived positive gradients in settled and well-behaved disks.
This could imply that the observed positive gradients are also instantaneous
states in the galaxies’ lifetimes. However, the ubiquity of the observed
positive gradients could potentially imply that the simulations are missing
some key physical process. It is worth noting that the metallicity
observations themselves potentially suffer from limited spatial resolution,
signal-to-noise issues, scatter between different metallicity diagnostics
used, and uncertainties in the underlying photoionization models. Although
these uncertainties can lead to artificially flatter gradients (e.g., Yuan et
al., 2013a; Acharyya et al., 2020), none of these effects have yet been
demonstrated to bias the observed gradients towards positive values.

Our results unambiguously show significant scatter in the metallicity gradient
evolution, particularly at high-redshift. Not only is there galaxy-to-galaxy
scatter of $\sim 0.1$–0.4 dex but there also is considerable stochasticity
(variations of $\sim 0.1$–0.2 dex in slope over a $\sim$ few Myr) in the
evolution of an individual galaxy. Therefore, the short timescale variations
that FOGGIE’s high time cadence data outputs allow us to study are not
negligible. Lack of a high time-cadence would lead to sparser sampling of this
stochastic behavior, and therefore such low-cadence studies would be less
likely to capture the large excursions the metallicity gradient makes along
its evolutionary course. We discuss this further in Section 5.5. Given the
stochastic behavior of the models, it would be challenging to interpret the
metallicity gradient evolution via forthcoming high-$z$ observations with JWST
and large ground based telescopes (e.g., the Giant Magellan Telescope; GMT).

####  3.1.1 Evolution of the radial gradient

The top panel of Figure 4 shows the metallicity gradient evolution of one
FOGGIE halo, Hurricane, as a function of time. The circles mark the times
corresponding to whole number redshift values i.e., $z=0$, 1, 2, and so on,
from right to left. The faint line in the background shows the evolution
sampled at a lower cadence of $\sim$ 500 Myr, whereas the thicker line in the
foreground represents the intrinsic cadence ($\sim 5$ Myr) of FOGGIE
simulations. We see an overall flattening of the metallicity gradient after
$z\sim 1$. We provide figures similar to the top panel of Figure 4, but
corresponding to the other FOGGIE halos, in Appendix A.

While the lower-cadence curve captures the overall general trend in
metallicity gradient, it misses out on the short timescale variations,
particularly at high-$z$ ($z\geq 1$). We therefore quantify the fraction of
evolution that would be otherwise “missed” by studies with such lower cadence.
The shaded area represents a typical observational uncertainty of $\pm 0.03$
dex/kpc around the lower-cadence evolution. We emphasise that the shaded area
does not represent the uncertainty in the simulated gradient, but the usual
amount of uncertainty seen in observational studies. Every time the solid line
exceeds the shaded region, we consider that as a significant excursion of the
metallicity gradient evolution away from the “general” (smoothed) behavior.
This behavior suggests that lower-cadence simulations or static models are
missing a significant degree of the stochastic variation. But observations are
snapshots in time and can easily catch a real galaxy in one of these
excursions. We find that this galaxy is more than the typical observational
uncertainty away from its general behavior for $\sim 50$ % of the time at
$z>1$. In other words, a cadence of $\sim 500$ Myr fails to capture almost
half of this galaxy’s evolution at high-$z$. Observations, by definition,
capture a single snapshot during a galaxy’s lifetime. Our work demonstrates
that, for a Milky-Way type galaxy, it is challenging to interpret an
instantaneous high-$z$ metallicity gradient observation because any snapshot
may be a large and/or short-lived excursion away from the general pattern.
This further highlights the challenges to interpreting forthcoming high-$z$
metallicity gradient measurements with JWST or GMT.

####  3.1.2 Rapid variations in radial gradient

An example of rapid variation of the radial metallicity gradient at high-
redshift is demonstrated in the lower two panels of Figure 4. The middle and
bottom panels of Figure 4 allow a closer look at specific time stamps during
the evolution of Hurricane, denoted by the brown and green star-symbols in the
top panel, “A” and “B”, respectively. “A” displays a strong negative
metallicity gradient, whereas “B”—only $\sim 50$ Myr later—shows a
significantly flatter (but still negative) gradient. Each row of these two
bottom panels shows, from left to right, the projected mass-weighted
metallicity map, the mass-weighted line-of-sight velocity map, the radial
metallicity profile, and the mass-weighted metallicity distribution. The
projected metallicity map, radial profile and distribution plots follow the
same template as in Figure 1. The line-of-sight velocity map depicts the
kinematics of the gas in the disk at each snapshot. Our goal is to investigate
the cause of the rapid flattening of the gradient in just $\sim 50$ Myr by
examining the kinematics of the gas. We choose the line-of-sight velocity as a
metric because it is a readily accessible quantity to spectroscopic
observations.

By comparing the snapshots “A” and “B” it is clear that the flatter
metallicity gradient at the later time is due to the presence of metal-rich
clumps on the outskirts of the central disk in “B”. This is clearly not a
resolution-effect, although resolution can play a crucial role in such
measurements (see Section 4). In this case, the extent of our analysis—dense
gas within $10h^{-1}$ ckpc of the disk center—serendipitously includes
interacting systems which leads to high metallicity on the outskirts, thereby
flattening the overall radial profile222We have repeated our analysis by
considering a smaller disk, down to $3h^{-1}$ ckpc, and arrived at the same
qualitative conclusions.. It appears that these interacting systems are
impacting the line-of-sight velocity map too, making it less smooth and
organised compared to that in “A”. Therefore, disturbed kinematics could
potentially be correlated with rapid interactions and therefore flatter
gradients. Such interactions are more common at high-$z$ when the galaxies are
undergoing more frequent mergers and violent starbursts and only last for a
short time. Indeed, $\lesssim$100 Myr after “B” the gradient is strongly
negative again, as seen in the top panel of Figure 4.

Figure 4: Top: The bold purple line shows the full evolution of the
metallicity gradient of Hurricane—one of the FOGGIE galaxies—with a cadence of
$\approx$ 5 Myr. The thinner purple line in the background denotes the same
evolution with a lower cadence of 500 Myr. This thin line is a proxy for the
relatively low-stochasticity, average evolution of the gradient over time. The
shaded region of $\pm$ 0.03 dex/kpc around the lower-cadence evolution is
intended to depict the typical observational uncertainty in gradient
measurements. The purple circles denote whole number redshifts $z=$ 0, 1, 2, 3
and 4 (from right to left). This instantaneous gradient for this galaxy is
outside this typical 1$\sigma$ uncertainty 50% of the time up to $z>=1$,
denoted by the black dashed line.
Middle: From left to right, the panels denote the projected metallicity map,
line of sight velocity map, radial profile and full distribution of mass-
weighted metallicity of the time stamp corresponding to A, as denoted in the
top panel. Bottom: Same as above but corresponding to the time stamp B,
denoted in the top panel. The time stamps A & B are only $\sim 50$ Myr apart,
although they show a stark contrast ($\sim 0.25$ dex) in the gradients. This
demonstrates the short timescale variations in metallicity distribution,
leading to the large scatter in gradients, particularly at high-$z$.

###  3.2 Non-parametric characterization

We investigate a non-parametric approach to quantify the distribution of
metallicity within a certain distance from the center, which we describe in
Section 2.5. Our new approach focuses on the distribution of different
metallicity values present within the disk, rather than where in the disk the
different metallicities occur.

The top panel of Figure 5 shows the time evolution of the metallicity
gradient. This is the same as the top panel of Figure 4. We reproduce it here
again for ease of comparison against the time evolution of other quantities
derived from the non-parametric approach, as discussed in the following
paragraphs. In each panel, the circles represent whole number redshifts, i.e.,
$z=0$, 1, 2, 3 and 4 from right to left.

In the middle panel of Figure 5 we show the evolution of the median (brown)
and inter-quartile range (IQR, blue) of the full metallicity distribution (see
Figure 1). Both of these quantities are independent of geometric assumptions
or fitting the data to a model. Therefore, the median and the IQR represent a
non-parametric, quantitative description of the distribution of the mass-
weighted gas-phase metallicity. We do not see any significant overall
evolution from $z=4$ to $z=0$ in either the median or the IQR. However, both
exhibit significant scatter, particularly downward detours at several epochs
throughout the lifetime of the galaxy.

The bottom panel of Figure 5 shows the star formation rate of Hurricane as a
function of time (Lochhaas et al., 2021). We provide figures similar to Figure
5, but corresponding to the the other FOGGIE halos, in Appendix A.

Our goal is to investigate the changes in the metallicity distribution in
response to star formation. Comparing the middle and the bottom panel, it
appears that the width (IQR) of the metallicity distribution tends to decrease
following prominent starburst events, as highlighted in the figure. Rapid star
formation epochs are closely ($\sim$ few Myr) followed by supernovae events,
which leads to deposition of metals as well as momentum on to the gas. However
the momentum feedback expels metal-rich gas (from the vicinity of the
starburst) before it can mix with the ISM, leading to a decrease in median
metallicity (brown curve in the middle panel) and a narrower distribution
(blue curve). It takes the IQR $\sim 2$ Gyr to “recover” to its initial value
at the start of the starburst, hinting at the timescale for mixing of metals
throughout the disk.

It is worth noting that the radial gradient measurement (top panel) does not
exhibit any strong response to the star formation history. A slight flattening
(upward excursion) is occasionally seen in radial gradient due to mixing of
gas following starbursts, but other factors, such as the rapid stochasticity
in the gradient and interactions with nearby galaxies, wash out any reliable
signature in response to starbursts, particularly at high redshift.

Figure 5: Each panel shows the time evolution of various quantities of
Hurricane. Circles in each panel denote whole number redshifts $z=$ 0, 1, 2, 3
and 4 (from right to left). The first panel shows metallicity gradient, the
same as top panel of Figure 4. The second panel shows the median and inter-
quartile range (in log-space) of the full distribution of the metallicity
within 10 $h^{-1}$ ckpc, in brown and blue respectively. The bottom panel
denotes the star formation rate history. The highlighted regions in gray show
that the inter-quartile range of the metallicity distribution decreases in
response to star formation episodes.

##  4 What flattens or inverts the metallicity gradients?

In this section we elaborate on the ubiquity of positive and flat metallicity
gradients in observations, and their scarcity in simulations. We demonstrate
that although FOGGIE galaxies sometimes exhibit positive gradients, it is a
short-lived phenomenon—and does not occur sufficiently often to be able to
explain the observations. We emphasise that the inability to reproduce long-
lived positive gradients is common to all simulations, and therefore indicates
a major gap in our understanding of chemical evolution of galaxies.

###  4.1 Comparison with literature

Comparison with existing observational studies shows tensions in cases where
positive metallicity gradients have been observed. Observations of flat or
positive metallicity gradients, particularly at $1<z<3$, is a robust finding
given it has been reported by various studies using a multitude of techniques
(Wuyts et al., 2016; Wang et al., 2017, 2019, 2020; Carton et al., 2018; Curti
et al., 2020; Simons et al., 2021). Recently Wang et al. (2022) and Venturi et
al. (2024) pushed the redshift-boundary to $z\gtrsim 3$ and $z\sim 8$,
respectively, using JWST observations. Vallini et al. (2024) measured
metallicity gradients at $z\sim 7$ based on ALMA observations. All of these
recent measurements report flat/positive gradients all the way out to $z\sim
8$. This is something that theoretical models need to be able to explain.
Simulation-based studies (including this work) have attempted to do so, but
without any conclusive results yet.

The evolution of metallicity gradients of FOGGIE galaxies over long timescales
broadly agrees with existing simulations. Multiple theory-based
studies—including MUGS (Gibson et al., 2013), FIRE (Ma et al., 2017), and
Illustris (Hemler et al., 2021)—have addressed this problem using a variety of
different input physics. All these simulations have reported high-$z$
metallicity gradients to be usually negative, only sometimes flat, and almost
never positive. This is consistent with FOGGIE as well. This broad agreement,
in spite of the diversity among these simulation suites in numerical
treatments of hydrodynamics, input physics, and analysis methods is
significant. It might imply that the simulations are underpredicting the
timescale of certain processes, such as mixing of metal rich gas accreted on
the disk outskirts, leading to very short-lived positive gradients.
Alternatively, it is possible that poor numerical resolution leads to
inaccurate mixing of metals through numerical viscosity (Wadsley et al.,
2008).

While the instantaneous metallicity gradients of FOGGIE galaxies do sometimes
achieve positive values, these are the result of the highly stochastic
behavior of the gradient during the evolution of a given galaxy (Section 4.2).
The gradient quickly settles back to being negative on short timescales
($\sim$ 50 Myr). The FOGGIE galaxies do not develop a stable, long-lived
positive metallicity gradient at any point during their lifetime.

###  4.2 Positive gradients from interactions

Figure 6: Top: From left to right, the panels denote the projected metallicity
map, line of sight velocity map, radial profile and full distribution of mass-
weighted metallicity at a specific time stamp during the evolution of one of
the FOGGIE halos, Hurricane. These panels are shown at the native resolution
of the simulation. Bottom: Same as above but after binning the quantities at
$0.1^{\prime\prime}$/pixel resolution. The resolution comparison demonstrates
the detectability of such positive gradients with new observatories e.g.,
JWST.

Figure 6 depicts one of the rare instances of positive metallicity being seen
in the FOGGIE (Hurricane) simulations, and how it would appear at a lower
spatial resolution. The quantities depicted are the same as in Figure 4. While
the top row depicts the quantities at the native FOGGIE resolution ($\sim 77$
pc at $z\sim 2.55$), the bottom row shows the same quantities at 0.1′′
resolution ($\sim 0.8$ kpc at $z\sim 2.55$). We achieved the coarsening simply
by re-binning the intrinsic 2-D projected map of each relevant quantity,
without incorporating a full point spread function (PSF). This provides an
approximate expectation of how this simulated galaxy might appear to
JWST/NIRSpec.

The main takeaway from comparing the two rows of Figure 6 is that the
occurrence of positive gradients is independent of spatial resolution.
Examining the higher resolution counterparts—particularly the top left
panel—it is clear that the positive gradient is due to a close passage of a
metal-rich satellite on the top left of the central disk. In the lower
resolution version, however, it would be extremely challenging to distinguish
between the two separate disks in the 2-D metallicity map. Hence, data at such
a resolution would hint that the observed positive gradient is associated with
a single system, without accounting for any mergers/interactions. Turning to
2-D kinematic maps, even at low spatial resolution, might prove useful in
distinguishing such systems. As an example, let us consider just the low
resolution data, as shown in the bottom row of Figure 6. We see a positive
metallicity gradient based on a 2-D metallicity map that does not give us any
obvious hints about an interacting system. However, the departure from a
smooth line-of-sight velocity map (dark purple in the second panel of bottom
row) hints at disturbed kinematics, thus pointing us to the cause of the
measured positive gradient.

These minor interactions at high-$z$ occur intermittently and only last for a
short time ($\lesssim 100$ Myr), following the generic behavior of minor
mergers in $\Lambda$CDM structure formation. Although positive gradients are
rare in FOGGIE simulations ($<5\%$ of the time within $2<z<3$ and $<1\%$
throughout cosmic time), whenever they do occur, it is due to a close passage
of a merger, or interaction with a metal-rich companion. While it might be
difficult to identify interacting systems in metallicity-space at the
resolution of current space-based observations, the line-of-sight kinematics
might provide useful data.

Finding flat/positive metallicity gradients in association with mergers and
interactions is not unexpected. In some simulations, merger events can lead
initially to steepening of gradients, owing to delivery of less enriched gas
to the disk outskirts (Buck et al., 2023) with subsequent mixing later
settling the disk to flat gradients (Tissera et al., 2016). Indeed Tissera et
al. (2016) reported both positive as well as strongly negative gradients for
galaxies with disturbed morphology and/or undergoing interactions with a close
neighbour. The evolutionary stage of the merging event dictates whether the
gradient would be negative or positive.

Rosas-Guevara et al. (2022) showed that galaxies in the voids (as opposed to
close to the large scale filaments) in the EAGLE simulations are less likely
to harbour positive metallicity gradients. Our work with FOGGIE emphasises
that environment and interactions play a key role in flattening or inverting
the metallicity gradients.

##  5 What drives the general evolution of the metallicity gradients?

Now we discuss the general evolution of the FOGGIE metallicity gradients in
the context of other observational as well as simulation studies. In
particular, we emphasise the lack of strong correlations between gradients and
other global galaxy properties such as SFR, size, and circular velocity. We
highlight the challenges in understanding the physics driving the gradients if
we reduce each property to a single globally averaged value. We therefore
propose to utilise the entire information embedded in the full 2D maps of
metallicity and velocity, via a non-parametric approach.

###  5.1 Redshift evolution of gradients

In terms of tracking the gradients of individual simulated galaxies across
time, Ma et al. (2017) find opposite trends compared to this work (FOGGIE) as
well as Gibson et al. (2013) and Hemler et al. (2021). Figure 9 of Ma et al.
(2017), which shows the evolution of an individual halo of the FIRE
simulations, displays a rapid initial steepening of the metallicity gradient
early on, followed by little evolution in the past $\sim 6$ Gyr owing to a
“quiet” phase of evolution, while both FOGGIE (Figure 4) and MUGS (Figure 1 of
Gibson et al., 2013) display a general flattening with time. Hemler et al.
(2021) reported predominantly negative metallicity gradients, with significant
scatter, and an overall flattening of gradients at low-$z$ in TNG50
simulations. This contrast can not entirely be explained by contrasting merger
histories in the different simulations. While FOGGIE galaxies have been chosen
to have fewer major mergers at $z\lesssim 2$, Ma et al. (2017) also note the
naturally decreasing merger rates in FIRE simulations at low-$z$. The
similarity in merger history leads to the gradients not evolving significantly
in the last few billion years. Therefore, having steepened at earlier times
the FIRE gradients stay steep at low-$z$ due to lack of mergers, whereas
FOGGIE gradients generally stay weakly negative at low-$z$, having already
flattened by $z\sim 1$. Such flattening can be attributed to an accretion-
dominated phase at later times (see below). Hemler et al. (2021) note that
mergers flatten the gradients at high-$z$, whereas AGN feedback drives the
flattening of gradients at low-$z$ in their TNG50 simulations.

In the analytic models of Sharda et al. (2021a) the gradient is primarily
governed by four factors: (a) in-situ star formation, (b) accretion of
pristine gas, (c) advection of gas, and (d) diffusion of gas. The FOGGIE
gradient evolution qualitatively agrees with the Sharda et al. (2021a) models,
such that both works hint at a steepening of gradient with time, followed by a
“turnaround” towards flatter metallicity gradients at late times. Though we
have higher absolute metallicity values than predicted by Sharda et al.
(2021a), the predicted scatter in the gradients matches that of FOGGIE. They
attributed the scatter to the “yield reduction factor” $\phi_{y}$, which
determines the fraction of newly produced metals that are mixed with the ISM
rather than being ejected as galactic winds. However, their models hint that
the scatter in the gradient is low at high-$z$ and increases around $z\sim
0.5$ before decreasing again at lower-$z$. This is in contrast to what we
observe with FOGGIE. The scatter in FOGGIE is largest at high-$z$ when the
galaxy mass is lowest, and the scatter reduces with time and increasing
stellar mass. This discrepancy could be attributed to the lack of merger
events in the analytic prescription of Sharda et al. (2021a), because mergers
seem to be a key driver of the scatter and stochastic variations of the FOGGIE
gradients at high-$z$. Moreover, the Sharda et al. (2021a) model assumes the
ISM to be in steady-state equilibrium, which the FOGGIE simulations suggest,
via high stochasticity, is not really being achieved.

The short timescale variations in gradients enabled by the unique high time-
cadence of FOGGIE are seen in other cosmological simulations as well. Ma et
al. (2017) found that the scatter in gradient of a single galaxy in the FIRE
simulations is statistically the same as that in the full ensemble, indicating
that at high-$z$ gradients allow us just an instantaneous view rather than the
accumulated history of the galaxy evolution. They also find that bursty star
formation can lead to flatter gradients on short timescales.

Semi-analytic models of Fu et al. (2013) show that radial gas flows have
little effect on setting the gradient at late-times. Instead, the gradient is
set by the fraction of metals ejected into the hot halo before being mixed
with the cool ISM because it is this pre-enriched gas that accretes on to the
outskirts of the disk at later times, flattening the gradient at low-redshift
(Werk et al., 2011). Indeed, Simons et al. (2021) show with the help of a toy
model in their Figure 11, that in the absence of metal re-distribution, the
metallicity gradient of an isolated galaxy will rapidly ($\sim 10-100$ Myr
depending on mass) become strongly negative. Therefore, simulations must be
able to be analyzed on timescales significantly shorter than this in order to
capture the short time scale variations in metallicity distribution. With an
output cadence of $\sim 5$ Myr, our result demonstrates that in the presence
of a realistic environment and feedback to re-distribute metals, gradients
become flatter with time. This is in line with Gibson et al. (2013), who found
that enhanced feedback was important for reproducing “flat” metallicity
gradients. The feedback in the FOGGIE simulations is strong enough to launch
metal-rich outflows (Lochhaas et al., 2023), and we do see an overall
flattening of metallicity gradients.

###  5.2 Correlation with other galaxy properties

A significant observational finding is that the metallicity gradients show
little-to-no correlation with other global galaxy properties. Figure 3 of
Belfiore et al. (2017) compared the observed radial gradients in low-$z$
galaxies with galaxy stellar mass and found a correlation, particularly at
stellar masses below $\log{\rm M_{*}/M_{\odot}}<10$. However, at high-$z$
(Figure 7 of Simons et al., 2021) the correlation is weak. Moreover, no
meaningful correlation has been observed with other galaxy properties such as
SFR, sSFR, velocity dispersion or galaxy size (Figure 9 of Simons et al.,
2021). The lack of clues as to what might be governing the positive
metallicity gradients makes it a challenging task to ascertain its causes.
This is perhaps an indication that collapsing the 2-D spatial information of
metallicity distribution to one value—the metallicity gradient—has already
taught us all it can about galaxy evolution on its own, and going forward we
need to combine it with the diagnostic power of the full spatial distribution.

Simulation-based studies have reported a diverse degree of correlation between
metallicity gradient and stellar mass. Gradients in Illustris TNG50 galaxies
show little correlation with mass (for an ensemble of galaxies at various
cosmic times—not the individual galaxies tracked over cosmic time) in Figure 8
of Hemler et al. (2021), while Figure 6 of Ma et al. (2017) shows a weak
correlation with mass where lower mass galaxies exhibit flatter gradients.
This result of Ma et al. (2017) is perhaps the closest available simulation-
based counterpart to the high-redshift observations of (Simons et al., 2021,
Figure 9).

Both Hemler et al. (2021) and Ma et al. (2017) find weak-to-no correlation
between metallicity gradient ($\nabla Z$) and galaxy size, SFR or circular
velocity ($v_{c}$). Such poor correlations have been reported by observational
studies as well (Simons et al., 2021; Sharda et al., 2021b). A common link
between all such theoretical and observational studies is that they use
globally averaged quantities that reduce complexity of chemical evolution in
the galaxy disk. For instance, in a plot of $\nabla Z$ vs $v_{c}$ both
quantities have been derived by reducing the spatial information to a global
average, discarding a lot of potentially useful information along the way. We
show in Figure 4 that considering the 2-D map of metallicity and velocity can
reveal correlations (e.g., disturbed kinematics correlated with positive
gradients) that are lost when each galaxy is reduced to one value of $\nabla
Z$ and $v_{c}$.

###  5.3 Implications of a non-parametric approach

Given the difficulties with fitted symmetric gradients, we consider a non-
parametric approach to quantifying the full distribution of metallicity.
Unlike the radial gradient approach that assumes the existence of a well-
behaved disk, which may not hold at high-$z$, the non-parametric approach is
independent of geometric assumptions. Even if an observed galaxy’s center is
not well-constrained because of complex morphology, a non-parametric treatment
is still feasible in such a scenario because all one needs is a distribution
of observed metallicity values.

Moreover, spatially resolved metallicity studies at high-$z$ often involve
lensed galaxies and therefore employ lens model reconstruction. Many high-$z$
galaxies have been observed via gravitational lensing, tying in the lens model
uncertainties in the gradients (e.g., Jones et al., 2013a, 2015; Wang et al.,
2019, 2020). Our non-parametric approach is less susceptible to the systematic
uncertainties of lens models because the metallicity distribution is
independent of the physical geometry of the galactic disk. In the context of
an observational study of a lensed galaxy, accurately estimating the center of
the disk in the source-place is sensitive to the lensing model and associated
uncertainties. The metallicity distribution however, can be recovered without
an accurate estimate of the disk center, therefore being less sensitive to
lensing models.

The new approach, however, is not intended to replace the conventional
approach, given that there are scenarios that can be explored by only the
latter. For instance, radial transport of metals on short timescales will have
its imprint on the metallicity gradient (by changing the location of the
‘high’ and ‘low’ metallicity pockets along the radius) but will not affect the
full distribution of metallicity values. Therefore, the new non-parametric
approach may not be ideal for studying radial transport of metals. However,
the non-parametric approach can serve as a good “complement” to the
conventional gradient method, because of its ability to study metallicity
distribution of galaxies with disturbed morphology.

###  5.4 Breaks in radial metallicity profiles

We observe a “break” in the radial metallicity profile in the FOGGIE galaxies
throughout most of their lifetimes. The FOGGIE galaxies have a relatively
steep, high-metallicity radial profile in the central part, out to $\approx$
4–5 kpc, followed by a relatively flat metallicity profile. At times
(typically $z\gtrsim 0.5$), we observe a sharp drop in metallicity between the
inner steep profile and the outer flat profile (see the middle panel of Figure
1 as an example). This flat metallicity profile on the outskirts arises from a
cold, outer gas disk that is not forming stars. This is in agreement with
Garcia et al. (2023) who studied the evolution of the location of the break in
Illustris-TNG simulations and concluded that a typical disk in their
simulations is formed of two parts—an inner star formation dominated part with
a steep gradient, and an outer mixing dominated part with a flat metallicity
profile. Multiple observational studies have reported a similar break in the
radial gas-phase metallicity profile (e.g., Sánchez-Menguiano et al., 2018).
Although limited spatial resolution could be a potential contributor to the
“flattening” of observed metallicity profiles (e.g. Yuan et al., 2013b; Carton
et al., 2018; Acharyya et al., 2020), the “flattened” outskirts have recently
been observed, even with good spatial resolution. This calls for an
explanation of possible reasons behind the flat, low-metallicity outskirts
such as those observed in the FOGGIE simulations.

In FOGGIE galaxies the sharp drop in metallicity from the inner to the outer
region, when it occurs, is at least partially due to misalignment of the
angular momentum vector of the inner and outer gas disks (Simons et al. in
prep). Such a misalignment acts as a barrier to efficient mixing of gas across
the inner and outer disks and therefore hinders metal transport. The inner
disk continues to be metal enriched owing to local star formation whereas the
non-star-forming outer disk remains metal poor. With time, the inner and outer
disk orientations align, which eventually leads to one continuous, co-rotating
gas disk. Based on animations of gas density projections of the disks of
FOGGIE galaxies, this is a slow process and takes $\sim$ 1 Gyr. But once the
inner and outer disks are merged into a contiguous structure, it facilitates
metal mixing across the inter-disk boundary and makes the radial metallicity
profile smoother, and consequently, the gradient flatter.

###  5.5 Galaxy-to-galaxy scatter vs time evolution scatter

Individual FOGGIE galaxies exhibit significant scatter in the metallicity
gradient, particularly at high-$z$. Following Simons et al. (2021) we compared
the scatter in the gradient as a function of galaxy stellar mass. The
resulting relation was considerably noisy, though, overall we found the
scatter in the gradient to be decreasing with stellar mass, in agreement with
Simons et al. (2021).

The inter-galaxy variation of metallicity gradient between the six FOGGIE
galaxies is of the same order as the variation seen in a given halo’s gradient
during its lifetime. Note that the substantial diversity in simulated
gradients is not unique to FOGGIE; Ma et al. (2017) and Hemler et al. (2021)
predict large diversity in gradients as well. With such large scatter in the
metallicity gradient evolution, it is challenging to interpret the
‘instantaneous’ observations of gas phase metallicity gradients. This is
particularly true for high-$z$ observations with JWST because the
stochasticity is larger at high-$z$, when the ‘disk’ of the galaxy has not
settled. FOGGIE galaxies spend $\sim$30-50% of the time up to $z>1$ more than
the typical observational uncertainty ($\sigma\approx$0.03 dex/kpc) away from
the mean behavior (see Figure 4). As observations become more precise in the
future, and the typical observational uncertainty decreases, the fraction
described above will only increase. We conclude that it would be even less
likely for observations to capture the overall evolution in metallicity
gradient by measuring an instantaneous gradient.

##  6 Caveats and challenges

We now briefly highlight the several challenges in studying metallicity
gradient at high-$z$, both from an observational as well as a simulation
standpoint, and discuss the relevant caveats specific to the FOGGIE
simulations.

##### Uncertainties in metallicity diagnostics:

Simulations and observations measure metallicity gradients in inherently
different ways—one uses intrinsic values accessible in the simulations while
the other relies on accurately diagnosing the starlight processed by the
(often dusty) ISM via empirical and theoretical models. Therefore, direct
comparison between simulated and observed metallicity gradients should be
interpreted with caution. It is worth noting that unlike the simulations, the
observed gradients involve varied metallicity diagnostics, which in turn
depend on diverse photoionization models, leading to a scatter in the observed
gradients that may be larger than in the true population (Kewley & Ellison,
2008).

##### Challenges in defining the disk center:

Another reason to be cautious about interpreting fitted linear gradients is
that they refer to a center that is not known precisely. At high-$z$, the
morphology is more likely to be irregular and clumpy rather than smooth and
disk-like, making it difficult to reliably measure a disk center, and
consequently, a reliable radial gradient. In the simulations we have access to
the full dark matter distribution which can help find the center.
Observations, however, typically do not contain such information and are
therefore reliant on the surface brightness map to deduce the center, which is
affected by the chaotic morphology of the disk. Metallicity distributions that
routinely depart from clean linear gradients, and high time variability of the
gradient at high-$z$, warrant a search for a different approach to quantify
the spatial distribution of metallicity, as discussed in Section 3.2.

##### FOGGIE merger history:

The initial bursty star formation episodes followed by quieter phases of star
formation history, as seen in the FOGGIE galaxies, could potentially be a
selection effect. The FOGGIE galaxies are selected to not have major mergers
after $z\sim 2$, and therefore have fewer external stimuli for rejuvenating
star formation at lower redshifts. This works in favour of our current
analysis because this gives us an ideal laboratory for studying the evolution
of metallicity gradients both during bursty (high-$z$) as well as at quiet
(low-$z$) star formation episodes.

##### Star formation & supernova feedback in FOGGIE:

The supernova feedback prescription in current FOGGIE simulations is overall
underpowered, and therefore unable to expel enough metal-rich gas. Moreover,
the star formation scheme employed is slightly over- and under-efficient in
more and less dense regions respectively. For a detailed description of the
central concentration of star formation in FOGGIE, see Section 2.7 of Wright
et al. (2023). The combination of these effects leads to too many metals being
locked up in the ISM, resulting in too high absolute values of metallicity. We
consider that the centrally concentrated star formation may lead to a
systematic impact on the metallicity gradient. However, the gradients
reproduced in our simulations are consistent with that of other state-of-the-
art cosmological simulations, wherein we all fail to reproduce the
observations. Therefore the issues with the star formation prescription does
not impact our key results. The upcoming generation of FOGGIE simulations will
have improved feedback models. Although the absolute values of metallicities
in the current FOGGIE simulations may not be perfect, the current simulations
are sufficient for comparing relative metallicity quantities, such as the
gradient.

##### Numerical mixing:

Achieving accurate “mixing” of gas in simulations through numerical methods is
a challenging problem (Wadsley et al., 2008), and yet it impacts the metal
mixing and therefore, metallicity distributions. Particle-based simulations
often employ a sub-grid “mixing model”. This necessitates calibration with one
or more parameters and therefore introduces potential sources of
uncertainties. On the other hand, grid-based simulations, such as FOGGIE,
often suffer from the issue of over-mixing owing to the numerical recipes
used. However, the speed at which metals can diffuse from their point of
origin to larger radii is still limited by physical timescales and processes.
Therefore we believe that the mixing problem does not impact our main
conclusions.

##### Challenges of forward modeling:

A potential next step to better understand this discrepancy is to produce
synthetic Integral Field Spectroscopy (IFS) datacubes from cosmological
simulations, with realistic beam smearing and noise. Producing synthetic IFS
datacubes is generally a powerful tool, but requires a complete understanding
of the underlying ionization conditions that power the H II regions and
therefore govern the emitted light. In the absence of a detailed ionization
model, which is challenging to do at the resolutions achieved in our
cosmological simulations, this approach would require significant assumptions
regarding the sub-grid physics. Therefore, even after accounting for beam
smearing, noise, and metallicity diagnositcs, some systematic uncertainties
stemming from the internal ionization structure of H II regions, and
reliability of the diagnostics themselves, will continue to remain and are
very challenging to remedy entirely. The fact that no galaxy simulation is yet
able to reproduce the occurrence of observed positive and flat metallicity
gradients (discussed in detail in Section 3.1) hints that we might still be
missing something fundamental in our models. Therefore, we believe that the
forward modeling of existing inadequate simulations, compounded by sub-grid
physics assumptions, will not help us understand the discrepancy between
observed and simulated gradients. We therefore consider it of higher priority
to investigate what underlying physics in the simulations can reproduce the
observed positive gradients.

##  7 Summary & Conclusions

Spatially resolved gas-phase metallicity studies will imminently be
ubiquitous, especially now that JWST is performing beyond expectations. We
expect several forthcoming observational studies to target gravitationally
lensed high-redshift galaxies in order to perform spatially resolved studies
(e.g. Wang et al., 2022). As and when more such observational measurements
become available we will need physical models to help us interpret
observations.

We use the FOGGIE simulations—hydrodynamic, zoom-in cosmological simulations
with high spatial resolution ($\sim 270$ comoving pc in the ISM) and high
output cadence ($\sim 5$ Myr)—to study the evolution of the spatial
distribution of gas phase metallicity. Capitalizing on our higher cadence of
outputs, we show that gas-phase metallicity gradients of FOGGIE galaxies
display substantial variations that look stochastic and can exhibit big swings
over short time intervals. Consistent with other simulations, including MUGS
(Gibson et al., 2013), FIRE (Ma et al., 2017), and Illustris TNG50 (Hemler et
al., 2021), we find that gradients are usually negative, sometimes flat, and
rarely positive. The short-lived positive gradients in FOGGIE result from
minor interactions leading to a higher-metallicity patch on the outskirts of
the disk, not from a “real” positive gradient in the disk itself. We find that
non-parametric measures of metal distribution, and 2-D maps of velocity, are
more responsive to short-time evolution than simple fitted gradients. This
approach has other advantages for lensed galaxy studies (see Section 5.3).

Our main conclusions are as follows.

  * •

We corroborate the tension between observed and simulated metallicity
gradients in terms of the lack of positive gradients in simulations (Figure
3).

  * •

We find that the conventional radial metallicity gradient measurement exhibits
large ($\sim 1$1–2 dex) scatter on short timescales, particularly at high-
redshift, when the galaxy disk might not have completely settled. The
magnitude of the scatter is of the same order as the galaxy-to-galaxy
variation, thereby making it hard to disentangle time variability and inter-
galaxy variation as the driver of the scatter in global relations such as the
mass-metallicity relation (MZR) (Figure 4).

  * •

However, the overall gradient evolution of all FOGGIE halos exhibit a general
flattening over large timescales. This general behavior can be significantly
different from the instantaneous gradient at any given time.

  * •

In fact, FOGGIE halos spend $\sim 30$–50 % ($\sim 10$–25 %) of their lifetime
more than the typical observational uncertainty away from the mean behavior in
gradient evolution up to $z=1$ ($z=0$). Studies with lower cadence outputs
would likely miss such small scale variability and therefore be unable to
capture $\sim 30$–50 % of the gradient evolution at high-redshift (Figure 4).

  * •

We propose a non-parametric approach of quantifying the metallicity
distribution—by characterizing the full distribution of metallicities in the
galaxy disk, without any geometric assumptions. We demonstrate that the median
and the IQR of such a distribution is more stable at high-redshift, even when
the disk may not have fully formed yet. These metrics, particularly the IQR,
exhibit a strong response to star formation and stellar feedback (Figure 5).

  * •

We demonstrate that, although rarely occurring, positive gradients in the
FOGGIE simulations result from closely interacting or merging systems, which
can be challenging to distinguish at the spatial resolution of current space-
based observations (Figure 6).

Our work will pave the way for, and help interpret, upcoming JWST high-
redshift metallicity distribution studies.

## Acknowledgments

AA’s efforts for this work were supported by NSF-AST 1910414 and HST AR
#16151. RA, CL, AA, and MSP were supported in part by NASA via an Astrophysics
Theory Program grant 80NSSC18K1105. RA and CL also acknowledge financial
support from the STScI Director’s Discretionary Research Fund (DDRF). RA
acknowledges funding by the European Research Council through ERC-AdG SPECMAP-
CGM, GA 101020943. BWO acknowledges support from NSF grants #1908109 and
#2106575 and NASA ATP grants NNX15AP39G and 80NSSC18K1105. JT and ACW
acknowledge support from the Nancy Grace Roman Space Telescope Project, under
the Milky Way Science Investigation Team. RA’s efforts for this work were
additionally supported by HST AR #15012 and HST GO #16730. EHL acknowledges
support from the 2023 STScI Space Astronomy Summer Program. Computations
described in this work were performed using the publicly-available Enzo code
(http://enzo-project.org), which is the product of a collaborative effort of
many independent scientists from numerous institutions around the world. Their
commitment to open science has helped make this work possible. The python
packages matplotlib (Hunter, 2007), numpy (Walt et al., 2011), scipy (Virtanen
et al., 2020), yt (Turk et al., 2011), datashader (Bednar et al., 2022), and
Astropy (Astropy Collaboration et al., 2013, 2018, 2022) were all used in
parts of this analysis.

##  Appendix A Appendix

In the main body of the paper we focused our analysis on one FOGGIE
halo—Hurricane. In this section we present our results for the other five
FOGGIE halos—Tempest, Maelstrom, Squall, Blizzard and Cyclone. Our main
conclusions qualitatively hold true for all these halos.

In Figures A1-A6 we present snapshots of mass-weighted metallicity map, radial
profile and distribution at redshifts $z=$2, 1, and 0 for each halo. These
figures are similar to Figure 1 in the main text.

In Figures A7-A9 we present the time evolution of the metallicity gradient for
each halo, accompanied with the quantification for what fraction of time a
given galaxy spends significantly away from its general evolution (see Section
3.1.2). These figures are similar to the top panel of Figure 4 in the main
text.

In Figures A10-A14 we present evolution of the metallicity gradient as well as
the non-parametric quantities—median and IQR of the metallicity
distribution—in comparison with the star formation history of each halo. These
figures are similar to Figure 5 in the main text.

Figure A1: Similar to Figure 1, depicting the mass-weighted projected
metallicity (top row), radial metallicity profile (middle row) and metallicity
distribution (bottom row) of all gas cells of the halo Tempest, within 10
$h^{-1}$ comoving kpc (ckpc) of the center, and above a certain density
threshold (see Section 2.2). The columns denote three different redshifts:
$z=2$ (left), $z=1$ (middle) and $z=0$ (right). Limiting our analysis to
within 10 $h^{-1}$ ckpc results in a different physical extent of our analysis
at each epoch. This is demonstrated by (1) the scale bar in the top row
(depicting 1 physical kpc) being a different length in each panel, and (2) the
extent of the radial fit in the middle row being different in each panel.
Figure A2: Similar to Figure A1 but for another FOGGIE halo—Maelstrom. Figure
A3: Similar to Figure A1 but for another FOGGIE halo—Squall. Figure A4:
Similar to Figure A1 but for another FOGGIE halo—Blizzard. Figure A5: Similar
to Figure A1 but for another FOGGIE halo—Cyclone. This FOGGIE halo had been
simulated only up to $z=0.5$ at the time of this analysis. In this case the
rightmost column corresponds to $z=0.6$ (instead of $z=0$). Figure A6: Similar
to Figure A1 but for another FOGGIE halo—Hurricane.

Figure A7: Similar to Figure 4 but for the other FOGGIE halos—Tempest (top)
and Maelstrom (bottom).

Figure A8: Similar to Figure 4 but for the other FOGGIE halos—Squall (top) and
Blizzard (bottom). Figure A9: Similar to Figure 4 but for another FOGGIE
halo—Cyclone. This FOGGIE halo had been simulated only up to $z=0.5$ at the
time of this analysis. So the rightmost purple circle corresponds to $z=0.5$
instead of $z=0$. Figure A10: Similar to Figure 5 but for a different
halo—Tempest. Figure A11: Similar to Figure 5 but for a different
halo—Maelstrom. Figure A12: Similar to Figure 5 but for a different
halo—Squall. Figure A13: Similar to Figure 5 but for a different
halo—Blizzard. Figure A14: Similar to Figure 5 but for a different
halo—Cyclone. This FOGGIE halo had been simulated only up to $z=0.5$ at the
time of this analysis. So the rightmost circle in each panel corresponds to
$z=0.5$ instead of $z=0$.

## References

  * Acharyya et al. (2020) Acharyya, A., Krumholz, M. R., Federrath, C., et al. 2020, MNRAS, 495, 3819, doi: 10.1093/mnras/staa1100
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bednar et al. (2022) Bednar, J. A., Crail, J., Crist-Harif, J., et al. 2022, holoviz/datashader: Version 0.14.3, v0.14.3, Zenodo, doi: 10.5281/zenodo.7331952
  * Belfiore et al. (2017) Belfiore, F., Maiolino, R., Tremonti, C., et al. 2017, MNRAS, 469, 151, doi: 10.1093/mnras/stx789
  * Bresolin (2007) Bresolin, F. 2007, ApJ, 656, 186, doi: 10.1086/510380
  * Bresolin et al. (2002) Bresolin, F., Gieren, W., Kudritzki, R.-P., Pietrzyński, G., & Przybilla, N. 2002, ApJ, 567, 277, doi: 10.1086/338505
  * Brummel-Smith et al. (2019) Brummel-Smith, C., Bryan, G., Butsky, I., et al. 2019, The Journal of Open Source Software, 4, 1636, doi: 10.21105/joss.01636
  * Bryan et al. (2014) Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19, doi: 10.1088/0067-0049/211/2/19
  * Buck et al. (2023) Buck, T., Obreja, A., Ratcliffe, B., et al. 2023, MNRAS, 523, 1565, doi: 10.1093/mnras/stad1503
  * Bundy et al. (2015) Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7, doi: 10.1088/0004-637X/798/1/7
  * Carton et al. (2018) Carton, D., Brinchmann, J., Contini, T., et al. 2018, MNRAS, 478, 4293, doi: 10.1093/mnras/sty1343
  * Cecil et al. (2014) Cecil, G. N., Croom, S., & The SAMI Galaxy Survey Team. 2014, in American Astronomical Society Meeting Abstracts, Vol. 223, American Astronomical Society Meeting Abstracts #223, 246.03
  * Chiti et al. (2020) Chiti, A., Frebel, A., Jerjen, H., Kim, D., & Norris, J. E. 2020, ApJ, 891, 8, doi: 10.3847/1538-4357/ab6d72
  * Corlies et al. (2020) Corlies, L., Peeples, M. S., Tumlinson, J., et al. 2020, ApJ, 896, 125, doi: 10.3847/1538-4357/ab9310
  * Curti et al. (2020) Curti, M., Maiolino, R., Cirasuolo, M., et al. 2020, MNRAS, 492, 821, doi: 10.1093/mnras/stz3379
  * Durrell et al. (2001) Durrell, P. R., Harris, W. E., & Pritchet, C. J. 2001, AJ, 121, 2557, doi: 10.1086/320403
  * Förster Schreiber et al. (2018) Förster Schreiber, N. M., Renzini, A., Mancini, C., et al. 2018, ApJS, 238, 21, doi: 10.3847/1538-4365/aadd49
  * Fu et al. (2013) Fu, J., Kauffmann, G., Huang, M.-l., et al. 2013, MNRAS, 434, 1531, doi: 10.1093/mnras/stt1117
  * Garcia et al. (2023) Garcia, A. M., Torrey, P., Hemler, Z. S., et al. 2023, MNRAS, 519, 4716, doi: 10.1093/mnras/stac3749
  * Garnett et al. (1997) Garnett, D. R., Shields, G. A., Skillman, E. D., Sagan, S. P., & Dufour, R. J. 1997, ApJ, 489, 63, doi: 10.1086/304775
  * Gibson et al. (2013) Gibson, B. K., Pilkington, K., Brook, C. B., Stinson, G. S., & Bailin, J. 2013, A&A, 554, A47, doi: 10.1051/0004-6361/201321239
  * Hayden et al. (2015) Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132, doi: 10.1088/0004-637X/808/2/132
  * Hemler et al. (2021) Hemler, Z. S., Torrey, P., Qi, J., et al. 2021, MNRAS, 506, 3024, doi: 10.1093/mnras/stab1803
  * Hummels et al. (2019) Hummels, C. B., Smith, B. D., Hopkins, P. F., et al. 2019, ApJ, 882, 156, doi: 10.3847/1538-4357/ab378f
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jones et al. (2013a) Jones, T., Ellis, R. S., Richard, J., & Jullo, E. 2013a, ApJ, 765, 48, doi: 10.1088/0004-637X/765/1/48
  * Jones et al. (2015) Jones, T., Wang, X., Schmidt, K. B., et al. 2015, AJ, 149, 107, doi: 10.1088/0004-6256/149/3/107
  * Jones et al. (2013b) Jones, T. A., Ellis, R. S., Schenker, M. A., & Stark, D. P. 2013b, ApJ, 779, 52, doi: 10.1088/0004-637X/779/1/52
  * Kennicutt et al. (2003) Kennicutt, Jr., R. C., Bresolin, F., & Garnett, D. R. 2003, ApJ, 591, 801, doi: 10.1086/375398
  * Kewley & Ellison (2008) Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183, doi: 10.1086/587500
  * Kewley et al. (2010) Kewley, L. J., Rupke, D., Zahid, H. J., Geller, M. J., & Barton, E. J. 2010, ApJ, 721, L48, doi: 10.1088/2041-8205/721/1/L48
  * Kirby et al. (2011) Kirby, E. N., Lanfranchi, G. A., Simon, J. D., Cohen, J. G., & Guhathakurta, P. 2011, ApJ, 727, 78, doi: 10.1088/0004-637X/727/2/78
  * Krabbe et al. (2008) Krabbe, A. C., Pastoriza, M. G., Winge, C., Rodrigues, I., & Ferreiro, D. L. 2008, MNRAS, 389, 1593, doi: 10.1111/j.1365-2966.2008.13701.x
  * Leethochawalit et al. (2016) Leethochawalit, N., Jones, T. A., Ellis, R. S., Stark, D. P., & Zitrin, A. 2016, ApJ, 831, 152, doi: 10.3847/0004-637X/831/2/152
  * Lochhaas et al. (2021) Lochhaas, C., Tumlinson, J., O’Shea, B. W., et al. 2021, ApJ, 922, 121, doi: 10.3847/1538-4357/ac2496
  * Lochhaas et al. (2023) Lochhaas, C., Tumlinson, J., Peeples, M. S., et al. 2023, ApJ, 948, 43, doi: 10.3847/1538-4357/acbb06
  * Loebman et al. (2016) Loebman, S. R., Debattista, V. P., Nidever, D. L., et al. 2016, ApJ, 818, L6, doi: 10.3847/2041-8205/818/1/L6
  * Ma et al. (2017) Ma, X., Hopkins, P. F., Feldmann, R., et al. 2017, MNRAS, 466, 4780, doi: 10.1093/mnras/stx034
  * Marino et al. (2016) Marino, R. A., Gil de Paz, A., Sánchez, S. F., et al. 2016, A&A, 585, A47, doi: 10.1051/0004-6361/201526986
  * Mingozzi et al. (2020) Mingozzi, M., Belfiore, F., Cresci, G., et al. 2020, A&A, 636, A42, doi: 10.1051/0004-6361/201937203
  * Minniti & Zoccali (2008) Minniti, D., & Zoccali, M. 2008, in Formation and Evolution of Galaxy Bulges, ed. M. Bureau, E. Athanassoula, & B. Barbuy, Vol. 245, 323–332, doi: 10.1017/S1743921308018048
  * Miralles-Caballero et al. (2014) Miralles-Caballero, D., Díaz, A. I., Rosales-Ortega, F. F., Pérez-Montero, E., & Sánchez, S. F. 2014, MNRAS, 440, 2265, doi: 10.1093/mnras/stu435
  * Molina et al. (2017) Molina, J., Ibar, E., Swinbank, A. M., et al. 2017, MNRAS, 466, 892, doi: 10.1093/mnras/stw3120
  * Moustakas et al. (2010) Moustakas, J., Kennicutt, Robert C., J., Tremonti, C. A., et al. 2010, The Astrophysical Journal Supplement Series, 190, 233, doi: 10.1088/0067-0049/190/2/233
  * Muñoz-Elgueta et al. (2018) Muñoz-Elgueta, N., Torres-Flores, S., Amram, P., et al. 2018, MNRAS, 480, 3257, doi: 10.1093/mnras/sty1934
  * Peeples et al. (2019) Peeples, M. S., Corlies, L., Tumlinson, J., et al. 2019, ApJ, 873, 129, doi: 10.3847/1538-4357/ab0654
  * Poetrodjojo et al. (2018) Poetrodjojo, H., Groves, B., Kewley, L. J., et al. 2018, MNRAS, 479, 5235, doi: 10.1093/mnras/sty1782
  * Rich et al. (2012) Rich, J. A., Torrey, P., Kewley, L. J., Dopita, M. A., & Rupke, D. S. N. 2012, ApJ, 753, 5, doi: 10.1088/0004-637X/753/1/5
  * Rosas-Guevara et al. (2022) Rosas-Guevara, Y., Tissera, P., Lagos, C. d. P., Paillas, E., & Padilla, N. 2022, MNRAS, 517, 712, doi: 10.1093/mnras/stac2583
  * Sánchez et al. (2014) Sánchez, S. F., Rosales-Ortega, F. F., Iglesias-Páramo, J., et al. 2014, A&A, 563, A49, doi: 10.1051/0004-6361/201322343
  * Sánchez-Menguiano et al. (2016) Sánchez-Menguiano, L., Sánchez, S. F., Pérez, I., et al. 2016, A&A, 587, A70, doi: 10.1051/0004-6361/201527450
  * Sánchez-Menguiano et al. (2018) —. 2018, A&A, 609, A119, doi: 10.1051/0004-6361/201731486
  * Sharda et al. (2021a) Sharda, P., Krumholz, M. R., Wisnioski, E., et al. 2021a, MNRAS, 502, 5935, doi: 10.1093/mnras/stab252
  * Sharda et al. (2021b) Sharda, P., Wisnioski, E., Krumholz, M. R., & Federrath, C. 2021b, MNRAS, 506, 1295, doi: 10.1093/mnras/stab1836
  * Simons et al. (2020) Simons, R. C., Peeples, M. S., Tumlinson, J., et al. 2020, ApJ, 905, 167, doi: 10.3847/1538-4357/abc5b8
  * Simons et al. (2021) Simons, R. C., Papovich, C., Momcheva, I., et al. 2021, ApJ, 923, 203, doi: 10.3847/1538-4357/ac28f4
  * Swinbank et al. (2012) Swinbank, A. M., Smail, I., Sobral, D., et al. 2012, ApJ, 760, 130, doi: 10.1088/0004-637X/760/2/130
  * Tissera et al. (2016) Tissera, P. B., Pedrosa, S. E., Sillero, E., & Vilchez, J. M. 2016, MNRAS, 456, 2982, doi: 10.1093/mnras/stv2736
  * Turk et al. (2011) Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, The Astrophysical Journal Supplement Series, 192, 9, doi: 10.1088/0067-0049/192/1/9
  * Vallini et al. (2024) Vallini, L., Witstok, J., Sommovigo, L., et al. 2024, MNRAS, 527, 10, doi: 10.1093/mnras/stad3150
  * van de Voort et al. (2019) van de Voort, F., Springel, V., Mandelker, N., van den Bosch, F. C., & Pakmor, R. 2019, MNRAS, 482, L85, doi: 10.1093/mnrasl/sly190
  * Venturi et al. (2024) Venturi, G., Carniani, S., Parlanti, E., et al. 2024, arXiv e-prints, arXiv:2403.03977.  https://arxiv.org/abs/2403.03977
  * Vila-Costas & Edmunds (1992) Vila-Costas, M. B., & Edmunds, M. G. 1992, MNRAS, 259, 121, doi: 10.1093/mnras/259.1.121
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vogt et al. (2015) Vogt, F. P. A., Dopita, M. A., Borthakur, S., et al. 2015, MNRAS, 450, 2593, doi: 10.1093/mnras/stv749
  * Wadsley et al. (2008) Wadsley, J. W., Veeravalli, G., & Couchman, H. M. P. 2008, MNRAS, 387, 427, doi: 10.1111/j.1365-2966.2008.13260.x
  * Walt et al. (2011) Walt, S. v. d., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Wang et al. (2017) Wang, X., Jones, T. A., Treu, T., et al. 2017, ApJ, 837, 89, doi: 10.3847/1538-4357/aa603c
  * Wang et al. (2019) —. 2019, ApJ, 882, 94, doi: 10.3847/1538-4357/ab3861
  * Wang et al. (2020) Wang, X., Jones, T. A., Treu, T., et al. 2020, The Astrophysical Journal, 900, 183, doi: 10.3847/1538-4357/abacce
  * Wang et al. (2022) Wang, X., Jones, T., Vulcani, B., et al. 2022, ApJ, 938, L16, doi: 10.3847/2041-8213/ac959e
  * Werk et al. (2011) Werk, J. K., Putman, M. E., Meurer, G. R., & Santiago-Figueroa, N. 2011, ApJ, 735, 71, doi: 10.1088/0004-637X/735/2/71
  * Wright et al. (2023) Wright, A. C., Tumlinson, J., Peeples, M. S., et al. 2023, arXiv e-prints, arXiv:2309.10039, doi: 10.48550/arXiv.2309.10039
  * Wuyts et al. (2016) Wuyts, E., Wisnioski, E., Fossati, M., et al. 2016, ApJ, 827, 74, doi: 10.3847/0004-637X/827/1/74
  * Yuan et al. (2013a) Yuan, T.-T., Kewley, L. J., & Rich, J. 2013a, ApJ, 767, 106, doi: 10.1088/0004-637X/767/2/106
  * Yuan et al. (2013b) Yuan, T.-T., Kewley, L. J., & Richard, J. 2013b, ApJ, 763, 9, doi: 10.1088/0004-637X/763/1/9
  * Zheng et al. (2020) Zheng, Y., Peeples, M. S., O’Shea, B. W., et al. 2020, ApJ, 896, 143, doi: 10.3847/1538-4357/ab960a