Explaining the scatter in the galaxy mass-metallicity relation with gas flows
MMon. Not. R. Astron. Soc. , 1–10 (2021) Printed 28 January 2021 (MN L A TEX style file v2.2)
Explaining the scatter in the galaxy mass-metallicity relation withgas flows
Maria L. van Loon (cid:63) , Peter D. Mitchell and Joop Schaye Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, the Netherlands
28 January 2021
ABSTRACT
The physical origin of the scatter in the relation between galaxy stellar mass and the metallic-ity of the interstellar medium, i.e. the Mass-Metallicity Relation (MZR), reflects the relativeimportance of key processes in galaxy evolution. The
EAGLE cosmological hydrodynamicalsimulation is used to investigate the correlations between the residuals of the MZR and theresiduals of the relations between stellar mass and, respectively, specific inflow, outflow andstar formation rate as well as the gas fraction for central galaxies. At low redshift, all theseresiduals are found to be anti-correlated with the residuals of the MZR for M (cid:63) / M (cid:12) (cid:46) .The correlations between the residuals of the MZR and the residuals of the other relationswith mass are interrelated, but we find that gas fraction, specific inflow rate and specific out-flow rate all have at least some independent influence on the scatter of the MZR. We find that,while for M (cid:63) / M (cid:12) > . the specific mass of the nuclear black hole is most important, for M (cid:63) / M (cid:12) (cid:46) . gas fraction and specific inflow rate are the variables that correlate moststrongly with the MZR scatter. The timescales involved in the residual correlations and thetime that galaxies stay above the MZR are revealed to be a few Gyr. However, most galaxiesthat are below the MZR at z = 0 have been below the MZR throughout their lifetimes. Key words: galaxies: formation – galaxies: evolution – metallicity
The processes involved in galaxy formation and evolution are cou-pled and take place over a range of length and time scales, whichmakes it challenging to develop a full understanding. Hydrodynam-ical simulations are an established tool to assist in this, and mod-ern simulations of cosmological volumes now reproduce many ob-servational constraints well enough (e.g. Vogelsberger et al. 2014;Schaye et al. 2015; Dav´e et al. 2016; Dubois et al. 2016) that theycan be used with confidence to study aspects of galaxy evolutionthat are not readily observable, such as galactic gas accretion.Observationally, such processes can be linked indirectly toquantities that are observable, such as the metallicity of the inster-stellar medium (ISM). Gas-phase metallicity is influenced by pro-cesses that create or deplete heavy elements within galaxies, and byprocesses that increase or decrease the total gas content of galaxies(see Maiolino & Mannucci 2019 for a recent review). The metal-licity of the ISM is therefore influenced by stellar evolution via thedistribution of the metals created inside stars, by galactic outflowsthat carry metals and gas away from the ISM, and by gas accre-tion onto the ISM, which dilutes its metal abundance assuming that (cid:63)
E-mail: [email protected] the surrounding circum-galactic medium (CGM) is comparativelymetal poor.It is well established observationally that galaxy gas-phasemetallicity correlates with galaxy stellar mass (e.g. Tremonti et al.2004; Zahid et al. 2014; Curti et al. 2020), with a scaling relation-ship that is generally referred to as the Mass-Metallicity Relation(MZR). The MZR has a positive slope that is approximately con-stant for M (cid:63) / M (cid:12) (cid:46) . and flattens for higher masses. It isusually accepted that the positive slope of the MZR for lower stel-lar mass galaxies is determined by the higher efficiency of galacticoutflows at lower masses, due to their shallower gravitational po-tentials (e.g. Larson 1974).The scatter in the MZR is 0.1-0.2 dex, which is more than ex-pected from measurement errors (Maiolino & Mannucci 2019), andunderstanding this scatter offers the prospect of advancing in turnour understanding of the underlying network of inflows, outflows,and star formation within galaxies (e.g. Forbes et al. 2014; Finlator& Dav´e 2008; Guo et al. 2016; Lara-L´opez et al. 2019; Wang &Lilly 2020a; De Lucia et al. 2020a). Observationally, it has beenreported that metallicity anti-correlates with star formation rate orspecific star formation rate at fixed stellar mass (e.g. Ellison et al.2008; Mannucci et al. 2010; Curti et al. 2020), possibly forming a“fundamental” relation between the three quantities that is invari- © 2021 RAS a r X i v : . [ a s t r o - ph . GA ] J a n M. L. van Loon et al. ant with redshift (Lara-L´opez et al. 2010; Mannucci et al. 2010).Other observational studies do not find a residual correlation be-tween star formation rate and metallicity at fixed stellar mass how-ever (S´anchez et al. 2013 2019). Correlated residuals of the MZRwith atomic gas mass or gas fraction (Bothwell et al. 2013) andfor molecular hydrogen mass (Bothwell et al. 2016) have also beenreported. Generally speaking, theoretical models and cosmologicalsimulations predict qualitatively similar anti-correlations of metal-licity with star formation rate or gas mass at fixed stellar mass (e.g.Yates et al. 2012; Lilly et al. 2013; Lagos et al. 2016; De Rossi et al.2017; Torrey et al. 2019; De Lucia et al. 2020b).Theoretical work has also been used to study the connectionbetwen the MZR and variables that are not readily observable. Forexample, using simple idealised galaxy evolution models, Forbeset al. (2014) and Wang & Lilly (2020a) study how fluctuations ingas inflows rates can drive the MZR scatter, emphasing the im-portance of the relationship between the fluctuation timescale forinflows and the gas consumption timescale in the ISM. De Luciaet al. (2020b) use a more complex semi-analytic model of galaxyformation to study the physical drivers of the MZR relation, andarrive at the conclusion that gas accretion drives the MZR scatter,with galaxies below the MZR relation being driven by sustainedperiods of gaseous inflow. Using hydrodynamical simulations of acosmological volume, Torrey et al. (2019) arrive at similar conclu-sions by studying how galaxies evolve in metallicity, gas mass, andstar formation rate over finite time intervals.Here, we extend this body of work by using cosmological hy-drodynamical simulations of a representative volume to directlystudy the connection between the MZR and inflow rates, outflowrates, and star formation rates. Specifically, we use the
EAGLE sim-ulations, building on previous analyses of the MZR in these sim-ulations by Schaye et al. (2015), Crain et al. (2015), Lagos et al.(2016), De Rossi et al. (2017), Almeida & Dalla Vecchia (2018),Trayford & Schaye (2019) and Zenocratti et al. (2020). Our workcomplements that of more idealised modelling that generally infersthe properties of inflows and outflows based on observables, andalso complements semi-analytic modelling by self-consistently al-lowing for potential coupling between outflows and inflows (i.e.“preventative feedback”, e.g. Dav´e et al. 2012; Mo & Mao 2002;van de Voort et al. 2011).This study is organised as follows. Section 2 briefly describesthe
EAGLE simulations and how the various variables used in thisstudy have been measured and defined. Section 3 presents the re-sults of the study, including some individual examples, and the cor-relations between the residuals of the MZR and residuals of the re-lations between stellar mass and the other variables (specific inflow,outflow and star formation rates and the gas fraction). We explorethe stellar mass dependence (section 3.1), inter-variable coupling(section 3.2), the contribution to the MZR scatter (section 3.3), andthe relevant timescales (section 3.4). Section 4 provides a summaryof the main results.
We use data from the
EAGLE simulation project, which consistsof a suite of cosmological hydrodynamical simulations (Schayeet al. 2015; Crain et al. 2015), which have been made public(McAlpine et al. 2016).
EAGLE simulates representative cosmolog-ical volumes, including gas, stars, supermassive black holes anddark matter, using smoothed particle hydrodynamics (SPH) to solvethe equations of hydrodynamics with a modified version of the
GADGET -3 code (last described in Springel 2005; see Schaller et al.2015 for the employed formulation of SPH.) The simulation as-sumes a Λ CDM cosmological model and uses ’subgrid’ models forrelevant unresolved physics. The parameters for the subgrid mod-els governing stellar and AGN feedback were calibrated to fit theobserved galaxy stellar mass function at z ≈ , as well as the lo-cal disk galaxy size - stellar mass relation and the local black holemass - stellar mass relation.The EAGLE simulation suite consists of a set of “Reference”simulations with a common set of model parameters, as well asa number of variations. All
EAGLE measurements in this arti-cle are taken from the largest Reference model simulation, Ref-L100N1504, which simulates a (100 Mpc) box volume with dark matter and baryonic particles. The (initial) par-ticle masses in this simulation are . × M (cid:12) for gas and stars,and . × M (cid:12) for dark matter.Since EAGLE employs SPH, it is possible to track simulationparticles across time. We use this information to measure galac-tic inflow and outflow rates by determining how much gas entersor leaves the ISM of a given galaxy between two simulation snap-shots. The precise method used to measure inflow and outflow ratesis described in detail in Mitchell et al. (2020ab). We define the ISMas the combination of star-forming gas (meaning its density exceedsthe metallicity-dependent threshold introduced in Schaye 2004),and also non-star forming gas with density n H > .
01 cm − that islocated within 20% of the halo virial radius, approximately match-ing a selection of neutral atomic hydrogen (e.g. Rahmati et al.2013).To account for particles fluctuating across the ISM boundary,outflowing gas particles are selected by their radial velocity, ensur-ing that they have moved a significant distance. The instantaneousradial velocity and the difference in radius between two subsequentsnapshots, divided by the time elapsed between them, must bothexceed 0.25 times the maximum circular velocity of the halo. Par-ticles that leave the ISM without meeting these criteria are trackedthrough subsequent snapshots until they: a) meet the outflow cri-teria, or b) have re-entered the ISM (at which point they are notcounted towards the inflow rate), or c) until three halo dynamicaltimes have passed.Inflowing gas particles can be distinguished into three differ-ent types, as presented by Mitchell et al. (2020a). First, particles canbe accreted onto a galaxy for the first time. Second, particles can berecycled from a progenitor of the current galaxy. Last, particles canbe transferred from the ISM of another galaxy. These categories areadded together to create the total inflow onto a galaxy, which is thevariable used in this study. Note that a distinction is made betweeninflows and gas accretion due to galaxy mergers; we choose not toinclude the latter in our analysis (defining mergers as satellites witha maximum past halo mass greater than . × M (cid:12) , correspond-ing to 100 dark matter particles at fiducial EAGLE resolution).Inflow, outflow, and star formation rates are measured usingadjacent snapshots in a 200 snapshot grid. Inflowing (outflowing)gas mass, or formed stellar mass between snapshots is divided bythe time elapsed between snapshots, which is ≈ Myr at lowredshift and becomes smaller at higher redshifts (for the exact tem-poral spacing, see figure A1 in Mitchell et al. 2020b). In parts ofour analysis, the temporal spacing is adjusted by tracking galaxiesacross snapshots and integrating their inflow, outflow or star for-mation rates over longer timescales, which we fix to be a constantfraction of the halo dynamical time.This study focusses only on central galaxies, thereby simpli-fying the interpretation of the results by excluding environmen- © 2021 RAS, MNRAS , 1–10 catter in MZR and gas flows t (Gyr) ∆ l og ( x ) M / M fl = 10 . t (Gyr) ∆ l og ( x ) M / M fl = 10 . x = ˙ M in ( M ) x = ˙ M out ( M ) x = ˙ M ( M ) x = Z ISM ( M ) Figure 1.
Examples of the residual inflow rate ( ˙ M in ), outflow rate ( ˙ M out ),star formation rate ( ˙ M (cid:63) ) and metallicity of the interstellar medium ( Z ISM )of two individual galaxies tracked through time, with present-day masses M (cid:63) / M (cid:12) = 10 . (top) and M (cid:63) / M (cid:12) = 10 . (bottom). At z = 0 ,samples of galaxies are chosen based on stellar mass (Top: M (cid:63) / M (cid:12) =10 . ± . ; Bottom: M (cid:63) / M (cid:12) = 10 . ± . ). Residuals are taken withrespect to the medians for these same samples at each snapshot. The inflow,outflow and star formation rate show fluctuations across a wide range oftimescales. Metallicity anti-correlates with these variations but it fluctuatesless strongly. tal effects. Furthermore, galaxies with zero ISM gas particles areexcluded due to their undefined metallicity. Similarly, for figuresshowing data concerning inflows, galaxies with zero inflow ratehave been excluded (the same holds for outflows and star formationrate (SFR)). Apart from this selection, we select galaxies only bystellar mass and redshift. However, most figures show (with trans-parent and dashed lines) where the interpretation becomes uncer-tain due to a majority of galaxies ( > ) containing <
20 inflow-ing, outflowing, newly formed star particles, or <
200 star forminggas particles. Furthermore, we discard data points including fewerthan 10 galaxies, or where > of galaxies have <
10 inflow-ing, outflowing, newly formed star particles, or <
100 star forminggas particles, to reduce the noise in the plots. None of these choicesstrongly influence the results.
Fig. 1 shows two examples of typical star forming galaxies, whoseinflow rate, outflow rate, star formation rate and metallicity havebeen tracked though time. At z = 0 samples of the galaxies arecreated based on their stellar mass (bin width = 0.3 dex). Thesesamples are then used to determine the median values in each snap-shot, with respect to which the residuals of the variables are calcu-lated. Fig. 1 plots the evolution of there residuals for two individ-ual galaxies. It illustrates clearly that when the inflow, outflow andstar formation rates are higher than the median of the sample (i.e.the residual is positive), the metallicity tends to be lower (i.e. theresidual is negative) and vice versa. The variables fluctuate over arange of timescales spanning ≈ . Gyr to ≈ Gyr, where thelonger timescale fluctuations show a clear anti-correlation betweenmetallicity and the other variables. Inflow rates (blue line) gener-ally seem to be the first to rise or fall, after which other variablesreact, if these (anti-) correlations are interpreted causally. The shorttimescale fluctuations of inflow, outflow and star formation rate arelarge compared to those in the metallicity.Fig. 2 shows the median residuals of the relations betweenstellar mass and, respectively, the specific inflow rate, specificoutflow rate, gas fraction and specific star formation rate plottedagainst the residuals of the mass-metallicity relation for differ-ent stellar mass bins ( M (cid:63) / M (cid:12) = 10 . ± . (blue), M (cid:63) / M (cid:12) =10 . ± . (green), M (cid:63) / M (cid:12) = 10 . ± . (magenta)), at z < . .The ’specific’ in specific inflow/outflow/star formation rate meansthe variables are divided by the stellar mass of the galaxy, thus fac-toring out any residual stellar mass dependence remaining withinthe finite stellar mass bins.Focusing first on the stellar mass bin M (cid:63) / M (cid:12) = 10 . ± . ,the residual specific inflow rate, specific outflow rate, gas fractionand specific star formation rate all anti-correlate with the resid-ual metallicity (see Fig. 2, green lines in the top left, top right,bottom left, bottom right panels, respectively). When interpretedcausally, this suggests that the processes responsible for increasingthese variables generally lower the gas metallicity in a galaxy, witha possible lower limit to the effect (the correlations flatten at ex-tremely low values). Due to the coupling of the variables however(see sections 1 and 3.2), the correlations need to be investigatedfurther before understanding which physical processes cause them.A starting point for uncovering the physical origin of the MZRscatter can be a comparison between the correlations shown in Fig.2. A comparison of slopes reveals that the correlations of the resid-ual metallicity with the residual specific inflow and outflow ratematch closely in steepness and shape, whereas the residual correla-tions with the gas fraction and sSFR residual correlations are muchsteeper, the gas fraction residual correlation being the steepest. Thissupports the idea that the correlation with the gas fraction may becaused by a combination of the influence of multiple other vari-ables. Fig. 2 depicts the residual correlations for three stellar mass bins,of M (cid:63) / M (cid:12) = 10 . ± . , . ± . and . ± . , at z < . .The redshift evolution of the residual correlations has also been in-vestigated and is found to be small, as can be viewed in AppendixA. Considering the stellar mass dependence of the residual correla-tions, Fig. 2 shows that the lower ( M (cid:63) / M (cid:12) = 10 . ± . , blue linesin Fig. 2) and intermediate ( M (cid:63) / M (cid:12) = 10 . ± . , green lines in © 2021 RAS, MNRAS , 1–10 M. L. van Loon et al. ∆ log ( ˙ M in ( M ) /M ) ∆ l og ( Z I S M ( M )) z < . M / M fl = 10 . ± . M / M fl = 10 . ± . M / M fl = 10 . ± . ∆ log ( ˙ M out ( M ) /M ) ∆log ( f SF , gas ( M )) ∆ l og ( Z I S M ( M )) ∆ log ( ˙ M ( M ) /M ) Figure 2.
The relations between the median residual metallicity of the interstellar medium and residual specific inflow rate (top left), residual specific outflowrate (top right), residual gas fraction (bottom left) and residual sSFR (bottom right) for central galaxies with M (cid:63) / M (cid:12) = 10 . ± . (blue lines), M (cid:63) / M (cid:12) =10 . ± . (green lines), M (cid:63) / M (cid:12) = 10 . ± . (magenta lines), at z < . . The residuals are taken w.r.t. the medians as a function of stellar mass. Adashed and transparent line style indicates regions where more than 50% of galaxies in the bin contain <
20 inflowing (top left panel), outflowing (top rightpanel), newly formed star particles (bottom right panel) or <
200 star forming gas particles (bottom left panel). The lines are cut off when more than 50% ofgalaxies in the bin contain <
10 inflowing (top left panel), outflowing (top right panel), newly formed star particles (bottom right panel), or <
100 star forminggas particles (bottom left panel). or when bins contain fewer than 10 galaxies. For stellar masses < . M (cid:12) , the shapes and slopes are comparable acrossthe residual correlations. For stellar masses > . M (cid:12) , they deviate, forming two groups: specific inflow and outflow rate (top), and gas fraction and sSFR(bottom). Fig. 2) stellar mass bins of all four residual correlations show acomparable progression, varying only slightly in shape and slope.The higher stellar mass bin ( M (cid:63) / M (cid:12) = 10 . ± . ) deviatesin shape from the lower two, especially for the residual trends withgas fraction and star formation rate, which change in sign to clearpositive correlations, consistent with the previous analysis of EA - GLE by De Rossi et al. (2017) (see also Yates et al. 2012 for a sim-ilar result from a semi-analytic model). The residual correlationswith inflow and outflow rate deviate in a more complex manner, ap-pearing to show a positive correlation for negative x-axis residuals,but retaining an overall anti-correlation for positive x-axis residu-als. By comparing
EAGLE simulations with and without AGN feedback, De Rossi et al. (2017) find that the inversion of the depen-dence on the gas fraction and star formation rate residuals at highstellar mass can be attributed to the AGN feedback. We explore thisaspect further in Fig. 3, which shows the correlation between resid-ual ISM metallicity and residual specific black hole mass. Fig. 3supports the idea that for M (cid:63) / M (cid:12) > . the shapes of theresidual correlations (Fig. 2) are likely caused by the influence ofAGN feedback on larger mass galaxies. For the highest mass bin,the correlation between residual metallicity and residual specificblack hole mass is steeper and negative, suggesting that high blackhole masses correlate strongly with a low metallicity. For the lowerstellar mass bins, the residual correlations are almost flat until ex-treme positive residuals in black hole mass, indicating a possible © 2021 RAS, MNRAS , 1–10 catter in MZR and gas flows ∆ log ( M BH ( M ) /M ) ∆ l og ( Z I S M ( M )) z < . M / M fl = 10 . ± . M / M fl = 10 . ± . M / M fl = 10 . ± . Figure 3.
The median residual of the relation between ISM metallicity andstellar mass as a function of the residual between specific black hole massand stellar mass for central galaxies with M (cid:63) / M (cid:12) = 10 . ± . (blueline), M (cid:63) / M (cid:12) = 10 . ± . (green line), M (cid:63) / M (cid:12) = 10 . ± . (ma-genta line), at z < . . Only bins containing more than 10 galaxies areshown. The anti-correlation between residual metallicity and residual spe-cific black hole mass is only significant for relatively high black hole massesand is steeper for the highest mass bin. causal connection whereby black holes above a threshold mass areable to influence the metallicity of the galaxy. While Fig. 2 shows clear correlations between metallicity and in-flow, outflow, and star formation rates, and also gas fraction, the lat-ter four variables can of course also be correlated with each other,rather than each directly driving a causal connection. One way toprobe the coupling of variables, is by correcting the residual corre-lations for others, thus making it possible to judge the influence ofeach variable on the other residual correlations. By comparing allpossibilities, insight is gained into which variable may be most fun-damental. In Fig. 4, each panel shows the residual correlation of oneof the four variables (black lines, M (cid:63) / M (cid:12) = 10 . ± . , z < . ;there curves are identical to the green lines in Fig. 2), and the resid-ual correlation after correcting for one of the other variables (blue:specific inflow rate, green: specific outflow rate, magenta: specificstar formation rate, cyan: gas fraction).The correction is calculated by taking the residuals of theresidual correlations. First, the median residuals of the relations be-tween stellar mass and, respectively, metallicity and (e.g.) specificoutflow rate are calculated (respectively ∆ log ( Z ISM ( M (cid:63) )) and ∆ log ( ˙ M out ( M (cid:63) ) /M (cid:63) ) , which are plotted in Fig. 2). Then, theresiduals of the relation between these two variables are calculated,creating a new variable: ∆ (cid:16) ∆ log ( Z ISM ( M (cid:63) )) (cid:104) ∆ log ( ˙ M out ( M (cid:63) ) /M (cid:63) ) (cid:105)(cid:17) which is then plotted against the residuals of the relation betweenstellar mass and another variable (e.g. specific inflow rate) to createFig. 4.Fig. 4 shows that each residual relation flattens consider-ably after correcting for another residual correlation (compare the coloured lines to the black lines), indicating that none of the vari-ables are fully independent. When studying the bottom right panelin Fig. 4, it becomes clear that the effect of the sSFR on the residualmetallicity is completely determined by the gas fraction (cyan line)and to a lesser extent by the specific inflow and outflow rates (blueand green lines, respectively), with the exception of extreme values( ∆ log (sSFR) < − . , ∆ log (sSFR) > . ). The relationsbetween residual metallicity and, respectively, residual specific in-flow rate, residual specific outflow rate and residual gas fractionflatten considerably when corrected for the other variables. How-ever, there is no significant difference between them.Combining all the information in Fig. 4, we can conclude thatthe gas fraction, specific inflow, outflow and star formation rate allplay an independent role in determining the metallicity of a galaxy,but that they are also coupled. The independent influence of starformation rate is least important. The gas fraction can be seen phys-ically as a time integrated aggregate of the three other processes,and will therefore naturally correlate strongly with the metallicityand affect other correlations. We now consider directly how much of the MZR scatter can be ex-plained by each of the variables considered, using a similar methodas used by Matthee & Schaye (2019) to explore the origin of thescatter in star formation rates at fixed stellar mass in the
EAGLE simulation. Fig. 5 shows how much of the original MZR scatter(black line) remains after correcting for the residual correlationwith another variable, using the procedure described in section 3.2.Fig. 5 then depicts the remaining (1 σ ) scatter. The residual correla-tions were recalculated for different stellar mass bins (width = 0.15dex) to account for the stellar mass dependencies of the residualcorrelations and include only central galaxies at z < . .Fig. 5 shows that up to ≈ . dex of the ≈ . dex scatterin the MZR can be explained by correcting for one of the variables(e.g. the gas fraction, cyan line). At low stellar masses ( M (cid:63) / M (cid:12) < . ) correcting for the gas fraction (cyan line) or the specific in-flow rate (blue line) reduces the scatter the most, marking these twovariables as most influential. When correcting for the sSFR (ma-genta line) or specific outflow rate (green line) residual correlation,only ≈ . dex in scatter is explained at M (cid:63) / M (cid:12) < . , com-pared to the ≈ . dex at M (cid:63) / M (cid:12) ≈ . . At M (cid:63) / M (cid:12) > the amount of explained scatter diminishes for all variables, withthe exception of specific black hole mass (yellow line), which onlyaccounts for scatter at higher stellar masses and is the most impor-tant variable for M (cid:63) / M (cid:12) > . . This can be explained by theeffects of AGN feedback on the metallicity of high-mass galaxies.Fig. 5 provides strong evidence supporting the idea of inflowsof gas having a larger influence on the scatter in the MZR thanoutflows or star formation for M (cid:63) / M (cid:12) (cid:46) . The gas fractionreduces the scatter in the MZR the most, which supports the ideathat the physical origin of the scatter lies in more than one variable,under the assumption that the gas fraction reflects the influence ofmultiple variables (inflow, outflow and star formation rates). Notethat this does not exclude the possibility of a direct influence of thegas fraction on the metallicity (higher gas surface densities couldchange the effectiveness of feedback, for example). Fig. 5 shows that not all of the scatter in the MZR can be explainedby just one of the individual variables considered. We remind the © 2021 RAS, MNRAS , 1–10
M. L. van Loon et al. ∆ log ( ˙ M in ( M ) /M ) ∆ l og ( Z I S M ( x )) x = ˙ M in /M [ Gyr − ] x = ˙ M in ( ˙ M out ) /Mx = ˙ M in ( ˙ M ) /Mx = ˙ M in ( f SF , gas ) /M ∆ log ( ˙ M out ( M ) /M ) M / M fl = 10 . ± . , z < . x = ˙ M out ( M ) /M [ Gyr − ] x = ˙ M out ( M , ˙ M in ) /Mx = ˙ M out ( M , ˙ M /M ) x = ˙ M out ( M , f SF , gas ) /M ∆log ( f SF , gas ( M )) ∆ l og ( Z I S M ( x )) x = f SF , gas ( M ) x = f SF , gas ( M , ˙ M in ) x = f SF , gas ( M , ˙ M out ) x = f SF , gas ( M , ˙ M ) ∆ log ( ˙ M ( M ) /M ) x = ˙ M ( M ) /M [ Gyr − ] x = ˙ M ( M , ˙ M in ) /Mx = ˙ M ( M , ˙ M out ) /Mx = ˙ M ( M , f SF , gas ) /M Figure 4.
The median residual of the gas mass-metallicity relation as a function of the residual specific inflow rate (top left), specific outflow rate (top right),gas fraction (top left), and sSFR (bottom right) for central galaxies with M (cid:63) / M (cid:12) = 10 . ± . , at z < . . Black:
The residuals of the relations with stellarmass.
Blue:
The residual correlation corrected for the relation between residual metallicity and residual specific inflow rate.
Green:
The residual correlationcorrected for the relation between residual metallicity and residual specific outflow rate.
Magenta:
The residual correlation corrected for the relation betweenresidual metallicity and residual specific star formation rate.
Cyan:
The residual correlation corrected for the relation between residual metallicity and residualthe gas fraction. Line styles are as in Fig.2. Each residual correlation flattens significantly when corrected for the dependence of metallicity on one of the othervariables implying that all variables are coupled. reader however that until now we have presented inflow and out-flow rates that are measured over a ≈
120 Myr time interval. Re-calling from Fig. 1 that the metallicity shows long timescale ( (cid:29) Gyr) fluctuations, we now investigate the timescales involved in theresidual correlations by averaging the variables over longer time in-tervals.Fig. 6 illustrates that the correlations with the residual spe-cific inflow, outflow and star formation rate are likely dominatedby fluctuations on long timescales. The figure shows the resid-ual correlations for measurements of the variables over differenttime intervals, which smooth the fluctuations that occur on shortertimescales. Galaxies selected at z < . are tracked back in time, across the desired time interval, and their inflows/outflows/formedstars are added together.When the time interval increases, smoothing the shortertimescale fluctuations, the relations with the residual specific in-flow and outflow rates steepen somewhat (top panels). This slightincrease of steepness is due to the decrease of ’noise’ (shortertimescale fluctuations in the variables), indicating that fluctuationson longer timescales are likely the main cause of the correlations.The residual correlation for specific outflow rate saturates at a timeinterval of > Gyr, suggesting that fluctuations on this timescalecan no longer be considered to be noise. The inflow and outflowpanels indicate that short timescales are not dominating the resid- © 2021 RAS, MNRAS , 1–10 catter in MZR and gas flows log ( M / M fl ) σ ( ∆ l og ( Z I S M ( M ))) z < . σ (∆ log ( Z ISM ( M ))) σ (∆ log ( Z ISM ( M , ˙ M in /M ))) σ (∆ log ( Z ISM ( M , ˙ M out /M ))) σ (∆ log ( Z ISM ( M , ˙ M /M ))) σ (∆ log ( Z ISM ( M , f SF , gas )) σ (∆ log ( Z ISM ( M , M BH /M )) Figure 5.
The mass dependence of the (1 σ ) scatter in the median MZR(black), when corrected for the correlations between residual metallic-ity and residual specific inflow rate (blue), residual specific outflow rate(green), residual sSFR (magenta), residual gas fraction (cyan) and residualspecific black hole mass (yellow) for central galaxies, at z < . (where theresiduals are taken w.r.t. the median relations with stellar mass). At low stel-lar masses ( M (cid:63) / M (cid:12) < . ), the gas fraction and specific inflow ratereduce the scatter the most. At high stellar masses ( M (cid:63) / M (cid:12) > . ),the specific black hole mass reduces the scatter the most. ual correlations, as the correlations would have flattened when thefluctuations over these timescales are smoothed.The relation with the residual specific star formation rate doesnot change when measured over different time intervals, indicatingthat for this correlation, short timescales are not the cause. Whythe residual correlation does not become steeper like the others isunclear.We next consider the timescales over which galaxies remainabove/below the MZR, independently of the other variables con-sidered. The top panel of Fig. 7 shows the ’memory’ of the scatterin the MZR by tracking a sample of galaxies that lie above (or be-low) the MZR at z = 0 to higher redshifts. At each cosmic timethe fraction of galaxies above (or below) the MZR at that time isdetermined, until a fraction smaller than 0.2 of the original galax-ies remain. The bottom panel of Fig. 7 shows instead the medianresidual metallicities of these galaxies at the different redshifts.If the metallicity in the ISM of galaxies had no ’memory’,a fraction of 0.5 (grey line, top panel) and a median residual of0 (grey line, bottom panel) would be expected, as the distributionwould be random at z > . Fig. 7 shows clearly that the metallic-ities of the galaxies do have a memory, as the scenario mentionedabove does not hold.Galaxies above the MZR have on average been above theMZR for (cid:38) , if their current mass is M (cid:63) / M (cid:12) (cid:38) .Galaxies with present-day mass M (cid:63) / M (cid:12) = 10 . ± . have beenabove the MZR for on average ≈
12 Gyr , i.e., their entire lifetime.Galaxies below the MZR have also been below the MZR for onaverage ≈
12 Gyr . This idea is consistent with our hypothesis thatprocesses like star formation, inflows and outflows determine themetallicity, since these also have trends over longer time scales (seefigures 1 and 6). Furthermore, the long timescales involved suggestthat the metallicity is driven by the age of the halo formation it- self. When comparing to the work by Matthee & Schaye (2019)on the memory of the SFR - stellar mass sequence in
EAGLE , wesee similar timescales ( ≈
10 Gyr ) appear, which they showed areconnected to halo formation. Our results are also broadly speakingconsistent with the analysis of the relevant timescales for MZR andstar formation rate evolution in the Illustris-TNG simulations, aspresented in Torrey et al. (2018), who find that both fluctuate overlong timescales.The asymmetry between above and below the MZR may bedue to the fact that galaxies below the MZR have a higher SFR,resulting in a more rapid decrease in stellar mass as we travel backin time, creating a difference between the two samples. Why thisdifference would result in a longer timescale is not immediatelyobvious however.Together, Fig. 6 and Fig. 7 show that the correlations betweenresidual metallicity and residual specific inflow, outflow and starformation rate are consistent with being driven by fluctuations overlonger timescales, coinciding with the typical timescales that galax-ies have been above the MZR. An extension of our analysis wouldbe to attempt to decompose the fluctuation spectrum explicitly be-tween different timescales, either using Fourier analysis (e.g. Wang& Lilly 2020b; Iyer et al. 2020) or principal component analysis(Matthee & Schaye 2019). We leave this as a possible avenue forfuture work.
We have presented our findings concerning the physical originof the scatter in the Mass-Metallicity Relation (MZR) of centralgalaxies in the
EAGLE simulation and the roles of the gas fraction,specific inflow, outflow and star formation rates. By studying theresidual correlations between ISM metallicity and these variables,comparing them, correcting them for one another and studying thetimescales involved, it has become clear that although all of thesevariables are coupled, they also have an independent impact on thescatter in the MZR. The gas fraction, which is most likely an aggre-gate of the other variables, is the variable correlating best with thescatter in the MZR and the inflow rate seems to play a larger rolethan the outflow or star formation rates. This section summarisesthe main results.(i) There are anti-correlations between the residuals of the MZR,and the residuals of the relations between stellar mass and, re-spectively, the gas fraction, specific inflow, outflow and star for-mation rate, at M (cid:63) / M (cid:12) (cid:46) (see Fig. 2). These residualcorrelations are comparable in shape for M (cid:63) / M (cid:12) (cid:46) . For M (cid:63) / M (cid:12) (cid:38) the shapes deviate compared to lower mass bins:the correlations between residual metallicity and residual sSFR andgas fraction change in sign, showing a positive correlation; the cor-relations between residual metallicity and residual specific inflowand outflow rate also deviate in shape, showing a positive corre-lation at negative x-axis residuals, i.e. for relatively low flow rates,and an anti-correlation at positive x-axis residuals, i.e. for relativelyhigh flow rates.(ii) The residual specific black hole mass (i.e. log ( M BH ( M (cid:63) ) /M (cid:63) ) ) anti-correlates with residual metal-licity at M (cid:63) / M (cid:12) (cid:38) and for lower stellar masses also forrelatively high black hole masses (see Fig. 3). This suggests that athigh black hole masses, AGN feedback strongly influences the gasmetallicity of galaxies.(iii) The relations between the residuals of the MZR and theresidual gas fraction, residual specific inflow, outflow and star for- © 2021 RAS, MNRAS , 1–10 M. L. van Loon et al. ∆ log ( ˙ M in ( M ) /M ) ∆ l og ( Z I S M ( M )) M / M fl = 10 . ± . , z < . ∆ log ( ˙ M out ( M ) /M ) ∆ log ( ˙ M ( M ) /M ) ∆ l og ( Z I S M ( M )) δt = 0 .
12 Gyr δt = 0 .
24 Gyr δt = 0 . δt = 1 .
08 Gyr δt = 2 . δt = 3 .
02 Gyr δt = 5 . Figure 6.
The relations between the median residual of the median MZR and the residuals of the relations with stellar mass of, respectively, specific inflow rate(top left), specific outflow rate (top right) and sSFR (bottom left) for central galaxies with M (cid:63) / M (cid:12) = 10 . ± . , at z < . . Different colours correspondto variables averaged over different time intervals. Line styles are as in Fig.2. The specific inflow and outflow rate residual relations steepen when the timeinterval increases, indicating the importance of longer timescales. mation rates flatten considerably when corrected for the depen-dence of metallicity on any of the other variables (see Fig. 4). Therelation between residual metallicity and residual specific star for-mation rate (sSFR) is completely determined by the gas fraction,and almost completely by specific inflow rate and specific outflowrate for non-extreme values ( − . < ∆ log (sSFR) < . ). Thissuggests that all variables are coupled, but that they have some in-dependent influence on the metallicity, of which the influence ofthe star formation rate is least strong.(iv) When correcting the scatter in the MZR for the correlationsbetween residual metallicity and residual specific inflow rate, resid-ual specific outflow rate, residual specific star formation rate orthe residual gas fraction, and measuring the remaining scatter, thescatter is reduced most (from ≈ . dex down to ≈ . dex)by the gas fraction, and similarly by the specific inflow rate, at M (cid:63) / M (cid:12) < . (see Fig. 5). Accounting for the specific out-flow or star formation rate reduces the scatter by ≈ . dex at M (cid:63) / M (cid:12) < . , compared to ≈ . dex at M (cid:63) / M (cid:12) ≈ . .At high stellar masses ( M (cid:63) / M (cid:12) > . ) correcting for the spe-cific black hole mass reduces the scatter whereas the contributionof the other variables diminishes.(v) When averaged over longer time intervals, the anti-correlations of the residual metallicity with the residual specific in-flow and outflow rates steepen ( M (cid:63) / M (cid:12) = 10 . ± . , z < . ),indicating that the fluctuations on longer timescales are responsiblefor the residual correlations (see Fig. 6). The relation with the resid-ual specific star formation rate ( M (cid:63) / M (cid:12) = 10 . ± . , z < . )does not change when measured over longer time intervals, indicat-ing that fluctuations on short timescales are not responsible for thisresidual correlation. © 2021 RAS, MNRAS , 1–10 catter in MZR and gas flows t / Gyr N u m b e r a b o v e ( b e l o w ) M Z R ( z ) N u m b e r a b o v e ( b e l o w ) M Z R ( z = ) M / M fl = 10 . ± . M / M fl = 10 . ± . M / M fl = 10 . ± . AboveBelow t /
Gyr ∆ l og ( Z ( M )) z Figure 7.
Top:
The fraction of galaxies above (solid line style) or below(dashed line style) the MZR as a function of cosmic time, tracking the mainprogenitors of the sample of galaxies that are above or below the MZR at z = 0 ( t = 13 . ). The grey line at y = 0 . indicates what the frac-tion would be if the distribution were random within the samples at eachredshift. Bottom:
The median residual metallicity of galaxies in a sampleabove (solid line style) or below (dashed line style) the MZR at z = 0 as a function of cosmic times. The grey line at y = 0 indicates a ran-dom median residual at each redshift. In both panels the stellar mass bins(blue: M (cid:63) / M (cid:12) = 10 . ± . , green: M (cid:63) / M (cid:12) = 10 . ± . , magenta: M (cid:63) / M (cid:12) = 10 . ± . ) refer to masses at z = 0 . The sample above theMZR has been above the MZR for (cid:38) . The sample below the MZRhas been below the MZR for ≈
12 Gyr . (vi) When galaxies with M (cid:63) / M (cid:12) (cid:38) above (below) theMZR at z = 0 are tracked back in time, these have generally beenabove (below) the MZR for ≈ or ≈
12 Gyr , respectively (seeFig. 7). Lower-mass galaxies ( M (cid:63) / M (cid:12) = 10 . ± . ) have gener-ally been on one side of the MZR for ≈
12 Gyr .Summarising, star formation, inflows and outflows of gas arecoupled processes that all seem to have an independent impact onthe metallicity of the ISM, of which the inflows have the most influ-ence. Consistent with the analyis of the Illustris-TNG simulationspresented in Torrey et al. (2018), our analysis of the
EAGLE sim-ulations indicates that important timescales for these processes are of the order of a few
Gyr or longer, which is reflected in how longgalaxies stay above/below the MZR, on average. There timescalesare also consistent with the recent inference presented in Wang& Lilly (2020a), based on the residuals measured in observationaldata.Overall, we find that the gas fraction is still the variable thatcorrelates best with the metallicity (see also Lagos et al. 2016;De Rossi et al. 2017). We expect that this is primarily because thegas fraction reflects the inflows, outflows and star formation overa range of timescales. This reinforces the general notion that theMZR and its scatter represents an important diagnostic for a rangeof processes that are critical for galaxy evolution, but which aredifficult to observe directly.
ACKNOWLEDGEMENTS
DATA AVAILABILITY
The EAGLE simulation data are publicly available fromhttp://icc.dur.ac.uk/Eagle/database.php, as described in McAlpineet al. (2016). The processed data underlying this work will beshared on reasonable request to the corresponding author.
REFERENCES
Almeida J. S., Dalla Vecchia C., 2018, The Astrophysical Journal,859, 109Bothwell M. S., Maiolino R., Cicone C., Peng Y., Wagg J., 2016,A&A, 595, A48Bothwell M. S., Maiolino R., Kennicutt R., Cresci G., MannucciF., Marconi A., Cicone C., 2013, MNRAS, 433, 1425Crain R. A. et al., 2015, Monthly Notices of the Royal Astronom-ical Society, 450, 1937Curti M., Mannucci F., Cresci G., Maiolino R., 2020, MNRAS,491, 944Dav´e R., Finlator K., Oppenheimer B. D., 2012, MNRAS, 421, 98Dav´e R., Thompson R., Hopkins P. F., 2016, MNRAS, 462, 3265De Lucia G., Xie L., Fontanot F., Hirschmann M., 2020a, MN-RAS, 498, 3215De Lucia G., Xie L., Fontanot F., Hirschmann M., 2020b, MN-RAS, 498, 3215De Rossi M. E., Bower R. G., Font A. S., Schaye J., Theuns T.,2017, Monthly Notices of the Royal Astronomical Society, 472,3354–3377Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., WelkerC., Volonteri M., 2016, MNRAS, 463, 3948Ellison S. L., Patton D. R., Simard L., McConnachie A. W., 2008,ApJ, 672, L107Finlator K., Dav´e R., 2008, Monthly Notices of the Royal Astro-nomical Society, 385, 2181 © 2021 RAS, MNRAS , 1–10 M. L. van Loon et al.
Forbes J. C., Krumholz M. R., Burkert A., Dekel A., 2014, MN-RAS, 443, 168Guo Y. et al., 2016, The Astrophysical Journal, 822, 103Iyer K. G. et al., 2020, MNRAS, 498, 430Lagos C. d. P. et al., 2016, MNRAS, 459, 2632Lara-L´opez M. A. et al., 2010, A&A, 521, L53Lara-L´opez M. A., De Rossi M. E., Pilyugin L. S., Gallazzi A.,Hughes T. M., Zinchenko I. A., 2019, MNRAS, 490, 868Larson R. B., 1974, MNRAS, 169, 229Lilly S. J., Carollo C. M., Pipino A., Renzini A., Peng Y., 2013,ApJ, 772, 119Maiolino R., Mannucci F., 2019, The Astronomy and Astro-physics Review, 27Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A.,2010, MNRAS, 408, 2115Matthee J., Schaye J., 2019, Monthly Notices of the Royal Astro-nomical Society, 484, 915–932McAlpine S. et al., 2016, Astronomy and Computing, 15, 72Mitchell P. D., Schaye J., Bower R. G., 2020a, Galactic inflow andwind recycling rates in the eagle simulationsMitchell P. D., Schaye J., Bower R. G., Crain R. A., 2020b,Monthly Notices of the Royal Astronomical Society, 494,3971–3997Mo H. J., Mao S., 2002, MNRAS, 333, 768Rahmati A., Pawlik A. H., Raiˇcevi´c M., Schaye J., 2013, MNRAS,430, 2427S´anchez S. F. et al., 2019, MNRAS, 484, 3042S´anchez S. F. et al., 2013, A&A, 554, A58Schaller M., Dalla Vecchia C., Schaye J., Bower R. G., Theuns T.,Crain R. A., Furlong M., McCarthy I. G., 2015, MNRAS, 454,2277Schaye J., 2004, The Astrophysical Journal, 609, 667Schaye J. et al., 2015, MNRAS, 446, 521Springel V., 2005, MNRAS, 364, 1105Torrey P. et al., 2018, MNRAS, 477, L16Torrey P. et al., 2019, MNRAS, 484, 5587Trayford J. W., Schaye J., 2019, Monthly Notices of the RoyalAstronomical Society, 485, 5715Tremonti C. A. et al., 2004, ApJ, 613, 898van de Voort F., Schaye J., Booth C. M., Haas M. R., Dalla Vec-chia C., 2011, MNRAS, 414, 2458Vogelsberger M. et al., 2014, MNRAS, 444, 1518Wang E., Lilly S. J., 2020a, arXiv e-prints, arXiv:2009.01935Wang E., Lilly S. J., 2020b, ApJ, 895, 25Yates R. M., Kauffmann G., Guo Q., 2012, MNRAS, 422, 215Zahid H. J., Dima G. I., Kudritzki R.-P., Kewley L. J., Geller M. J.,Hwang H. S., Silverman J. D., Kashino D., 2014, ApJ, 791, 130Zenocratti L. J., De Rossi M. E., Lara-L´opez M. A., Theuns T.,2020, MNRAS, 496, L33
APPENDIX A: REDSHIFT EVOLUTION
The redshift evolution varies between the specific inflow rate, spe-cific outflow rate, sSFR and gas fraction residual correlations (fig-ure A1). For the calculations of the specific inflow, outflow andstar formation rates different time intervals were used for differentredshifts, in approximation of 0.2 times the halo dynamical time,to account for the expansion of the universe. Both the specific in-flow rate and gas fraction residual correlations show an increasein steepness of the correlation at higher redshifts for galaxies with M (cid:63) / M (cid:12) = 10 . ± . . Contrarily, the specific outflow residual correlation shows no clear redshift evolution trend and the sSFRresidual correlation seems to flatten at higher redshift althoughthere is no smooth progression. However, overal the redshift evolu-tion is mild. © 2021 RAS, MNRAS , 1–10 catter in MZR and gas flows ∆ log ( ˙ M in ( M ) /M ) ∆ l og ( Z I S M ( M )) M / M fl = 10 . ± . ∆ log ( ˙ M out ( M ) /M ) z = [0 . , . z = [0 . , . z = [0 . , . z = [1 . , . z = [3 . , . ∆log ( f SF , gas ( M )) ∆ l og ( Z I S M ( M )) ∆ log ( ˙ M ( M ) /M ) Figure A1.
The correlations between residual metallicity and residual specific inflow rate (top left), residual specific outflow rate (top right) residual gasfraction (bottom left) and residual sSFR (bottom right) for galaxies with M (cid:63) / M (cid:12) = 10 . ± . , at different redshifts. A time interval of 0.24 Gyr for z < . ; 0.24 Gyr for z = [0 . , . ; 0.12 Gyr for z = [0 . , . , 0.085 Gyr for z = [1 . , . , 0.043 Gyr for z = [3 . , . between measurements isused, in approximation of 0.2 times the halo dynamical time at these redshifts. Line styles are as in Fig.2. The gas fraction and specific inflow rate residualcorrelations steepen with redshift while the sSFR residual correlation flattens with redshift. At z > . the specific outflow rate residual correlation shows noclear redshift evolution. Overall, the evolution is weak.© 2021 RAS, MNRAS000