Figuring Out Gas and Galaxies In Enzo (FOGGIE). VIII. Complex and Stochastic Metallicity Gradients at z > 2

We use cosmological zoom-in simulations from the Figuring Out Gas and Galaxies In Enzo (FOGGIE) project¹¹ (M. S. Peeples et al. 2019; L. Corlies et al. 2020; R. C. Simons et al. 2020; Y. Zheng et al. 2020; C. Lochhaas et al. 2021, 2023; A. C. Wright et al. 2024). The simulations have been described in detail in previous FOGGIE papers; however, we summarize them briefly here.

First, a large cosmological box of 100h⁻¹ comoving Mpc was run with Enzo at relatively low resolution to provide halo catalogs from which to select for zooms. A flat cold dark matter (ΛCDM) cosmology (Planck Collaboration et al. 2014) was used for these simulations, with 1 − Ω_λ = Ω_m = 0.285, Ω_b = 0.0461, and h = 0.695. Six halos were picked such that they are Milky Way mass at z = 0. Then, these six halos were re-simulated with Enzo (G. L. Bryan et al. 2014; last described in C. 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, 200h⁻¹ 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 (R. C. Simons et al. 2020), the current production FOGGIE simulations impose a minimum level of refinement n_ref = 9 (1100 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_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 (≳​​​​​​95%) are resolved at the highest level, i.e., 91 pc at z = 2 (R. C. Simons et al. 2020), and the median cell mass in the CGM and disk–halo interface is <200 M_⊙ (C. 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 C. B. Hummels et al. 2019; V. van de Voort 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 ≳ 2) and without (z ≲ 2) major mergers. We refer the reader to Table 1 for a detailed list of global properties of all of 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, and 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 of 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 the Appendix that our main conclusions remain qualitatively unchanged for the other halos as well.

Observed metallicity gradients are limited to the gaseous, star-forming regions of galaxies and their immediate environment. We therefore limit our analysis to within 10h⁻¹ comoving kpc of the galaxy center. We choose this extent because it is equivalent to ​​​​​​≈2.5 × r_e, where r_e ∼ 4 kpc is a typical effective radius of the gas disk of star-forming galaxies at z ∼ 0. Additionally, we impose a simple redshift-dependent density criterion such that ρ_cut = 2 × 10⁻²⁶ g cm⁻³ if z > 0.5, ρ_cut = 2 × 10⁻²⁷ g cm⁻³ if z < 0.25, and a linear interpolation with respect to physical time (not redshift) between 0.25 ≤ z ≤ 0.5. Note that this density criterion is based on proper density, and not comoving density. We select only the material with density above ρ_cut, and apply our analysis methods on this dense material. The choice of this density criterion was used by C. 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 (C. Lochhaas et al. 2023; R. C. Simons et al. 2024).

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 of the panels correspond to the z = 0.07 snapshot of the Hurricane halo, and show metallicities within a 10h⁻¹ 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 other FOGGIE halos and several redshift epochs, in Figures A1–A6 in the Appendix.

Zoom In Zoom Out Reset image size

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 2D 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 kiloparsecs 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⁻¹ 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 (D. Minniti & M. 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.

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 3D distribution of the physical property concerned. Therefore, it is important to compare the full 3D 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 3D 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 3D metallicity profile. The top row of Figure 2 corresponds to the projected 2D 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 3D metallicity distribution, the 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 3D distribution all agree with each other within 1σ uncertainties. Therefore, we conclude the projection effects do not impact the main science results of this paper. Henceforth, we proceed with the full 3D metal distribution for the rest of the paper, disregarding projection effects.

Zoom In Zoom Out Reset image size

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.

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 2D 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 nonparametric 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 has already been utilized by the stellar metallicity community to great scientific success, ranging from the Milky Way (M. R. Hayden et al. 2015; S. R. Loebman et al. 2016) to stellar halos (P. R. Durrell et al. 2001) and satellites (E. N. Kirby et al. 2011; A. Chiti et al. 2020).

We characterize the spatial distribution of metallicity within 10h⁻¹ 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 interquartile range (Z_IQR) as the difference between the 75th (Z₇₅) and 25th (Z₂₅) percentiles, and use it, along with the median, to trace the evolution of the metallicity distribution across cosmic time (Section 3.2). All of the percentiles are in $\mathrm{log}({Z}_{\odot })$ units, and therefore, Z_IQR is the difference between Z₇₅ and Z₂₅ 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 ≲ 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 multipeaked distributions with skewed Gaussian models have proven to be nontrivial, 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 3D distribution of metals in the galaxy disk by parameterizing it with simple 1D models.

In Table 1 we tabulate the fitted and other measured quantities for all FOGGIE halos, including the total stellar mass and total metallicity. The stellar mass is defined as the total mass in stars within 10h⁻¹ ckpc. The total metallicity is the ratio of metal mass to gas mass within 10h⁻¹ ckpc, after applying the density criteria to select the disk (see Section 2.2). We then normalize it using solar metallicity of 0.01295, which is the default solar metallicity in yt (M. J. Turk et al. 2011), and is consistent with the solar metallicity used in Cloudy photoionization models. In Table 1 we report this normalized total metallicity in log-scale. The total metallicities of FOGGIE galaxies are higher than those observed for galaxies of similar mass and type, particularly at high redshift. FOGGIE galaxies exhibit supersolar metallicity ∼50% of the time. This is due to a combination of the current feedback and star formation prescription in FOGGIE, which leads to centrally concentrated star formation and too many metals being locked up in the stars. For a detailed explanation, see A. C. Wright et al. (2024). Although the absolute metallicity values are somewhat higher than observed, the relative variation in metallicity, i.e., metallicity gradients, is robust and therefore renders our primary science conclusions unaffected.

Table 1. Measured Quantities within 10h⁻¹ ckpc for Each Snapshot of Each FOGGIE Halo: Stellar Mass, Total Metallicity, Median and Interquartile Range (IQR) of the Metallicity Distribution, and the Fitted Radial Gradient (See Section 2.3)

Halo	z	Time	M⋆/M⊙	$\mathrm{log}{Z}_{{\rm{total}}}$/Z⊙	$\mathrm{log}{Z}_{{\rm{median}}}$/Z⊙	$\mathrm{log}{Z}_{{\rm{IQR}}}$/Z⊙	∇ Z
		(Gyr)					(dex/kpc)
Tempest	0.0	13.62	10.72	−0.13	−0.18	0.36	−​​​​​​0.064 ± 0.000
Tempest	1.0	5.89	10.56	0.17	0.17	0.12	−0.096 ± 0.000
Tempest	2.0	3.32	10.19	−​​​​​​0.10	−0.31	0.41	−0.699 ± 0.001
Maelstrom	0.0	13.62	11.02	0.13	0.12	0.20	−0.032 ± 0.000
Maelstrom	1.0	5.89	10.79	−0.05	−0.09	0.47	−0.143 ± 0.000
Maelstrom	2.0	3.32	10.48	−0.20	−0.42	0.61	−0.153 ± 0.000
Squall	0.0	13.62	11.07	−0.08	−0.38	0.64	−0.110 ± 0.000
Squall	1.0	5.89	10.45	0.09	0.09	0.21	−0.080 ± 0.000
Squall	2.0	3.32	9.18	0.07	0.00	0.18	−0.241 ± 0.000
Blizzard	0.0	13.62	11.14	−0.13	−0.20	0.42	−0.068 ± 0.000
Blizzard	1.0	5.89	10.91	0.16	0.13	0.23	−0.059 ± 0.000
Cyclone	1.0	5.90	11.12	−0.04	−0.09	0.09	−0.036 ± 0.000
Cyclone	2.0	3.32	10.49	0.16	0.19	0.26	−0.220 ± 0.000
Hurricane	0.0	13.62	11.37	−0.20	−0.45	0.54	−0.080 ± 0.000
Hurricane	1.0	5.89	11.04	0.10	0.06	0.28	−0.124 ± 0.000
Hurricane	2.0	3.32	10.42	0.08	0.08	0.65	−0.301 ± 0.001

Note. All of the quantities are quoted in log₁₀ scale, with the exception of redshift, time (in gigayears) and metallicity gradient (in dex kpc⁻¹).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  Machine-readable (MRT)Typeset image

In the table presented here, we only include a few snapshots, intended to serve as examples. The full table is available in machine-readable format, while the original copy can be found at: https://ayanacharyya.github.io/data/master_table_Zpaper_upto10.0ckpchinv_wtby_mass_fitmultiple_wdencut_islog.txt.

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 X. 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.

Zoom In Zoom Out Reset image size

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 (∼−0.5 to −0.1 dex kpc⁻¹) and evolve toward shallower negative gradients (∼−0.3 to −0.1 dex kpc⁻¹). The bottom panel of Figure 3 depicts the spread in metallicity gradients predicted by different simulations—MUGS (B. K. Gibson et al. 2013), FIRE-1 (X. Ma et al. 2017), FIRE-2 (M. A. Bellardini et al. 2022), and Illustris TNG (Z. S. Hemler et al. 2021) in comparison to that of FOGGIE (see also Figure 3 of X. 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 ∼5.3 Myr, which makes the rapid stochastic fluctuations of the gradients readily apparent. FIRE-2 (M. A. Bellardini et al. 2021, 2022; R. L. Graf et al. 2024) simulations have an output cadence of 20–25 Myr, making these the most suitable comparison sample to FOGGIE (detailed in Section 4.1). However, both these simulations contrast with the 50–100 Myr cadence of most previous simulations (e.g., B. K. Gibson et al. 2013; X. Ma et al. 2017; Z. S. 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 gigayears 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., T.-T. Yuan et al. 2013a; A. Acharyya et al. 2020), none of these effects have yet been demonstrated to bias the observed gradients toward 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 ∼​​​​​​0.1–0.4 dex but there also is considerable stochasticity (variations of ∼0.1–0.2 dex in slope over a ​​​​​​∼few megayears) 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 ∼​​​​​​500 Myr, whereas the thicker line in the foreground represents the intrinsic cadence (∼5 Myr) of FOGGIE simulations. We see an overall flattening of the metallicity gradient after z ∼ 1. We provide figures similar to the top panel of Figure 4, but corresponding to the other FOGGIE halos, in the Appendix.

Zoom In Zoom Out Reset image size

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 ≥ 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 ±0.03 dex kpc⁻¹ around the lower-cadence evolution. We emphasize 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 ∼50% of the time at z > 1. In other words, a cadence of ∼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

We investigate a nonparametric 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 nonparametric 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.

Zoom In Zoom Out Reset image size

In the middle panel of Figure 5 we show the evolution of the median (brown) and interquartile 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 nonparametric, 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 the IQR. The median metallicity shows a slight decline over this redshift range, discussed in detail in Section 3.2.1. 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 (C. Lochhaas et al. 2021). We provide figures similar to Figure 5, but corresponding to the other FOGGIE halos, in the Appendix.

3.2.1. How Does Star Formation Shape the Metallicity Distribution?

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 (E. Wuyts et al. 2016; X. Wang et al. 2017, 2019, 2020; D. Carton et al. 2018; M. Curti et al. 2020; R. C. Simons et al. 2021; Z. Li et al. 2022). Recently X. Wang et al. (2022) and G. Venturi et al. (2024) pushed the redshift-boundary to z ≳ 3 and z ∼ 8, respectively, using JWST observations. L. Vallini et al. (2024) measured metallicity gradients at z ∼ 7 based on ALMA observations. All of these recent measurements report flat/positive gradients all the way out to z ∼ 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. Note that directly comparing observed and simulated gradients implicitly assumes a flat mass-to-light ratio profile, which is a caveat we discuss further in Section 6.

The evolution of metallicity gradients of FOGGIE galaxies over long timescales broadly agrees with existing simulations. Multiple theory-based studies—including MUGS (B. K. Gibson et al. 2013), FIRE-1 (X. Ma et al. 2017), FIRE-2 (M. A. Bellardini et al. 2021, 2022; R. L. Graf et al. 2024), and Illustris (Z. S. Hemler et al. 2021)—have addressed this problem using a variety of different input physics. All of these simulations have reported high-z metallicity gradients to be usually negative, only sometimes flat, and almost never positive. Moreover, there is a clear broad agreement between all of these simulations at low redshift. 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 (J. W. 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 (∼50 Myr). The FOGGIE galaxies do not develop a stable, long-lived positive metallicity gradient at any point during their lifetime.

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 (∼77 pc at z ∼ 2.55), the bottom row shows the same quantities at 0$\mathop{.}\limits^{\unicode{x02033}}$1 resolution (∼0.8 kpc at z ∼ 2.55). We achieved the coarsening simply by re-binning the intrinsic 2D projected map of each relevant quantity, without incorporating a full point-spread function. This provides an approximate expectation of how this simulated galaxy might appear to JWST/NIRSpec.

Zoom In Zoom Out Reset image size

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 2D 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 2D 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 2D 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 (≲100 Myr), following the generic behavior of minor mergers in Λ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 (T. Buck et al. 2023) with subsequent mixing later settling the disk to flat gradients (P. B. Tissera et al. 2016). Indeed P. B. Tissera et al. (2016) and P. B. Tissera et al. (2022) reported both positive as well as strongly negative gradients for galaxies with disturbed morphology and/or undergoing interactions with a close neighbor. The evolutionary stage of the merging event dictates whether the gradient would be negative or positive.

Y. 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 harbor positive metallicity gradients. Our work with FOGGIE emphasizes that environment and interactions play a key role in flattening or inverting the metallicity gradients.

In terms of tracking the gradients of individual simulated galaxies across time, X. Ma et al. (2017) found opposite trends compared to this work (FOGGIE) as well as B. K. Gibson et al. (2013) and Z. S. Hemler et al. (2021). Figure 9 of X. Ma et al. (2017), which shows the evolution of an individual halo of the FIRE-1 simulations, displays a rapid initial steepening of the metallicity gradient early on, followed by little evolution in the past ∼6 Gyr owing to a "quiet" phase of evolution, while both FOGGIE (Figure 4) and MUGS (Figure 1 of B. K. Gibson et al. 2013) display a general flattening with time. Z. S. 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 cannot entirely be explained by contrasting merger histories in the different simulations. While FOGGIE galaxies have been chosen to have fewer major mergers at z ≲ 2, X. Ma et al. (2017) also noted 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 ∼ 1. Such flattening can be attributed to an accretion-dominated phase at later times (see below). Z. S. Hemler et al. (2021) noted that mergers flatten the gradients at high-z, whereas AGN feedback drives the flattening of gradients at low-z in their TNG50 simulations. However, in the improved FIRE-2 simulations M. A. Bellardini et al. (2021; see also M. A. Bellardini et al. 2022; R. L. Graf et al. 2024) found a systematic steepening of the gradient with time, which is at odds with the behavior of FOGGIE galaxies, although they do agree at low redshift. The steepening in FIRE-2 is due to the disk settling down at lower redshift, and consequently the ISM becoming less turbulent, hindering the radial mixing of metals. This qualitative difference between FIRE-2 and FOGGIE is likely due to different numerical techniques (AMR in FOGGIE versus Lagrangian in FIRE) and feedback prescriptions.

In the analytic models of P. 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 P. Sharda et al. (2021a) models, such that both works hint at an initial steepening of gradient with time, followed by a "turnaround" toward flatter metallicity gradients at late times. This "turnaround" is attributed to accretion of pristine gas as the galaxies become more massive, and consequent mixing of this pristine gas within the galaxy. We note that without this accretion-dominated phase, the inside-out quenching and inside-out disk assembly that is typical in MW-type galaxies (e.g., S. Tacchella et al. 2016; V. Avila-Reese et al. 2018) would lead to a build-up for negative metallicity gradient in the galaxy. R. C. Simons et al. (2021) confirmed this behavior with the help of closed-box, toy models (see below).

Although we have higher absolute metallicity values than predicted by P. Sharda et al. (2021a), the predicted scatter in the gradients matches that of FOGGIE. They attributed the scatter to the "yield reduction factor" ϕ_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 ∼ 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 P. 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 P. 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. X. 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.

Semianalytic models of J. 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 onto the outskirts of the disk at later times, flattening the gradient at low redshift (J. K. Werk et al. 2011). Indeed, R. C. Simons et al. (2021) showed 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 (∼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 timescale variations in metallicity distribution. With an output cadence of ∼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 B. K. 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 (C. Lochhaas et al. 2023), and we do see an overall flattening of metallicity gradients.

A significant observational finding is that the metallicity gradients show little-to-no correlation with other global galaxy properties. Figure 3 of F. 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 $\mathrm{log}{M}_{* }/{M}_{\odot }\lt 10$. However, at high-z (Figure 7 of R. C. 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 R. C. 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 2D 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 Z. S. Hemler et al. (2021), while Figure 6 of X. Ma et al. (2017) shows a weak correlation with mass where lower mass galaxies exhibit flatter gradients. This result of X. Ma et al. (2017) is perhaps the closest available simulation-based counterpart to the high-redshift observations of R. C. Simons et al. (2021, Figure 9).

X. Ma et al. (2017) and M. A. Bellardini et al. (2022) reported an anticorrelation between radial gradient and the relative strength of circular velocity to velocity dispersion (v_c/v_σ), such that more rotationally supported galaxies tend to have steeper gradients, due to inefficient radial mixing. Z. S. Hemler et al. (2021), on the other hand, found weak-to-no correlation between metallicity gradient (∇Z) and galaxy size, SFR, or v_c. Such poor correlations have been reported by observational studies as well (P. Sharda et al. 2021b; R. C. Simons et al. 2021). 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 ∇Z versus 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 2D 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 ∇Z and v_c.

Given the difficulties with fitted symmetric gradients, we consider a nonparametric 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 nonparametric approach is independent of geometric assumptions. Even if an observed galaxy's center is not well constrained because of complex morphology, a nonparametric 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., T. Jones et al. 2013, 2015; X. Wang et al. 2019, 2020). Our nonparametric 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 nonparametric approach may not be ideal for studying radial transport of metals. However, the nonparametric 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.

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 ≈​​​​​​4–5 kpc, followed by a relatively flat metallicity profile. At times (typically z ≳ 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 A. M. 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., L. Sánchez-Menguiano et al. 2018). Although limited spatial resolution could be a potential contributor to the "flattening" of observed metallicity profiles (e.g., T.-T. Yuan et al. 2013b; D. Carton et al. 2018; A. 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 (R. C. Simons et al. 2024). 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, corotating gas disk. Based on animations of gas density projections of the disks of FOGGIE galaxies, this is a slow process and takes ∼​​​​​​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.

Individual FOGGIE galaxies exhibit significant scatter in the metallicity gradient, particularly at high-z. Following R. C. 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 R. C. 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; X. Ma et al. (2017) and Z. S. Hemler et al. (2021) predicted 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 ∼30%–50% of the time up to z > 1 more than the typical observational uncertainty (σ ≈ 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.

The scatter in metallicity gradients is highly dependent on the method of measuring the gradient, the metallicity diagnostics used (L. J. Kewley & S. L. Ellison 2008; H. Poetrodjojo et al. 2021), and whether the gradients have been normalized to a scale-radius. Indeed, L. Sánchez-Menguiano et al. (2018) demonstrated that upon normalizing the gradient to the effective radii, local galaxies have similar gradients, i.e., the scatter reduces. We have attempted to compute normalized gradients for the FOGGIE galaxies, by using the stellar half-mass–radius as the effective radius. This led the absolute gradient values in dex/r_e units to become steeper than the pre-normalized values in dex kpc⁻¹ units (as expected, given the multiplication with r_e). However, the overall galaxy-to-galaxy scatter or the inter-galaxy scatter did not change significantly upon normalization.