Interplay of Stellar and Gas-Phase Metallicities: Unveiling Insights for Stellar Feedback Modeling with Illustris, IllustrisTNG, and EAGLE

The metal content of galaxies provides a window into their formation in the full context of the cosmic baryon cycle. In this study, we examine the relationship between stellar mass and stellar metallicity (${\rm MZ}_*{\rm R}$) in the hydrodynamic simulations Illustris, TNG, and EAGLE to understand the global properties of stellar metallicities within the feedback paradigm employed by these simulations. Interestingly, we observe significant variations in the overall normalization and redshift evolution of the ${\rm MZ}_*{\rm R}$ across the three simulations. However, all simulations consistently demonstrate a tertiary dependence on the specific star formation rate (sSFR) of galaxies. This finding parallels the relationship seen in both simulations and observations between stellar mass, gas-phase metallicity, and some proxy of galaxy gas content (e.g., SFR, gas fraction, atomic gas mass). Since we find this correlation exists in all three simulations, each employing a sub-grid treatment of the dense, star-forming interstellar medium (ISM) to simulate smooth stellar feedback, we interpret this result as a fairly general feature of simulations of this kind. Furthermore, with a toy analytic model, we propose that the tertiary correlation in the stellar component is sensitive to the extent of the ``burstiness'' of feedback within galaxies.


INTRODUCTION
Heavy elements are valuable tracers of the gas flows that are an integral part of galaxy formation.As aging stellar populations return metals to the interstellar medium (ISM), these metals are then dispersed and advected with any subsequent gas flows, and are incorporated into future generations of stars (Kobayashi et al. 2020).Thus, the metal content of galaxies is a balance that is sensitive not only to the amount of heavy elements produced, but also on the efficacy of gas mixing within the ISM (e.g., Elmegreen 1999;Veilleux et al. 2005), prevalence of outflows (e.g., Veilleux et al. 2020), and rate of fresh/pristine gas accretion (Kereš et al. 2005).Since the metals are intrinsically tied to these physical processes within a galaxy, they may act as an observable diagnostic for the physics driving galaxy evolution (Dalcanton 2007;Kewley et al. 2019).For example, since outflows depend strongly on stellar feedback within a galaxy and out-★ E-mail: alexgarcia@virginia.edu † ARC DECRA Fellow flows have a distinct impact on galactic metallicity, the metal content of galaxies can thus be used as a probe of stellar feedback.
On large scales, trends between physical properties and metallicity within galaxy populations have been observed.For example, the stellar mass-gas-phase metallicity relation (MZ g R) describes a trend of increasing metallicity with increasing galaxy mass (Tremonti et al. 2004; Lee et al. 2006) and is regarded as a key scaling relation between the stellar masses of galaxies and the metallicity of the gas component.The reason that high mass galaxies correspond to elevated metallicities has been postulated to be a result of the high mass galaxies having either higher effective yields (i.e.retaining a higher fraction of their produced metals in the gas-phase; Tremonti et al. 2004), lower gas fractions (which results in higher metallicities at a fixed amount of metal mass; Pasquali et al. 2012), or perhaps both.The slope of the MZ g R has been noted to change as a function of mass with a steep power-law relation at low and intermediate masses that flattens at the highest masses (e.g., Blanc et al. 2019).These changes in slope are attributed to the efficiency (or lack thereof) of a number of different physical processes, including (but not limited to) stellar feedback, ISM turbulence, and the mixing of enriched and pristine gas in the circumgalactic medium (CGM).Moreover, there is a marked evolution in the MZ g R with time, whereby higher redshift galaxies form a clear MZ g R but with a lower overall normalisation (e.g., Savaglio et al. 2005;Maiolino et al. 2008;Zahid et al. 2011;Langeroodi et al. 2022).The redshift evolution is attributed to the lack of chemical maturity at high redshift wherein fewer generations of stars have formed to synthesize heavier elements.The MZ g R itself helps to reveal the net trends that must be present in terms of the gas mass evolution -coupled to the outflow properties -with galaxy mass and redshift for the metal content of galaxies to appropriately evolve.
In addition to the MZ g R, there is also a stellar mass-stellar metallicity relation (MZ * R) whereby the stellar metallicity increases with increasing stellar mass much in the same way as the MZ g R (e.g.Gallazzi et al. 2005;Panter et al. 2008).This relationship, much like its gas-phase counterpart, follows a trend of increasing metals with increasing stellar mass before flattening at high masses (e.g., Sommariva et al. 2012;Zahid et al. 2017).Just as the overall normalisation of the MZ g R decreases as a function of redshift, so too does the normalisation of the MZ * R (e.g., Gallazzi et al. 2014;Cullen et al. 2019).
Many studies tend to focus on just the gas-phase or stellar metallicities, yet recent works (Topping et al. 2020;Cullen et al. 2021;Fraser-McKelvie et al. 2022;Greener et al. 2022, etc) have shown that there is potentially a large amount of power in examining them in conjunction.Physically, the link between stellar and gas-phase metals follows from the idea that newly formed stellar populations simply inherit the metallicity from the gas from which they were formed.Thus, the stellar metallicity of the galaxy reflects the mass (or luminosity) weighted average metallicity integrated over the galaxy's formation history.Conversely, the gas-phase traces recent evolutionary processes.As such, the stellar metallicity is not directly influenced by gas inflows and/or metal production from aging stellar populations, but is instead constantly producing new stellar populations that reflect the gas phase metallicity value at the time of formation.Stellar metallicities tend to be systematically lower than their gasphase counterparts across stellar mass (e.g., Gallazzi et al. 2005;Lian et al. 2018).In more detail, González Delgado et al. (2014) find that the largest divergence of the two metallicities is in the lowest mass galaxies.Fraser-McKelvie et al. (2022) attribute this divergence to the efficacy of inflows and outflows at different masses as well as a strong dependence on a galaxy's star formation history (Lian et al. 2018 indicate that initial mass function variations are another possibility).It should be noted that the absolute values of gas-phase metallicites are sensitive to metallicity diagnostic (Kewley & Ellison 2008).Therefore trends between gas-phase and stellar metallicities are likely also sensitive to the chosen diagnostic.The magnitude of the discrepancy between the two metallicities in simulations has been seen to be both larger (Yates et al. 2012) and smaller (De Rossi et al. 2017) than observations.Comparisons between simulated metallicities and observed metallicities is not always straightforward, however, without the use of mock observational techniques (i.e., radiative transfer codes to generate spectra; see discussion in Nelson et al. 2018).
Very high redshift ( ∼ 8) measurements of gas-phase metallicities have already been attained (e.g., Langeroodi et al. 2022;Shapley et al. 2023;Sanders et al. 2023), but stellar metallicities at these high redshifts are significantly more difficult to measure.The difference between the viability of obtaining gas-phase and stellar metallicities at high  can be attributed to the different methods of deriving the two metallicities.Briefly, gas-phase metallicities are derived from emission diagnostics, more specifically the direct and strong line methods (see Kewley & Ellison 2008;Kewley et al. 2019, for com-plete reviews).Stellar metallicities, on the other hand, are typically determined by absorption features in the stellar continuum.The stellar diagnostics require high continuum signal-to-noise and as such, at present, the majority of stellar metallicity measurements come from low redshift systems (e.g., Chisholm et al. 2019), stacking galaxies (e.g., Halliday et al. 2008;Cullen et al. 2019;Sextl et al. 2023), observing lensed galaxies (e.g., Patrício et al. 2016), or very long exposure times (Kriek et al. 2019).
It is now recognized that the scatter of galaxies about the massmetallicity relations is not random.Instead, systems with metallicities above (below) the MZ g R tend to have low (high) gas fractions and low (high) star formation rates (Ellison et al. 2008;Mannucci et al. 2010;Lara-López et al. 2010;Bothwell et al. 2013Bothwell et al. , 2016;;Yang et al. 2022, etc).This "correlated scatter" is naturally driven by the relationship between gas inflows, star formation, and metal production (though this is only true on the global scale; Baker et al. 2023).More specifically, gas inflows tend to drive down galactic metallicities while driving up gas content and star formation rates.Conversely, galaxies with low gas inflow rates steadily consume their gas, driving up metallicities while the gas fractions and star formation rates decay (Torrey et al. 2018;hereafter T18).This correlated scatter can therefore be thought of as a three-parameter relation between stellar mass, gas-phase metallicity, and some proxy for the gas richness of the galaxy (gas mass, SFR, etc.).Perhaps even more remarkably, the same three parameter fit that relates stellar mass, gas-phase metallicity, and gas richness can simultaneously fit galaxies at low-and high-redshift (Cresci et al. 2019;Huang et al. 2019, however recent JWST observations suggest that galaxies at  ≳ 6 deviate from the three parameter fit of lower redshift systems, see Curti et al. 2023).This so-called Fundamental Metallicity Relation (FMR) hints that the redshift evolution of the MZ g R is not disjoint from its own scatter.Instead, the same physical processes (e.g., inflows, gas consumption) that drive galaxies about the MZ g R at any given redshift, also cause an evolution in the MZ g R with redshift.
However, since the stellar metallicities, and therefore the MZ * R, are not directly influenced by gas inflows, it is not immediately obvious that there should be an equivalent FMR for stars.Yet, observationally (e.g., Zahid et al. 2017;Sextl et al. 2023) and in simulations (De Rossi et al. 2018), there is evidence for this, insofar as a relationship between the MZ * R and star formation exists.As of yet, there have been no observational studies systematically studying or quantifying the existence of an "FMR for stars".In simulations, however, there is some evidence for this relationship not being "fundamental" as redshift evolution has been seen (Fontanot et al. 2021).
However, within simulations as a whole, significant uncertainty still exists in the preferred method to model both star formation and outflows in galaxy formation models -on which the metallicity critically depends.While empirical relations such as the Kennicutt-Schmidt relation (Schmidt 1959;Kennicutt 1998; KS relation) provide guidance or constraints for setting star formation rates (SFRs) based on the current gas content within galaxies, the implications for galactic star formation rates in very low mass or high redshift galaxies remains unclear.Indeed, many cosmological hydrodynamic simulations (e.g., Illustris, TNG, EAGLE) employ sub-grid ISM pressurization that manifestly enforce the KS relation (Vogelsberger et al. 2014b;Pillepich et al. 2018b;Crain et al. 2015), giving rise to SFRs that evolve reasonably smoothly out to the highest redshifts.Yet, models that attempt to treat both the multi-phase nature of the ISM and locality of star formation and feedback generally predict SFRs more explicitly with a significantly increased level of time variability.(i.e., burstiness; Faucher-Giguère 2018).
Notably, both "smooth" and "bursty" galaxy formation models are able to match a wide range of galaxy scaling relations (e.g., the galaxy stellar mass function, star formation main sequence, mass-metallicity relation; see, e.g., Hopkins et al. 2014;Vogelsberger et al. 2014a).Yet, the burstiness or time-variability of star formation (and its associated feedback) is of somewhat critical importance to our understanding of galaxy formation.In particular, bursty feedback -and the rapid, periodic, expulsion of material from the galactic center -has been invoked as a preferred mechanism that converts dark matter cuspy profiles into cored profiles -alleviating an important tension between ΛCDM and observations (e.g., Lazar et al. 2020).
In this paper, we examine this relationship between the MZ * R and star formation using the Illustris, IllustrisTNG, and EAGLE hydrodynamical galaxy formation simulations.All three of these simulations are sufficiently large to capture galaxy populations where not only can the MZ * R be defined, but the scatter can also be characterized.Notably, while these simulations share several similarities in how the ISM, star formation, and metal enrichment are treated, they also model the feedback processes that shape both the MZ g R and MZ * R in detail differently.As such, these simulations can be used to gain insight into variations resulting from changes in adopted galaxy formation models.Moreover, we develop an analytical model to probe not only how and why the residual correlation emerges, but also how its properties change with redshift.
The structure of this paper is as follows.In Section 2 we describe our methods including a brief description of the simulations employed, and our definitions of stellar metallicity.In Section 3 we present our results including a comparison of the MZ * R's in Illustris, IllustrisTNG, and EAGLE, the residual correlations in the scatter, and its origin.In Section 4 we discuss our results, including a presentation of a simple toy model that can be used to capture the essential properties of the residual correlations, including the redshift evolution.Finally, in Section 5 we provide conclusions.

METHODS
In this work, we analyse the gas-phase and stellar metallicities in Illustris, IllustrisTNG, and EAGLE galaxies to understand their dependence on star formation and stellar mass (i.e., fundamental metallicity relations).Importantly, as will be detailed in the following sections, these three models have a commonality: sub-grid ISM pressurisation and smooth stellar feedback.Thus, though the detailed implementations of all the models are appreciably different, we believe that any results found in all of the models should comprise a fairly generic result of simulations with this ISM treatment; however, more testing is required within (and without, as we discuss in Section 4.2.1) the sub-grid ISM paradigm to confirm this.
In this section, we briefly describe each of the simulations and our selection criteria for galaxies within the sample analysed.All measurements are in physical units, except for the simulation box sizes, which are reported in co-moving units.We summarize the details of this section in Table 1.

Illustris
For our analysis, we utilize the Illustris (Vogelsberger et al. 2013(Vogelsberger et al. , 2014a,b;,b;Genel et al. 2014;Torrey et al. 2014) cosmological simulations.These runs were performed with the moving-mesh hydrodynamical code arepo (Springel 2010).The Illustris model includes several important astrophysical processes including gravity, hydrodynamics, star formation, stellar evolution, chemical enrichment, radiative cooling and heating of the ISM, stellar feedback, black hole growth, and active galactic nuclei (AGN) feedback.Most important for the purposes of our work are self consistent implementations of star formation, stellar feedback, and chemical enrichment.Illustris treats the dense, star forming ISM with the effective equation of state of Springel & Hernquist (2003;hereafter, SH03).In this model, gas beyond a threshold density ( H > 0.13 cm −3 ) forms stars probabilistically assuming a Chabrier (2003) initial mass function (IMF).Crucially, the stars that are formed adopt their metallicity from the gas from which they form.As the stars evolve, they return both mass and metals to the ISM.Stellar mass return and yield tables are employed to allow for the direct simulation of time-dependent mass return and heavy element enrichment.Illustris explicitly tracks nine different chemical species (H, He, C, N, O, Ne, Mg, Si, and Fe).
The main Illustris simulation included a single volume of size (106.5 Mpc) 3 at three different resolutions, which are denoted as: Illustris-1 for the highest resolution run (with 2 × 1820 3 particles) and Illustris-3 for the lowest resolution run (with 2 × 455 3 particles).In this work, we use the highest resolution run, Illustris-1, which we will henceforth use synonymously with Illustris itself.

IllustrisTNG
IllustrisTNG (Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2018;Pillepich et al. 2018b;Springel et al. 2018;Pillepich et al. 2019;Nelson et al. 2019a,b, hereafter TNG) is the successor to the original Illustris simulations.TNG includes updates to the original Illustris model in addition to alleviating some deficiencies.Importantly, TNG also models the dense, star forming ISM with the SH03 equation of state.As such, these two simulations are related, yet have appreciably different physical implementations (for a complete list of differences see Weinberger et al. 2017;Pillepich et al. 2018a).Critical for the discussion later in this work, TNG implements redshift dependent winds.A wind velocity floor is enforced in TNG to ensure that low mass haloes do not have unphysically large mass loading factors.The effect of this is that winds are more efficient at suppressing lowredshift star formation and, consequently, feedback is more effective at higher redshifts.TNG tracks the same nine elements as Illustris, but it additionally adds an "other metals" item as a proxy for other metals not explicitly tracked.
Whereas the original Illustris suite consisted of only one volume at several resolution levels, TNG is comprised of three different volumes, each with their own set of resolutions.The simulations within TNG are designated by their approximate box size -TNG50 (51.7 Mpc) 3 , TNG100 (110.7 Mpc) 3 , and TNG300 (302.6 Mpc) 3 , with the different resolution levels indicated as in the Illustris suite.Here, we employ the highest resolution run of TNG100 (TNG100-1; 2×1820 3 particles, hereafter simply TNG) to make a fair comparison with the original Illustris simulation.

EAGLE
Finally, we use the "Evolution and Assembly of GaLaxies and their Environment" (EAGLE, Crain et al. 2015;Schaye et al. 2015;McAlpine et al. 2016) cosmological simulations.Unlike Illustris and TNG, EAGLE employs a heavily modified version of gadget-3 (Springel 2005) employing the anarchy formulation of smoothed particle hydrodynamics (SPH; see Schaye et al. 2015 Appendix A), which alleviates problems with the original version of SPH (Sĳacki et al. 2012;Vogelsberger et al. 2012;Kereš et al. 2012;Torrey et al. 2012).EAGLE is a full-physics simulation that includes many of the same baryonic processes with hydrodynamics (star-formation, chemical enrichment, radiative cooling and heating, etc).The dense (unresolved) ISM is treated with a sub-grid equation of state (Schaye & Dalla Vecchia 2008;hereafter, SDV08), much like that of SH03.The SDV08 prescription forms stars according to a Chabrier (2003) IMF from the dense ISM gas.The density threshold for star formation is given by the metallicity-dependent transition from atomic to molecular gas computed by Schaye (2004), with an additional temperature-dependent criterion (see Schaye et al. 2015, their Section 4.3).Stellar populations evolve according to the Wiersma et al. (2009) evolutionary model and eventually return their mass and metals back into the ISM.While the SH03 equation of state describes the multiphase ISM with the average gas density, SDV08 use a polytropic equation of state.EAGLE tracks eleven different chemical species (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe).
EAGLE is comprised of several simulations ranging from size (12 Mpc) 3 to (100 Mpc) 3 .For this work, as an even-handed comparison to the selected Illustris and TNG runs, we employ data products of an intermediate resolution (2 × 1504 3 particles) run with a boxsize of (100 Mpc) 3 referred to as RefL0100N1504 (hereafter simply EAGLE).

Galaxy Selection
In all of the aforementioned simulations, gravitationally-bound substructures are identified using a subfind algorithm (Springel et al. 2001;Dolag et al. 2009) which relies on a friends-of-friends (FoF; Davis et al. 1985) algorithm to identify parent groups.In this paper, the total mass identified by subfind is used for calculating global galaxy properties whereas observational works typically estimate quantities within a given aperture.We limit the sample from each simulation to central galaxies with stellar masses 8.0 < log( * [ ⊙ ]) < 12.0 and (total) gas mass1 limits of log( gas [ ⊙ ]) > 8.5.This ensures that we have ∼ 10 2 star particles and ∼ 5 × 10 2 gas particles per selected galaxy, which we consider well-resolved (see Table 1 for initial baryon mass resolutions of each simulation).We note that while this resolution is too coarse to make statements about spatially resolved properties within the galaxies (e.g., metallicity gradients), it is sufficient for global properties of the systems.Further, following from Donnari et al. (2019), Pillepich et al. (2019), Nelson et al. (2021), Hemler et al. (2021), and Garcia et al. (2023), we define a specific star formation main sequence (sSFMS) in order to select only star forming galaxies.The sSFMS is defined through a linear-least squares fit to the median relation of galaxies with stellar mass log( * [ ⊙ ]) < 10.2 in mass bins of 0.2 dex.Beyond log( * [ ⊙ ]) > 10.2, we extrapolate the linear fit.From this sSFMS, we define (and omit) quiescent galaxies as those whose sSFRs lie greater than 0.5 dex below this relation.We note that our key results are not sensitive to these selection criteria for the sSFMS, which we discuss further in Appendix B.
With all of these choices, at  = 0 we obtain 44,965 galaxies in Illustris, 22,133 in TNG, and 14,720 in EAGLE.While all of these simulations have very similar halo mass functions, we note that the overall number of galaxies in the sample is quite different between the three.The difference between the number of galaxies in EAGLE and TNG can be attributed to: (i) the larger box of TNG simply containing more galaxies and (ii) lower gas fractions in EAGLE making it more sensitive to our gas mass cut at the low mass end.Illustris, on the other hand, has an elevated stellar mass function compared to TNG, which is due to the different feedback and wind implementations (Pillepich et al. 2018a;Pillepich et al. 2018b), as such there are a factor of ∼2 more galaxies at  = 0 than in TNG.
Throughout all of this analysis gas-phase metallicities are from star-forming gas particles only as a fair comparison with observations, whose gas-phase metallicities are limited to emission diagnosistics from star forming regions of galaxies (see Kewley & Ellison 2008;Kewley et al. 2019).Gas is determined to be star-forming based on criteria outlined in Section 2.1.1 for Illustris/TNG and Section 2.1.3for EAGLE.
We aim to offer predictions for the high redshift ( ≲ 10) MZ * R in each of these simulations.However, due to the aforementioned mass limits of our sample, beyond  = 8 there are not enough galaxies to create a statistically significant sample.We therefore offer predictions for galaxies only out to  = 8.

Stellar MZR
Critically, the analysis done in this work relies on the underlying relationship between the mass of a system and the total metal content locked in stars within galaxies.To that end, we first quantify the MZ * R in Illustris, TNG, and EAGLE as a baseline for the rest of the results in this work.In these three simulations we find that a relationship between stellar mass and stellar metallicity exists, as illustrated in Figure 1 for each simulation at  = 0. We include the MZ * Rs from Gallazzi et al. (2005) and Panter et al. (2008) as a point of comparison in Figure 1; however, we caution the reader that comparing metallicities from observations and simulations is not even-handed (see Nelson et al. 2018).Specifically, computing metallicities by averaging the metal content of star particles weighted by their mass (as is done in all three simulations analysed here) produces systematically elevated stellar metallicities compared to using radiative transfer codes more akin to mock observations of galaxies (Guidi et al. 2016).In that light, it is perhaps unsurprising that the  = 0 MZ * Rs presented in Figure 1 and higher redshift versions in Figure 2 are all elevated compared to the observations.Nevertheless, it should be noted that while Illustris and TNG appear to follow much closer to the Gallazzi et al. (2005) and Panter et al. (2008) MZ * Rs compared to EAGLE, a higher resolution run of EAGLE with modified sub-grid parameters (recal-L025N0752) more closely reproduces the observed MZ * R (see Schaye et al. 2015, their Figure 13).
In terms of the shape of the relation in the simulations, we find qualitative agreement between the three: roughly power-law with constant slope at low masses, leveling out at higher masses, qualitatively similar to the observations.In detail, however, the three simulations exhibit different behaviours.Specifically, we find that the MZ * R normalisation at  = 0 for the lowest mass galaxies vary Stellar mass versus stellar metallicity (both scaled to solar;  ⊙ = 0.0127) for star forming galaxies in TNG, Illustris, and EAGLE (left, centre, and right panels, respectively) at  = 0.Each panel is colour-coded by the total number of galaxies within each pixel, with lighter colours corresponding to fewer galaxies.We overplot the MZ * R from local observations (green, Gallazzi et al. 2005;orange, Panter et al. 2008).It is important to note that comparisons between stellar metallicities from simulations to stellar metallicities from observations is not necessarily one-to-one (see text for more details), thus we caution the reader in any comparisons against observed MZ * R trends.
Figure 2. MZ * R, as defined in Section 3.1, for star-forming galaxies in the TNG, Illustris, and EAGLE simulations (left, centre, and right panels, respectively) for  = 0 − 8. Galaxies are split into several different mass bins, and the median metallicity within each bin is measured.We require that there must be at least 10 galaxies in each mass bin shown (i.e., higher masses are less populated at higher redshift).The shaded regions correspond to the standard deviation about the median mass bins.The black lines and points correspond to observations from Gallazzi et al. (2005) at  = 0 (dashed line), Kashino et al. (2022) ranging from 1.6 <  < 3.0 (dotted line), and Cullen et al. ( 2019) from 2.5 <  < 5.0 (points).
by nearly 1.0 dex between the simulation with the EAGLE galaxies being more enriched than those in TNG, which, in turn, are more enriched than the galaxies in Illustris.At higher masses,  * ∼ 10 10  ⊙ , however, the normalisation of the metallicities appears more similar between the different simulations.We also find that the slope of the power-law relation at low masses differs between the simulations: Illustris galaxies increase in metallicity more rapidly as a function of stellar mass, than those in TNG and EAGLE.At  > 0, we find that differences in normalisation persist, and, in addition, the normalisations appear to evolve differently between the various simulations.In Figure 2, we trace the MZ * R back to  = 8 for each simulation, which we compare to observations from Gallazzi et al. (2005) at  = 0, Kashino et al. ( 2022) at 1.6 <  < 3.0, and Cullen et al. ( 2019) at 2.5 <  < 5.0 (though, see above discussion about comparison of stellar metallicities with observations).For simplicity, we define the MZ * R simply as the median metallicity within fixed mass bins2 .The shaded regions about each MZ * R represents the scatter about the relation, which we define as the standard deviation of metallicities within each mass bin.Consistent with the  = 0 MZ * R, the higher redshift Illustris galaxies are less enriched at low masses and follow a steeper trend than their TNG and EA-GLE counterparts.However, around  ≳ 4, the EAGLE galaxies have changed normalisation enough to be more consistent with Illustris and there is a significant increase in the slope of the relation.
To that end, there is a marked difference in the redshift evolution of the MZ * R between the three simulations: TNG has some modest normalisation evolution with lower redshift galaxies being more enriched than high redshift, EAGLE has a similar, albeit more extreme, evolution, and Illustris shows virtually no redshift evolution.This lack of evolution in the stellar metallicities over cosmic time in Illustris has been noted previously (D'Souza & Bell 2018).These authors attribute the modest redshift evolution of the MZ * R shown in Gallazzi et al. (2014) at  ≈ 0.7 as evidence that, in the absence of robust predictions at higher redshifts, the lack of redshift evolution in the MZ * R is physical.However, the predictions made here in TNG and EAGLE of an MZ * R that evolves with redshift is in agreement with other simulation models (e.g., Ma et al. 2016;Yates et al. 2021) and more recent observations (e.g., Choi et al. 2014;Leethochawalit et al. 2018;Cullen et al. 2019;Beverage et al. 2021;Kashino et al. 2022).

Residual correlations in the scatter
Within the scatter of the MZ g R residual correlations have been shown to exist.In particular, a secondary correlation with the star formation rate (SFR) or specific star formation rate (sSFR; SFR normalized by galaxy mass) has been seen both observationally (e.g., Ellison et al. 2008;Mannucci et al. 2010;Andrews & Martini 2013) and in simulations (e.g., De Rossi et al. 2017 in EAGLE;T18 in TNG).We find that galaxies' stellar metallicities also follow this trend: at a fixed stellar mass, galaxies with lower stellar metallicities generally have higher sSFRs (Figures 3 and 4), echoing a similar result at  = 0 for a higher resolution run of EAGLE by De Rossi et al. (2018; see Section 4.2 for a more in-depth comparison) .Figure 3 shows the MZ * R for each simulation at  = 0 colour coded by their sSFR.Qualitatively, we find that for all simulations, save for TNG at  = 0, there is a clear residual correlation of the MZ * R with sSFR.At  = 0 in TNG, lower mass galaxies ( * ≲ 10 9.5  ⊙ ) the residual correlation appears weaker than their Illustris or EAGLE counterparts.However, as we show in Appendix A, at  > 0, the galaxies fall into the more familiar relation like that of Illustris and EAGLE at  = 0. Given the overall model similarities -and differences -between Illustris and TNG, we attribute the lack of a clear residual correlation at  = 0 in TNG to the redshift-dependent winds in the model (which was not present in Illustris).
More quantitatively, Figure 4 shows the stellar metallicities as a function of sSFR in fixed mass bins (coded by their colour) for each simulation at  = 0.The overplotted lines are the best fit within the corresponding mass bin.We find that the majority of these best fits (at all redshifts) have negative slopes, indicating that, at a fixed stellar mass, galaxies with elevated stellar metallicities generally have lower specific star formation rates.As previously mentioned, one notable exception to this trend is the lowest mass galaxies in TNG at  = 0.This complements the qualitative statements above: the slope of the correlations in TNG at  = 0 here in the lowest mass bins (10 8.0 − 10 9.0  ⊙ ) is inverted.We show in Appendix A that at higher redshift, the residual correlation is much more apparent in this mass range for TNG.Another notable exception is in the highest mass bins.In all three simulations, the highest mass bins (purple in Figure 4), appear much flatter than, or inverted with respect to, their lower-mass counterparts out to  = 2 for EAGLE and  = 1 for Illustris and TNG.Interestingly, this flattening and inversion of the residual correlation within higher stellar masses has been noted previously in both simulations (Yates et al. 2012, l-galaxies;De Rossi et al. 2017, EAGLE) and observations (Yates et al. 2012, Sloan Digitial Sky Survey DR7) within the gas-phase metallicities and in simulations (De Rossi et al. 2018, EAGLE) in the stellar metallicities.De Rossi et al. (2017Rossi et al. ( , 2018) ) show, using modified AGN parameters within EAGLE, that this inversion is due to AGN feedback, which is almost certainly the case in the highest mass galaxies analysed here.Mannucci et al. (2010;henceforth M10) posit that if the metallicity, SFR, and stellar mass are all correlated, then perhaps there exists a tighter parameterisation than a traditional MZR.They fit the mass-(gas-phase) metallicity-SFR ( * −  gas − SFR) relation with a linear combination of the stellar mass and star formation rates,

Projection of least scatter
where  is a free parameter3 .This free parameter encodes the relative importances of stellar mass and SFR.For example, taking  = 0.0, the original form of the MZR is recovered, which would suggest that the SFR is unimportant in minimzing the scatter of the relation.
Conversely, with  = 1.0 the primary driving factor of the scatter is the specific star-formation rates (sSFR) of the galaxies.We employ this same parameterisation for the stellar mass-stellar metallicity-SFR ( * −  * − SFR) relation, later in this section.
To determine the best   we fit a linear regression of log( * [ ⊙ ]) versus   for values of  ranging from 0 to 1. Whichever value of  corresponds to the smallest standard deviation about fixed   bins from the linear fit is deemed the best fit (demonstrated in Figure 5 for EAGLE at  = 2).Further, we define an uncertainty of this parameter by taking any  that only varies by 1% of the minimum dispersion (see dashed gray lines in top panel of Figure 5).
Figure 6 shows the evolution of the  parameter as a function of redshift in each of the simulations with the errorbars corresponding to the aforementioned 1% uncertainty of .We compare against observational values determined in M10 ( = 0.33), Andrews & Martini (2013, henceforth AM13;  = 0.66), and Curti et al. (2020, henceforth C20;  = 0.55) for the gas-phase.Since these observational values specifically refer to the gas-phase and not stellar metallicities, we caution the reader to take comparisons against observations lightly.Another important consideration is that the relative weighting of the mass and SFR has been shown to depend strongly on the metallicity diagnostic used (AM13).
Regardless, we ultimately find a non-zero evolution in the best fit  as a function of redshift.At  = 0, we find relatively weak dependence on the SFR; in fact, in TNG there is virtually no dependence on the SFR.This lack of importance placed on the SFR is perhaps unsurprising in light of the lack of a significant  * −  * − SFR relation at that redshift in TNG (mentioned in the previous section and discussed further in Appendix A).In Illustris and TNG,  increases slightly with increasing redshift.In Eagle, on the other hand,  decreases weakly with increasing redshift.While the  * −  gas − SFR relation is oftentimes referred to as the fundamental metallicity relation (FMR; M10; AM13; C20) and has been shown to have no evolution in the relation out to  ∼ 2.5 − 3.5 (e.g., M10; Lara-López et al. 2010), we avoid interpreting the  * −  * − SFR as an "FMR for stars" due to the redshift evolution seen in each of the simulations presented here.Interestingly, the redshift dependence on this relationship is also noted in Fontanot et al. (2021) with a semianalytic model (SAM; see Section 4.2 for a more in-depth comparison , total star formation rate per unit stellar mass of the galaxy; sSFR).Lighter colours correspond to higher star forming galaxies while darker colours are lower star forming galaxies.At fixed stellar mass lower metallicity galaxies tend to have higher sSFRs.The inset on the right panel is a comparison with De Rossi et al. ( 2018) who use a higher resolution run of EAGLE at  = 0; full discussion of this comparison can be found in Section 4.2.Note that the left panel of TNG galaxies, at lower masses, does not clearly display this relationship at  = 0, we discuss this further in Section 3.3 and Appendix A. ).Linear least-squares fits are provided by the outlined lines for each mass bin.We note that most of these slopes, across redshift, are negative, indicating a decrease in metallicity with increasing sSFR for a given mass.
with this model).Evidently, as we have shown, this  * −  * − SFR relationship is not "fundamental" insofar as it evolves with time.

Origin of 𝑀 * − 𝑍 * − SFR Relation
The  * −  gas −SFR relation (i.e., FMR) follows from fairly straightforward arguments: galactic inflows of pristine gas drive down the metallicity, but increase the available gas reservoir for star formation.Conversely, systems with limited or no gas inflows steadily consume their gas reservoirs while creating new metals.Stellar metallicities, however, are not directly impacted by inflows.Newly formed stars inherit the metallicity of the gas in which they form.As such, the stellar metallicity represents the average gas-phase metallicity of the galaxy integrated over its entire lifetime.It follows, then, that the gasphase and stellar metallicity of a system may be correlated, where the gas phase metallicity has the ability to more rapidly change in the presence of, e.g., rapid pristine gas inflows or significant enrichment events.In this section we investigate the extent of the correlation between gas and stellar metallicities in Illustris, TNG, and EAGLE.By defining the MZ g R in the same fashion as outlined in Section 3.1 for the MZ * R and interpolating between the mass bin centers, we obtain an offset from the median relation in both the gas-phase (Δ gas ) and stellar (Δ * ) components for each galaxy.We find that there is a linear relationship between these offsets in all three simulations (see Figure 7 for  = 0 offsets).This linear relationship implies that the gas and stellar phase metallicities are tightly coupled.We fit these offsets, Δ gas versus Δ * , with a linear least squares fit requiring that the intercept pass through the origin4 .As we discuss further in Section 4.1, slopes approaching a value of unity (like that of Illusthus we expect the distributions of Δ * and Δ gas to be centered around 0 (which, indeed, can be seen is roughly true from the histograms in the top panels of Figure 7).
tris' at  = 0) indicate a more direct relationship between the offsets of the stellar and gas phase metallicities with respect to their respective mass-metallicity relationship.Flatter slopes (e.g., EAGLE at  = 0) indicate a weaker relationship.Similar to the slope, understanding the strength of the correlation also depends on the magnitude of scatter within the offsets.A very tight correlation implies that stellar metallicities are more directly linked to the gas-phase, while broader distributions imply that stellar metallicities are less directly linked to the gas-phase.Thus, both the slope of the offsets and scatter within the relation inform how quickly the stars can respond to perturbations within the gas-phase metallicity.We find the slope of the offsets to be a more natural measure of the relationship here and thus use it as our primary diagnostic in this discussion.
It is important to note that while Δ gas and Δ * are related to the absolute  gas and  * of a system, they are different quantities.We therefore caution the reader against taking variations in the bulk offsets from the MZRs across redshifts as representative of the difference between a single galaxy's metal evolution through two cosmic times (we discuss this subtle difference with our treatment of a toy model in Section 4.1).
Performing the same offset analysis on at all redshifts, we obtain the time evolution shown in Figure 8.We find general agreement between the TNG and EAGLE simulations at  ≲ 3: a shallower slope that steepens with redshift.Illustris, however, follows the opposite trend in this same redshift range, a very steep relation that weakens going back in time.Beyond  ≳ 4, all three simulations have slopes that increases with increasing redshift.
In spite of the  * −  * − SFR having non-negligible redshift evolution in the strength of the residual correlation across redshift, the  * −  * −SFR is rooted in the close relationship of gas-phase and stellar metallicities across time which is consistently present across time, albeit with some change in strength.

Toy Model
Stellar metallicities are naturally linked to gas-phase metallicities.When stars are formed, they inherit the metallicity of the gas from which they were born.Thus, changes in the gas-phase metallicity will eventually impact stellar metallicities, albeit with a time delay that accounts for the fact that it takes time for a sufficient mass of new stars to be formed with the "new" inherited gas-phase metallicity value.We model this behaviour with a simple analytic framework to better understand the origins of the MZ * R and the scatter about it.
We start by defining Δ * as the stellar metallicity of an individual galaxy,  * , subtracted from the median relation (i.e., MZ * R), ⟨ * ⟩, at that galaxy's particular mass.Thus, the rate of change of a galaxy's offset from the MZ * R can be written as Equation 2 tells us that the changes Δ * are sensitive to not only changes in the absolute metallicity of a galaxy, but, critically, on the evolution of the normalisation of the MZ * R, as well.From this we obtain two interesting regimes: (i) in the limit that the normalisation of the MZ * R does not evolve, or evolves much more slowly compared to the time-scale of stellar formation ( SF ), Δ * 's are driven only by changes in the absolute metallicity and, conversely, (ii) a galaxy with no metal evolution over some time could experience a change in Δ * via the overall normalisation of the MZ * R changing.We begin by considering the first of the two regimes -a steady state model in as defined in Section 3.2.1),versus redshift in EAGLE (blue cirlces), Illustris (orange triangles), and TNG (green stars).The dashed line represents the relative importance found by M10 ( = 0.33), the dot-dashed line is that found by AM13 ( = 0.66), and the solid line is C20 ( = 0.55; note that all these works present  with the gas-phase metals).The errorbars correspond to the uncertainty associated with 1% deviations from the minimum dispersion of the best fit   (see Figure 5 and Section 3.2.1 for more details).  .Bottom Panels: Offsets of both the gas-phase (abscissa) and stellar (ordinate) metallicities of star forming galaxies in EAGLE, Illustris, and TNG (left to right) at  = 0 from the MZRs.In each panel, a linear least squares best fit line to this relation is overplotted (black line; fixing the intercept to the origin) as well as a one-to-one line (gray solid line).The slope of the linear least squares fit is the relevant quantity in this relation and allows us to quantify how correlated the two metallicities are to each other; a slope of unity would suggest that the gas-phase and stellar metals are exactly correlated, while a slope of 0 suggests that they are uncorrelated.Top Panels: Normalised distributions of the offsets from each simulation's MZ g R. We fit both a Gaussian distribution (red) and skewed Gaussian (blue) to the distribution.The typical mean and standard deviation of these distributions (across redshift) is ∼0.0 dex and ∼0.12 dex, respectively.We find that the distributions have some non-zero skew, which we discuss in more detail in Section 4.1.which there is no MZ * R evolution -to better understand how a toy galaxy's gas-phase and stellar metallicities interact.Then we explore the second regime and discuss the implications of the co-evolution of the MZ * R as a whole underneath galaxy-by-galaxy evolution.

Steady State Model
In principle, the rate of change of the stellar metallicity of an individual system should be a function of (i) the amount of stars being formed, (ii) the present gas-phase metallicity of gas forming new stars, and (iii) the present stellar metallicity.To attain an analytic relation that is characterised by all three, we first define the stellar metallicity of the galaxy as the total mass of metals locked within stars, normalized by the total stellar mass: By differentiating equation 3, and understanding that the rate of change of metals in stars is simply the gas-phase metallicity scaled by the amount (by mass) of stars that are being formed5 , we obtain a straight-forward expression for the rate of change of the stellar metallicity that depends on exactly the three parameters described above, From equation 4, it follows that the stellar metallicity should tend towards the gas-phase metallicity value (of the star forming gas) over time.Moreover, the rate at which this occurs is set only by the amount of new stars being formed normalized by the total stellar mass of the system, i.e. the specific star formation rate.However, the quantity that we are comparing to in this model is the offsets in stellar metallicity, Δ * .We must therefore use equation 2 to model how the offsets in stellar metallicity respond.Substituting the absolute metallicity with its offset and MZR terms, as well as rearranging, we are left with For now we will assume that the MZ * R does not evolve (d⟨ * ⟩/d = 0) and that the stellar and gas-phase MZRs are at the same value (⟨ gas ⟩ = ⟨ * ⟩) such that the first and third term of the left-hand side are zero.We will revisit these assumptions in Section 4.1.2.Applying these assumptions to equation 5 we are left with With this representation of the evolution of galaxy metal content, we model the response of a toy galaxy in this paradigm.For simplicity, we assume that our toy galaxy initially lies directly on both the MZ g R and MZ * R, so that its offsets from either are zero 6 .Next, we draw a deviation for the gas-phase metallicity (Δ gas ) from a normal distribution with a standard deviation of 0.12 dex, corresponding to injected back into the ISM from the stars from supernovae and AGB winds (see, e.g., Peng & Maiolino 2014).For simplicity, and to isolate its effect, we neglect this term here.We discuss the implications of this in Section 4.1.3. 6We note that non-zero offsets in either the gas-phase or stellar metallicities (or both) do not significantly impact the results of this toy model.a typical standard deviation of the distribution of gas-phase metallicity offsets in the simulations (see, e.g., top panels of Figure 7), centered at an offset of 0 dex.The toy galaxy's gas-phase metallicity tends towards the randomly drawn value over a time-scale drawn from an exponential distribution with some coherence time-scale,  c , nominally set to one-tenth a Hubble time (as per T18; see discussion later in Section 4.1.3).Physically these perturbations are meant to mimic a potential pathway that a galaxy's gas-phase metallicity will undertake as it evolves.As its gas-phase metallicity changes, we model the response of the toy galaxy's stellar metallicities by equation 6.Specifically, we model its formation of new stars with a stellar formation time-scale,  SF , as the inverse of sSFRs consistent with observations (one over the hubble time).
Simply put, the coherence time-scale represents how quickly the gas-phase metallicity is perturbed, while the stellar formation timescale informs how quickly the stars can respond to these changes.Together the coherence and star formation time-scales fully parameterise the behavior of our toy model.Thus, it is informative to define a (dimensionless) ratio of the two: Since, nominally,  c is one-tenth a Hubble time and  SF is a Hubble time, Γ takes a fiducial value of 0.1.It is worth noting that this Γ holds all of the predictive power of our model.By increasing Γ, the gasphase metallicity perturbations take longer (or identically, the system forms stars more quickly) and thus the stellar metallicity might track the gas-phase metallicity closely.Conversely, by decreasing Γ, we expect that the stellar metallicities take longer to equilibrate to the gas-phase metallicities, and therefore the two metallicities are not likely to remain matched as the gas phase metallicity changes.Thus, the strength of the relation of the offsets of the metallicities from their respective MZR's is a direct effect of Γ.
By scaling the initial conditions for the sSFR7 we find the strength (i.e., slope) of the offsets from the MZR values can be uniquely described by a value of Γ (see Figure 9).Lower (higher) Γ corresponds to a weaker (stronger) correlation between the two metallicities.This result agrees with naïve expectations: if the gas-phase metallicity evolves very quickly compared to the time it takes to form new stars, then the stellar metallicities will significantly lag behind the gas-phase metallicities.
We do note that the distribution of gas-phase metallicity offsets in the simulations more closely follows a skew-normal distribution with skewness parameters ranging from -1 to -3.5 (with the exceptions of TNG  = 2, Illustris  = 0, and EAGLE  = 7 and 8, which have skews of 0.0, -0.59, +0.36, and +0.70, respectively).Despite the magnitude of skew > 1 indicating that the distribution is non-Gaussian, when sampling gas-phase metallicity perturbations from a skew-normal distribution in our analytic model, we find that the overall slope of the correlation remains unchanged.Further, as mentioned in Section 3.3, we forced the intercept of the offset relation to the origin; however, when sampling from a skewed normal distribution, we remove this requirement.This is sensible since the skewed normal distribution potentially shifts the 'typical' value of the distribution away from an offset of 0. We note, however, that even if we continue to require the regression to pass through the origin, the quantitative change in the slope is modest.Thus, since the key results found here are not significantly impacted by using a skew normal distribution, we opt to use a regular Gaussian distribution.
With the predictions offered by our analytic model and the slopes found in Section 3.3, we predict the range of values for Γ in each simulation (Figure 9).For EAGLE, we find that the log Γ values range from ∼ −0.5 to 0.2, in TNG we find log Γ ranges from ∼ −0.55 to 0.8, and in Illustris they range from ∼ 0.2 to 0.5.

Impact of the underlying MZRs
It is important to note that the slopes quoted in Figures 7 and 8 are based on Δ gas and Δ * -offsets from median relations.As we noted at the beginning of this discussion, changes in Δ are sensitive to not only changes in a system's absolute  but to changes in the underlying MZRs as well.Yet, in the previous section, we assume that the median relation does not evolve with time and that the gasphase median relation is equal to that of the stellar relation.We relax that assumption in this section and consider how the evolution and interplay of the underlying relations impact the interpretation of our toy model.
Both the first and third terms on the right-hand side of equation 5 imply that a galaxy's Δ * changes owing to bulk changes in the underlying median relation.While these terms are related, each represents slightly different physical mechanisms.The first term expresses the difference between the median MZ g R and median MZ * R.This term is qualitatively similar to that of the offset term presented in Section 4.1.1:the entire MZ * R tends towards the MZ g R over some time-scale.The larger the difference in the two median relations, the faster the MZ * R will respond.Physically (both globally and individually), this links to the idea that the stellar metallicity is adopted from the gas where stars form.The next generation of stars will have a higher metallicity when the gas-phase metals are higher than that of the stars.On the other hand, the third term models the evolution of the normalisation of the MZ * R. Individual Δ * s will experience changes owing to the evolution in the MZ * R, regardless of what the absolute stellar metallicity is doing.Changes in the MZ * R are set partially by the whole relation tending towards the MZ g R, but also include newly formed metals within the stars of the galaxies themselves.Importantly, both the first and third terms are sensitive to changes in the MZRs.
The evolution of MZ * R plays a more explicit role in interpreting our toy model, we therefore isolate its effect here by only considering the second and third terms on the right-hand side of equation 5.The evolution of the normalisation only becomes important when it is significantly faster than stellar formation time-scales, however.Once MZ * R evolution becomes significant enough, we can no longer neglect these additional terms in equation 5.By taking the difference in metallicity within all the mass bins within the MZ * R at each redshift and dividing by difference in Hubble times between two given redshits, we find that the median evolution of the normalisation of the MZ * R is 0.08 dex/Gyr in EAGLE, 0.02 dex/Gyr in TNG, and 0.0 dex/Gyr in Illustris.We note that these values are a crude approximation that do not account for possible mass and redshift dependencies on the evolution rate of the MZ * R; nevertheless, they are useful for considering the typical redshift evolution within each simulation.Furthermore, these evolution rates reiterate the qualitative interpretations of Figure 2 in Section 3.1 -the most redshift evolution of the MZ * R is found within EAGLE, with modest evolution in TNG, and virtually no evolution in Illustris.
When adding the MZ * R evolutionary term with both the Illustris and TNG rates, we find no significant effect on the overall obtained slope of the relation of Δ gas versus Δ * compared to the slope from the steady state model.For EAGLE, since the MZ * R normalisation changes more quickly, the evolutionary term becomes more important.Specifically, when using the EAGLE MZ * R evolution rate, log Γ ≲ −0.25 (corresponding to only  = 0) is where the model begins to deviate significantly from the steady state.In this regime, the MZ * R is evolving so rapidly compared to the stellar formation time-scales that the second term in Equation 2 dominates the change in Δ * over time.The significant MZ * R evolution has the net effect of systematically moving the average Δ * away from 0. Yet, by definition, these Δ * values should be centered around ∼ 0. The reason behind the systematic shift is that our model only runs a single toy galaxy's metallicity track; the additional evolutionary term just represents the typical amount the MZ * R evolves -the model does not actually recalculate an MZ * R based on a population of galaxies' stellar metallicities at each time step.As such, fitting the relation with a line that passes through the origin (as we did in Section 3.3) produces ill-defined results.In order to rectify the toy model in this regime, a more sophisticated approach would be required that traces evolutionary tracks of multiple systems at once all while evolving the overall normalisation of the MZ * R by a fixed amount at each time step, from which Δ * 's from the median relation could be calculated.While it is possible to construct such a model, we favor the simplicity of the toy model detailed above in light of the relatively good job that it does reproducing the simulations' offsets in conjunction with the limited range in which such a careful tracing of MZ * R evolution is necessary.
The first term on the right hand-side of equation 5 explicitly encodes information about the difference in absolute stellar and gasphase MZRs but also implicitly encodes information about their respective evolution.In the limit that the neither MZR evolves, this term is just a constant that changes the rate at which Δ * evolves.A constant, non-zero offset should steepen the slope of Δ gas versus Δ * for ⟨ gas ⟩ > ⟨ * ⟩ as it increases the rate of change of Δ * .It should be noted that for large differences in ⟨ gas ⟩ and ⟨ * ⟩ this term may play a significant role.Both the MZ * R and MZ g R do have some evolution associated with them, however.The median gas-phase metallicity evolutionary rates are 0.09 dex/Gyr in EAGLE, 0.08 dex/Gyr in TNG, and 0.05 dex/Gyr in Illustris (all faster than their MZ * R counterparts).A complete treatment of this model would need to explicitly track the absolute metallicities of the MZRs, as well as their evolution during the toy galaxy's evolution.Since the MZ g R evolutionary values are all significant, the current model presented here does not have the capabilities to perform such a task (for reasons detailed above).
Besides the impact on the change on ⟨ gas ⟩, MZ g R evolution has no impact on our toy model.The reason that evolution of the MZ g R is otherwise unimportant in our model is that we draw the gas-phase metallicity perturbations about the MZ g R itself (Δ gas values), centered at an offset of 0.Moreover, we find that the scatter about the MZ g R does not vary significantly across time in the simulations.Therefore, our perturbations to the gas-phase metallicity are not sensitive to changes of the MZ g R.
Finally, we note that the underlying MZRs have some structure as a function of mass.We make the simplifying assumption that the stellar mass of the galaxy does not evolve significantly with time.In reality, there would likely be some contribution to Δ * owing to the galaxy increasing in stellar mass as it assembles its stellar populations.
In summary, the additional MZ * R evolutionary term may appear necessary to properly compare our toy model's predictions to the obtained slopes of the offsets from Section 3.3; however, with the exception of  = 0 in EAGLE, there is no difference in the interpreta-tion by including the secondary term.Regardless, we caution that the significant MZR evolution in EAGLE makes the previously derived values of log Γ more uncertain compared to, e.g., Illustris which has very little MZR evolution.Additionally, we find that the difference in the median MZRs increases the slope of Δ gas versus Δ * when ⟨ gas ⟩ > ⟨ * ⟩ and is a constant offset.

Limitations of the model
We note that there appears to be some tension with the T18 model, which (effectively) predicts that log Γ = −1 for TNG.Recall that we neglect the return of mass and metals back into the ISM from stars dying.This removes a (1 − ) term on the sSFR, where  is the mass return fraction from stars (identically, a 1/(1 − ) scaling on  SF ).If we instead included this return fraction, assuming that is a constant over  c , it is equivalent to just scaling Γ by (1 − ).This does indeed shift us closer to agreement with T18, the extent of the agreement, however, depends on the adopted value of .In order to produce a slope that corresponds to a log Γ = −1, we must assume stellar mass return values ranging from ∼ 66 − 97%.This value is much larger than current understanding predicts (i.e., ≲40%; Pillepich et al. 2018a;Hopkins et al. 2023).It is also possible that T18 underestimates the amount of time that a galaxy will persist above or below the MZR.The equivalent coherence time-scale from T18 (called metallicity evolution time-scale in that work) uses a Pearson correlation coefficient to describe the strength of the correlations in the offsets from the MZ g R.After one coherence (metallicity evolution) time-scale the strength of the correlation is weaker, but not gone entirely.Indeed, van Loon et al. (2021) find that the time-scale for evolution of gas-phase metallicities is ≫ 1 Gyr in EAGLE, which (particularly at high redshift) is significantly greater than a tenth of the Hubble time.Thus, the subtle difference in interpretation of this time-scale may push the coherence time-scale (as defined in this work) to larger values, shifting the T18 predictions closer towards agreement with this work.Furthermore, the stellar formation timescale could also be dependent on the mass of the system, as Matthee & Schaye (2019) show, with less massive galaxies having more short time-scale fluctuations.Another possibility is that our method of randomizing metallicities isn't really a fair match to what's happening in the simulations.For example, we model the gas-phase metallicity changing at a constant rate over the entire coherence time-scale.The rate of change of gas-phase metallicity might change as a function of time (from, e.g., infalling material with a non-uniform density distribution).Moreover, it is possible that after a single coherence time-scale the system's gas-phase metallicity is not immediately perturbed again, giving additional time for the stellar metallicity value to tend towards the gas-phase value.
It is also worth noting that this model does not include every physical mechanism with which a system's metal content could evolve.For one, the role that mergers play is underestimated by our treatment here.While the role that infalling gas is accounted for (i.e., a possible mechanism for the gas coherence time-scale), merging galaxies also bring in stars that could change the stellar metallicity completely disjoint from new star formation.Dry mergers in particular would bypass both the gas-phase coherence and stellar formation time-scales entirely by changing the stellar metallicity of a system very quickly without forming new stars or perturbing the gas-phase metallicity.The fraction of stars that are born ex-situ (outside of the galaxy) depends on the mass of the galaxy (Rodriguez-Gomez et al. 2016;Davison et al. 2020).The ex-situ fraction is negligible for a majority of the sample analysed in this work (see, e.g., Pillepich et al. 2015) and therefore its impact on the overall stellar metallicity is likely negligible.Past log  * ≳ 10.2 log  ⊙ this fraction becomes more significant, however (> 25%).Additionally, AGN feedback plays a vital role in the evolution of metal content of high mass galaxies, yet does not fit into our toy model here.De Rossi et al. (2017, 2018) and van Loon et al. (2021), using the EAGLE model, show that AGN have the effect of suppressing both the gas-phase and stellar metal content in the highest mass galaxies via quenching star formation and ejecting enriched material with outflows.These effects are thought to cause the inversion of the  * −  gas − SFR relation at the highest masses (as in Yates et al. 2012), and are likely contributing to the inversion of the  * −  * − SFR relation in the highest mass galaxies shown in Figure 4. Regardless, such quenching or outflow mechanisms are simply not taken into account in this toy model, yet it is clear that they are important to some extent in the metal evolution of the highest mass systems.

Comparisons with other models
Just as in observations, there is limited work quantifying the  * −  * − SFR relation in simulations.Fontanot et al. (2021) used the semi-analytic model GAlaxy Evolution and Assembly (gaea) and found a residual trend between  * −  * − SFR.These authors investigate the evolution of both the gas-phase and stellar metallicities through time and recover both the typical gas-phase FMR and a similar relation for stars.Similar to the results shown in this work (see Section 3.2; Figure 6), they find the  * −  * − SFR correlation evolves with redshift.Differently than this work, however, they quantify this change in the relation not in terms of a three component relation, but instead by measuring an offset from the  = 0 defined  * −  * − SFR relation8 .As such, this comparison is limited to systems with overlapping masses or sSFRs at  = 0 and the desired redshift ( = 2.24, 4.18, and 4.88).In these systems, they find that offsets from the  = 0  * −  * − SFR relation are a factor of ≳2 more than in the gas-phase.These authors consider this increased scatter evidence for significant redshift evolution of the  * −  * − SFR and attribute this evolution to the increase of baryonic mass locked in the stars as galaxies evolve.
De Rossi et al. (2018) also investigated the correlation of gas content within the MZ * R within EAGLE (though a smaller box, higherresolution run; recal-L025N0752,  baryon = 0.226 × 10 6  ⊙ ) at  = 0.Those authors found that the scatter of the MZ * R at  = 0 is correlated with not just the sSFR, but also the gas fraction and massweighted age of a galaxy's stellar population, specifically that higher metallicity systems: (i) are older, (ii) have a lower gas fractions, and (iii) have lower sSFRs (and vice-versa).The inset in the right-most panel of Figure 3 shows the quantitative comparison between the high resolution EAGLE and our analysis here.Qualitatively, these three findings corroborate what we find in Section 3.2 with the sSFR of galaxies.While the EAGLE run we analyse has slightly elevated metalicities compared to the higher resolution run of De Rossi et al. (2018), the overall relationship is remarkably similar between the two runs, suggesting the  * −  * − SFR relationship is not dependent on the resolution of the simulation and rather a feature of the model at least at this low redshift.However, more robust comparisons at varying resolutions and at higher redshifts are required to confirm this.

Implications for stellar feedback modeling
As mentioned in Section 2.1, all of the cosmological simulations presented in this work are of the sub-grid ISM pressurization type (i.e., SH03 for Illustris and TNG, and SDV08 for EAGLE) which are characterized by gradual (i.e.relatively non-bursty) star formation and stellar feedback.Thus, it appears that the existence of a  * −  * − SFR correlation in all of these models -which is naturally explained by our toy model with smoothly evolving star formation rates and regularly evolving metallicities -constitutes a reasonably generic prediction of these smooth (i.e.non-bursty) stellar feedback models.These sub-grid implementations of the star forming ISM are not the only ways to model the ISM, however.The Feedback in Realistic Environments (FIRE; Hopkins et al. 2014) model, for one, attempts to more explicitly model the multi-phase ISM.As a consequence, the star formation (and, as a result, stellar feedback) in the FIRE model can occur in large bursts -especially in low-mass and high-redshift systems.These large bursts of feedback can quickly remove a significant fraction of the gas from the galaxies.
The analytic model presented in Section 4.1 is implicitly smooth star formation: gas is accreted onto the galaxy and, over some stellar formation time-scale  SF , forms into stars.While this model is certainly simplistic, it produces a sensible result in the stellar metallicity of the system lags behind the gas-phase metallicity as subsequent generations of stars form from the enriched ISM, evolve, and return the metals to form new stars.Our model breaks down in a bursty star formation paradigm.Rapid ejections of gas from the discs reduces the time that stars have to form as the star forming gas is removed from the system.It is likely that these bursts would significantly disrupt the process of stellar metallicites catching up to the gas-phase values.Thus, the strength of the residual correlation within the MZ * R, even if such a thing exists, is likely to be sensitive to the "burstiness" of stellar feedback within galaxies.Further, these interruptions would likely leave their mark on the correlation between the gas-phase and stellar metallicities, as well, potentially weakening the overall correlation.
Interestingly, the feedback prescription used in the SAM employed in Fontanot et al. (2021;gaea, Hirschmann et al. 2016) takes some inspiration from the FIRE model (Hopkins et al. 2014) insofar as it attempts to model bursty feedback, yet, as previously mentioned, Fontanot et al. (2021) still recover a similar residual correlation within the MZ * R. While quantitative comparisons on the strength of the correlation are difficult to make, owing to the largely different methodology in the two works, it is worth noting that the relationship indeed does exist to some extent in more bursty feedback models.However, SAMs reconstruct feedback via approximations to observations combined with theoretical models.Thus, it is possible that these approximations might curtail the effectiveness of the feedback compared to hydrodynamic simulations and a residual correlation may or may not appear in full hydrodynamic simulations, e.g., FIRE.
Further examination, in these hydrodynamic simulations with bursty feedback is required to fully understand the effect of burstiness on this residual trend.However, we suggest that the (lack of) correlation between gas and stellar metallicities for galaxy populations may provide a promising opportunity for constraining or confirming the role of bursty feedback in galaxy assembly.

CONCLUSIONS
In this work, we select star-forming, central galaxies with stellar masses 8.0 < log( * [ ⊙ ]) < 12.0 and gas masses log( gas [ ⊙ ]) > 8.5 from the Illustris, TNG, and EAGLE cosmological simulations.Taking the global metallicity of these galaxies, we create a stellar mass-stellar metallicity relation (MZ * R) and stellar mass-gas-phase metallicity relation (MZ g R) for each at  = 0 − 8.With these, we can find each galaxy's offset from each MZR and examine residual trends.We examine these residual trends by using a simple analytic model to track the co-evolution of the gas-phase and stellar metallicities.
Our conclusions are as follows: • The general trend of each of the MZ * R in each simulation is similar: less metals in low mass galaxies which increases with stellar mass before plateauing at the highest masses.The overall normalisation of the MZ * R (Figure 1) and the redshift evolution (Figure 2), however, vary significantly in between Illustris, TNG, and EAGLE.
• The scatter about the MZ * R is correlated with the sSFR of the galaxies.Higher (lower) metallicity galaxies tend to have lower (higher) specific star formation rates (Figure 3 and 4).This result follows suit from observations of the scatter about the MZ g R (e.g., Ellison et al. 2008;M10;Lara-López et al. 2010).However, we find that the dependence on SFR varies as a function of redshift (Figure 6) in contrast to observations and results from simulations of an "FMR" in the gas-phase.
• The origin of both this correlated scatter and the MZ * R can be traced back to the intercorrelation of stellar and gas-phase metallicities, which we have fairly strong correlations within these simulations (Figures 7 and 8).
• We present a simple toy model for the causal evolution of the stellar metallicities from the gas-phase metallicities.We determine that the driving factor between the correlation of the gas-phase and stellar metallicities depends on the interplay of the coherence timescale (how quickly gas-phase metallicities change) and the stellar formation time-scale (how quickly new stars form).By definition a dimensionless parameter, Γ, as the ratio of these two time-scales, we can quantify the correlation between the gas-phase and stellar metallicities in the simulations we analyse (Figure 9).The broad-strokes agreement in the residual correlations between Illustris, TNG, and EAGLE suggest that this  * −  * − SFR may be a fairly generic prediction of galaxy formation models with this ISM pressure support creating smooth star formation histories.As the gas-phase metallicities are perturbed, the stellar metallicities have sufficient time to respond and "catch-up" to the gas-phase metals.This, however, is not obviously true a priori in models with more bursty assembly histories, thus the existence of this residual correlation could potentially provide an observable diagnostic for ISM feedback models.1, for EAGLE, TNG, and Illustris (top, middle, bottom, respectively) as a function of redshift for the four sSFMS selection variations.Right: Value for the slope of the offset relations (as in Figure 7) for EAGLE, TNG, and Illustris (top, middle, bottom, respectively) as a function of redshift for the four sSFMS selection variations.
Figure1.Stellar mass versus stellar metallicity (both scaled to solar;  ⊙ = 0.0127) for star forming galaxies in TNG, Illustris, and EAGLE (left, centre, and right panels, respectively) at  = 0.Each panel is colour-coded by the total number of galaxies within each pixel, with lighter colours corresponding to fewer galaxies.We overplot the MZ * R from local observations (green,Gallazzi et al. 2005; orange, Panter et al. 2008).It is important to note that comparisons between stellar metallicities from simulations to stellar metallicities from observations is not necessarily one-to-one (see text for more details), thus we caution the reader in any comparisons against observed MZ * R trends.

Figure 3 .
Figure 3. Same as Figure1colour coded by specific star formation rate (i.e., total star formation rate per unit stellar mass of the galaxy; sSFR).Lighter colours correspond to higher star forming galaxies while darker colours are lower star forming galaxies.At fixed stellar mass lower metallicity galaxies tend to have higher sSFRs.The inset on the right panel is a comparison with DeRossi et al. (2018) who use a higher resolution run of EAGLE at  = 0; full discussion of this comparison can be found in Section 4.2.Note that the left panel of TNG galaxies, at lower masses, does not clearly display this relationship at  = 0, we discuss this further in Section 3.3 and Appendix A.

Figure 4 .
Figure 4. Stellar metallicity as a function of sSFR for star-forming galaxies in TNG, Illustris, and EAGLE (left, centre, and right, respectively) at  = 0 colour-coded by stellar mass bins.Galaxies are in mass bins of width 1.0 dex, centered at 8.5, 9.5, 10.5, and 11.5 log(  * [  ⊙ ] ).Linear least-squares fits are provided by the outlined lines for each mass bin.We note that most of these slopes, across redshift, are negative, indicating a decrease in metallicity with increasing sSFR for a given mass.

Figure 5 .
Figure 5. Demonstration of parameterisation of  * −  * − SFR correlation in terms of   (Equation 1) following from M10. Top: The dispersion (i.e., standard deviation of the residuals) from the linear fit as a function of , the free parameter in   .Shown here is the dispersions for fits of the EAGLE  = 2  * −  * − SFR relation.The red dashed lines indicate where the dispersion is minimized and at what value of .The gray dashed lines represent the uncertainty in  values corresponding to 1% deviations from the minimum dispersion.Bottom: The  * −  * − SFR relation parameterised by the   that minimizes the dispersion about the linear relation.For EAGLE at  = 2, we find that the best  = 0.37.

Figure 6 .
Figure 6., the relative weighting of star formation to stellar mass in the three component relation ( * −  * − SFR; as defined in Section 3.2.1),versus redshift in EAGLE (blue cirlces), Illustris (orange triangles), and TNG (green stars).The dashed line represents the relative importance found by M10 ( = 0.33), the dot-dashed line is that found by AM13 ( = 0.66), and the solid line is C20 ( = 0.55; note that all these works present  with the gas-phase metals).The errorbars correspond to the uncertainty associated with 1% deviations from the minimum dispersion of the best fit   (see Figure5and Section 3.2.1 for more details).

Figure 7
Figure7.Bottom Panels: Offsets of both the gas-phase (abscissa) and stellar (ordinate) metallicities of star forming galaxies in EAGLE, Illustris, and TNG (left to right) at  = 0 from the MZRs.In each panel, a linear least squares best fit line to this relation is overplotted (black line; fixing the intercept to the origin) as well as a one-to-one line (gray solid line).The slope of the linear least squares fit is the relevant quantity in this relation and allows us to quantify how correlated the two metallicities are to each other; a slope of unity would suggest that the gas-phase and stellar metals are exactly correlated, while a slope of 0 suggests that they are uncorrelated.Top Panels: Normalised distributions of the offsets from each simulation's MZ g R. We fit both a Gaussian distribution (red) and skewed Gaussian (blue) to the distribution.The typical mean and standard deviation of these distributions (across redshift) is ∼0.0 dex and ∼0.12 dex, respectively.We find that the distributions have some non-zero skew, which we discuss in more detail in Section 4.1.

Figure 8 .Figure 9 .
Figure 8. Slope of the offsets from the MZRs for EAGLE (blue cirlces), Illustris (orange triangles), and TNG (green stars) as a function of redshift with horizontal lines corresponding slopes obtained by running our analytic model (Section 4.1) with certain values of Γ (equation 7).The gradient background indicates that slopes closer to unity indicate a stronger relationship between Δ gas and Δ * .We note that the Γ = 0.1 line corresponds to predictions by T18.
Figure B1.Left: Determination of , the free parameter in Equation 1, for EAGLE, TNG, and Illustris (top, middle, bottom, respectively) as a function of redshift for the four sSFMS selection variations.Right: Value for the slope of the offset relations (as in Figure7) for EAGLE, TNG, and Illustris (top, middle, bottom, respectively) as a function of redshift for the four sSFMS selection variations.

Table 1 .
Summary of relevant properties of the three simulation models used in this work (Illustris, TNG, and EAGLE). baryon is the initial baryonic matter mass resolution in each simulation.