## The mass-metallicity relation as a ruler for galaxy evolution: insights from the James Webb Space Telescope

A. Pallottini 1, ⋆ , A. Ferrara 1 , S. Gallerani 1 , L. Sommovigo 2 , S. Carniani 1 , L. Vallini 3 , M. Kohandel 1 , and G. Venturi 1

- 1 Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy
- 2 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
- 3 INAF-Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93 3, I-40129, Bologna, Italy /

Received July 31, 2024; accepted XX XX, XXXX

## ABSTRACT

Context. Galaxy evolution emerges from the balance between cosmic gas accretion, fueling star formation, and supernova feedback regulating metal enrichment of the interstellar medium. Hence, the relation between stellar mass ( M ⋆ ) and gas metallicity ( Zg ) is fundamental to understand the physics of galaxies. High-quality spectroscopic JWST data enable accurate measurements of both M ⋆ and Zg up to redshift z ≃ 10.

Aims. Our aims are to understand (i) the nature of the observed mass-metallicity relation (MZR), (ii) its connection with the star formation rate (SFR), (iii) the role played by SFR stochasticity (flickering), and (iv) how it is regulated by stellar feedback.

Methods. We compare the MZR obtained by the JADES, CEERS, and UNCOVER surveys, which comprise about 180 galaxies at z = ≃ 3 -10 with 10 M 6 ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 10 M , with ⊙ ≃ 200 simulated galaxies in the same mass range from the serra high-resolution ( ≃ 20 pc) suite of cosmological radiation-hydrodynamic simulations. To interpret the MZR, we develop a minimal, physically motivated, model for galaxy evolution that includes: cosmic accretion, possibly modulated with an amplitude A 100 on 100 Myr timescales; a time delay t d between SFR and supernova feedback; SN-driven outflows with a varying mass loading factor ϵ SN, that is normalized to the fire simulations predictions for ϵ SN = 1.

Results. Using our minimal model, we find the observed mean MZR is reproduced for relatively ine ffi cient outflows ( ϵ SN = 1 4), / in line with findings from JADES. Matching the observed MZR dispersion across the full stellar mass range requires a delay time t d = 20 Myr, in addition to a significant modulation ( A 100 = 1 3) of the accretion rate. Successful models are characterized by a / relatively low flickering ( σ SFR ≃ 0 2), corresponding to a metallicity dispersion of . σ Z ≃ 0 2. Such values are close but slightly lower . than predicted from serra ( σ SFR ≃ 0 24, . σ Z ≃ 0 3), clarifying why simulations show no clear trend for the MZR, and show some . tension with the observations especially at M ⋆ ≃ 10 10 M . ⊙

Conclusions. As the MZR is very sensitive to SFR stochasticity, models predicting high r.m.s. values ( σ SFR ≃ 0 5) result . in a 'chemical chaos' (i.e. σ Z ≃ 1 4), . virtually destroying the observed MZR. As a consequence, invoking a highly stochastic SFR ( σ SFR ≃ 0 8) to explain the overabundance of bright, super-early galaxies leads to inconsistencies with the observed MZR. .

Key words. Galaxies: star formation - evolution - high-redshift

## 1. Introduction

The baryon cycle plays a pivotal role in regulating the galaxy formation and evolution process through cosmic time. As cosmic gas accretes in the dark matter halo potential well, it can cool and eventually form stars; massive stars explode as supernovae, enriching the surrounding interstellar medium with metals, that can potentially be ejected from the galaxy. Therefore key information regarding the star formation and enrichment history is encoded in the relation between stellar mass ( M ⋆ ) and gas metallicity ( Zg ), the so called mass-metallicity relation (MZR), making it e ff ectively a standard ruler for galaxy evolution (see Maiolino &amp; Mannucci 2019, for a review).

Over the past 20 years, it has become possible to measure the MZR from z ≃ 0 (Tremonti et al. 2004; Kirby et al. 2013; Andrews &amp; Martini 2013) up to the beginning of cosmic noon ( z ≃ 3, Erb et al. 2006; Maiolino et al. 2008; Zahid et al. 2011; Sanders et al. 2018; Curti et al. 2020), mostly by using the oxygen to hydrogen abundance ratio [O H] as a proxy for / Zg . Lo-

⋆ andrea.pallottini@sns.it

cally, these observations have shown that Zg increases with stellar mass up to M ⋆ ≃ 10 10 M and beyond that point it flattens. ⊙ This can be understood as due to the fact that low mass galaxies are prone to metal ejection due to their shallower potential well, while massive galaxies can retain most of the metals they have produced (Ferrara 2008). In this view, galaxy evolution is guided by a 'bathtub' equilibrium (Bouché et al. 2010; Dekel &amp; Mandelker 2014) which determines the metallicity by balancing infall and outflows (Lilly et al. 2013).

Further, observations have highlighted that Zg and M ⋆ are tightly connected with the star formation rate (SFR) in the socalled fundamental mass-metallicity relation (FMR, Mannucci et al. 2010), which seems to be redshift-independent up to z ≃ 3 (Curti et al. 2020); other factors as galaxy size (Ellison et al. 2008) and molecular gas content (Bothwell et al. 2013) might also play some secondary role in connecting Zg to M ⋆ .

The MZR evolution at z &gt; ∼ 4 has been predicted via semianalytical models (Dayal et al. 2013; Zahid et al. 2014; Ucci et al. 2023) or cosmological simulations (Torrey et al. 2019; Liu &amp; Bromm 2020; Langan et al. 2020; Wilkins et al. 2023).

Most models predict almost no evolution or a very weak redshift dependence for z &gt; 3 (Ma et al. 2016; Kannan et al. 2022; Marszewski et al. 2024); some others suggest a break down of the relation at z &gt; ∼ 9 -12 (Pallottini et al. 2014; Sarmento et al. 2018).

Thanks to the exquisite spectroscopic capabilities of JWST, we can accurately infer Zg for galaxies well within the Epoch of Reionization (EoR, z &gt; 6). At present, surveys like JADES (Bunker et al. 2023) and CEERS (Finkelstein et al. 2023), have provided metallicity measurements at z ≃ 3 -10 (Curti et al. 2024a) and z ≃ 4 -10 (Nakajima et al. 2023), for a combined sample of about 170 galaxies in the 5 × 10 6 M ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 10 M ⊙ stellar mass range.

These JWST observations suggest that the MZR is already in place in the EoR, albeit downshifted by 0.5 dex with respect to the local one. The evolution in the z ≃ 3 -10 range is very mild, as relatively low mass ( M ⋆ ≃ 10 9 -10 M ) galaxies can be ⊙ chemically mature ( Zg ≃ 0 3 Z . ⊙ ) already at z ≃ 6. The MZR shows a r.m.s. dispersion σ Z = 0 3 dex, and possibly a hint of . flattening at z &gt; ∼ 6 (Curti et al. 2024a).

Furthermore, multiple line detections in the same target can be used to constrain abundance patterns of individual elements (Kobayashi &amp; Ferrara 2024; Curti et al. 2024b). Resolved observations are beginning to probe the presence of metallicity gradients up to z ≃ 8 (Venturi et al. 2024), which should provide deeper insights into the formation process of the galaxy and the history of the stellar mass buildup (see Vallini et al. 2024, for an ALMA perspective). Observations targeting lensed fields (Bezanson et al. 2022; Chemerynska et al. 2024b) are uncovering the metallicities of even the fainter galaxies (Chemerynska et al. 2024a), which is key in order to explore the infall outflow / interplay due to stellar feedback in the less massive objects. Importantly, using observations of lensed targets at z ≃ 3 -9, Morishita et al. (2024) suggest that the 0.5 dex shift of normalization and the higher scattering (0.3 dex) w.r.t. the local MZR can be caused by the high level of burstiness of the SFR expected in highz galaxies (Pallottini &amp; Ferrara 2023; Sun et al. 2023b).

Indeed, as we move from the local Universe to high z (see Madau &amp; Dickinson 2014; Dayal &amp; Ferrara 2018; Förster Schreiber &amp; Wuyts 2020, for reviews), galaxies are expected to become more bursty, as a consequence of the increase of specific SFR (sSFR, González et al. 2010; Stark et al. 2013; Smit et al. 2014) combined with a decrease in size (Shibuya et al. 2015), which can make the stellar feedback more e ff ective (Krumholz &amp;Burkhart 2016) in maintaining a higher level of turbulence (Simons et al. 2017; Genzel et al. 2017). With the recent availability of JWST observations (Roberts-Borsani et al. 2022; Castellano et al. 2022; Finkelstein &amp; Bagley 2022; Naidu et al. 2022; Treu et al. 2022; Adams et al. 2023; Atek et al. 2023; Donnan et al. 2023; Harikane et al. 2023; Santini et al. 2023), the time varying stochastic SFR behavior (in short, flickering or burstiness), has been invoked as a possible mechanism (Mason et al. 2023; Mirocha &amp; Furlanetto 2023; Shen et al. 2023; Sun et al. 2023b; Kravtsov &amp; Belokurov 2024) to explain the overabundance of bright galaxies at z &gt; ∼ 10 (see Dekel et al. 2023; Ferrara et al. 2023, for alternatives). However, some works predict that r.m.s. amplitude of the SFR variability falls short of explaining the said overabundance (Pallottini &amp; Ferrara 2023). Thus, it is unclear if flickering alone can explain the phenomenon, particularly since the required high level of variability (Muñoz et al. 2023) is apparently not seen in the data (Ciesla et al. 2024).

Studies of the SFR flickering have experienced a recent surge, as stochastic SFR variations can have an impact on our ability to observe distant galaxies (Sun et al. 2023a) and a ff ect the possibility of detecting Population III stars (Vikaeus et al. 2022; Katz et al. 2023). In addition, the flickering can modify the UV spectral slopes of highz galaxies and the escape fraction of ionizing photons (Gelli et al. 2024a), ultimately controlling the cosmic reionization history (Davies &amp; Furlanetto 2016; Nikoli´ c et al. 2024).

Analyzing the amplitude and temporal variation of SFR flickering provides unique insights into the various feedback processes that regulate the formation and evolution of early galaxies (Pallottini &amp; Ferrara 2023). On the one hand, high-amplitude and high-frequency flickering, on timescales shorter than the delay between SFR and SN explosions ( ≃ 20 -40 Myr) are invoked (Gelli et al. 2023, 2024b) to explain the observed abrupt quenching of SFR in ≃ 10 9 M galaxies already at ⊙ z ≃ 7 (Dome et al. 2024; Looser et al. 2024). On the other hand, the presence of disks observed in M ⋆ ≃ 10 9 -10 M galaxies up to ⊙ z ≃ 7 (Rowland et al. 2024; Fujimoto et al. 2024) implies a limited level of stochasticity, as otherwise the disk would be disrupted (Kohandel et al. 2024).

In the present work, we aim to clarify the e ff ects of SFR flickering and feedback regulation of highz galaxies by using the MZR as a standard ruler. The paper is structured as follows; in Sec. 2 we present the serra simulation suite and compare it with MZR observations from JADES, CEERS, and UNCOVER; in Sec. 3 we develop a minimal model to give a physical interpretation to the MZR; in Sec. 4 we discuss the discrepancies found between simulations and data in the light of the analytical model; conclusions are given in Sec. 5.

## 2. The mass-metallicity relation at high redshift

To study the MZR at high- , firstly we adopt the data from z serra , a suite of cosmological zoom-in simulations that follows the evolution of z &gt; ∼ 6 galaxies (Pallottini et al. 2022). In serra we use the code ramses (Teyssier 2002) to evolve DM, gas, and stars. We enable the module ramses rt -(Rosdahl et al. 2013) to track radiation on-the-fly by coupling it (Pallottini et al. 2019; Decataldo et al. 2019) with krome (Grassi et al. 2014), in order to account for non-equilibrium chemistry up to H2 formation (Pallottini et al. 2017a). Stars are formed using a Schmidt (1959)Kennicutt (1998) like relation based on H2. Once formed, stars produce UV radiation that ionizes the gas, photo-dissociates H2, and builds a radiation pressure in the gas. Massive stars generate stellar winds and explode as SNe. Depending on the kind of feedback, energy injection in the gas can be of thermal (subject to radiative losses) and or turbulent (later dissipated) nature / (Pallottini et al. 2017b). Adopting a ≃ 1 2 . × 10 4 M ( ⊙ ≃ 20 pc) gas mass (spatial) resolution, with music (Hahn &amp; Abel 2011) we initialize cosmological initial conditions 1 generated at z = 100, and follow the evolution of galaxies down to z = 6 (and z = 4 for a sub-sample, see Kohandel et al. 2024).

For this work, we select ≃ 200 central serra galaxies at z = 7 7 with stellar mass 10 M . 6 ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 10 M ; we compute ⊙ the gas metallicity as Zg = M Z g / Mg by considering the total gas ( Mg ) and metal ( MZ ) mass contained within 3 stellar e ff ective radii from the galaxy center, which is close to the gas half-mass radius . 2 Note that, for a more direct comparison with observation, we should compute line and continuum emission from

1 We adopt a Λ CDMmodel with vacuum, matter, and baryon densities in units of the critical density ΩΛ = 0 692, . Ω m = 0 308, . Ω = b 0 0481, . normalized Hubble constant h = H0 / (100 km s -1 Mpc -1 ) = 0 678, spec-. tral index n = 0 967, and . σ 8 = 0 826 (Planck Collaboration et al. 2014). . 2 Di ff erent assumptions on the integration radius produce only small changes in the inferred metallicity.

Fig. 1: The mass-metallicity relation (MZR) at high-redshift ( ). z Weplot the gas metallicity ( Zg ) as a function of stellar mass ( M ⋆ ) for the serra simulated galaxies (Pallottini et al. 2022, z ≃ 7) and the JWST observations from CEERS Nakajima et al. (2023, z ≃ 4 -10), JADES (Curti et al. 2024a, z ≃ 3 -6 and 6 -10), and UNCOVER (Chemerynska et al. 2024a, z ≃ 6 -8). To guide the eye, the observational fits from Nakajima et al. (2023) and Curti et al. (2024a) are also reported as lines with the same colors of corresponding data. As a reference, the log( Zg / Z ) ⊙ = 0 37 log( . M ⋆/ M ) ⊙ -4 3 fit from . fire -2 simulations (Marszewski et al. 2024, z = 5 -12) is shown, along with their binned data and scatter at z = 8. The right axis reports the logarithmic oxygen abundance [O H] assuming a solar composition (Asplund et al. / 2009).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-67de279d-9177-4057-be8d-61b94006c5cc.jpg)

serra galaxies; then, we could infer M ⋆ from synthetic SED fitting and combine line emission to prepare metallicity calibrators. Such an approach would reduce potential mismatches due to uncertainties in using di ff erent calibrators (Curti et al. 2017; Chemerynska et al. 2024a), e.g. to relax the single zone model assumption (Marconi et al. 2024) and correct for biases in the comparison (Cameron et al. 2023), also avoiding potential issues in the M ⋆ determinations (Narayanan et al. 2024). While this kind of forward modeling approach is very powerful, also allowing to prepare novel observational strategies (Zanella et al. 2021; Rizzo et al. 2022), here it hinders the possibility to compare with most of the other models, as the MZR is usually computed integrating metal, gas, and stellar masses; thus we avoid it in the present work.

The prediction for serra are plotted in Fig. 1, along with JWST data from Curti et al. (2024a, JADES sample, Bunker et al. 2023), Nakajima et al. (2023, CEERS sample, Finkelstein et al. 2023), and Chemerynska et al. (2024a, UNCOVER sample, Bezanson et al. 2022), for a total of about 180 galaxies at z ≃ 3 -10 with 10 6 M ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 10 M . The bulk of ⊙ serra galaxies is between 10 M 7 ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 9 M ; these galaxies have ⊙ 10 -1 5 . Z ⊙ ∼ &lt; Zg &lt; ∼ 10 -0 5 . Z , which is broadly consistent with the ⊙ observed JWST galaxies. However, at 10 7 M ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 8 M ⊙ some serra galaxies show Zg ≃ 0 25 dex higher than the ob-. served vales; moreover, for M ⋆ ≃ 10 10 M , ⊙ Zg is on average lower than observed, up to extreme cases of 0 5 dex. In general, . serra galaxies show no clear MZR trend and their dispersion appears to be somewhat larger than the observed one. Interestingly, with respect to Nakajima et al. (2023) and Curti et al. (2024a), the lensed galaxies from Chemerynska et al. (2024a) seem to indicate lower values of Zg at M ⋆ &lt; ∼ 10 7 M ; the lensed data seems ⊙ consistent with serra , but only a handful of simulated galaxies and observations are in that mass range.

As noted by Curti et al. (2024a), the lack of a clear trend and the large scatter in the serra relation resemble the behavior of their 6 &lt; z &lt; 10 sample. However clear di ff erences are present, such as the relatively metal rich galaxies at M ⋆ ≃ 10 8 M and ⊙ low metallicity galaxies at M ⋆ ≃ 10 10 M . Further, the lack of ⊙ a clear trend might be partially due galaxy sample selection, as most of the serra data has 10 8 M ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 9 M and only a ⊙ small fraction of galaxies has M ⋆ ≃ 10 7 M . Thus, the observed ⊙ galaxy mass range is not uniformly covered by serra galaxies, whose sample is drawn from a collection of zoom-in simulations focusing on M ⋆ ≃ 10 10 M galaxies and their close environment. ⊙ Indeed, some of the M ⋆ &lt; ∼ 10 9 M mass galaxies form in a highly ⊙ pre-enriched environment, which explains their relatively high metallicities ( Zg ≃ 10 -0 5 . , see Gelli et al. 2020). Note that considering di ff erent redshift intervals 3 yields qualitatively similar results.

Thus, in summary, the normalization of the MZR of serra is roughly similar to the observed galaxies, but no trend is clearly visible and some discrepancies are present, especially in the M ⋆ ≃ 10 10 M range. Similar (or larger) tensions are present ⊙ also when a comparison with other models is performed. These can be noted, e.g. from the fit obtained by the fire -2 simulations (Marszewski et al. 2024). fire -2 data are close to the observed Zg at M ⋆ ≃ 10 9 5 . M , but the predicted slope of the MZR is di ⊙ ff erent, thus the tension with Nakajima et al. (2023); Curti et al. (2024a) increase with decreasing M ⋆ . Finally, fire -2 galaxies present a much higher scatter, that increases with decreasing stellar mass up to ≃ 1 dex at M ⋆ ≃ 5 × 10 7 M . ⊙

An in depth comparison between observations and simulations is presented in Curti et al. (2024a), which, in addition to serra , considers the zoom-in simulations from firstlight (Langan et al. 2020, see Ceverino et al. 2017 for the main paper) and fire (Ma et al. 2016, see Hopkins et al. 2014), the IllustrisTNG cosmological simulations Torrey et al. (2019, see Pillepich et al. 2018), and the astraeus semi-analytical models (Ucci et al. 2023, see Hutter et al. 2021); as detailed in Curti et al. (2024a), most simulations seem to reasonably match the Zg data at ≃ 10 9 M but have steeper slopes and underpredict the ⊙ metallicities at lower M ⋆ .

Interestingly, some of the simulations seem to be consistent with what observed by Chemerynska et al. (2024a), which however have only eight targets, thus it is unclear if there is a change of slope of the MZR for M ⋆ &lt; ∼ 10 7 M . ⊙

## 3. A minimal physical model to explain the MZR

Given the very contrived match between simulations and observations, it seems worthwhile to step back and try to understand the trend and dispersion of the MZR using basic physical models. The advantage of the approach consists in simplifying the complexity of numerical simulation while retaining the most important physical processes shaping the MZR.

To this aim we devise a minimal, physical model describing the evolution of DM ( Mdm ), gas ( Mg ), star ( M ⋆ ), and gas-phase metal ( M Z g ) mass of a galaxy, following Dayal et al. (2013). As in Furlanetto &amp; Mirocha (2022), we introduce an explicit feedback delay resulting in a modulation of the star formation history (SFH).

3 We have checked the serra data at z = 7 , 8 , and 9.

## 3.1. Model setup

We assume that the DM halo increases as a result of cosmological accretion, which on average can be written as in Correa et al. (2015, see eq. 23 therein)

$$\dot { M } _ { d m } = 7 1. 6 \left ( \frac { M _ { d m } } { 1 0 ^ { 1 2 } M _ { \odot } } \right ) \left ( \frac { 0. 7 } { h } \right ) \, E ( z ) \, \mathsf M _ { \odot } \, y r ^ { - 1 } \,, \quad \quad \ \ ( 1 a ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where E z ( ) = [ -0 24 . + 0 75(1 . + z )] p Ω m (1 + z ) 3 + ΩΛ . We allow stars to form (eq. 2b) and include stellar feedback via SNe, exploding after a delay time t d . SN cause mass outflows with a rate ˙ Mout (eq. 2a), enrich the gas of metals with a yield y , and have a return fraction R . Thus, the stellar mass increases because of SFR and decreases after t d because of SN explosions:

$$\dot { M } _ { * } & = S F R ( t ) - R \, S F R ( t - t _ { d } ) \,, & ( 1 b ) & \quad \varepsilon$$

which for t d = 0 (no-SN-delay) would give an instantaneous recycling approximation as in Dayal et al. (2013).

The gas mass increases because of cosmic accretion and return from stars, and it decreases because of the SFR process and SN driven outflows:

$$\dot { M } _ { g } = f _ { b } \dot { M } _ { d m } ( t ) + R \, \text{SFR} ( t - t _ { d } ) - \text{SFR} ( t ) - \dot { M } _ { o u t } ( t - t _ { d } ) \,, \quad ( 1 c ) \quad \text{black} ), \, s$$

where f b = Ω b / Ω m is the cosmological baryion fraction. Finally, the gas metal content increases because of the SN yields while it decreases because of astration and outflows:

$$\dot { M } _ { g } ^ { Z } = f _ { b } Z _ { i n } \dot { M } _ { d m } ( t ) + y \, \text{SFR} ( t - t _ { d } ) - Z _ { g } [ \text{SFR} ( t ) - \dot { M } _ { o u t } ( t - t _ { d } ) ] \,, \, ( 1 d ) \quad \text{limits of} \quad.$$

where Zg = M Z g / Mg is the gas metallicity, and we include the possibility for the cosmic gas to be already enriched at Zin , which we set to Zin = 10 -7 Z , i.e. virtually metal free. ⊙

To solve the system in eq.s 1, we need to specify a functional form for ˙ Mout and SFR. We take the outflow rate as in Muratov et al. (2015):

$$\dot { M } _ { o u t } = 3. 6 \, \epsilon _ { \text{SN} } \, S F R \left ( \frac { M _ { * } } { 1 0 ^ { 1 0 } M _ { \odot } } \right ) ^ { - 0. 3 5 } M _ { \odot } \, y r ^ { - 1 } \,, \quad \quad ( 2 a ) \quad \text{$which gi$} \\ \quad \text{ing an e}.$$

with ϵ SN an e ffi ciency parameter that we use to calibrate the mass loading factor; ϵ SN = 1 correspond to the results from Muratov et al. ( fire simulations, Hopkins et al. 2014 2014). Note that ϵ SN = 1 yields a loading factor consistent with what HerreraCamus et al. (2021) observes for HZ4 (but see Parlanti et al. 2024). Further, using ϵ SN = 1 in eq. 2 would match the results from Pizzati et al. (2023) for z ≃ 5 -6 ALPINE galaxies (Le Fèvre et al. 2020) which feature a spatially extended gas component (Fujimoto et al. 2020). However, ϵ SN ≃ 1 3 is needed to / match the loading factor inferred from z ≃ 3 -9 JADES galaxies (Carniani et al. 2024).

We assume that the SFR is proportional to the gas mass:

$$S F R & = \frac { M _ { g } } { t _ { * } } \,, & ( 2 b ) & \text{outflow} \, \in \, \\ & \text{in Fi} \\ & \text{out to $...$}.$$

and we select a constant t ⋆ = 1 Gyr as our fiducial star formation time scale 4 . The choice of t ⋆ is roughly consistent with the upper

4 In the model, galaxy evolution is self-regulated by feedback, thus the sensitivity of the results from t ⋆ is limited in our case, and null in the t d = 0 case with constant loading factors for infall and outflow (Lilly et al. 2013). A posteriori , we note that the fiducial t ⋆ = 1 Gyr gives specific star formation rates sSFR = SFR / M ⋆ which are consistent with serra galaxies, i.e. ≃ sSFR ≃ 100 -1 Gyr -1 for galaxies with masses M = 10 6 -10 10 M . ⊙

⋆

Fig. 2: Example of the impact of SFR burstiness on the evolution of a galaxies hosted in a Mdm = 3 × 10 10 M DM halo at ⊙ z = 6 using the minimal physical model described in Sec. 3.1 (see in particular eq.s 1). The time evolution of the dark matter (DM, black), star (green), gas (brown), and metal mass ( M ) content is shown with solid (dashed) lines for a SN delay time t d = 10 Myr ( t d = 0, i.e. no delay). The upper axis shows the redshift z corresponding to cosmic time t .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-a1628c1b-a3c8-484d-977d-c6f21f66252a.jpg)

limits of the depletion times reported by Dessauges-Zavadsky et al. (2020) for ALPINE galaxies, and it is discussed in more detail in Sec. 4.2 (see also App. A).

To complete the model, we assume that halos can form stars after they reach a virial temperature for atomic cooling to be e ff ective (10 4 K); similarly to serra , we assume return fraction and yields appropriate for a Kroupa (2001) IMF and adopting Bertelli et al. (1994) tracks with a Z ⋆ = Z ⊙ stellar population 5 , which gives y ≃ 0 0228 and . R ≃ 0 3242. .

The system (eq.s 1 and 2) is solved along with z = z t ( ) using an embedded 5th 6th order Runge-Kutta method (Fehlberg / 1968) with adaptive timestep, which is selected in order to have a fractional (absolute) precision of at least ≤ 10 -5 ( ≤ 10 2 M ) for ⊙ Mdm , M ⋆ , Mg , and M Z g . The delayed SFR and stellar mass used to evaluate the stellar feedback are computed using a 7 points time stencil.

## 3.2. Overview of the minimal physical model

The delay between SFR and feedback t d and the amplitude of the outflow ( ϵ SN) regulate the burstiness of the galaxy in our single zone model (cfr. Furlanetto &amp; Mirocha 2022). To understand the impact of burstiness on galaxy evolution, it is instructive to consider a halo with mass Mdm = 3 × 10 10 M at ⊙ z = 6 and set the outflow e ffi ciency to the standard value of ϵ SN = 1.

In Fig. 2 we compare the evolution of the stellar, gas, and metal mass content for t d = 0 (no-SN-delay, dashed lines) and t d = 10 Myr (solid lines). By construction, the halo mass growth (eq. 1a) is una ff ected by t d , thus the two models follow the same DMevolution, which smoothly grows from Mdm ≃ 10 7 M at ⊙ z ≃ 18, when the universe was t ≃ 200 Myr old. At that time, for t d = 0 the gas and stars has ≃ 10 5 M and ⊙ ≃ 10 2 M , respectively, and ⊙

5 Changing the metallicity of the stellar population does modifies the yield and return fraction, but gives only minor modifications to the upcoming results.

Fig. 3: Example of the impact of SFR burstiness on the MZR for a system with Mdm = 3 × 10 10 M at ⊙ z = 6 evolved with our minimal physical model. Di ff erent lines indicates di ff erent delay times for the SN feedback, as indicated in the legend. Note that the model results have been resampled in 2 Myr steps and we do not consider the evolutionary phases when the galaxies have a negligible amount of gas ( Mg ≤ 10 -3 Mdm ).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-2f347be0-adb0-49d9-830d-5d741405e5d4.jpg)

both grow reaching approximately ≃ 5 × 10 7 M at ⊙ t ≃ 800 Myr, with M ⋆/ Mdm ≃ 5 × 10 -3 broadly consistent with what expected from abundance matching models (e.g. Behroozi et al. 2013).

In the case with t d = 10 Myr, the secular trends are qualitatively similar, but the system experiences a series of bursts on time scales of ≃ 20 Myr = 2 t d . Let us focus on a single burst cycle around z ≃ 10, when the halo has Mdm ≃ 10 9 M ⊙ and M ⋆ ≃ 10 6 M . In the first half period of the cycle (0 ⊙ -t d ), cosmic accretion quickly replenish the gas (gas infall rate is ≃ 0 5 M . ⊙ yr -1 from eq. 1a), reaching Mg ≃ 10 7 M ⊙ in a few Myr. Stars form at a modest rate SFR ≃ 10 -2 M yr ⊙ -1 (eq. 2b) under feedback-free conditions (cfr. Dekel et al. 2023). After t d , the feedback kicks in, first enriching the galaxy and eventually ( ≃ 2 t d ) becoming powerful enough to expel both gas (gas outflow is ≃ 1M yr ⊙ -1 using eq. 2a) and metals (cfr. Ferrara et al. 2023). However, in the example the impact of SFR flickering is limited already after z ≃ 15, when the galaxy has M ⋆ ≃ 10 5 M ⊙ or 10 × lower if no delay is considered. At z ≃ 6, the final stellar mass is M ⋆ ≃ 10 7 M within a factor of 2 for the two models. For ⊙ the gas and metals, the model with SN delay starts to converge to the no-SN-delay case at z ≃ 8, when M ⋆ ≃ 10 7 M ; oscillations ⊙ for the gas and metals are still present up to z ≃ 7, but they are milder.

As in Pallottini &amp; Ferrara (2023), it is convenient to quantify the burstiness of the SFR by defining the variable δ expressing the log of the stochastic SFR variation with respect to its mean (in short: flickering)

$$\delta \equiv \log \frac { S F R } { \langle S F R \rangle } \,, & & ( 3 a ) & & \text{Cherney}$$

where log SFR M ⟨ / ⊙ yr -1 ⟩ is a second order polynomial fit in the time variable t . The low order fit is needed to factor out the smooth behavior of the functions without removing time modulations (see Leja et al. 2019; Chaves-Montero &amp; Hearin 2021, for alternative methods).

Similarly to the flickering, we define the analog variable for the gas metallicity

$$\left \| _ { 7 } \quad \delta _ { Z } \equiv \log \frac { Z _ { g } } { \langle Z _ { g } \rangle } \,,$$

where log ⟨ Zg / Z ⊙⟩ is a second order polynomial fit in the variable log M ⋆/ M . For both ⊙ δ and δ Z we can define the typical variation as the root mean square (r.m.s.) deviation, i.e. σ SFR and σ Z , respectively.

Adopting eq.s 3, we find that the t d = 10 Myr case has a flickering σ SFR ≃ 0 3 before the two models start to converge, . which is slightly higher than the σ SFR ≃ 0 24 obtained for . serra galaxies (Pallottini &amp; Ferrara 2023). For di ff erent t d the qualitative behavior is very similar, as the model is o ff -balance and oscillates around the no-SN-delay case, that is similar to a bathtub solution (Bouché et al. 2010; Dekel &amp; Mandelker 2014). For shorter (longer) t d , σ SFR can be smaller (higher), as galaxy with a smaller (higher) M ⋆ can be o ff -balanced.

The balance between infall and outflow rates is critical, since the dynamical system (eq.s 1) modeling the galaxy evolution give rise to exponentially increasing decreasing trends; introduc-/ ing a t d , 0 can e ffi ciently break the balance, causing cycles of complete gas and metal depletion.

To explore the MZR variations induced by the flickering, in Fig. 3 we plot the tracks of the selected Mdm = 3 × 10 10 M at ⊙ z = 6 with t d from 5 to 25 Myr. Because of the longer t d , the SFR becomes increasingly bursty. Specifically, for the selected model the flickering ranges from σ SFR ≃ 0 1 to . σ SFR ≃ 0 6 for . t d = 5 Myr to t d = 25 Myr, respectively.

The the MZR modulations induced by the flickering are highly non linear. With t d = 5 -10 Myr, the MZR dispersion is negligible-small ( σ Z ≃ 0 05 . -0 2), as the flickering can a . ff ect galaxies only up to M ⋆ ≃ 10 5 -7 M . For ⊙ t d = 15 Myr, the full M ⋆ range experience flickering with a σ SFR ≃ 0 5, yielding an high . σ Z ≃ 0 8. Cases with . t d &gt; 15 Myr do not show a higher stochasticity, as an increasing delay only regulates which M ⋆ can be a ff ected by the feedback unbalance . Qualitatively, the behavior 6 is similar with higher Mdm , but the convergence to the bathtub solution is faster.

In summary, a feedback delay induces a stochastic, bursty star formation behavior that becomes progressively more enhanced as t d is increased. Such SFR flickering corresponds to analog Zg fluctuations which act to decrease the mean metallicity as gas and metals are e ff ectively ejected from the galaxy by outflows.

## 3.3. The impact of a stochastic SFR on the MZR

We run sets of models that at z = 6 have a Mdm from 10 9 M to ⊙ 10 13 M with a logarithmic binning of 0 1 dex, in order to cover ⊙ . the mass range 10 M 5 ⊙ ∼ &lt; M ⋆ &lt; ∼ 5 × 10 10 M . We analyze di ⊙ ff erent variations of the minimal models, reporting the results in di ff erent panels of Fig. 4, where they are compared with data from JWST observations (Nakajima et al. 2023; Curti et al. 2024a; Chemerynska et al. 2024a).

## Base models

Let us start by discussing the base models analyzed in Sec. 3, i.e. with ϵ SN = 1, consistent with fire simulations (Muratov et al.

6 Note that an higher lower level of flickering can be induced by chang-/ ing the infall outflow rates, as shown later in Sec. 3.3. /

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-9f32c762-83d8-4dee-8fd9-8f6529669d18.jpg)

Fig. 4: Predicted MZR relation for di ff erent sets of minimal models with 10 M 9 ⊙ ≤ Mdm ≤ 10 13 M halos at ⊙ z = 6. In the upper left panel we report the base models with ϵ S N and t d = 20 Myr (see Sec. 3.2). In the upper right panel we modify the base model by reducing the e ffi ciency of the feedback ( ϵ SN = 1 4 in eq. 2a). In the / lower left panel we consider both weak feedback and a cosmic accretion that is modulated on a time scales of 100 Myr ( A 100 = 1 3 in eq. 4). In the / lower left panel we show the e ff ect of turning o ff the SN-delay for the weak feedback models with modulated cosmic accretion. In each panel, pastel lines are the individual tracks sampled every 2 Myr; their median and 16% -84% dispersion is plotted with violet dashed and solid lines, respectively; the dashed black line is the control case with no-SN-delay ( t d = 0) and no modulated accretion ( A 100 = 0). For each model, we report the value values of the SFR flickering ( σ SFR) and the MZR dispersion ( σ Z , see Fig. 5). As in Fig. 1, we show JWST data and fits (Nakajima et al. 2023; Curti et al. 2024a; Chemerynska et al. 2024a).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-335de01d-9a1d-4826-96b3-28e3b8d5df59.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-e12d688a-36fd-4739-adf2-217212bad38a.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-9be2a052-def1-4d34-9be0-298f386b1867.jpg)

2015), and selecting t d = 20 Myr, so that most of the explored stellar mass range is expected to have a stochastic SFR.

extreme, as only part of the lensed galaxy data can be matched (Chemerynska et al. 2024a).

In the upper left panel of Fig. 4 we report both the individual tracks of the base models and their M ⋆ average trends. For the base models, the median case falls significantly below the observed average MZR. Also, the dispersion is much larger than seen in JWST data, that can be barely matched by the upper envelope of the predicted curves. The base models have a σ SFR ≃ 0 5 flickering, that induces a . σ Z &gt; ∼ 1 almost in the entire stellar mass range ( M ⋆ &lt; ∼ 10 9 M , Fig. 5). The flickering is due ⊙ to cycles of infall outflow which determine an oscillatory behav-/ ior of the evolution around the t d = 0 control case, that is largely independent of the hosting Mdm . The Zg -M ⋆ normalization is determined by the strength of the outflows, that completely devoids the galaxy of gas and metals, not allowing the metallicity to rise to the observed level. Even for t d = 0 the outflow seems too

In summary, the e ffi cient feedback ( ϵ SN = 1) combined with a t d = 20 Myr generates a SFR flickering of σ SFR ≃ 0 5, . inducing MZR variations of σ Z ≃ 1 4 resulting in a 'chemical . chaos' that is not observed. Barring the option that observations are tracing only the tip of the Zg -M ⋆ relation, in order to explain the observed MZR we modify the base models as follows.

## Weak feedback

First, we reduce the impact of stellar feedback by decreasing the outflow rate from Muratov et al. (2015) (eq. 2a), i.e. setting ϵ SN = 1 4; such a loading factor is closer to what Carniani et al. (2024) / infer for 10 7 5 . &lt; ∼ M ⋆ &lt; ∼ 10 9 M JADES galaxies. As shown in the ⊙ upper right panel of Fig. 4, this modification increases the average metallicity of the modeled galaxies, matching the observed

Fig. 5: SFR ( σ SFR, upper panel) and gas metallicity ( σ Z , lower panel) r.m.s. variations as a function of stellar mass for the sets of models presented in Fig. 4 (see eq.s 3 for the definitions). As a reference, we report the serra average σ SFR (Pallottini &amp; Ferrara 2023) and σ Z (this work), the σ SFR from Ciesla et al. (2024); Mason et al. (2023); Muñoz et al. (2023); Shen et al. (2023), and the average σ Z for JWST data. Note that ⟨ JWST is the r.m.s. of the ⟩ observations after subtraction of the fit to the MZR in the di ff erent data sets in Nakajima et al. (2023) and Curti et al. (2024a); the thickness of the shaded area encloses the min max of the / r.m.s. of the di ff erent data sets.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-948bd523-9dbe-47e2-a983-6dff5c7379db.jpg)

mean trends in the whole 10 6 M ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 10 M mass range. ⊙ Comparing to the t d = 0 control case, it is evident that the outflow e ffi ciency is the main driver of the mean Zg -M ⋆ trend.

Note that, reducing the outflow e ffi ciency while maintaining a fixed t d = 20 Myr results in σ Z ≃ 0 1 modulation of the MZR, . i.e. too low with respect to the observed ( σ Z ≃ 0 25, see Fig. 5). . This e ff ect is induced by the suppression of the amplitude of the SFR flickering, especially at the high mass end.

Indeed the models show a small flickering, with σ SFR ≃ 0 15 . for M ⋆ &lt; ∼ 10 8 M , and even slightly lower values for higher stel-⊙ lar masses (Fig. 5). As discussed in Pallottini &amp; Ferrara (2023), M ⋆ &lt; ∼ 10 9 M galaxies show SFR fluctuations with characteristic ⊙ periods of ≃ 20 -40 Myr, which correspond to the typical evolutionary times of SNe. Such flickering timescale increases to ≃ 80 -150 Myr in more massive galaxies which are more sensitive to mergers and accretion rate modulation. The latter e ff ects are not accounted for in our minimal model (eq. 1a). The next step is then to heuristically incorporate them in our treatment.

## Weak feedback, modulated accretion

To incorporate the extra stochastic behavior expected from variation in the cosmic accretion in a simple manner, we impose a sinusoidal modulation on a 100 Myr timescale

$$\dot { M } _ { d m } ^ { m o d } & = \dot { M } _ { d m } \left [ 1 + A _ { 1 0 0 } \sin \left ( \frac { t } { 1 0 0 \, M y r } \right ) \right ] \,, & \quad ( 4 ) \quad \text{snorer $s$}$$

with ˙ Mdm being the average cosmic accretion from eq. 1a and with A 100 scaling the amplitude of the modulation . 7

As feedback determines the mean trend of the MZR, we keep the reduced outflow e ffi ciency ( ϵ SN = 1 4 in eq. 2a) and compute / the evolution of a minimal model by setting the amplitude of the oscillation to A 100 = 1 3, which should yield a scatter below the / 0.3 dex dispersion expected for the distribution of growth rates of DM halos (Rodríguez-Puebla et al. 2016; Ren et al. 2019; Mirocha et al. 2021). The results from these modified models are shown in the lower left panel of Fig. 4: both the mean trend and scatter of the observed MZR are recovered (Fig. 5). Such a matching set of modified models yields a roughly constant σ SFR ≃ 0 2 and . σ Z ≃ 0 3, . satisfactorily matching the JWST data.

## Weak feedback, modulated accretion, no-SN-delay

As a final check, in the lower right panel of Fig. 4 we show the effect of removing the SN-delay from the best matching model just discussed. While the mean MZR trend is recovered, the metallicity dispersion ( σ Z ≃ 0 05) is . much smaller than observed ( σ Z ≃ 0 25). Also, the median of the models fails to reproduce . the relation for most of the low mass lensed galaxies (Chemerynska et al. 2024a).

With no-SN-delay, the SFR flickering is induced only by modulation of the cosmic accretion, which fails to unbalance the system through the ejection of a substantial fraction of the contained gas. With the selected A 100 = 1 3 and / t d = 0, the model reach σ SFR ≃ 0 1 for . M ⋆ &lt; ∼ 10 8 M , while the flickering gets to ⊙ σ SFR ≃ 0 2 only for . M ⋆ ≃ 5 × 10 8 M galaxies. ⊙

Compared to the control case ( A 100 = 0), the average Zg is slightly higher; practically, SFR flickering induced by the accretion modulation can increase Zg since a higher gas mass can be converted into stars without e ffi ciently ejecting the gas via outflows, i.e. the opposite outcome with respect to SFR flickering induced via SN-delay. Interestingly, despite the relatively high SFR flickering (average of σ SFR ≃ 0 2), the metallicity variation . is the lowest of all considered cases ( σ Z ≃ 0 05): only a delayed . SN feedback can e ffi ciently o ff -balance the system and induce modulations on both SFR and Zg .

We note that the serra cosmological simulations show a relatively modest SFR flickering ( σ SFR ≃ 0 24, Pallottini &amp; Ferrara . 2023), and close to the expectations for z ≃ 6 galaxies (Muñoz et al. 2023; Ciesla et al. 2024). Yet, given the sharp sensitivity of σ Z to SFR variations, such time variability results in a σ Z larger than observed.

To summarize, our analysis shows that: (i) a weak feedback, and (ii) a relatively mild stochasticity given by (iii) SN-delay combined with a long term accretion rate modulation are critical to simultaneously match both the mean and dispersion of the observed MZR relation in the entire stellar mass range of galaxies sampled by JWST at z = 3 -10.

Finally, we note that invoking stochasticity to explain the overabundance of bright JWST galaxies requires σ SFR &gt; 0 5 . (more specifically, σ SFR ≃ 0 8, . σ SFR ≃ 0 7, and . σ SFR ≃ 0 6 ac-. cording to Muñoz et al. 2023, Shen et al. 2023, and Mason et al. 2023, respectively; see Fig. 5). However, such extreme flickering amplitudes are (i) higher than that recovered from SED fitting al-

7 Note that small halos should be relatively less a ff ected by modulation on long timescales ( ≃ 100 Myr), since the hosted galaxies have shorter SFH. However, cosmological accretion variations can be coherently combined with the SN delay to give a higher σ SFR also in the low mass range.

Fig. 6: Deviation from the fundamental mass-metallicity relation (FMR) for highz galaxies. We adopt the µ ≡ log( M ⋆/ M ) ⊙ -0 65 log(SFR M . / ⊙ yr -1 ) FMR obtained by Curti et al. (2020) for the high sSFR galaxies sub-sample from the SDSS ( z &lt; ∼ 2 York et al. 2000). We report the results of the weak feedback ( ϵ SN = 1 4), modulated accretion ( / A 100 = 1 3) models along / with the data from JWST observations (Nakajima et al. 2023; Curti et al. 2024a; Chemerynska et al. 2024a) and serra simulations (Pallottini et al. 2022).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-db65f197-96ef-485e-96c8-d77cb275cc30.jpg)

lowing for SFH variability ( σ SFR &lt; ∼ 0 5, Ciesla et al. 2024), and . (ii) would induce a chemical chaos ( σ Z &gt; ∼ 1 4) that is not present . in the observed MZR at z = 3 -10 ( σ Z ≃ 0 25 Nakajima et al. . 2023; Curti et al. 2024a; Chemerynska et al. 2024a).

## 4. Discussion

The minimal models presented here can also be useful as a tool to guide the development and compare complex cosmological simulations, particularly when focusing on relations that - like the MZR - are very sensitive to small variations of the underlying physics, that is included via sub-grid models in galaxy simulations (Agertz et al. 2020).

## 4.1. Comparing numerical simulations

In galaxy simulations, ejective feedback, such as outflows caused by mechanical kinetic prescriptions, can e / ff ectively suppress the SFR activity as much as preventative feedback, such as turbulence delayed cooling prescriptions (Rosdahl et al. 2017), / leading to similar SFR histories when the same initial conditions are considered (Lupi et al. 2020). Despite similar SFR histories, other properties can di ff er significantly (Rosdahl et al. 2017), such as the dynamical state of a galaxy. For instance, Fujimoto et al. (2024) show that fire -2 galaxies (Wetzel et al. 2023) have a much lower incidence of disks compared to serra galaxies (Kohandel et al. 2024). Additionally, SERRA galaxies tend to have higher Zg and lower MZR scatter compared to fire -2 (Marszewski et al. 2024), especially in the M ⋆ ≃ 10 9 M mass ⊙ range.

Both e ff ects are likely caused by the higher SFR flickering of fire -2 (Sun et al. 2023a) with respect to serra (Pallottini &amp; Ferrara 2023), which is ultimately regulated by the e ffi ciency and kind (e.g. ejective vs preventative) of included stellar feedback.

Violent outflows are present in fire -2 (Hopkins et al. 2018), while the suppression of SFR due to the photo-dissociation of molecular hydrogen in serra (Pallottini et al. 2019) can more gently can modulate the SFR activity without drastically modifying the galaxy dynamics.

Similar considerations hold when comparing successive improvements of simulations suites. For instance, the loading factor from fire -2 (Pandya et al. 2021) are lower by a factor of ≃ 2 with respect to fire (Muratov et al. 2015); this di ff erence dictates that the MZR from fire -2 (Marszewski et al. 2024) is about 0.3 dex higher than the MZR from fire (Ma et al. 2016); in both cases the slope is steeper (normalization is lower) than what observed by Nakajima et al. (2023); Curti et al. (2024a) and is instead more consistent with Chemerynska et al. (2024a). Despite the simplicity of the minimal model presented, hints of such di ff erences in the simulations can be highlighted and explained.

However, it is di ffi cult to directly quantify the di ff erences in the simulations with the presented models. For instance, the loading factor in our simplified models (but also in relatively more complex ones as Thompson et al. 2016 and Pizzati et al. 2020) is a scalar quantity used as an input value for each galaxy, while in simulations is usually defined as an output that has a spatial dependence (Gallerani et al. 2018).

Additionally, when dealing with simulation snapshots, quantifying both the SFR event and the consequent outflows is tricky because of the time-delay between SFR and feedback (Pandya et al. 2021), particularly for galaxies far from a quasiequilibrium situation.

Further, the simplifications done in our minimal model impose we can only capture part of the complexity present in simulations. For instance, while the SN-feedback delay is uniquely defined, simulations finds a distribution of t d that depend on the underlying physical models included (Pallottini &amp; Ferrara 2023); thus, trying to find the set of parameters of a minimal model that gives the best match to a simulation will only give e ff ective parameters, a sort of summary statistic that can be used to compare di ff erent simulations.

Apart from these di ffi culties, a systematic parameter fitting procedure should give a better quantification of the comparison between di ff erent simulations, and can also be used to enable an even better match between the data and the minimal models. However, such models are lacking a few ingredients that prevents reproducing the full connection between metal build up and star formation history of galaxies.

## 4.2. Caveats of the minimal model

The limitations of the presented minimal models can be appreciated by considering the fundamental mass-metallicity relation (FMR), i.e. the link between Zg , M ⋆ , and SFR that is observed to be redshift independent up to z ∼ 3 (Mannucci et al. 2010; Curti et al. 2020) and seems to break down at z &gt; ∼ 5 (Nakajima et al. 2023; Curti et al. 2024a). We adopt the parameterization from Curti et al. (2020)

$$\underset { \sim 1 \dots \dots } { \text{to SERRA} } \ \mu \equiv \log ( M _ { * } / M _ { \odot } ) - 0. 6 5 \log ( S F R / M _ { \odot } \, y r ^ { - 1 } ) \,,$$

which is obtained by minimizing the Zg dispersion for the high sSFR sub-sample from the Sloan Digital Sky Survey (SDSS, York et al. 2000, z &lt; ∼ 2).

In Fig. 6 we show the FMR for the minimal models matching the mean MZR and its dispersion, i.e. the weak feedback ( ϵ SN = 1 4) set with modulated accretion ( / A 100 = 1 3). These models / are below the SDSS relation by about 0.5 dex at µ &lt; ∼ 7, with the

separation that decreases with increasing µ , down to 0.1 dex at µ &gt; ∼ 9, almost connecting with the FMR trend observed at low z . As expected (Nakajima et al. 2023; Curti et al. 2024a), JWST galaxies are o ff set below the SDSS relation. The best matching models for the MZR only partially match the JWST data in the FMR plane; the minimal models are closer to the JWST data when considering the lensed galaxies from Chemerynska et al. (2024a).

Further the minimal models only produce downward deviation with respect to the FMR from the SDSS. Instead, serra galaxies manage to better recover the trend for JWST galaxies around µ ≃ 8, reproducing both downward and upward deviations with respect to the FMR from the SDSS; however, serra galaxies show an overabundance of low Zg galaxies at µ ≃ 8 5 . -9 0, which is not seen in the JWST data. .

In an equilibrium model (Lilly et al. 2013), the shape of the MZR is independent from the SFR, which uniquely determines the speed a galaxy can climb up the Zg -M ⋆ curve. Note that the simplified model seems to deviate more from the data at low µ , i.e. a combination of relatively high SFR and low M ⋆ ; the adopted feedback increases in e ffi ciency for lower M ⋆ (eq. 2a). Thus, the observed FMR deviation might be recovered by the minimal model if we consider not only ejecting but also a preventing feedback, e.g. photodissociation of molecular hydrogen, which can suppress the star formation without a ff ecting the gas mass (e.g. see Alyssum, a serra galaxy described in Pallottini et al. 2022) and can play a role in explaining some of the bright blue monsters seen by JWST (Ferrara 2024).

As an alternative to keep matching the MZR and recovering the o ff set from the FMR, we could try to change the simplified SFR prescription (eq. 2b); this can be done for instance by adopting as the time scale for star formation a redshift dependent depletion time (Tacconi et al. 2010, 2020; Sommovigo et al. 2022), which are typically shorter than the selected t ⋆ = 1 Gyr by about a factor ≃ 10, similar to what Vallini et al. (2024) report for z ≃ 7 galaxies by using glam (Vallini et al. 2020, 2021) determinations, thus promoting a faster galaxy evolution at high- . z

Indeed, by adopting the depletion time from Sommovigo et al. (2022) as the time scale for SFR (see App. A, in particular Fig. A.1), the stellar mass build up is faster at M ⋆ &lt; ∼ 10 6 M , ⊙ while later saturates to values similar to those resulting from the t ⋆ = 1 Gyr case as both SFHs become feedback regulated, approaching a bathtub equilibrium.

However, on the one hand, the modifications a ff ect the predictions of the minimal model mostly in the M ⋆ &lt; ∼ 10 6 M range; ⊙ on the other hand, a z average depletion time is sensitive to the population of the sample, e.g. for main sequence is ≃ × 5 lower than for starburst (Tacconi et al. 2020), and the results from Sommovigo et al. (2022) are mostly relying on the data of M ⋆ &gt; ∼ 10 8 -10 M galaxies from the ALPINE (Le Fèvre et al. ⊙ 2020) and REBELS (Bouwens et al. 2022) surveys. To summarize, while trying to reconcile the FMR behavior is an interesting perspective and there are a few options, some care should be taken in extending modifying the prescriptions for the minimal / models; we leave such a possibility for a future work.

## 5. Conclusions

The exquisite spectroscopic data collected by JWST allow for the first time to derive the mass-metallicity relation (MZR) in galaxies up to redshift z ≃ 10.

In this work, we considered the MZR data for a combined sample of about 180 galaxies with stellar mass 10 6 M ⊙ ∼ &lt; M ⋆ &lt; ∼ 10 10 M ⊙ at z = 3 -10 given Nakajima et al.

(2023), Curti et al. (2024a), and Chemerynska et al. (2024a), and obtained from the CEERS (Finkelstein et al. 2023), JADES (Bunker et al. 2023), and UNCOVER (Bezanson et al. 2022) surveys, respectively.

We compared the observations with the predictions from the serra simulations (Pallottini et al. 2022), finding a broad agreement with the observed MZR. However, simulations show a lack of a clear metallicity trend, and some tension, particularly at M ≃ 10 10 M . ⊙

⋆

To better understand the observed highz MZR and to clarify the behaviour of the simulations, we have devised a minimal physical model for galaxy evolution in which star formation is fueled by cosmic gas accretion (Correa et al. 2015) and is regulated by supernova driven outflows (eq. 2a), with a loading factor from Muratov et al. (2015) that is controlled by an e ffi ciency ϵ SN; additionally, we incorporated an explicit delay ( t d ) between star formation and supernova feedback, inducing a stochastic SFR behavior. Some models also explore the possible modulation of cosmic accretion (eq. 4). s are:

- · To recover the average trend of the MZR observed at high- , z outflows must be less e ffi cient ( ϵ SN = 1 4) than predicted in / Muratov et al. (2015). However, this evidence is consistent with the mass loading factors inferred for JADES galaxies (Carniani et al. 2024).
- · To match the MZR dispersion in the full stellar mass range, both a delayed ( t d = 20 Myr) SN feedback, and a modulation of the cosmic accretion ( A 100 = 1 3) are required. /
- · The r.m.s. variation of the MZR ( σ Z ) is very sensitive to the SFR flickering ( σ SFR). The weak feedback necessary to reproduce the average MZR trend ( ϵ SN = 1 4) also re-/ sults in a low-amplitude SFR flickering ( σ SFR ≃ 0 2). This . prescription correctly reproduces the moderate MZR scatter ( σ Z ≃ 0 25) observed by JWST. .
- · serra galaxies have slightly higher σ SFR ≃ 0 24 (Pallottini &amp; . Ferrara 2023), resulting in a metallicity scatter of σ Z ≃ 0 3, . higher than observed. Di ff erently from the analytical model, serra galaxies also overpredict the number of low metallicity ( Zg ≃ 0 2Z . ⊙ ) galaxies at M ⋆ ≃ 10 10 M . ⊙
- · In general, any model predicting σ SFR &gt; ∼ 0 5 is likely to over-. shoot the observed metallicity r.m.s. scatter, possibly leading to a 'chemical chaos' ( σ Z &gt; ∼ 1 4) that is it not present in . JWST data, for which σ Z ≃ 0 25). .

The last point also entails that simultaneously explaining the observed MZR (Nakajima et al. 2023; Curti et al. 2024a) and the overabundance of bright galaxies observed by JWST (Finkelstein et al. 2022; Naidu et al. 2022) through high SFR flickering (Mason et al. 2023; Mirocha &amp; Furlanetto 2023; Shen et al. 2023; Sun et al. 2023b; Kravtsov &amp; Belokurov 2024) is very challenging, considering that a σ SFR ≃ 0 8 is required to explain . z &gt; 10 observations (Muñoz et al. 2023). We note that while both the SN-feedback delay and cosmic accretion modulation induce SFR flickering, the latter tends to give a gentler modulation of Zg , since high feedback e ffi ciency and a delay are needed in order to o ff -balance the MZR generated from a bathtub-like equilibrium (Lilly et al. 2013). However, variation of cosmic accretion can only go up to 0.3 dex (Rodríguez-Puebla et al. 2016; Ren et al. 2019; Mirocha et al. 2021). Thus, alternative explanations (Dekel et al. 2023; Ferrara et al. 2023) seem more favored to solve the overabundance problem.

While we showed that the presented minimal model can also be used to guide our intuition in comparing complex cosmological simulations, we recall that such a simplified model does not fully capture the physical complexity of the connection between

galaxy formation and metal enrichment, as highlighted by the poor match with the JWST data when analyzing the o ff -set from the fundamental mass-metallicity relation. Combining insights from such models with further analysis of ongoing observations, and more complex numerical simulations is crucial for understanding the galaxy formation and evolution process at high- . z

Acknowledgements. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support from the Class B project SERRA HP10BPUZ8F (PI: Pallottini). AF acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020 740120 / (PI: Ferrara). SC and GV acknowledge support support from the ERC Starting Grant WINGS H2020 101040227 (PI: Carniani). Any dissemination of results / must indicate that it reflects only the author's view and that the Commission is not responsible for any use that may be made of the information it contains. Partial support (AF) from the Carl Friedrich von Siemens-Forschungspreis der Alexander von Humboldt-Stiftung Research Award is kindly acknowledged. We gratefully acknowledge computational resources of the Center for High Performance Computing (CHPC) at SNS. We acknowledge usage of the Python programming language (Van Rossum &amp; de Boer 1991; Van Rossum &amp; Drake 2009), Astropy (Astropy Collaboration et al. 2013), Cython (Behnel et al. 2011), Matplotlib (Hunter 2007), Numba (Lam et al. 2015), NumPy (van der Walt et al. 2011), pynbody (Pontzen et al. 2013), and SciPy (Virtanen et al. 2020).

## References

Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, 518, 4755

Agertz, O., Pontzen, A., Read, J. I., et al. 2020, MNRAS, 491, 1656

- Andrews, B. H. &amp; Martini, P. 2013, ApJ, 765, 140
- Asplund, M., Grevesse, N., Sauval, A. J., &amp; Scott, P. 2009, ARA&amp;A, 47, 481
- Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&amp;A, 558, A33
- Atek, H., Shuntov, M., Furtak, L. J., et al. 2023, MNRAS, 519, 1201
- Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Computing in Science Engineering, 13, 31
- Behroozi, P. S., Wechsler, R. H., &amp; Wu, H.-Y. 2013, ApJ, 762, 109
- Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., &amp; Nasi, E. 1994, A&amp;A Supp., 106, 275
- Bezanson, R., Labbe, I., Whitaker, K. E., et al. 2022, arXiv e-prints, arXiv:2212.04026
- Bothwell, M. S., Maiolino, R., Kennicutt, R., et al. 2013, MNRAS, 433, 1425
- Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001
- Bouwens, R. J., Smit, R., Schouws, S., et al. 2022, ApJ, 931, 160
- Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2306.02467
- Cameron, A. J., Katz, H., &amp; Rey, M. P. 2023, MNRAS, 522, L89
- Carniani, S., Venturi, G., Parlanti, E., et al. 2024, A&amp;A, 685, A99
- Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJL, 938, L15
- Ceverino, D., Glover, S. C. O., &amp; Klessen, R. S. 2017, MNRAS, 470, 2791
- Chaves-Montero, J. &amp; Hearin, A. 2021, MNRAS, 506, 2373
- Chemerynska, I., Atek, H., Dayal, P., et al. 2024a, arXiv e-prints, arXiv:2407.17110
- Chemerynska, I., Atek, H., Furtak, L. J., et al. 2024b, MNRAS, 531, 2615
- Ciesla, L., Elbaz, D., Ilbert, O., et al. 2024, A&amp;A, 686, A128
- Correa, C. A., Wyithe, J. S. B., Schaye, J., &amp; Du ff y, A. R. 2015, MNRAS, 450, 1521
- Curti, M., Cresci, G., Mannucci, F., et al. 2017, MNRAS, 465, 1384
- Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2024a, A&amp;A, 684, A75
- Curti, M., Mannucci, F., Cresci, G., &amp; Maiolino, R. 2020, MNRAS, 491, 944
- Curti, M., Witstok, J., Jakobsen, P., et al. 2024b, arXiv e-prints, arXiv:2407.02575
- Davies, F. B. &amp; Furlanetto, S. R. 2016, MNRAS, 460, 1328
- Dayal, P. &amp; Ferrara, A. 2018, Phys. Rep., 780, 1
- Dayal, P., Ferrara, A., &amp; Dunlop, J. S. 2013, MNRAS, 430, 2891
- Decataldo, D., Pallottini, A., Ferrara, A., Vallini, L., &amp; Gallerani, S. 2019, MNRAS, 487, 3377
- Dekel, A. &amp; Mandelker, N. 2014, MNRAS, 444, 2071
- Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., &amp; Li, Z. 2023, MNRAS, 523, 3201
- Dessauges-Zavadsky, M., Ginolfi, M., Pozzi, F., et al. 2020, A&amp;A, 643, A5 Dome, T., Tacchella, S., Fialkov, A., et al. 2024, MNRAS, 527, 2139
- Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011
- Ellison, S. L., Patton, D. R., Simard, L., &amp; McConnachie, A. W. 2008, ApJL, 672, L107
- Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813
- Fehlberg, E. 1968, Classical Fifth-, Sixth-, Seventh-, and Eighth-Order RungeKutta Formulas with Stepsize Control
- Ferrara, A. 2008, in Low-Metallicity Star Formation: From the First Stars to Dwarf Galaxies, ed. L. K. Hunt, S. C. Madden, &amp; R. Schneider, Vol. 255, 86-99
- Ferrara, A. 2024, arXiv e-prints, arXiv:2405.20370
- Ferrara, A., Pallottini, A., &amp; Dayal, P. 2023, MNRAS, 522, 3986
- Finkelstein, S. L. &amp; Bagley, M. B. 2022, ApJ, 938, 25
- Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJL, 946, L13

Finkelstein, S. L., Bagley, M. B., Haro, P. A., et al. 2022, ApJL, 940, L55

- Förster Schreiber, N. M. &amp; Wuyts, S. 2020, ARA&amp;A, 58, 661
- Fujimoto, S., Ouchi, M., Kohno, K., et al. 2024, arXiv e-prints, arXiv:2402.18543
- Fujimoto, S., Silverman, J. D., Bethermin, M., et al. 2020, ApJ, 900, 1

Furlanetto, S. R. &amp; Mirocha, J. 2022, MNRAS, 511, 3895

- Gallerani, S., Pallottini, A., Feruglio, C., et al. 2018, MNRAS, 473, 1909
- Gelli, V., Mason, C., &amp; Hayward, C. C. 2024a, arXiv e-prints, arXiv:2405.13108
- Gelli, V., Salvadori, S., Ferrara, A., &amp; Pallottini, A. 2024b, ApJ, 964, 76
- Gelli, V., Salvadori, S., Ferrara, A., Pallottini, A., &amp; Carniani, S. 2023, ApJL, 954, L11
- Gelli, V., Salvadori, S., Pallottini, A., &amp; Ferrara, A. 2020, MNRAS, 498, 4134 Genzel, R., Förster Schreiber, N. M., Übler, H., et al. 2017, Nature, 543, 397 González, V., Labbé, I., Bouwens, R. J., et al. 2010, ApJ, 713, 115
- Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386

Hahn, O. &amp; Abel, T. 2011, MNRAS, 415, 2101

- Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5
- Herrera-Camus, R., Förster Schreiber, N., Genzel, R., et al. 2021, A&amp;A, 649, A31
- Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581
- Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800
- Hunter, J. D. 2007, Computing in Science Engineering, 9, 90
- Hutter, A., Dayal, P., Yepes, G., et al. 2021, MNRAS, 503, 3698
- Kannan, R., Garaldi, E., Smith, A., et al. 2022, MNRAS, 511, 4005
- Katz, H., Kimm, T., Ellis, R. S., Devriendt, J., &amp; Slyz, A. 2023, MNRAS, 524, 351
- Kennicutt, Jr., R. C. 1998, ApJ, 498, 541
- Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102
- Kobayashi, C. &amp; Ferrara, A. 2024, ApJL, 962, L6
- Kohandel, M., Pallottini, A., Ferrara, A., et al. 2024, A&amp;A, 685, A72

Kravtsov, A. &amp; Belokurov, V. 2024, arXiv e-prints, arXiv:2405.04578

- Kroupa, P. 2001, MNRAS, 322, 231
- Krumholz, M. R. &amp; Burkhart, B. 2016, MNRAS, 458, 1671
- Lam, S. K., Pitrou, A., &amp; Seibert, S. 2015, in Proc. Second Workshop on the LLVM Compiler Infrastructure in HPC, 1-6
- Langan, I., Ceverino, D., &amp; Finlator, K. 2020, MNRAS, 494, 1988
- Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&amp;A, 643, A1
- Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., &amp; Speagle, J. S. 2019, ApJ, 876, 3
- Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., &amp; Peng, Y. 2013, ApJ, 772, 119
- Liu, B. &amp; Bromm, V. 2020, MNRAS, 497, 2839
- Looser, T. J., D'Eugenio, F., Maiolino, R., et al. 2024, Nature, 629, 53
- Lupi, A., Pallottini, A., Ferrara, A., et al. 2020, MNRAS, 496, 5160
- Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140 Madau, P. &amp; Dickinson, M. 2014, ARA&amp;A, 52, 415
- Maiolino, R. &amp; Mannucci, F. 2019, A&amp;A Rev., 27, 3
- Maiolino, R., Nagao, T., Grazian, A., et al. 2008, A&amp;A, 488, 463
- Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., &amp; Gnerucci, A. 2010, MNRAS, 408, 2115
- Marconi, A., Amiri, A., Feltre, A., et al. 2024, arXiv e-prints, arXiv:2401.13028 Marszewski, A., Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., &amp; Feldmann, R. 2024, ApJL, 967, L41
- Mason, C. A., Trenti, M., &amp; Treu, T. 2023, MNRAS, 521, 497
- Mirocha, J. &amp; Furlanetto, S. R. 2023, MNRAS, 519, 843
- Mirocha, J., La Plante, P., &amp; Liu, A. 2021, MNRAS, 507, 3872
- Morishita, T., Stiavelli, M., Grillo, C., et al. 2024, arXiv e-prints, arXiv:2402.14084
- Muñoz, J. B., Mirocha, J., Furlanetto, S., &amp; Sabti, N. 2023, MNRAS, 526, L47 Muratov, A. L., Kereš, D., Faucher-Giguère, C.-A., et al. 2015, MNRAS, 454, 2691
- Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJL, 940, L14
- Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, ApJS, 269, 33
- Narayanan, D., Lower, S., Torrey, P., et al. 2024, ApJ, 961, 73
- Nikoli´ c, I., Mesinger, A., Davies, J. E., &amp; Prelogovi´ c, D. 2024, arXiv e-prints, arXiv:2406.15237
- Pallottini, A. &amp; Ferrara, A. 2023, A&amp;A, 677, L4
- Pallottini, A., Ferrara, A., Bovino, S., et al. 2017a, MNRAS, 471, 4128
- Pallottini, A., Ferrara, A., Decataldo, D., et al. 2019, MNRAS, 487, 1689
- Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621
- Pallottini, A., Ferrara, A., Gallerani, S., Salvadori, S., &amp; D'Odorico, V. 2014, MNRAS, 440, 2498
- Pallottini, A., Ferrara, A., Gallerani, S., et al. 2017b, MNRAS, 465, 2540

| Pandya, V., Fielding, D. B., Anglés-Alcázar, D., et al. 2021, MNRAS, 508, 2979 Parlanti, E., Carniani, S., Venturi, G., et al. 2024, arXiv e-prints, arXiv:2407.19008 A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077 E., Ferrara, A., Pallottini, A., et al. 2020, MNRAS, 495, 160 E., Ferrara, A., Pallottini, A., et al. 2023, MNRAS, 519, 4608 Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A16 A., Rovskar, R., Stinson, G. S., et al. 2013, pynbody: Astrophysics Sim- ulation Analysis for Python, astrophysics Source Code Library, ascl:1305.002   |
|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Pillepich, Pizzati,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
| Pizzati,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     |
| Planck                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
| Pontzen,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     |
| Ren, K., Trenti, M., &Mason, C. A. 2019, ApJ, 878, 114                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
| Rizzo, F., Kohandel, M., Pallottini, A., et al. 2022, A&A, 667, A5 Roberts-Borsani, G., Morishita, T., Treu, T., et al. 2022, ApJL, 938, L13                                                                                                                                                                                                                                                                                                                                                                                                                                                 |
| Rodríguez-Puebla, A., Primack, J. R., Behroozi, P., & Faber, S. M. 2016, MN- RAS, 455, 2592                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
| 436, 2188                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    |
| Rosdahl, J., Schaye, J., Dubois, Y., Kimm, T., & Teyssier, R. 2017, MNRAS, 466, 11                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           |
| Rowland, L. E., Hodge, J., Bouwens, R., et al. 2024, arXiv e-prints, arXiv:2405.06025 Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2018, ApJ, 858, 99                                                                                                                                                                                                                                                                                                                                                                                                                                   |
| Santini, P., Fontana, A., Castellano, M., et al. 2023, ApJL, 942, L27                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        |
| Sarmento, R., Scannapieco, E., &Cohen, S. 2018, ApJ, 854, 75                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 |
| Schmidt, M. 1959, ApJ, 129, 243                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              |
| Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   |
| R. 2023, MNRAS, 525, 3254                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    |
| S., & Kannan, Shibuya, T., Ouchi, M., &Harikane, Y. 2015, ApJS, 219, 15                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      |
| Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2017, ApJ, 843, 46                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       |
| Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
| Sommovigo, L., Ferrara, A., Carniani, S., et al. 2022, MNRAS, 517, Stark, D. P., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129                                                                                                                                                                                                                                                                                                                                                                                                                                                      |
| 5930 Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., & Shen, X. 2023a,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      |
| MNRAS, 526, 2665                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             |
| Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., et al. 2023b, ApJL, 955, L35 Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781                                                                                                                                                                                                                                                                                                                                                                                                                                            |
| Tacconi, L. J., Genzel, R., &Sternberg, A. 2020, ARA&A, 58, 157                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              |
| Teyssier, R. 2002, A&A, 385, 337 Thompson, T. A., Quataert, E., Zhang, D., & Weinberg, D. H. 2016,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           |
| 455, 1830                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    |
| Torrey, P., Vogelsberger, M., Marinacci, F., et al. 2019, MNRAS, 484,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        |
| MNRAS, 5587 898                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              |
| Tremonti, C. A., Heckman, T. M., Kau ff mann, G., et al. 2004, ApJ, 613,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     |
| Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110 Ucci, G., Dayal, P., Hutter, A., et al. 2023, MNRAS, 518, 3557                                                                                                                                                                                                                                                                                                                                                                                                                                                         |
| Vallini, L., Ferrara, A., Pallottini, A., Carniani, S., & Gallerani, S.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      |
| 2020, MN- RAS, 495, L22 Vallini, L., Ferrara, A., Pallottini, A., Carniani, S., & Gallerani, S.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              |
| 2021, MN- RAS, 505, 5543                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     |
| Vallini, L., Witstok, J., Sommovigo, L., et al. 2024, MNRAS, 527, 10 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Engineering, 13, 22                                                                                                                                                                                                                                                                                                                                                                                                                                |
| Science                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      |
| Van Rossum, G. &de Boer, J. 1991, CWI Quarterly, 4, 283 Van Rossum, G. &Drake, F. L. 2009, Python 3 Reference Manual (Scotts CA: CreateSpace)                                                                                                                                                                                                                                                                                                                                                                                                                                                |
| Valley, Venturi, G., Carniani, S., Parlanti, E., et al. 2024, arXiv                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |
| e-prints, arXiv:2403.03977 Vikaeus, A., Whalen, D. J., &Zackrisson, E. 2022, ApJL, 933, L8 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261                                                                                                                                                                                                                                                                                                                                                                                                                  |
| Wetzel, A., Hayward, C. C., Sanderson, R. E., et al. 2023, ApJS, 265, 44 Wilkins, S. M., Vijayan, A. P., Lovell, C. C., et al. 2023, MNRAS, 519, 3118                                                                                                                                                                                                                                                                                                                                                                                                                                        |
| York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  |
| Zahid, H. J., Dima, G. I., Kudritzki, R.-P., et al. 2014, ApJ, 791, 130                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      |
| Zahid, H. J., Kewley, L. J., &Bresolin, F. 2011, ApJ, 730, 137                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               |
| Zanella, A., Pallottini, A., Ferrara, A., et al. 2021, MNRAS, 500,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           |
| 118                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          |

## Appendix A: Modifications to the SFR time scales

In bathtub models (Lilly et al. 2013), the shape of the MZR is mostly independent from the functional form of the SFR, as the equilibrium solution is determined by the balance between infall and outflow, while the SFR uniquely determines the speed a model climbs the MZR from low to high M ⋆ . Instead, the SFR prescription can be important if a stochastic behavior is present in the system and when analyzing the FMR.

In the fiducial model presented in the paper we assumed that the SFR is proportional to the gas mass and is regulated by a fixed time scale for star formation, i.e. t ⋆ = 1 Gyr. As a natural alternatively for a fixed t ⋆ , a model can consider to adopt t ⋆ via empirical determinations (Tacconi et al. 2010) or modelbased predictions (Sommovigo et al. 2022) of the depletion time, which is expected to decrease with redshift, thus giving a faster star formation rate at higher z .

As a test, we modify our models by changing t ⋆ in our SFR prescription (eq. 2b) by adopting the scaling suggested by Sommovigo et al. (2022, see eq. 16 therein)

$$t _ { * } ^ { m o d } & \simeq \frac { 8. 9 5 \, G y r } { E ( z ) } & & ( \text{A.1} ) \\ & = \frac { 8. 9 5 \, G y r } { [ - 0. 2 4 + 0. 7 5 ( 1 + z ) ] \sqrt { \Omega _ { m } ( 1 + z ) ^ { 3 } + \Omega _ { \Lambda } } } \,, \\. &. &. &. &. &. &. & -$$

and setting the proportionality constant in order to match Tacconi et al. (2010) determinations at z = 0. As a reference, t ⋆ ≃ 20 t mod ⋆ at z ≃ 10, thus a faster stellar mass build is expected for the modified model.

As a benchmark, we select the halo with Mdm = 3 × 10 10 M ⊙ at z = 6 analyzed in Sec. 3.2 (see Fig. 2), evolve it with t ⋆ = 1 Gyr and t mod ⋆ from (eq. A.1), and plot the results in Fig. A.1 for both t d = 10 Myr and a reference no-SN-delay case.

The largest di ff erences between the t ⋆ = 1 Gyr and t mod ⋆ models are present for the stellar mass build up at z &gt; ∼ 10, when the galaxy has M ⋆ &lt; ∼ 10 7 M , The ⊙ t mod ⋆ model promotes a faster build up, particularly when a SN-delay is present, as the lack of feedback combined with the short SFR timescales make he galaxy reach M ⋆ &lt; ∼ 10 6 M ⊙ during the first burst. As the time progresses, di ff erences between the models tend to wash out, as the galaxy approaches M ⋆ ≃ 10 8 M , the SFH starts to be-⊙ come feedback-regulated, and at z ≃ 6, stellar, gas, and metal masses are roughly similar in all cases. Interestingly, while the amplitude of the flickering is similar ( σ SFR ≃ 0 5), the modula-. tion acts on slightly shorter time scales for the t mod ⋆ model, as a consequence of the faster SFR rise and the heavier o ff -balance caused by the stellar feedback. This implies that the t mod ⋆ model at z &lt; ∼ 7 5 gives an even enhanced dispersion for the MZR (from . σ Z ≃ 0 7 to . σ Z ≃ 1 7). Note that when no-SN-delay is present, . the SFH evolution is feedback regulated from the start, and the main di ff erence is an increased gas fraction for the t ⋆ = 1 Gyr model.

Note that adopting Sommovigo et al. (2022) implies a depletion time with a ∝ z -2 5 . redshift dependence, i.e. steeper with respect to the ∝ z -1 5 . dependence from Tacconi et al. (2010). Thus, by assuming Tacconi et al. (2010) instead of Sommovigo et al. (2022) to modify the SFR timescale, we expect smaller variation with respect to the fiducial t ⋆ = 1 Gyr model.

Fig. A.1: Di ff erence in the evolution of a Mdm = 3 × 10 10 M computed by prescribing the fiducial ( ⊙ t ⋆ = 1 Gry, left panel ) and faster ( t mod ⋆ , eq. A.1, right panel ) star formation timescale (eq. 2b). For both panels, Dark matter, stellar, gas, and metal mass evolution is reported for both the no-SN-delay and t d = 10 Myr model with the same notation of as in Fig. 2. Note that the left panel is exactly Fig. 2, it is reported here to have a clearer visual comparison.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-5a21a78a-9ccc-4a11-bf72-a3e1a85b7d99.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00061v1-figure-6d953d91-d62a-4d87-8254-ba08838189ef.jpg)
A. Pallottini A. Ferrara S. Gallerani L. Sommovigo S. Carniani L. Vallini M. Kohandel G. Venturi
The mass-metallicity relation as a ruler for galaxy evolution: insights from the James Webb Space Telescope
Galaxy evolution emerges from the balance between cosmic gas accretion,
fueling star formation, and supernova (SN) feedback, regulating the metal
enrichment. Hence, the stellar mass ($M_*$) - gas metallicity relation (MZR) is
key to understand the physics of galaxies. High-quality JWST data enable
accurate measurements of the MZR up to redshift z=10. Our aims are to
understand the observed MZR, its connection with the star formation rate (SFR),
the role played by SFR stochasticity, and how it is regulated by SN feedback.
We compare the MZR from the JADES, CEERS, and UNCOVER surveys, which comprise
about 180 galaxies at $z=3-10$ with $10^6<M_*/M_\odot<10^{10}$, with 200
galaxies from the SERRA cosmological simulations. To interpret the MZR, we
develop a minimal model for galaxy evolution that includes: cosmic accretion
modulated with an amplitude $A_{100}$ on 100 Myr; a time delay $t_d$ between
SFR and SN; SN-driven outflows with a varying mass loading factor
$\epsilon_{SN}$. Using our minimal model, we find the observed mean MZR is
reproduced by weak outflows ($\epsilon_{SN}=1/4$), in line with findings from
JADES. Matching the observed MZR dispersion requires $t_d=20$ Myr and a
$A_{100}=1/3$ modulation of the accretion rate. Successful models have low
stochasticity ($\sigma_{SFR}=0.2$), yielding a MZR dispersion of
$\sigma_{Z}=0.2$. Such values are close but lower than SERRA predictions
($\sigma_{SFR}=0.24$, $\sigma_{Z}=0.3$), clarifying why SERRA show no clear MZR
trend and some tension with the observations. As the MZR is very sensitive to
SFR stochasticity, models predicting high r.m.s. values ($\sigma_{SFR}=0.5$)
result in a ``chemical chaos'' (i.e. $\sigma_{Z}=1.4$), virtually destroying
the MZR. As a consequence, invoking a highly stochastic SFR
($\sigma_{SFR}=0.8$) to explain the overabundance of bright, super-early
galaxies leads to inconsistencies with the observed MZR.