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

For our analysis, we utilize the Illustris (Vogelsberger et al., 2013, 2014a, 2014b; 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 ($n_{\rm 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\times 1820^{3}$ particles) and Illustris-3 for the lowest resolution run (with $2\times 455^{3}$ particles). In this work, we use the highest resolution run, Illustris-1, which we will henceforth use synonymously with Illustris itself.

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 low-redshift 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\times 1820^{3}$ particles, hereafter simply TNG) to make a fair comparison with the original Illustris simulation.

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 (Sijacki 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~{}{\rm Mpc})^{3}$ to $(100~{}{\rm 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\times 1504^{3}$ particles) run with a box-size of ($100~{}{\rm Mpc})^{3}$ referred to as RefL0100N1504 (hereafter simply EAGLE).

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 ($M_{*}-Z_{\rm gas}-{\rm SFR}$) relation with a linear combination of the stellar mass and star formation rates,

where $\alpha$ is a free parameter³³3We note that more complex relationships parameterising this residual correlation are available (e.g., Lara-López et al., 2010, 2013; Curti et al., 2020). We opt to use this specific parameterisation in view of its simplicity.. This free parameter encodes the relative importances of stellar mass and SFR. For example, taking $\alpha=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 $\alpha=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 ($M_{*}-Z_{*}-{\rm SFR}$) relation, later in this section.

To determine the best $\mu_{\alpha}$ we fit a linear regression of $\log(Z_{*}~{}[Z_{\odot}])$ versus $\mu_{\alpha}$ for values of $\alpha$ ranging from 0 to 1. Whichever value of $\alpha$ corresponds to the smallest standard deviation about fixed $\mu_{\alpha}$ bins from the linear fit is deemed the best fit (demonstrated in Figure 5 for EAGLE at $z=2$). Further, we define an uncertainty of this parameter by taking any $\alpha$ 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 $\alpha$ parameter as a function of redshift in each of the simulations with the errorbars corresponding to the aforementioned 1% uncertainty of $\alpha$. We compare against observational values determined in M10 ($\alpha=0.33$), Andrews & Martini (2013, henceforth AM13; $\alpha=0.66$), and Curti et al. (2020, henceforth C20; $\alpha=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 $\alpha$ as a function of redshift. At $z=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 $M_{*}-Z_{*}-{\rm SFR}$ relation at that redshift in TNG (mentioned in the previous section and discussed further in Appendix A). In Illustris and TNG, $\alpha$ increases slightly with increasing redshift. In Eagle, on the other hand, $\alpha$ decreases weakly with increasing redshift. While the $M_{*}-Z_{\rm gas}-{\rm 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 $z\sim 2.5-3.5$ (e.g., M10; Lara-López et al. 2010), we avoid interpreting the $M_{*}-Z_{*}-{\rm 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 semi-analytic model (SAM; see Section 4.2 for a more in-depth comparison with this model). Evidently, as we have shown, this $M_{*}-Z_{*}-{\rm SFR}$ relationship is not “fundamental” insofar as it evolves with time.

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 formed⁵⁵5We note that to be complete there should be some exchange of metals 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., 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, $\Delta Z_{*}$. 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 (${\rm d}\langle Z_{*}\rangle/{\rm d}t=0$) and that the stellar and gas-phase MZRs are at the same value ($\langle Z_{\rm gas}\rangle=\langle Z_{*}\rangle$) 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${}_{\rm g}$R and MZ${}_{*}$R, so that its offsets from either are zero⁶⁶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.. Next, we draw a deviation for the gas-phase metallicity ($\Delta Z_{\rm gas}$) from a normal distribution with a standard deviation of 0.12 dex, corresponding to 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, $\tau_{\rm 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, $\tau_{\rm 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 time-scale 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, $\tau_{\rm c}$ is one-tenth a Hubble time and $\tau_{\rm SF}$ is a Hubble time, $\Gamma$ takes a fiducial value of 0.1. It is worth noting that this $\Gamma$ holds all of the predictive power of our model. By increasing $\Gamma$, the gas-phase 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 $\Gamma$, 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 $\Gamma$.

By scaling the initial conditions for the sSFR⁷⁷7Though we choose to vary sSFR, this is, in principle, exactly equivalent to decreasing the coherence time-scale. We choose the former for the sake of convenience. we find the strength (i.e., slope) of the offsets from the MZR values can be uniquely described by a value of $\Gamma$ (see Figure 9). Lower (higher) $\Gamma$ 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 $z=2$, Illustris $z=0$, and EAGLE $z=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 $\Gamma$ in each simulation (Figure 9). For EAGLE, we find that the $\log\Gamma$ values range from $\sim-0.5$ to $0.2$, in TNG we find $\log\Gamma$ ranges from $\sim-0.55$ to $0.8$, and in Illustris they range from $\sim 0.2$ to $0.5$.

It is important to note that the slopes quoted in Figures 7 and 8 are based on $\Delta Z_{\rm gas}$ and $\Delta Z_{*}$ – offsets from median relations. As we noted at the beginning of this discussion, changes in $\Delta Z$ are sensitive to not only changes in a system’s absolute $Z$ 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 gas-phase 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 $\Delta Z_{*}$ 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${}_{\rm 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${}_{\rm 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 $\Delta Z_{*}$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${}_{\rm 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 $\Delta Z_{\rm gas}$ versus $\Delta Z_{*}$ 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\Gamma\lesssim-0.25$ (corresponding to only $z=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 $\Delta Z_{*}$ over time. The significant MZ${}_{*}$R evolution has the net effect of systematically moving the average $\Delta Z_{*}$ away from 0. Yet, by definition, these $\Delta Z_{*}$ values should be centered around $\sim 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 $\Delta Z_{*}$’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 gas-phase 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 $\Delta Z_{*}$ evolves. A constant, non-zero offset should steepen the slope of $\Delta Z_{\rm gas}$ versus $\Delta Z_{*}$ for $\langle Z_{\rm gas}\rangle>\langle Z_{*}\rangle$ as it increases the rate of change of $\Delta Z_{*}$. It should be noted that for large differences in $\langle Z_{\rm gas}\rangle$ and $\langle Z_{*}\rangle$ this term may play a significant role. Both the MZ${}_{*}$R and MZ${}_{\rm 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${}_{\rm 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 $\langle Z_{\rm gas}\rangle$, MZ${}_{\rm g}$R evolution has no impact on our toy model. The reason that evolution of the MZ${}_{\rm g}$R is otherwise unimportant in our model is that we draw the gas-phase metallicity perturbations about the MZ${}_{\rm g}$R itself ($\Delta Z_{\rm gas}$ values), centered at an offset of 0. Moreover, we find that the scatter about the MZ${}_{\rm 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${}_{\rm 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 $\Delta Z_{\rm*}$ 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 $z=0$ in EAGLE, there is no difference in the interpretation by including the secondary term. Regardless, we caution that the significant MZR evolution in EAGLE makes the previously derived values of $\log\Gamma$ 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 $\Delta Z_{\rm gas}$ versus $\Delta Z_{*}$ when $\langle Z_{\rm gas}\rangle>\langle Z_{*}\rangle$ and is a constant offset.

We note that there appears to be some tension with the T18 model, which (effectively) predicts that $\log\Gamma=-1$ for TNG. Recall that we neglect the return of mass and metals back into the ISM from stars dying. This removes a $(1-R)$ term on the sSFR, where $R$ is the mass return fraction from stars (identically, a $1/(1-R)$ scaling on $\tau_{\rm SF}$). If we instead included this return fraction, assuming that is a constant over $\tau_{\rm c}$, it is equivalent to just scaling $\Gamma$ by $(1-R)$. This does indeed shift us closer to agreement with T18, the extent of the agreement, however, depends on the adopted value of $R$. In order to produce a slope that corresponds to a $\log\Gamma=-1$, we must assume stellar mass return values ranging from $\sim 66-97$%. This value is much larger than current understanding predicts (i.e., $\lesssim$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${}_{\rm 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 $\gg$ 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 time-scale 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 M_{*}\gtrsim 10.2\log M_{\odot}$ 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 $M_{*}-Z_{\rm gas}-{\rm SFR}$ relation at the highest masses (as in Yates et al. 2012), and are likely contributing to the inversion of the $M_{*}-Z_{*}-{\rm 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.

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 $M_{*}-Z_{*}-{\rm 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 $\tau_{\rm 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.