Mini-Quenching of z = 4 − 8 Galaxies by Bursty Star Formation

The recent reported discovery of a low-mass z = 5 . 2 and an intermediate-mass z = 7 . 3 quenched galaxy with JWST / NIRSpec is the first evidence of halted star formation above z ≈ 5. Here we show how bursty star formation at z = 4 − 8 gives rise to temporarily quenched, or mini-quenched galaxies in the mass range M ⋆ = 10 7 − 10 9 M ⊙ using four models of galaxy formation: the periodic box simulation I llustris TNG, the zoom-in simulations VELA and F irst L ight and an empirical halo model. The main causes for mini-quenching are stellar feedback, lack of gas accretion onto galaxies and galaxy-galaxy interactions. The abundance of (mini-)quenched galaxies agrees across the models: the population first appears below z ≈ 8, after which their proportion increases with cosmic time, from ∼ 0 . 5 − 1 . 0% at z = 7 to ∼ 2 − 4% at z = 4, corresponding to comoving number densities of ∼ 10 − 5 Mpc − 3 and ∼ 10 − 3 Mpc − 3 , respectively. These numbers are consistent with star formation rate duty cycles inferred for VELA and F irst L ight galaxies. Their star formation histories (SFHs) suggest that mini-quenching at z = 4 − 8 is short-lived with a duration of ∼ 20 − 40 Myr, which is close to the free-fall timescale of the inner halo. However, mock spectral energy distributions of mini-quenched galaxies in I llustris TNG and VELA do not match JADES-GS-z7-01-QU photometry, unless their SFHs are artificially altered to be more bursty on timescales of ∼ 40 Myr. Studying mini-quenched galaxies might aid in calibrating sub-grid models governing galaxy formation, as these may not generate su ffi cient burstiness at high redshift to explain the SFH inferred for JADES-GS-z7-01-QU.


INTRODUCTION
The tight relation between stellar mass and star formation rate (SFR) of galaxies, known as the star-forming main sequence (MS), persists out to high redshift (Speagle et al. 2014;D'Silva et al. 2023;Popesso et al. 2023), even if the shape, scatter and normalisation have yet to be accurately determined above redshift z ≈ 5.For star-forming galaxies, the 'regulator' or 'bathtub' model posits that the SFR in a galaxy is controlled by the self-regulation of gas inflow, gas outflow, and gas consumption through star formation, providing a framework for understanding continuous star formation over long timescales on the MS (Bouché et al. 2010;Dekel et al. 2013;Lilly et al. 2013;Tacchella et al. 2016).Massive galaxies at and above the knee of the galaxy stellar mass function, then reduce their star-formation activity and leave the MS.Those galaxies are typically called 'quiescent' (or 'quenched') and increase significantly in abundance from redshift z ≈ 2 to today (Faber et al. 2007;Peng et al. 2010).
A plethora of mechanisms has been proposed to explain the ⋆ E-mail: td448@cam.ac.uk quenching of galaxies, each acting on different timescales and classes of galaxies.To match the observed galaxy number densities, internal mechanisms (Croton et al. 2006;Merlin et al. 2012;Sherman et al. 2020;Zinger et al. 2020;Tacchella et al. 2022b) are invoked at both the low-mass and high-mass end of the galaxy stellar mass function.Its suppression below the knee (M ⋆ ∼ 10 10.5 M ⊙ at z ≈ 4, Grazian et al. 2015) is typically accomplished via stellar feedback and above the knee via feedback from active galactic nuclei (AGN, Curtis & Sijacki 2016;Henden et al. 2018;Nelson et al. 2021), both of which expel gas and heat the circumgalactic medium.Some of these internal mechanisms, such as bulge formation, could lead to morphological quenching (Martig et al. 2009;Gensior et al. 2020;Lu et al. 2021;Shin et al. 2022).Since star formation preferentially occurs in gravitationally unstable gas discs, the stabilization thereof by the presence of the bulge (its growth and evolution regulated by AGN feedback) effectively reduces the star-formation activity.Especially at low redshift of z < 2 where galaxies can fall into dense clusters, environmental quenching effects (Dekel & Birnboim 2006;Peng et al. 2012;Ji et al. 2018;Contini et al. 2020;Whitaker et al. 2021;Williams et al. 2021) such as the loss or removal of gas from a galaxy due to ram pressure, tidal interactions and virial shock heating (all three processes on timescales of 100 Myr −2 Gyr) are believed to play a major role in regulating SFRs of low-mass (M ⋆ < 10 10 M ⊙ ) satellite galaxies.Since density contrasts between (proto-)clusters and the field are relatively small at high redshift of z > 2 (Overzier 2016), environmental quenching is less efficient.Nonetheless, quiescent galaxies have been detected at much higher redshift with their number density predicted to correlate inversely with redshift, down to ∼ 10 −5.5 Mpc −3 (CANDELS, Merlin et al. 2019) or ∼ 10 −4.5 Mpc −3 (JWST CEERS, Carnall et al. 2023a) at z ≈ 5. Most of the early quiescent galaxies that have been identified to date (Glazebrook et al. 2017;Carnall et al. 2020;Forrest et al. 2020;Valentino et al. 2020;Santini et al. 2021;Nanayakkara et al. 2022;Long et al. 2023) are massive (M ⋆ > 10 10 M ⊙ ), possibly due to observational limitations.Specifically, the observability of galaxies near a survey's limiting flux (typically bursty low-mass and/or high-redshift galaxies) can be highly time-dependent due to the SFR variability (Sun et al. 2023).
Can the population of low-mass high-redshift quiescent galaxies be explained with simple models of bursty star formation?At high redshift and lower stellar masses, an equilibrium between stellar feedback and gravity cannot be sustained, leading to bursty star formation (Anglés-Alcázar et al. 2017;Faucher-Giguère 2018).Therefore, although the galaxy population might still follow the MS at high redshift, their trajectory about the MS could be dominated by short-term bursts of star formation triggered by various processes such as the merging of galaxies or the inflow of gas (Rodríguez-Puebla et al. 2016;Tacchella et al. 2016Tacchella et al. , 2020)).These bursts are separated by periods of relative quiescence in which galaxies fall off the star-forming MS.It is therefore of great interest to infer SFRs that probe different timescales in order to observationally constrain the burstiness of low-mass and high redshift galaxies (Weisz et al. 2012;Emami et al. 2019;Faisst et al. 2019).
The discoveries by Looser et al. (2023) and Strait et al. (2023) motivated us to investigate the properties of high-redshift quiescent galaxies with low stellar masses (M ⋆ = 10 7 − 10 9 M ⊙ ) using theoretical models.Specifically, we use four galaxy formation models: the periodic box simulation IllustrisTNG (Pillepich et al. 2018a(Pillepich et al. , 2019;;Nelson et al. 2019), the zoom-in simulation VELA (Ceverino et al. 2014;Zolotov et al. 2015) and an empirical halo model (EHM; Tacchella et al. 2018).These four models allow us to study this galaxy population from complementary viewpoints since they provide access to different dynamic ranges and timescales.The combination of those four models allows us to cross-check the consistency of predic-tions on the theory side.All four models produce quiescent galaxies with M ⋆ = 10 7 − 10 9 M ⊙ at z > 4.However, these galaxies are only temporarily quiescent and rejuvenate back onto the MS.In order to differentiate this process from permanent 1 star formation quenching at higher stellar masses (M ⋆ > 10 10 M ⊙ ), we refer to this temporary quenching as mini-quenching.We show in this work that the main causes for mini-quenching are stellar feedback, lack of gas accretion onto galaxies, mergers and tidal interactions.The abundance of these mini-quenched galaxies (MQGs) is in good agreement between the four galaxy formation models.However, we find that the level of burstiness inferred from the four models on short timescales ∼ 40 Myr is lower than what is observed.
The organisation of the paper is as follows: In Secs.2.1, 2.2, 2.3 and 2.4, we describe the four galaxy formation models: IllustrisTNG, EHM, FirstLight and VELA.While the periodic box simulation IllustrisTNG provides a statistical sample for estimating abundances, FirstLight and VELA zoom-ins of galaxies have more resolved star-formation histories (SFHs).We strive to treat all four models on an equal footing and present their differences in Sec.2.5.Details on how we calculate stellar masses and star-formation rates can be found in Secs.2.6, 2.7 and 2.8.Methods to generate mock spectral energy distributions (SEDs) are described in Sec.2.9.In Sec. 3, we highlight that the properties, abundances and timescales of mini-quenching events at high redshift are roughly consistent across all four models, which we show by providing quantitative estimates.We compare simulated SEDs to JADES-GS-z7-01-QU photometry in Sec. 4. We discuss our main findings in Sec. 5.

IllustrisTNG Simulations
In order to study a representative sample of galaxies, we use both IllustrisTNG50 and IllustrisTNG100.These simulations were performed with the state-of-the-art code Arepo described by Springel (2010) and Weinberger et al. (2020).The hydrodynamical equations are solved on a moving Voronoi mesh using a finite volume method.Various astrophysical processes such as metalline cooling, star formation and feedback remain unresolved in IllustrisTNG and are approximated by subgrid models (Pillepich et al. 2018b).Gas above a density threshold of n H ∼ 0.1 cm −3 forms stars stochastically following the empirical Kennicutt-Schmidt relation and assuming a Chabrier (Chabrier 2003) initial mass function (IMF), see Table 1.
The dark matter particle mass resolution is m DM = 7.5×10 6 M ⊙ in TNG100 and m DM = 4.5×10 5 M ⊙ in TNG50.Dark matter haloes are identified using the friends-of-friends (FoF) algorithm with a standard linking length of b = 0.2 × (mean inter-particle separation) (Springel et al. 2001a).Within each FoF halo, sub-haloes identified by the SUBFIND (Springel et al. 2001b) algorithm are made up of all the resolution elements (gas, stars, dark matter, and black holes) which are gravitationally bound to the subhalo.In our framework, galaxies are subhaloes with at least ∼ 200 stellar particles, i.e.M ⋆,min = 2 × 10 7 M ⊙ for TNG50 and M ⋆,min = 2 × 10 8 M ⊙ for TNG100.The sample includes both central and satellite galaxies.We also perform a combined simulation study of TNG100+50, with appropriate weights when calculating comoving number densities  2023) is given in the form of a common parametrization (dashed curves) below z ≤ 6 for M ⋆ ≥ 10 8.5 M ⊙ .The MS ridges in the four models are consistent to within ∼ 30%, with TNG100+50 and FirstLight10+20 predictions closest to observational constraints. of MQGs.We disregard TNG300 since sufficiently resolved galaxies (M ⋆ > 2 × 10 9 M ⊙ ) are not expected to exhibit mini-quenching events in the simulation (see Sect. 3.5), even though some do exhibit AGN-induced quiescence by z ≈ 4.2 (Hartley et al. 2023).

Empirical Halo Model
Empirical models are galaxy formation models whose physical prescription is largely motivated and calibrated based on observations.The workings for semi-analytical and empirical models have become increasingly similar, though the former solves physical equations, while the later focuses on effective prescriptions.Empirical models are now successful in describing the galaxy population over a wide range of redshifts (Rodríguez-Puebla et al. 2016;Behroozi et al. 2019;Tacchella et al. 2013Tacchella et al. , 2018)).We will focus on the empirical model introduced in Tacchella et al. (2018), based on halo merger trees extracted from the COLOR simulations (Hellwing et al. 2016) in a (70.4 Mpc/h) 3 box with dark matter particle mass resolution m DM = 8.8×10 6 M ⊙ .Therein, the SFR of a galaxy is assumed to be proportional to the gas accretion rate of its parent halo, Ṁgas , normalised by a redshift-independent efficiency, ϵ(M h ), of converting gas into stars, (1) The shortest timescale the model can probe, ≈ 0.05 t H , is 2 − 3 times shorter than the dynamical timescale of the halo (∼ 0.15 t H , Peebles 1980).The model clearly misses many effects such as stellar feedback which can happen on even shorter timescales.At z = 7, the shortest probed timescale is thus 0.05 t H ∼ 42 Myr.We impose M ⋆,min = 2 × 10 7 M ⊙ in accordance with TNG50.

FirstLight Simulations
The FirstLight simulations (Ceverino et al. 2017) were run with the Adaptive Mesh Refinement (AMR) code ART (Kravtsov et al. 1997;Kravtsov 2003;Ceverino & Klypin 2009).Besides gravity and hydrodynamics, on sub-grid level the code incorporates gas cooling due to atomic hydrogen and helium, metal and molecular hydrogen cooling, photoionization heating by a constant cosmological UV background with partial self-shielding, star formation and feedback (thermal + kinetic + radiative), as described in Ceverino et al. (2018).The parent haloes were selected at z ∼ 5 from N-body simulations with box sizes 10 Mpc/h and 20 Mpc/h such that their maximum circular velocity V max lies between 50 and 250 km/s.This range excludes very massive and rare haloes with number densities lower than ∼ 3 × 10 −4 h 3 Mpc −3 , as well as small haloes in which galaxy formation is inefficient.The dark matter particle mass resolution is m DM = 10 4 M ⊙ .In this work, we perform a joint study of zoom-ins based on both box sizes 10 Mpc/h and 20 Mpc/h, to which we will refer to as FirstLight10+20.

VELA Simulations
This set of hydrodynamic simulations was likewise performed with the ART code, referred to as the VELA runs.The main features can be summarised as follows.The Eulerian gas dynamics is followed using an adaptive mesh refinement (AMR) approach.The dark matter particle mass resolution is m DM = 8 × 10 4 M ⊙ while the AMR maximum resolution is 17 − 35 pc at all times.In the circumgalactic medium (at the virial radius of the dark matter halo), the median resolution amounts to ∼ 500 pc.The virial masses of the 29 galaxies are chosen to be in the range M vir = 2 × 10 11 − 2 × 10 12 M ⊙ at z = 1 about a median of 4.6 × 10 11 M ⊙ .Beside gravity and hydrodynamics, the code includes many physical processes relevant for galaxy formation: gas cooling by atomic hydrogen and helium, metal and molecular hydrogen cooling, photoionization heating by the UV background with partial self-shielding, star formation, stellar mass loss, metal enrichment of the ISM and stellar feedback.Supernovae and stellar winds are implemented by local injection of thermal energy as described in Ceverino & Klypin (2009); Ceverino et al. (2010Ceverino et al. ( , 2012)).Radiative stellar feedback is implemented at a moderate level Ceverino et al. (2014).AGN feedback and feedback associated with cosmic rays and magnetic fields are not implemented.Note that we do not impose a M ⋆ = 10 7 − 10 9 M ⊙ selection on VELA and FirstLight galaxies when calculating SFR duty cycles and mini-quenching timescales in Secs.3.3 and 3.4.

Model Differences
While IllustrisTNG aims to realistically capture the formation and evolution of the large-scale structure and the galaxies embedded therein, the simulations required significant computational resources (TNG100 alone took 18.62 million CPU core hours, Nelson et al. 2018).In contrast, EHM is much simpler and provides a flexible framework for modelling SFRs which can be calibrated to fit observational data.However, the evolution of individual galaxies cannot be followed with fine spatio-temporal resolution, a shortcoming that the zoom-in simulations VELA and FirstLight address.In Table 1, we compare the cosmologies and IMFs assumed in the four galaxy formation models.The value of Ω m in the Planck 2015 cosmology (Ω m = 0.3089) is significantly higher than the WMAP5 (Ω m = 0.258) value, while WMAP7 (Ω m = 0.273) is in between the two.We estimate the global abundance of (mini-)quenched galaxies from IllustrisTNG and EHM, finding consistent results in Sec.
3.2 despite the differences in Ω m .Duty cycle estimates of individual galaxies in VELA and FirstLight (both adopting WMAP5) depend more on subgrid modelling than cosmology and are found to be consistent with these global abundance estimates as well, see Sec. 3.3.While the Chabrier IMF provides a better fit to observations of lowmass stars and brown dwarfs in the Galactic disc than the Salpeter IMF (Chabrier 2003), neither is well-motivated at high redshift of z > 4 where background radiation fields, gas temperatures and densities are considerably different (Riaz et al. 2021).

Averaging Timescales
Both simulations and observations consistently demonstrate that using shorter SFR averaging timescales results in a higher normalization and increased scatter of the main sequence, particularly at the low-mass end (Schaerer et al. 2013;Hayward et al. 2014;Speagle et al. 2014;Sparre et al. 2015;Caplar & Tacchella 2019;Donnari et al. 2019).In observational studies, part of this effect can be attributed to a sampling bias towards stars in their young evolution-ary stages or those located in regions with intense star formation activity.Our choice of ∆t = 10 Myr for IllustrisTNG, VELA and FirstLightis on the low end of timescales, yet is necessary to resolve SFHs of galaxies with high enough fidelity to infer mini-quenching timescales.For EHM, we choose the shortest timescale we can probe in this model, ∆t = 0.05 t H , as the natural SFR averaging timescale.

Aperture Choice
Stellar mass and SFR of a galaxy depend on the radius within which they are computed.We aim to establish consistency with Donnari et al. (2019), who have also demonstrated that too small apertures lead to an underestimate of the galaxy SFR, and Merlin et al. ( 2019) as well as to mimic as accurately as possible the observational approach.To that end, in IllustrisTNG and VELA, we use the values estimated within the 3D spherical galactocentric distance that corresponds to twice the stellar half mass radius 2 × R star,h of each galaxy.Recall that for IllustrisTNG, there is an additional gravitational boundedness criterion imposed under the hood by the SUB-FIND algorithm.For FirstLight, galaxy stellar mass and SFRs are estimated within 0.15 R vir from the center of the parent halo, where R vir is the virial radius of the halo as per the spherical collapse result (Bryan & Norman 1998).

Identification of the Star-Forming Main Sequence
To define the star-forming MS in IllustrisTNG and EHM, we follow Donnari et al. (2019) and iteratively remove quiescent galaxies until the median SFR in a given mass bin converges.The MS is thus allowed to deviate from a log-linear trend ('bending MS', Whitaker et al. 2014;Donnari et al. 2019), as in fact it does by exhibiting a turnover at about M ⋆ ≳ 10 10.5 M ⊙ (Lee et al. 2018;Tomczak et al. 2016;Popesso et al. 2023).Depending on how a star-forming galaxy is defined observationally, this trend may persist up to z ≈ 6.The bending MS estimation is performed at z = 4, 5, 6, 7, 8 using stellar mass bins of width 0.25 dex with results shown in Fig. 1.We recover the well-known fact that the normalization increases with redshift at fixed stellar mass.This can be traced to the stellar mass growth time-scale (the ratio of stellar mass to star formation rate) which is expected to be comparable to the Hubble time (Ma et al. 2018;Ceverino et al. 2018).As seen in Fig. 1, this trend becomes weaker towards higher redshift.To obtain SFHs in between these redshifts, we interpolate the MS.In FirstLight10+20 and VELA, the MS ridge is parametrized as (2) For VELA, the best-fit parameters are s b = 0.046, β = 0.14 and µ = 5/2 (Tacchella et al. 2016) while for FirstLight10+20 we find s b = 0.064, β = 0.05 and µ = 2.45.Across the mass range M ⋆ = 10 7 − 10 10 M ⊙ , the MS ridges are only consistent to within ∼ 30% between the four galaxy formation models at any given redshift, see Fig. 1.Compared to observations of the star-forming MS as compiled by Popesso et al. (2023), while TNG100+50 and FirstLight10+20 are in good agreement therewith, EHM and VELA typically have smaller normalizations.This deviation is evident for VELA galaxies whose stellar-to-virial mass ratios are higher than deduced from observations (Ceverino et al. 2014).Quiescent (star-forming) galaxies are those whose SFR falls below (above) a certain relative distance from the median SFR at the corresponding mass.One popular choice is to define as quiescent those galaxies whose logarithmic specific star-formation rate ).The mini-quenching event of V5 (cf.Fig. 5) is highlighted with a triangle symbol.The abundance of (mini-)quenched galaxies agrees across the models, with an overall fraction of f QG = 0.5 − 1.0% at z = 7.
is sSFR< 10 −11 yr −1 (Donnari et al. 2019;Merlin et al. 2019).Alternatively, some authors adopt a threshold such as 2 × σ MS ∼ 0.6 dex (Tacchella et al. 2016) from the MS.However, here we choose a selection criterion of 0.75 dex from the MS at all redshifts unless explicitly noted.This criterion should be seen as a necessary though not sufficient condition for mini-quenching.For VELA and FirstLight galaxies whose SFHs we can resolve well, we further require that the quenching is only temporary.Mini-quenching timescales as well as duty cycles are quantified in Secs.3.3 and 3.4, respectively.

Calculating Spectral Energy Distributions
To calculate a dust-free SED of a galaxy, we treat each stellar particle in the respective simulation as a simple stellar population (SSP) using a stellar population synthesis method.To correct for dust attenuation and extinction, we follow Nelson et al. (2018).When applied to IllustrisTNG, the model accurately reproduces the observed distribution of optical (g−r) colors from the Sloan Digital Sky Survey.In addition to adopting a simple powerlaw extinction model (Charlot & Fall 2000) for the attenuation by finite-lifetime birth clouds surrounding young stellar populations as well as the ambient diffuse ISM, we follow the distribution of metals and (neutral) hydrogen gas in and around each simulated galaxy.
To obtain neutral hydrogen fraction estimates for star-forming cells in IllustrisTNG, we post-process the outputs and recalculate (following Villaescusa-Navarro et al. 2018) the equilibrium fractions according to the Springel & Hernquist (2003) model to account for the multiphase interstellar medium, including the presence of molecular hydrogen, H 2 .For VELA, we instead employ the total hydrogen column density to estimate the optical depth.The Nelson et al. (2018) resolved dust model then attributes a neighborhoodand viewing angle-dependent attenuation to each stellar particle.We choose the viewing angle randomly as one of the vertices of the N s = 1 HEALPIX sphere (Górski et al. 2005) oriented in simulation coordinates.

WHAT IS MINI-QUENCHING?
By mini-quenching we refer to galaxies in the mass range M ⋆ = 10 7 − 10 9 M ⊙ (such as the one reported by Strait et al. 2023;Looser et al. 2023) that are temporarily quiescent.In this stellar mass range, star formation is bursty and regulated mainly by fluctuations in the gas inflow, stellar feedback and environmental effects such as gas-rich mergers.In IllustrisTNG, black holes are seeded when M DM = 5 × 10 10 h −1 M ⊙ and AGN feedback is thus designed to be inefficient at the low-mass end (Weinberger et al. 2017).However, overmassive black haloes in the centers of dwarfs can give rise to efficient AGN feedback at high redshift without violating observed HI gas mass constraints (Koudmani et al. 2022), though direct observational evidence of such feedback channels is still needed.There is some contribution from photoheating in the presence of an ionizing background (Okamoto et al. 2008;Pawlik & Schaye 2009;Brown et al. 2014), but Wu et al. (2019) report that in self-consistent radia-Table 2. MS scatter σ MS at various redshifts z = 4 − 8, averaged across M ⋆ = 10 7 − 10 9 M ⊙ and decomposed into stellar mass bins of width 0.5 dex.All estimates are given at redshifts z = 4, 5, 6, 7, 8 for TNG100+50, EHM and FirstLight10+20.Results at z ≤ 5 for FirstLight10+20 are not available.In accordance with quenched galaxy number densities in Fig. 2 (first two panels), we calculate a weighted standard deviation to avoid overcounting in the mass range M ⋆ > 2 × 10 8 M ⊙ into which both TNG100 and TNG50 galaxies fall.2018) find that M ⋆ ∼ 10 8 M ⊙ is the (weakly redshift dependent) mass transition threshold where SFHs begin to transition from bursty to stable.Note that the physical mechanisms underlying the transition from bursty to steady star formation are complex (Sparre et al. 2017;Hopkins et al. 2023;Gurvich et al. 2023) and might be related to the virialization of the inner circumgalactic medium (Stern et al. 2021).
As mentioned in the introduction, mini-quenching does not refer to long-term quenching of galaxies with M ⋆ > 10 10 M ⊙ .This form of quenching of high-mass galaxies is commonly observed at lower redshift, but see Glazebrook et et al. (2023b).The primary factor governing star formation within high-mass systems is likely radio mode feedback from supermassive black holes, observable via radio lobes, X-ray cavities and radiooptical correlations (Kormendy & Ho 2013;Terrazas et al. 2020;Houston et al. 2023).
In the following, we compare various statistics between IllustrisTNG, EHM, FirstLight and VELA to demonstrate that their predictions for high-redshift mini-quenching are in good agreement with each other.

Scatter around the MS
The distribution of galaxies around the star-forming MS (see Sec. 2.8) in the various galaxy formation models is shown in Fig. 2. Quiescent galaxies suffer from a higher SFR uncertainty, which we accommodate by adding error estimates from bootstrapping to the TNG100+50 results.Specifically, for each galaxy we bootstrap over the initial masses of stellar particles that formed in the last 10 Myr (averaging timescale) before z = 7.We assume a scatter of m ⋆ /2 if only 1 stellar particle of mass m ⋆ formed over the last 10 Myr.Note that upper and lower error bars are asymmetric.
We find that the scatter σ MS around the MS (as quantified in Ta-ble 2) decreases towards higher stellar mass2 .For instance, at z = 7 TNG100+50 galaxies in the stellar mass range M ⋆ = 10 7 − 10 7.5 M ⊙ have a scatter of σ MS = 0.18 while those of mass M ⋆ = 10 8.5 − 10 9 M ⊙ have σ MS = 0.15.This trend is in agreement with results from FIRE-2 (Ma et al. 2018) and is a result of bursty star formation regulated mainly by stellar feedback and environmental effects.However, the decrease of the scatter towards higher mass is less pronounced in Flares simulations (Lovell et al. 2022), which resolve galaxies above M ⋆ ∼ 10 9 M ⊙ .Flares employs a physically motivated model for AGN feedback that takes into account the dynamics of the accretion disk and the surrounding gas as well as the radiation emitted by the AGN, which better reproduces the observed properties of massive galaxies at high redshifts such as sizes, masses and stellar populations (Vijayan et al. 2020;Roper et al. 2022).
The most common approach to infer the evolution of the starforming MS observationally is by combining galaxy samples at different redshifts for which the SFRs are estimated for every galaxy.More robust determinations of total SFRs and the portions associated with the unobscured, SFR FUV , and obscured, SFR IR , regimes bridge the gap with other constraints such the galaxy stellar mass function of star-forming galaxies and the FUV and IR luminosity functions (Rodríguez-Puebla et al. 2020).UV-derived SFRs tend to indicate a scatter of ∼ 0.25 − 0.30 dex (Elbaz et al. 2007;Whitaker et al. 2012;Speagle et al. 2014;Fang et al. 2018).Driven by stochastic, bursty SFHs at the low-mass end and the presence of bulges, bars but also AGN at the high-mass end, the dispersion can increase to ∼ 0.4 dex (Santini et al. 2017;Popesso et al. 2019;Guo et al. 2015).According to Guo et al. (2015); Willett et al. (2015); Davies et al. (2019), the scatter might follow a minimum vertex parabolic 'U'-shape decreasing with stellar mass from log(M ⋆ /M ⊙ ) ∼ 8 − 10 and then increasing at log(M ⋆ /M ⊙ ) ∼ 10 − 11.5.Note though that few studies go above z ≈ 5.
IllustrisTNG likely underestimates the scatter around the MS at higher mass compared to observations given that the AGN feedback is modeled as a simplified subgrid process, where the (fully isotropic) energy and momentum injection from the AGN is regu-Table 3. Fraction and abundance of quenched galaxies (QGs) in the stellar mass range M ⋆ = 10 7 − 10 9 M ⊙ : (1) fraction of quenched galaxies in percent, with error estimates from bootstrap resampling; (2) fraction of environmentally mini-quenched galaxies (EMQGs) among all quenched galaxies in percent; (3) comoving number density of quenched galaxies in Mpc −3 .All estimates are given at redshifts z = 4, 5, 6, 7, 8 for TNG100, TNG50, a joint analysis of TNG100 and TNG50 as well as EHM.The fraction of EMQGs cannot be calculated within EHM and is thus not shown.lated by a set of rules based on the accretion rate and black hole mass (Weinberger et al. 2017).The EHM, VELA, FIRE-2 and FirstLight models do not include feedback from AGN and likewise cannot fully capture the complex feedback processes in massive galaxies.
For EHM (Fig. 2, second panel), we in addition see that the overall number of galaxies is larger than in the case of TNG100+50 (17226 vs 2369 at z = 7 in the range M ⋆ = 10 7 − 10 9 M ⊙ ).EHM might indeed overpredict the galaxy stellar mass function at the low-mass end as hinted at by Tacchella et al. (2018) whilst the NIR band luminosity function and by extension the galaxy stellar mass function predicted by IllustrisTNG are largely consistent with observations at high redshift (Shen et al. 2022).Despite that, the total EHM scatter σ MS = 0.21 dex around the MS agrees well with the scatter σ MS = 0.17 dex found for TNG100+50, see Table 2.The MS scatter in EHM comes directly from the scatter in the halo mass accretion rate.FirstLight10+20 exhibits a higher MS scatter σ MS = 0.38 dex than IllustrisTNG and EHM, which can be traced to FirstLight being able to resolve intense SF bursts on timescales smaller than ∼ 10 Myr (also see Ceverino et al. 2018).
In Table 2, we also quantify the redshift evolution of the total and stellar mass bin-decomposed MS scatter.TNG100+50, EHM and FirstLight10+20 galaxies all exhibit a larger scatter towards lower redshift.This is in agreement with Ma et al. (2018); Caplar & Tacchella (2019) and is a result of hierarchical structure formation which leads to a more diverse range of galaxy properties at later epochs.At lower redshift, we also probe a wider range of largescale environments such as galaxy clusters, which contributes to the increased scatter in the star-forming MS.

Abundance of (Mini-)Quenching Events
How common are (mini-)quenching events at high redshift z = 4−8?While each of the four galaxy formation models probes a slightly different stellar mass range, we find that estimates from all four are largely consistent with a picture in which the population first appears below z ≈ 8, after which the fraction of (mini-)quenched galaxies increases with cosmic time, from ∼ 0.5 − 1.0% at z = 7 to ∼ 2 − 4% at z = 4, corresponding to comoving number densities of 10 −5 Mpc −3 and 10 −3 Mpc −3 , respectively.Fig. 3 compares the fraction of quenched galaxies in IllustrisTNG and EHM.A combined analysis of TNG100+50 yields an overall quenched galaxy fraction of f QG = 0.6% at z = 7, corresponding to a comoving number density of 8.7 × 10 −6 Mpc −3 .For a comparison of redshifts z = 4, 5, 6, 7, 8 see Table 3.For  (dashed).The population of quenched galaxies first appears below z ≈ 8, after which their fraction increases with cosmic time, from ∼ 0.5 − 1.0% in the mini-quenching mass range M ⋆ = 10 7 − 10 9 M ⊙ at z = 7 to ∼ 2 − 4% at z = 4.The error bar on each estimate denotes the standard deviation as obtained with bootstrap resampling.The overall mass-agnostic quenched galaxy fraction is shown in Table 3.There is qualitative agreement between IllustrisTNG and the simplistic EHM approach, even though the latter does not model stellar feedback.TNG100, we find that at z = 7 an overall f QG = 1.0% of galaxies are mini-quenched (8 out of 826 galaxies), translating 3 to a comoving number density of 5.9 × 10 −6 Mpc −3 .In TNG50, the fraction of lowest-mass quenched galaxies around M ⋆ ∼ 2 × 10 7 M ⊙ at z = 7 reads 0.7% while the overall one is f QG = 0.3% (4 out of 1543 galaxies).The comoving number density is thus 3.6 × 10 −5 Mpc −3 , which is higher than the TNG100 value since the minimum stellar mass M ⋆,min = 2 × 10 7 M ⊙ is lower.
Since the dependence on stellar mass is less pronounced in EHM and TNG100 than in TNG50, Fig. 3 also demonstrates that EHM (based on COLOR merger trees) and TNG100 explore a different dynamic range than TNG50, and better sample the intermediate-mass range of the galaxy stellar mass function.Since EHM might overpredict the galaxy stellar mass function at the low-mass end (cf.Sec.3.1), quenched galaxy number densities are higher by about 1/2 dex than in the combined TNG100+50 analysis, see Table 3.

Duty Cycles in FirstLight and VELA
It is straightforward to look at the evolution of selected zoom-in galaxies around the ridge of the star-forming MS.In Fig. 2 (right panel), we show SFHs of two VELA galaxies with bursty SFHs (V5 and V10) and one galaxy with a more stable SFH (V32).However, extracting MQG number densities from zoom-ins is ill-defined4 .How do we know whether the abundance of MQGs in FirstLight and VELA is consistent with IllustrisTNG and EHM?To this end, it is illustrative to investigate the SFR duty cycle.It denotes the fraction of time spent in an active phase, i.e., the ratio between the actively star forming (∆ MS > [−0.75, −0.6] dex) time interval ∆t on , and the time elapsed between the first star formation event t form and the epoch of observation t obs (see Gelli et al. 2023), In contrast, FirstLight10+20 galaxies display systematically lower values of f duty compared to VELA.At z = 5, the duty cycle Figure 4. Bursty SFHs of FirstLight10+20 galaxies around 12 selected mini-quenching events in the redshift range z = 4 − 8.Each trajectory shows the deviation from the MS from 100 Myr before until 100 Myr after the miniquenching event.One galaxy undergoes two mini-quenching events separated by ∼ 30 Myr.The selection threshold of 0.75 dex below the MS is shown in dashed black.Mini-quenching timescales τ MQ in Table 5 are obtained by summing up the time spent below this threshold.We color-code trajectories by the value of τ MQ : the higher τ MQ , the darker the blue shade.
drops as low as ∼ 91%.Some of this discrepancy can be attributed to FirstLight galaxies benefiting from an eight-fold increase in resolution compared to VELA galaxies.This enhanced resolution allows FirstLight to resolve many lower-mass systems, which tend to exhibit more bursty star formation behavior (see Fig. 2).However, even when focusing the f duty analysis on galaxies within a specific mass range, the lower tail for FirstLight galaxies remains substantial.For instance, at z = 6, when restricting the analysis to galaxies with M ⋆ ∈ [10 8 − 10 9 ] M ⊙ , the duty cycle is estimated to be 98.5 +0.5  −0.8 % for VELA and 99.9 +0.1 −15.5 % for FirstLight10+20.This indicates that the difference in duty cycle between FirstLight10+20 and VELA remains significant, even within the confines of a specific mass range where both VELA and FirstLight10+20 galaxies are well resolved.While the duty cycle of a galaxy also depends on large-scale environments (see Sec. 3.5), we thus conclude that the stronger feedback in FirstLight (see Ceverino et al. 2018) compared to VELA is also reflected in lower duty cycles.
The SFR duty cycle estimates are in good agreement with MQG number densities obtained for IllustrisTNG (see Table 3).At z = 7, 0.6% of galaxies are mini-quenched in a joint analysis of TNG100 and TNG50.Assuming this fraction remains constant at all z > 7, this would correspond to a duty cycle of f duty = 99.4% compared to the value f duty = 98.6 +1.4 −4.1 % we infer for VELA and f duty = 99.9 +0.1 −37.9 % for FirstLight10+20.However, the f duty values inferred are not consistent with those from the SERRA simulations ( f duty ∼ 0.60 − 0.99, Pallottini et al. 2022;Gelli et al. 2023), a discrepancy which we will elaborate on in Sec.4.4.

How Long is a Mini-Quenching Event?
In the compaction-triggered quenching model (Dekel & Burkert 2014;Zolotov et al. 2015;Lapiner et al. 2023), galaxies at z ≈ 2 − 3 undergo three evolutionary phases: cold gas accretion, compaction and post-compaction, and quenching.Before this final successful quenching attempt, galaxies oscillate about the MS ridgeline on timescales of ∼ 0.4 t H , as found by Tacchella et al. (2016) in the VELA simulation suite at z = 2 − 4. Houston et al. (2023) found that the oscillation period increases towards higher stellar mass galaxies to about 1.98 ± 2.27 Gyr for galaxies with M ⋆ = 10 11.5 − 10 12 M ⊙ .
Here we are interested in a different but related timescale.Given a SFH, we estimate the accumulated time each MQG spends 0.75 dex or further below the MS, which we define as the mini-quenching timescale τ MQ .
Many SFHs of MQGs in IllustrisTNG at high redshift trace out an irregular pattern (especially) around the mini-quenching event (±20 Myr), suggesting a bursty mode of star formation.However, SFHs are poorly resolved, with many displaying a distinctive Vshaped feature around the mini-quenching event.We thus only capture a glimpse into said event and cannot put meaningful constraints on τ MQ .A similar challenge is encountered with EHM which is based on halo merger trees extracted from the COLOR simulations in a (70.4 Mpc/h) 3 box.By tracking the main progenitor and main descendant branch of halos, Tacchella et al. (2018) showed that in this empirical framework many galaxies oscillate around the MS at z > 4, hinting at bursty star formation as a result of halo mergers.However, due to limitations in resolution, EHM does not provide a robust framework for estimating τ MQ either.
For FirstLight10+20 galaxies, we show SFHs around 12 selected mini-quenching events in Fig. 4. The well-resolved trajectories demonstrate the bursty, stochastic nature of star formation at high redshift.Mini-quenching events typically last longer than in VELA as expected from the lower SFR duty cycles.The distribution is centered around τ MQ ∼ 40 Myr with a heavy upper tail.In fact, the longest mini-quenching events can last up to τ MQ ∼ 100 Myr before the galaxy rejuvenates.
Expressed in terms of the Hubble time t H (z = 7) ∼ 760 Myr, the typical timescale at z = 7 is τ MQ ∼ 0.02 t H .This timescale is more than an order of magnitude smaller than the oscillation timescale ∼ 0.4 t H about the ridge of the MS, indicating that the SFR needs to change significantly on shorter timescales around the miniquenching event (with a more steady evolution closer to the MS).In the framework of correlated stochastic processes (Kelson 2014;Caplar & Tacchella 2019;Abramson & Kelson 2020;Tacchella et al. 2020;Iyer et al. 2020Iyer et al. , 2022)), this might suggest different amounts of power on different temporal scales.
The mini-quenching timescale τ MQ ∼ 20−40 Myr is several times (∼ 1/6) shorter than the overall free-fall time t ff = (3π/(32Gρ)) 1/2 ∼ 240 Myr of a halo of mass M h = 6 × 10 9 M ⊙ , hosting a galaxy such as V5 of mass M ⋆ = 2.1×10 7 M ⊙ (see Sec. 4) based on the observed stellar-to-halo-mass relation (Girelli et al. 2020).However, we argue that τ MQ is close to the local free-fall timescale t ff,loc of the inner halo where galaxies typically reside.In this picture, after the expulsion of gas in the wake of stellar feedback the galaxy finds itself in a state of mini-quenching until gas gets reaccreted/falls onto the galaxy.The mini-quenching timescale is influenced by the mass distribution of the halo, the density and temperature of the gas, and feedback effects such as radiative winds (Stern et al. 2021;Gelli et al. 2023).At high densities and low metallicities, t ff,loc can be shorter than the stellar feedback timescales (Dekel et al. 2023), in which case τ MQ will be the closest to t ff,loc .

What Drives Mini-Quenching?
The fluctuation of SFRs from a state of mini-quenching to an episodic burst is typically mediated by the interplay between gasrich mergers and the steady influx of cold gas streams, periodically inhibited by the feedback from evolved stars.EHM captures mergers insofar that an increased dark matter accretion rate is assumed to result in enhanced gas accretion.IllustrisTNG, FirstLight and VELA take account of stellar feedback and thus capture bursty star formation processes on smaller temporal scales.
Even though the effects of ram pressure and virial shocks are less efficient at high redshift (Fujita 2001;Birnboim & Dekel 2003;Maier 2021), in IllustrisTNG we find that tidal interactions during close galaxy-galaxy encounters (which are more common at high redshift than today) play an important role in determining SFRs of (especially) low-mass galaxies.We follow the main progenitor and main descendant branches of each galaxy 5 around the mini-quenching event, and search for discontinuities in the stellar mass evolution which are not matched by a corresponding change in SFRs.If said discontinuity occurs right before or after the miniquenching event, we label the galaxy an environmentally miniquenched galaxy (EMQG).We refrain from adopting a more rigorous merger-only selection based on merger trees since we find that several quenched galaxies are tidally distorted (sometimes with an accompanying decrease of stellar material falling inside the aperture) without fully merging.
The fraction of EMQGs among all quenched galaxies (cf.Table 3) is highest for galaxies on the low-mass end (TNG50) and increases with redshift (up to # EMQG/# QG × 100 ∼ 66.7% at z = 7).For intermediate-mass galaxies sampled from TNG100, the redshift trend is reversed, and peaks at z = 4 with 25.5%.In a combined analysis of TNG100+50, we find an EMQG fraction of 27.3%.Tidally induced mini-quenching events such as complete and incomplete mergers thus play an important role in the context of bursty star formation.

MOCK SPECTRAL ENERGY DISTRIBUTIONS
COMPARED TO JADES-GS-Z7-01-QU To bridge the gap with observations, we search for the VELA miniquenching event that is most similar to JADES-GS-z7-01-QU (z = 7.3 and M ⋆ = 5 × 10 8 M ⊙ ).VELA galaxies at z ≈ 7 span the mass range M ⋆ ∼ 5.3 × 10 6 − 7.4 × 10 8 M ⊙ , yet the few galaxies with M ⋆ ≳ 10 8 M ⊙ are far from undergoing mini-quenching.We thus focus on one galaxy (labeled V5 and highlighted in Fig. 2) and its mini-quenching event at z = 6.7 and M ⋆ = 2.1 × 10 7 M ⊙ as the reference galaxy to compare to JADES-GS-z7-01-QU.
The account for differences in redshift and stellar mass between V5 and JADES-GS-z7-01-QU, we apply flux correction factors d L (7.3) 2 /d L (6.7) 2 (d L (z) being the luminosity distance at redshift z) and M ⋆,V5 /M ⋆,GS-z7-01-QU (assuming mass-to-light ratio scaling, see Schombert et al. 2019).However, incorporating these correction factors we find that the flux density of JADES-GS-z7-01-QU is about a factor of 5 stronger than V5 (see Fig. 5,right panel).Can this discrepancy be resolved by modifying metallicities and ages of the simple stellar populations entering the mock SED calculation?

Stellar Metallicity
Stellar populations in V5 are only moderately more metal-rich (Z ⋆ ∼ 0.04 Z ⊙ on average) than deduced for JADES-GS-z7-01-QU, Z ⋆ ∼ 0.01 Z ⊙ .Reducing Z ⋆ artificially for some stellar populations does not modify SEDs by more than ∼ 0.1 dex, in accordance with Gelli et al. (2023) (at most ∼ 0.2 dex for a reduction of Z ⋆ by 2 orders of magnitude from Z ⋆ ∼ Z ⊙ to Z ⋆ ∼ 0.01 Z ⊙ ), hence we do not explore this path further.

Ages of Stellar Populations
The SFH of V5 (see Fig. 5, left panel) reveals that the simulated galaxy underwent an extended burst of star formation at z ≈ 7.3, followed by a gradual suppression of its SFR.At the epoch of observation, z = 6.7, the average age of stellar populations is thus ∼ 230 Myr.While the shape of the SED (see Fig. 5, right panel) agrees moderately well with JADES-GS-z7-01-QU photometry, the overall normalization does not.
Using BAGPIPES, Looser et al. ( 2023) infer a top-hat-like SFH and the time elapsed between mini-quenching and the epoch of observation inferred by four different full spectral fitting codes is ∆t quench = 10 − 40 Myr.We thus evaluate the SED of an idealised top-hat SFH.We adopt ∆t quench = 15 Myr and vary the top-hat width between 10 Myr and 90 Myr.The corresponding top-hat heights are SFR = 2.02 M ⊙ yr −1 and SFR = 0.22 M ⊙ yr −1 , respectively, and the resulting SEDs are shown as the shaded area in Fig. 5 (right panel).The top-hat that gives rise to an SED closest to JADES-GS-z7-01-QU photometry is of width 60 Myr (cf.Fig. 5, left panel).The fluxes now reach ∼ 5 times higher in the UV and ∼ 2 times higher in the red part of the spectrum (observed frame), in good agreement with JADES-GS-z7-01-QU.Since the mock SED flux for wavelengths around the red filter F444W is too low compared to JADES-GS-z7-01-QU photometry, the agreement found for the SED shape is only modest.

Caveats
When artificially modifying SFHs and comparing to observations, it is important to mention two caveats.First, V5 is mini-quenching at z = 6.7 and M ⋆ = 2.1 × 10 7 M ⊙ as opposed to the inferred epoch of observation z = 7.3 and estimated mass M ⋆ ∼ 5 × 10 8 M ⊙ found for JADES-GS-z7-01-QU.The mass-to-light scaling that we adopt might be an invalid assumption.However, we have repeated the SED analysis for several mini-quenching events in VELA across z = 6.5− 7.5.For TNG100, we have likewise calculated SEDs of MQGs at z = 7, which is the closest redshift at which a resolved dust attenuated SED modelling can be performed6 , following Nelson et al. (2018).In all cases, we come to the same conclusion: The SEDs can only be reconciled with JADES-GS-z7-01-QU photometry when artificially modifying ages of stellar populations.
Secondly, there are different recipes for dust attenuation modelling, and while an extensive comparison is beyond the scope (cf.Nelson et al. 2018;Vogelsberger et al. 2020;Shen et al. 2020), we find that differences between e.g.adopting HI vs H column densities only lead to minute effects on the resulting mock SEDs.However, since JADES-GS-z7-01-QU photometry suggests a redder spectrum than we infer for the top-hat SFHs, we speculate that JADES-GS-z7-01-QU is more dust-obscured than V5 and/or has some older populations than the top-hat.

Discussion
We speculate that VELA and IllustrisTNG galaxies at high redshift are not bursty enough on small timescales ∼ 40 Myr to give rise to the high fluxes observed for JADES-GS-z7-01-QU.Only when allowing a maximum burst of star formation (top-hat SFH), matching ∆t quench = 15 Myr with the observed one, can the normalization of SEDs be reconciled.While the exact value observed in the red filter F444W is hard to reproduce in mocks, the overall SED shape inferred for a top-hat SFH is in good agreement with JADES-GS-z7-01-QU, including Balmer absorption lines, lack of emission lines, and UV continuum.There is a possibility that JADES-GS-z7-01-QU is in fact an obscured AGN, which would not only provide a mechanism for abrupt quenching but possibly explain the flux discrepancy in the F444W filter.The abundance of (dust-obscured) AGN at high redshift of z > 5 might be an order of magnitude higher than expected from extrapolating quasar UV luminosity functions (Matthee et al. 2023;Larson et al. 2023;Übler et al. 2023;Endsley et al. 2023).However, unlike JADES-GS-z7-01-QU many highredshift AGN appear to have a strong Balmer break (e.g.Kocevski et al. 2023).
Which sub-grid models governing galaxy formation could give rise to higher burstiness on small timescales (∼ 40 Myr) to match the observations?While for IllustrisTNG galaxies, resolution effects (e.g.M ⋆ of low-mass galaxies typically increases with resolution, see Pillepich et al. 2018a) might also be at play, IllustrisTNG has been shown to poorly predict some observables at high redshift, e.g. the scatter around the MS at the high-mass end (cf.Sec.3.1) or the abundance of dust-obscured, far-infrared galaxies and thus the obscured cosmic star formation rate density (Shen et al. 2022).
One possibility is that better modelling at high redshift resolves the discrepancies with observations found here.The fact that Gelli et al. ( 2023) also succeeds in reproducing the SED of JADES-GS-z7-01-QU when artificially modifying SFHs points in that direction.They analyse the SERRA suite of high-resolution zoom-in simulations that includes on-the-fly radiative transfer and a nonequilibrium chemical network.Even though these prescriptions are more suitably in the epoch of reionization (see Pallottini et al. 2022) than a spatially uniform UV background (Haardt & Madau 2012), the levels of burstiness they infer on small timescales ∼ 40 Myr are still too low.

CONCLUSIONS
Bursty star formation at high redshift gives rise to (likely only) temporarily quenched, or mini-quenched galaxies in the mass range M ⋆ = 10 7 − 10 9 M ⊙ .With the advent of JWST and the first observations of such galaxies (Looser et al. 2023;Strait et al. 2023), it is critical to gain a thorough understanding of the physical mechanisms and the timescales involved.Combining insights and leveraging periodic box simulations, zoom-in simulations and empirical models is an important first step in understanding the regulation of star formation in high redshift galaxies.
Methods: We employ four galaxy formation models, a periodic box simulation (IllustrisTNG), two zoom-in simulations (FirstLight and VELA) and an empirical halo model (EHM) to investigate the properties of high-redshift (mini-)quenched galaxies.We adopt an aperture of twice the stellar half mass radius 2 × R star,h (and 0.15 R vir for FirstLight), an averaging timescale of ∆t = 10 Myr (∆t = t dyn for EHM) and a selection threshold of 0.75 dex below the star-forming MS.
Abundance: We find that the abundance of quenched galaxies at high redshift inferred from IllustrisTNG, VELA and EHM is largely consistent with each other, implying that this galaxy population is rare.The quenched galaxy population first appears below z ≈ 8, after which their fraction increases with cosmic time, from ∼ 0.5 − 1.0% at z = 7 to ∼ 2 − 4% at z = 4 in the mass range M ⋆ = 10 7 − 10 9 M ⊙ .The corresponding comoving number densities read 10 −5 Mpc −3 at z = 7 and 10 −3 Mpc −3 at z = 4.The number of quenched galaxies decreases monotonically with increas-ing stellar mass M ⋆ , a dependence that is stronger in models which probe a smaller dynamic range (smaller box sizes) such as TNG50.Quenched galaxy fractions in IllustrisTNG and EHM are consistent with SFR duty cycle estimates ( f duty ∼ 98.6 − 99.9% at z = 7) inferred for FirstLight and VELA galaxies.
Duration: For MQGs, SFHs rapidly change before, during and after the mini-quenching event, consistent with the idea that miniquenching results from bursty star formation.The distribution of mini-quenching timescales (defined as the accumulated time a MQG spends ≥ 0.75 dex below the MS) averaged across z = 4 − 8 in FirstLight and VELA peaks around τ MQ ∼ 20 − 40 Myr.While the upper tail of the distribution is highly sensitive to the threshold adopted, only one simulated galaxy stays quiescent for an extended period of time (V28 for ∼ 220 Myr at redshifts z ≈ 4.9 − 5.7).This mini-quenching timescale is close to the local free-fall timescale t ff of the inner halo.
Cause: In EHM, quenching is by construction caused by a lack of gas inflow, which itself is tied to dark matter accretion rates.In IllustrisTNG, FirstLight and VELA, we in addition find that the periodic injection of energy and momentum into the circum-and intergalactic medium via stellar feedback in the context of bursty star formation contributes to the regulation of star formation and thus the phenomenon of mini-quenching.However, by following the main progenitor branch and main descendant branch of quenched galaxies in TNG100 and TNG50, we show that many quenched galaxies (∼ 27% at z = 7) are gravitationally interacting with other galaxies, and even when not fully merging are tidally disrupted.
Consistency with Observations: Simulated SEDs can only be reconciled with JADES-GS-z7-01-QU photometry in both IllustrisTNG and VELA when artificially modifying ages of simulated stellar populations.In particular, a top-hat SFH of width 60 Myr shows best agreement with JADES-GS-z7-01-QU, consistent with observationally inferred SFHs.While simulated SED shapes agree moderately well including Balmer absorption lines, the flux density in the red F444W filter is lower in the top-hat SEDs than observed by a factor of ∼ 1.5.This is likely caused by higher levels of dust obscuration for JADES-GS-z7-01-QU compared to simulated galaxies, though some older populations would be needed for even better agreement, disallowed by the top-hat.Alternatively, JADES-GS-z7-01-QU could be an obscured AGN.
Outlook: The fact that we need to artificially modify ages of stellar populations to find agreement with the observed SED lets us conclude that sub-grid models governing galaxy formation at high redshift have likely to be adapted, including in the higher-mass regime of M ⋆ = 10 8 − 10 9 M ⊙ in which star formation is expected to transition from bursty to stable.On-the-fly radiative transfer and a non-equilibrium chemical network (Gelli et al. 2023) adopted for SERRA simulations is not enough to remedy the discrepancies.In the context of bursty star formation, an improved understanding is needed of how much (stochastic) power exists on the temporal scales probed by observations.MQGs can thus be a useful probe for subgrid models, and will help close the gap between observations and theoretical models.Extending the concept of MQGs to lower redshift and studying the transition from bursty to steady star formation at z = 1 − 3 will be necessary to interpret the upcoming wealth of measurements on the low-mass quiescent population driven by deep JWST data.

Figure 1 .
Figure 1.The star-forming MS in the four galaxy formation models investigated in this work across redshifts z = 4 − 8: the combined simulation study TNG100+50 (first panel), EHM (second panel), FirstLight10+20 (third panel) and VELA (fourth panel).SFRs are averaged over 10 Myr.A compilation of observational studies by Popesso et al. (2023) is given in the form of a common parametrization (dashed curves) below z ≤ 6 for M ⋆ ≥ 10 8.5 M ⊙ .The MS ridges in the four models are consistent to within ∼ 30%, with TNG100+50 and FirstLight10+20 predictions closest to observational constraints.

Figure 2 .
Figure 2. Galaxies about the star-forming MS in a combined simulation study of TNG100+50 (first panel, shown for z = 7), EHM (second panel, shown for z = 7), FirstLight10+20 (third panel, shown for z = 7) and VELA (fourth panel, shown for z = 5 − 10).Quiescent galaxies are highlighted with blue circles and lie below the selection threshold of 0.75 dex below the MS, shown in dashed black.For TNG100+50, error estimates from bootstrapping over masses of newly formed stellar particles are added to each quiescent galaxy.Galaxies with SFR = 0 are formally assigned a MS deviation of ∆ MS = −∞.The colorbar indicates the comoving number density of galaxies.The total MS scatter averaged across M ⋆ = 10 7 − 10 9 M ⊙ is indicated in the bottom right of the panels.VELA: Two galaxies with bursty SFHs (V5 and V10) and one galaxy with a more stable SFH (V32) are shown in different colours.SFRs are averaged over 10 Myr, corresponding to the distance between two dots in the SFHs.Star symbols indicate the first star formation event (z = 10 for V5 and V32 and z = 8 for V10).The mini-quenching event of V5 (cf.Fig.5) is highlighted with a triangle symbol.The abundance of (mini-)quenched galaxies agrees across the models, with an overall fraction of f QG = 0.5 − 1.0% at z = 7.

Figure 3 .
Figure3.Quenched galaxy fraction f QG vs stellar mass M ⋆ at various redshifts (see legend) in a combined simulation study of TNG100+50 (solid), in the individual boxes TNG100 (dash-dotted), TNG50 (dotted) and EHM (dashed).The population of quenched galaxies first appears below z ≈ 8, after which their fraction increases with cosmic time, from ∼ 0.5 − 1.0% in the mini-quenching mass range M ⋆ = 10 7 − 10 9 M ⊙ at z = 7 to ∼ 2 − 4% at z = 4.The error bar on each estimate denotes the standard deviation as obtained with bootstrap resampling.The overall mass-agnostic quenched galaxy fraction is shown in Table3.There is qualitative agreement between IllustrisTNG and the simplistic EHM approach, even though the latter does not model stellar feedback.

)
Table 4 presents a comparison of duty cycle estimates for redshifts z = 4 − 8, with two different selection thresholds ∆ MS = [−0.75,−0.6] dex.Notably, VELA galaxies demonstrate an active star formation phase for ∼ 96.1 − 98.4% of their time at z = 4.The duty cycle increases with the epoch of observation z and at z = 8 attains values ∼ 99%.

Figure 5 .
Figure 5. Original and modified SFH of the VELA galaxy V5 at z > z MQ = 6.6 (left panel) and corresponding SEDs compared to JADES-GS-z7-01-QU (right panel).Shown are: (a) Original SFH and dust attenuated SED of V5 (blue), (b) Idealised top-hat SFH model of width 60 Myr (orange), with corresponding dust attenuated SED (orange), (c) JADES-GS-z7-01-QU photometry (green dots) for filters F150W, F200W, F277W, F335M and F444W.Horizontal bars on the photometry indicate the wavelength range probed by each passband.All SEDs are flux-corrected to account for redshift and stellar-mass offsets between V5 and JADES-GS-z7-01-QU.The orange shade on the right panel is obtained by varying the width of the top-hat from 10 to 90 Myr.The level of burstiness inferred from VELA on short timescales ∼ 40 Myr is lower than what is observed.

Table 1 .
Galaxy (unlike our four galaxy formation models), photoheating due to reionization can suppress SFRs by more than 50% only in low-mass haloes, specifically M h < 10 8.4 M ⊙ at z = 6.In this work, we do not investigate such reionization quenching at very low masses.Ma et al. (

Table 4 .
SFR duty cycles in VELA and FirstLight10+20: (1) epoch of observation z; (2) selection threshold in dex below the MS; (3) median and 16/84 percentile of the SFR duty cycle f duty in percent.