Metal and dust evolution in ALMA REBELS galaxies: insights for future JWST observations

ALMA observations revealed the presence of significant amounts of dust in the first Gyr of Cosmic time. However, the metal and dust buildup picture remains very uncertain due to the lack of constraints on metallicity. JWST has started to reveal the metal content of high-redshift targets, which may lead to firmer constraints on high-redshift dusty galaxies evolution. In this work, we use detailed chemical and dust evolution models to explore the evolution of galaxies within the ALMA REBELS survey, testing different metallicity scenarios that could be inferred from JWST observations. In the models, we track the buildup of stellar mass by using non-parametric SFHs for REBELS galaxies. Different scenarios for metal and dust evolution are simulated by allowing different prescriptions for gas flows and dust processes. The model outputs are compared with measured dust scaling relations, by employing metallicity-dependent calibrations for the gas mass based on the [CII]158micron line. Independently of the galaxies metal content, we found no need for extreme dust prescriptions to explain the dust masses revealed by ALMA. However, different levels of metal enrichment will lead to different dominant dust production mechanisms, with stardust production dominant over other ISM dust processes only in the metal-poor case. This points out how metallicity measurements from JWST will significantly improve our understanding of the dust buildup in high-redshift galaxies. We also show that models struggle to reproduce observables such as dust-to-gas and dust-to-stellar ratios simultaneously, possibly indicating an overestimation of the gas mass through current calibrations, especially at high metallicities.


INTRODUCTION
The understanding of the evolution of galaxies across cosmic time is one of the fundamental goals of present-day astrophysics.Over the ★ E-mail: marco.palla@inaf.itlast decades, a constantly increasing number of galaxies in the Epoch of Reionization (EoR,  ≳ 6) have been observed through their restframe UV window (e.g.Bouwens et al. 2010;McLure et al. 2013;Finkelstein et al. 2015;Oesch et al. 2016;Stefanon et al. 2019), which is now opening up to rest-frame near-IR (NIR) wavelengths thanks to the launch of the James Webb Space Telescope (JWST, e.g., Harikane et al. 2022;Naidu et al. 2022;Atek et al. 2023;Castellano et al. 2023).
However, observations of galaxies across redshifts have long demonstrated the severe impact of dust on galaxy detectability and information that can be extracted from UV-to-NIR emission (Draine 1989(Draine , 2003;;Calzetti et al. 2000;Casey et al. 2014).Only in the last several years, the Atacama Large Millimeter/submillimeter Array (ALMA) has enabled the study of the dust content in distant galaxies at rest-frame far-IR (FIR) wavelengths (see Hodge & da Cunha 2020 for a review).These observations have allowed, despite the uncertainties, fairly robust estimates of the dust mass budget in galaxies at  ≳ 6 (e.g.Riechers et al. 2013;Venemans et al. 2018;Tamura et al. 2019;Schouws et al. 2022).In particular, they reveal dust-to-stellar mass ratios ranging between ≃ 0.01−3% for "normal" (i.e. with conditions representative of the galaxy population at these early epochs) star forming galaxies (e.g. Watson et al. 2015;Laporte et al. 2017;Bakx et al. 2020;Fudamoto et al. 2021).At the same time, a range of theoretical models were tested to understand the role of different processes in the buildup of the observed dust masses, i.e. dust formation from supernovae (SNe) and asymptotic giant branch (AGB) stars, grain growth in the interstellar medium (ISM), and dust destruction from SN shocks.Either analytic (e.g.Valiante et al. 2009Valiante et al. , 2017;;Dayal et al. 2010;de Bennassuti et al. 2014;Gioannini et al. 2017), semi-analytic (e.g.Mancini et al. 2015;Popping et al. 2017;Ginolfi et al. 2018;Triani et al. 2020) or hydrodynamical (e.g.McKinnon et al. 2018;Graziani et al. 2020) schemes have been used to study the impact of these different processes on the dust content of high-redshift galaxies.However, these studies present quite different outcomes depending on their physical treatment of the ISM and dust physics (see, e.g., Popping & Péroux 2022).Furthermore, the low number of sources with dust detections at  ≳ 6 to compare with has strongly limited our understanding of dust evolution in the EoR.
The ALMA Reionization Era Bright Emission Line Survey (REBELS, Bouwens et al. 2022) program has significantly increased the number of star-forming galaxies with robust dust detections during the EoR.By searching for bright ISM cooling lines ([C ii]158m and [O ii]88m) in 40 UV-bright ( UV < −21.5) galaxies at  ∼ 6.5−9.5, this survey has yielded a sample of 16 galaxies detected in the dust continuum, among which 15 galaxies are also detected in the [C ii]158m emission line (from now on referred to as [C ii] line (Inami et al. 2022, Schouws et al., in prep).The fairly large number of sources detected in REBELS, together with the diversity probed in stellar masses and star formation rates (e.g.Topping et al. 2022), represents a new benchmark for the study of galaxies with dust mass measurements at high-redshift (Ferrara et al. 2022;Sommovigo et al. 2022a, see also Sommovigo et al. 2022b).Dayal et al. (2022) provided the first comparison between this set of sources and galaxy evolution models.By using predictions from the DELPHI semi-analytical code, they showed that their model struggles to reproduce the observed dust masses in Sommovigo et al. (2022a), especially when dealing with galaxies in the low-mass end of the sample (log( * /M ⊙ ) ∼ 9,  * from Stefanon et al., in prep.).In particular, for such galaxies, dust scaling relations were only reproduced with their "maximal dust mass" model, which relies on extreme assumptions concerning the dust evolution (e.g., grain growth timescale below 1 Myr, no dust destruction by SN shocks).Different conclusions are reached by Di Cesare et al. (2023), who investigated the same galaxies using dustygadget simulations (Graziani et al. 2020).In this case, predictions instead show an overestimation of the dust-to-stellar mass ratios relative to the data for galaxies at the highmass end (log( * /M ⊙ ) > 9.5) of the sample.In turn, this suggests either an overproduction of dust in the simulated galaxies or an overestimation of the dust temperature (and hence an underestimation of the dust mass) in the data.However, the discrepancy in the results of Dayal et al. (2022) and Di Cesare et al. (2023) can (at least) partly be attributed to the use of different sets of stellar and dust masses by Di Cesare et al. (2023) (from Topping et al. 2022 andSommovigo et al. 2022b, respectively).
Despite the progress in the characterization and understanding of dusty,  ≳ 6 sources, there is still a large amount of information that is lacking for these objects.Until now, our comprehension of the dust evolution in the EoR has been strictly limited to the analysis of the dust mass and dust-to-stellar mass ratio as a function of the stellar mass in different galaxies (e.g.Dayal et al. 2022, see also Algera et al. 2023).Fundamental diagnostic quantities for dust evolution, such as the dust-to-gas and dust-to-metal ratios, have been inaccessible and this prevented us from strictly pinning down the role of the different dust formation and destruction mechanisms in the high-redshift Universe (see Schneider & Maiolino 2023 for a recent review).However, the advent of JWST has opened up new and exciting possibilities in the investigation of these galaxies.Aside from the improvement in the photometric sampling towards the rest-frame optical and NIR, which can provide a much better characterization of the galactic SED and improve the stellar mass and SFR determination, the JWST/NIRSpec spectrograph (Jakobsen et al. 2022) allows intermediate resolution spectroscopy (R ∼ 1000 − 2700) in a wavelength regime comprising several of the most prominent optical nebular emission lines (e.g., [OII]3726, 3729, [OIII]4960, 5008, H, H).In turn, this allows robust log(O/H) + 12 metallicity determinations through strong-line diagnostics (see Kewley & Ellison 2008;Kewley et al. 2019;Maiolino & Mannucci 2019) and possibly also through the direct-T e method, in the case of detection of the auroral [OIII]4363 line.Despite of the significant intrinsic uncertainty in these pioneering high-redshift metallicity measurements (especially in the case of strong-line diagnostics), these are fundamental for many reasons.First, they allow to explore the upper mass end of the galactic scaling relations, such as the mass-metallicity relation (MZR), which is extremely poorly sampled above 10 9−9.5 M ⊙ at  ≳ 6 (Langeroodi et al. 2022;Nakajima et al. 2023;Curti et al. 2023).Moreover, ISM abundances from JWST spectroscopic observations allow us to use dust diagnostics other than the dust-to-stellar mass ratio (hereafter, DtS) to study the dust evolution in high-redshift sources.In fact, metallicity-dependent calibrations of the galactic gas mass from the [C ii] line luminosity (Heintz et al. 2021(Heintz et al. , 2022;;Vizgan et al. 2022) open up the possibility to estimate the ISM gas and metal mass in galaxies and therefore probe other fundamental dust relations, such as the dust-to-gas ratio (DtG) and dust-to-metal ratio (DtM) as function of the galactic metallicity.These dust-related diagrams have been crucial in assessing the role of mechanisms contributing to dust evolution in different galaxy types from the local universe up to redshift 4 ( Rémy-Ruyer et al. 2014;Schneider et al. 2016;De Vis et al. 2017;Ginolfi et al. 2018;Péroux & Howk 2020;Nanni et al. 2020;De Looze et al. 2020;Galliano et al. 2021, see also Galliano et al. 2018 for a review).Additionally, looking at the DtG and DtM in combination with the gas-independent DtS measurement may be crucial to directly test the goodness of the above-mentioned [C ii]-gas mass calibrations.In fact, the latter will play a crucial role in the study of the high-redshift Universe since they are usually the only way to estimate the gas mass in objects at such redshifts (due to the difficulties in getting other emission lines and/or the Rayleigh-Jeans end of the FIR SED; see, e.g.Heintz et al. 2021 and references therein).
Therefore, the aim of this paper is to provide an overview of the possible scenarios of metal and dust evolution that may arise from the analysis of REBELS galaxies with [C ii] and dust detections in the light of future JWST data.Starting from galactic chemical evolution models including the most recent stellar yields and dust prescriptions (SN and AGB dust production, dust grain growth, dust destruction) from the literature, we implement the galactic star formation histories (SFHs) as determined by Topping et al. (2022) for each specific REBELS source and test different setups of the gas flow history and dust evolution parameters to model various scenarios for the metal and dust evolution, respectively.In this way, we provide insights into the dominant dust production mechanisms and the efficiency of dust production of  ≳ 6 galaxies, currently under the assumption of different chemical enrichment scenarios, which will be constrained in the future by means of rest-frame optical ISM line spectroscopy.
The results of this paper therefore stress the urgent need for JWST spectroscopic follow-up of galaxies observed in the FIR with ALMA.JWST has started to acquire low-resolution (R ∼ 100) data for several REBELS sources (GO program ID: 1626, P.I.Stefanon), but we need larger samples and higher spectral resolution to robustly determine chemical abundances.Only by building statistically significant samples of JWST-observed galaxies, we can constrain the general metallic content of high-z dusty galaxies and therefore give tighter constraints on the dust regulation mechanisms during the EoR.
The paper is organised as follows: in Section 2 we briefly present the data products from the ALMA REBELS program used throughout this paper, focusing on the [C ii] luminosity-gas mass calibration adopted.In Section 3, we illustrate the model framework adopted for our predictions, specifying the details about the different setups used to produce different metal and dust scenarios.In Section 4, we show what the shortcomings of current data constraints are in deciphering dust evolution at high-redshift, while in Section 5 we show our model predictions for different metal and dust enrichment scenarios that may arise in REBELS targets.Finally, in Section 6 we draw our conclusions.

DATA
In this paper, we take advantage of a subsample of 13 galaxies for which both the [C ii]158m line and the 1900 GHz continuum were detected (at >5 and >3, respectively).For a full analysis of the [C ii]-detected sources, we refer the reader to Schouws et al. (in prep.), while the dust continuum detections are presented in Inami et al. (2022).
The stellar masses ( * ) adopted for these REBELS targets are presented in Topping et al. (2022).In particular, we make use of the estimations given by the SED-fitting code Prospector (Johnson et al. 2021) under the assumption of a non-parametric SFH.Such SFHs are particularly well-suited to model old stellar populations that may be present, avoiding the "outshining problem" caused by more recent bursts of star formation (see Topping et al. 2022 and references therein).Our selected subsample spans a stellar mass range of log(M * /M ⊙ ) ≃ 9.0 − 10.3, i.e. similar to the one encompassed in the full sample.
Among the different dust estimations available (e.g.Ferrara et al. 2022;Sommovigo et al. 2022a), we use the dust estimates presented in Sommovigo et al. (2022b), which are based on the methodology shown in Sommovigo et al. (2021Sommovigo et al. ( , 2022a) ) 1 .This method allows us to get fairly robust estimates of the dust temperature and mass by means of only one FIR continuum data point with Bayesian methods, by assuming the well established SFR-[C ii] (De Looze et al. 2014, see Schaerer et al. 2020;Ferrara et al. 2022 for high-z probes) and Schmidt-Kennicutt relations Kennicutt 1998) and imposing a metallicity-dependent DtG relation (e.g.Rémy-Ruyer et al. 2014).We refer to Sommovigo et al. (2022a) for a detailed explanation on the employed methodology.It is worth mentioning that for two REBELS galaxies at the high-mass end of the sample (REBELS-25 and REBELS-38) additional FIR continuum measurements around ∼ 90 m rest-frame are available (Algera et al. 2023), allowing for more robust constraints on the dust SED and thus on dust temperature and mass.We comment on the effects of these new measurements throughout the paper.
To make predictions on dust scaling relations including DtG or DtM we need to have estimates of the amount of gas available in the analyzed sources.However, measurements through usual indicators (e.g.CO lines, FIR continuum tail) are extremely hard to get at  ≳ 5, preventing analyses on high-z galaxy samples (see, e.g., Heintz et al. 2022).Therefore, the most promising approach to obtain a census of the gas content in EoR galaxies is to use a suitable tracer of the cold neutral gas phase, such as the [C ii] line.In this paper, we adopt the calibration by Heintz et al. (2021) to estimate the HI gas content relying on the [C ii] line luminosity and the metallicity: where / ⊙ is the relative solar metallicity with 12 + log(O/H) = 8.69 for log(/ ⊙ ) = 0 (Asplund et al. 2009) and  HI and  [C II] are in units of M ⊙ and L ⊙ , respectively.It is worth reminding that the metallicity  is a free parameter of Eq. ( 1), as metallicity measurements of any of the REBELS galaxies has not been published yet.For H 2 instead, we make use of the average [C ii]-to-H 2 ratio found in the EoR galaxy simulations by Vizgan et al. (2022): Other works in the literature indicate larger [C ii]-to-H 2 conversion factors (e.g.Zanella et al. 2018, see also Aravena et al. 2023).However, such calibrations are based on  ≲ 2 galaxies.Moreover, the difference between the conversion factors is around a factor of 1.5 which we found to have no implications for the studied trends presented in this work.It is worth noting that the adopted calibrations were already tested in the framework of REBELS galaxies by Heintz et al. (2022) by adopting the fundamental  * -SFR-metallicity relation by Curti et al. (2020).However, such a relation was established at  < 3.In this work, we also explore the effects of metallicity on the Heintz et al. (2021) calibration, testing different metal enrichment scenarios that could emerge from JWST observations of REBELS targets.
A summary of the derived properties for the REBELS galaxies subject of this study is shown in Tab. 1.For the gas masses and the gas fractions, we use the sum of the HI and H 2 masses inferred from Table 1.Summary of the derived properties of the REBELS galaxies within our sample. [C II] are from Schouws et al. (in prep). * are from non-parametric SED fitting presented in Topping et al. (2022). dust are from Sommovigo et al. (2022b)

MODELLING DUST EVOLUTION IN REBELS GALAXIES
To study the origin of trends in observed dust and metal scaling relations inferred from future JWST observations, we use detailed chemical evolution models for galaxies including interstellar dust evolution.To this aim, we adopt a significantly revised version of the Chemevol2 source code (De Vis et al. 2021), in which we implemented updated recipes for elemental nucleosynthesis in stars and dust processes in the ISM.
In the next subsections, we go over the details of the model, explaining the star formation histories (SFHs) implemented for the different targets and the adopted prescriptions to predict the metal and dust scaling relations in REBELS galaxies.

Chemical and dust evolution model
In our galactic evolution models we assume that galaxies start forming at redshift  = 25, by primordial gas accretion from the intergalactic medium (IGM) which progressively accumulates onto the system, and the galaxy evolves suffering galactic winds.As stated above, we assume gas accreted from the IGM to be devoid of both metals and dust.However, we also made tests in which reaccretion of the gas ejected through galactic winds also fuels subsequent gas infall onto galaxies, finding no significant differences in the model results.The models also relax the instantaneous recycling approximation (IRA), i.e. they take into account stellar lifetimes in the process of chemical and dust enrichment.Therefore, realistic predictions are possible for different chemical abundance ratios and dust enrichment.

Chemical evolution prescriptions
The evolution of a chemical element i in the ISM can be written as (see also Matteucci 2012): where () is the star formation rate (SFR) and  i () is the fraction of the element i in the ISM at the galactic evolutionary time , which in our framework starts from 0 at redshift  = 25.The first term on the right-hand side of Eq. (3) corresponds to the rate at which an element i is removed from the ISM due to the star formation process.
i () (see Palla et al. 2020b for the complete expression) takes into account the nucleosynthesis from low-intermediate mass stars (LIMS,  < 8M ⊙ ), core collapse (CC) SNe ( > 8M ⊙ ) and Type Ia SNe.The IMF assumed is the one from Chabrier (2003), in order to be consistent with the stellar masses and SFHs adopted from Topping et al. (2022).The stellar yields are taken from the most recent prescriptions of Ventura et al. (2013Ventura et al. ( , 2018Ventura et al. ( , 2020Ventura et al. ( , 2021) ) for LIMS and Limongi & Chieffi (2018) (set R300) for CC-SNe.For Type Ia SNe instead, we adopt the yields by Iwamoto et al. (1999), as extensively tested in chemical evolution works (see Palla 2021 and references therein).For the latter, we assume a  −1 delay time distribution (Totani et al. 2008;Maoz & Graur 2017), which is extensively used in models of chemical evolution (e.g.Yates et al. 2013;Côté et al. 2016;Vincenzo et al. 2017).
The third term in Eq. ( 3) is the gas infall rate, which in the framework of chemical evolution models is generally parametrised with an exponential decaying law, even for high-redshift galaxies (e.g.Romano et al. 2017;Palla et al. 2020b;Kobayashi & Ferrara 2023): where  inf is the normalization constant constrained to reproduce the total gas mass accreting onto the system  gas,tot ,  i,inf is the fraction of the element i in the infalling gas (with primordial composition) and  inf is the typical infall timescale, which is chosen according the shape of the galactic SFH adopted (see 3.2).However, during this work we also test multiple gas accretion episodes in case of SFHs showing strong peaks in SFR (see 3.2).The last term in Eq. 3 accounts for galactic winds.Here, we assume that the outflow is proportional to the SFR, as in this form: where  i is the wind mass loading factor for the element i.Since our main focus is to follow the general scaling relations in the modeled galaxies, we do not use differential winds, i.e. we assume the same loading factor  ≃ 3 for the different chemical elements in the gas.
The choice of such a loading factor allows to explain the extended [C ii] halo observed in several redshift  ≃ 6 galaxies in terms of a star-formation driven outflow, with similar [C ii] halo features relative to REBELS target stacks (see Pizzati et al. 2020;Fudamoto et al. 2022).
Despite of careful choice of model prescriptions, motivated either by consistency with the derivation of the physical quantities of REBELS galaxies or by other observational proofs, we want to highlight that a certain level of uncertainty persists in the modelling assumptions in high-redshift galaxies.For example, the treatment of stellar generations in the early stages of galaxy formation, with its assumptions on the IMF and the choice of stellar yields, have a non-negligible impact on the modelling results, leading to different conclusions on the evolution of the baryonic components in highredshift galaxies (e.g.Naab & Ostriker 2017 and references therein).
However, our choice of prescriptions remains robust accross our main goal, which is giving a first exploration of the parameter space in the explanation of dust-scaling relations, rather than drawing definitive conclusions on the evolution of EoR, UV-bright galaxies as those of the REBELS sample.In fact, the level of uncertainty in the treatment of dust evolution is even higher than the one for the rest of the baryonic physics, with severe limitation persisting even in the local Universe (e.g.De Looze et al. 2020;Galliano et al. 2021;Calura et al. 2023).

Dust evolution prescriptions
Here we apply the same formalism used in previous works which included the chemical evolution of dust in galaxies of various morphological types and at different redshifts (e.g., Calura et al. 2008;Gioannini et al. 2017;Gjergo et al. 2020;De Looze et al. 2020;Galliano et al. 2021).The equation used to determine the amount of dust at time  is: Similar to Eq. ( 3),  dust () represents the dust fraction in the ISM.
The first term on the right-hand side of Eq. ( 6) accounts for the amount of dust destroyed by astration, i.e. incorporated into new stars.
dust () accounts for stellar dust production, which in our model is from AGB stars (i.e. the final stage of a LIMS before the final white dwarf fate) and CC-SNe3 .For AGB stars, we use the dust yields of Dell 'Agli et al. (2017); Ventura et al. (2018Ventura et al. ( , 2020Ventura et al. ( , 2021)), while for CC-SNe we use the dust yields from Marassi et al. (2019) (set ROT FE).It is worth noting that both yield prescriptions are metallicity dependent.As detailed in 3.2, we allow the SN dust yields to be reduced by a factor  to simulate the uncertain effect of the SN reverse shock, which can destroy much of the dust produced in the same SN explosion (e.g.Bianchi & Schneider 2007;Bocchio et al. 2016;Kirchschlager et al. 2019Kirchschlager et al. , 2023;;Slavin et al. 2020).In fact, different models of this phenomenon lead to very different survival rates, from more than 50% to much less than 10% (see Micelotta et al. 2019;Schneider & Maiolino 2023 for a review).
The third term on the right side of Eq. ( 6) accounts for dust accretion (or growth) in the cold ISM component.In this process, pre-existing dust seeds grow in size due to accretion of refractory elements on their surface (e.g.Savage & Sembach 1996;Hirashita & Kuo 2011;Konstantopoulou et al. 2022), enabling a significant increase of the global dust mass (e.g.Dwek 1998;Asano et al. 2013;Mancini et al. 2015;De Vis et al. 2017;Di Cesare et al. 2023).As in Eq. ( 6) the dust growth rate is regulated by the accretion timescale  grow , which can be expressed as: where DtM() is the dust-to-metal ratio at a time , DtM max is the maximum dust-to-metal ratio (0.6, see Appendix A) and  0 is the characteristic dust growth timescale.We look at this term in detail in 3.2.The fourth term of Eq. ( 6) accounts for dust destruction due to exploding SNe (both CC and Type Ia), whose shocks efficiently destroy interstellar dust grains and return the heavy elements locked up in dust grains back to the gas phase (see, e.g.McKee 1989;Jones et al. 1994).Dust destruction is also expressed in terms of a grain destruction timescale  destr , which can be written as: where  rate is the sum of the CC and Type Ia SN rates, and  clear is the ISM mass completely cleared out of dust per SN.We discuss this last term in 3.2.Finally, the last term of the Eq. ( 6) accounts for the amount of dust lost due to galactic winds.As mentioned in 3.1.1,we do not assume differential winds.Therefore, we assume that gas and dust are well mixed in galactic outflows.

Model setup
As mentioned in Section 2, here we adopt the stellar masses from non-parametric SFHs presented in Topping et al. (2022).Therefore, the input SFHs for our models are taken from this set of data.
In Fig. 1 we show the SFHs (blue lines) from Topping et al. ( 2022) for REBELS-14, 29 and 25, i.e. three galaxies with rather different SFHs and stellar masses (from  * ≃ 10 9 M ⊙ to  * ≃ 2 × 10 10 M ⊙ ) within our sample.In this Figure, we also show the SFHs used as input for our chemical evolution models for each galaxy (orange lines).In fact, to avoid unrealistic spikes in the chemical enrichment, we smooth the Topping et al. (2022) SFH by adopting non-parametric regression from the PyQt-Fit package4 .The resulting smoothed SFH is then multiplied by a renormalization factor to keep the total stellar mass the same as the one from the original non-parametric SFH.To let the chemical evolution work, we also need ISM gas in our models to fuel star formation.Since variations in the gas mass are expected by varying the galactic metallicity (see Eq. ( 1)), we test multiple setups for the infall mass, i.e. the total amount of gas accreting the galactic system (see, e.g.Palla et al. 2020b).By having the galactic SFH fixed for each galaxy, the infall mass represents the main parameter with which we can vary the efficiency of star formation, i.e. the amount of gas converted into stars per unit time, which in turn has a direct effect on galactic metallicity, with lower efficiencies leading to lower levels of metal enrichment (e.g.Matteucci 2012, 2021 and references therein).As a further proof of this, despite the importance of galactic outflows, it was demonstrated by several works (e.g.Yin et al. 2011;Gioannini et al. 2017;De Vis et al. 2021) that galaxies with large gas mass fractions require models with low efficiencies of star formation to match their different observables.
For what concerns the gas infall parameters (see Eq. ( 4)), we set the typical infall timescale  inf as the evolutionary time at which half of the total mass of the galaxy has been assembled, to be in line with observed scaling relations between SFR and gas mass (see, e.g.Kennicutt 1998).Moreover, we also run additional models in which we allow for multiple gas accretion episodes in case the SFH shows distinct peaks in the SFR, i.e. when the SFR varies by at least a factor of 2 in a time interval less than Δ < 0.1 Gyr.When these conditions are met, an additional infall episode modeled as in Eq. ( 4) starts, with the fraction of gas accreted by this episode being equal to the fraction of stellar mass formed in the time between the previous and the subsequent burst of SF.This choice is made in order to mimic a burst of star formation triggered by a large amount of gas that becomes available after, e.g. a galaxy merger.It should be noted that we tested the viability of the above model assumption by varying the specific conditions assumed for the SFR change factor and the time interval Δ.
In Fig. 2, we show the ISM metallicity  evolution for the three REBELS galaxies with SFHs shown in Fig. 1 assuming different gas infall masses, namely 3, 10 and 15 times the stellar mass  * computed for that REBELS target (see Appendix B for a summary table on the prescriptions adopted).We stress that such an exploration of the parameter space is necessary due to the lack of metallicity measurements in REBELS galaxies (see also Section 1).In Fig. 2, the different model tracks and the colored areas in between want to highlight the large discrepancy in metallicity evolution by having different efficiencies of star formation, which in turn are caused by the different gas accretion histories assumed in the models.We indeed clearly see that by reducing the gas infall mass we monotonically increase the metallicity at a given evolutionary time.However, the amount of gas mass accreted to the galaxy has a minor effect on the shape of metallicity evolution in each of the panels in Fig. 2.This happens because of the fixed SFHs for the model galaxies (see Fig. 1).This is especially true for the enrichment of oxygen (adopted as global metallicity indicator due to its mass dominance over other metals), which closely traces the star formation rate in galaxies (see e.g.Maiolino & Mannucci 2019 and references therein), at variance with other elements for which the production is either largely dependent on the initial metallicity of a stellar population (e.g.nitrogen, Romano et al. 2010Romano et al. , 2017) ) or for which the chemical enrichment by a stellar population can be significantly delayed with respect to the time of its birth (e.g.iron, Matteucci 2012 and references therein).
We also note that the models with the largest amount of gas infall  gas,tot = 15  * ) show an average level of chemical enrichment around 1/10 solar ( ⊙ = 0.014, Asplund et al. 2009), which is indicative of a metal-poor scenario.On the other hand, models with the lowest amount of gas accreted onto the galaxies ( gas,tot = 3  * ) have metallicities around the solar value (see Fig. 2, central and right panels).Therefore, we can label such models as indicative of a metal-rich scenario.We also check that these model setups lead to "realistic" efficiencies of star formation in the framework of chemical evolution models.We find values in between 0.1 Gyr −1 and 10 Gyr −1 , which are typical in chemical evolution models from local up to high-redshift galaxies (e.g.Gioannini et al. 2017;Gjergo et al. 2020;Palla et al. 2020aPalla et al. ,b, 2022;;De Vis et al. 2021;Galliano et al. 2021).Therefore, we will consider these two scenarios for our model galaxies to test the impact of metallicity variations on galactic dust evolution.
We also show the metallicity evolution for models where multiple gas accretion events are included in case of a star formation burst added to the SFHs from Topping et al. (2022) (see Fig. 2, thin lines).The difference with the standard models with a single gas accretion episode is small, and therefore the adoption of different setups has negligible impact also on the evolution of the dust quantities.For this reason, in the rest of the paper we look at the results for the models with single gas infall episodes only.

Dust evolution scenarios
To probe different scenarios of dust evolution in EoR galaxies, we test various prescriptions for stardust production, dust growth, and SN dust destruction.
As described in 3.1.2,we adopt the dust yields from Dell'Agli et al. ( 2017); Ventura et al. (2018Ventura et al. ( , 2020Ventura et al. ( , 2021) ) for AGB stars, while for CC-SNe we use the dust yields from Marassi et al. (2019).For these latter, we adopt the set ROT FE, which uses the same chemical network and explosion setup as the adopted SN yields from Limongi & Chieffi (2018) (set R300).Therefore, our dust yields match perfectly with the adopted chemical yields.This is at variance with a large number of previous works on dust evolution, assuming dust and chemical yields from different stellar models or using scaling factors (e.g.Calura et al. 2008;Zhukovska 2014;Gioannini et al. 2017;De Vis et al. 2017, 2021, however see also e.g.Valiante et al. 2014;Mancini et al. 2015;Ginolfi et al. 2018;Graziani et al. 2020 for a consistent use of chemical and dust yields).However, SN dust yields are affected by severe uncertainties due to the largely unknown effect of the reverse shock (Gall & Hjorth 2018;Micelotta et al. 2019), which is not included in the SN dust formation calculations (Marassi et al. 2019, see however Bianchi & Schneider 2007 for previous yield grids).For this reason, we allow the reduction of the SN dust yields by a factor  = 1, 3, 10 to mimic the effect of the SN reverse shock and test its effect in the dust scaling relations for these galaxies.
For what concerns the dust accretion timescale, this is dependent on a characteristic timescale  0 .Among the several formulations in the literature, we use the one from Mattsson et al. (2012): where  is a free parameter, whose expected value is very uncertain, with results from local galaxies ranging from values of a few hundreds to greater than three thousand (e.g.Finally, the dust destruction rate in galaxies depends critically on the  clear parameter, which indicates the ISM mass completely cleared out of dust in the ISM.Throughout this paper, we assume the formulation of Priestley et al. (2022): where  clear and  0 are expressed in M ⊙ .Here, we allow  0 to vary to test different dust destruction efficiencies.In particular, we set  0 = 2 × 10 3 , 8 × 10 3 , 1. Thanks to these parametrisations, in this paper we concentrate on three different scenario for dust evolution in REBELS galaxies, which we summarise below (see also Appendix B): • dust-poor scenario: in this scenario, we reduce the SN dust yields by a factor  = 10.The dust growth parameter  is set to 200 (lowest value), while  0 for dust destruction to 1.5 × 10 4 M ⊙ (highest value); • intermediate-dust scenario: here, we reduce the SN dust yields by a factor  = 3.The dust growth parameter  and the dust destruction scaling  0 are set to the intermediate values among those tested, i.e. 1000 and 8 × 10 3 M ⊙ , respectively; • dust-rich scenario: in this case, the SN dust yields are the same as the ones presented in Marassi et al. (2019).The dust growth parameter  is set to 3000 (highest value), while  0 for dust destruction to 2 × 10 3 M ⊙ (lowest value).However, we also test "mixed" dust scenarios, i.e. models where dust process recipes favouring or disfavouring the dust mass buildup are mixed together (see Appendix D2).In turn, this allows us to probe in more detail the relative influence that the different dust processes have on shaping the dust scaling relations.

CURRENT CONSTRAINTS: DUST EVOLUTION WITHOUT METALLICITY INFORMATION
In this Section, we show what the current data products can reveal about the dust evolution in the REBELS galaxy sample.Until now, our understanding of the dust evolution in the EoR has been limited to few observational constraints, i.e. the galactic dust mass and stellar mass (e.g., Watson et al. 2015;Hashimoto et al. 2019;Reuter et al. 2020;Algera et al. 2023).In the following, we present how these constraints can be interpreted in light of the chemical and dust evolution models presented in Section 3.
In Fig. 3, we show a common diagnostic when looking at highredshift galaxies, i.e. the DtS evolution as a function of the stellar mass  * .In particular, we compare the results of the galaxy models shown in Fig. 2 (REBELS-14, 29, 25) with the REBELS data binned according to their stellar mass (REBELS-8, 14, 19, 39 for the low stellar mass bin, 12,18,25,27,29,32,38,40 for the high-mass one).The choice of comparing the three individual galaxy models with binned data is justified by the fact that REBELS-14 show a SFH that is quite typical of the low stellar mass bin, while REBELS-29 and 25 have two SFHs that are recurrent in the high stellar mass bin.In addition, this choice allows us to better show the typical evolution of an individual galaxy for different dust evolutionary setups.As mentioned in 3.2, in Fig. 3 we only display the results for models with the largest and smallest amounts of gas mass accretion, i.e.  gas,tot = 15  * , indicative of a metal-poor chemical evolution scenario (Fig. 3 upper panels), and  gas,tot = 3  * , indicative of a metal-rich scenario (Fig. 3

lower panels).
By looking at the Figure, it can be noted that the uncertainty related to the dust and stellar mass in REBELS galaxies (light orange and dark red errorbars indicate typical uncertainties in the low and high mass bins) prevents to draw definitive conclusions on the efficiency of dust buildup in our target sample.Only the dust models assuming relatively slow dust buildup (the dust-poor scenario in 3.2.1) in Fig. 3 left panels can be ruled out due to their very low DtS predictions at all stellar masses, in agreement with conclusions drawn in previous work (e.g.Dayal et al. 2022).The other setups (Fig. 3 central and right panels) show in general values that are compatible with current observations: galaxies in the lowest  * bin show a preference for very efficient dust buildup, while REBELS galaxies with large stellar masses seem to support a slower dust buildup scenario.In fact, as can be noted in Tab. 1, the dust masses show a rather flat trend compared to stellar mass, which instead can vary up to an order of magnitude: therefore, a decreasing DtS ratio with stellar mass is observed for REBELS galaxies.This is in reasonable agreement with the trend observed in local galaxies (from Casasola et al. 2020, small grey crosses), which shows, however, a shallower relation.
The DtS vs.  * diagnostic diagram (as well as the equivalent  dust vs.  * relation) also definitely struggles in discriminating the level of metal/chemical enrichment in galaxies.Despite a difference in metallicity of almost an order of magnitude between the different scenarios of metal enrichment (see Fig. 2), Fig. 3 upper and lower central panels show a rather small difference (less than a factor of 2) in the predicted dust masses.This happens because of the physical quantities affecting the SN dust destruction rate, which in such models is only secondary to stardust production.The destruction timescale is both proportional to the gas mass and The color code for the galaxies with different star formation history is the same as in Fig. 2. The thin, grey dashed line indicates a ratio equal to 1, corresponding to the equilibrium between the two dust processes.
the metallicity  destr ∝  gas (1 + ).Since models with larger gas masses have lower metallicity (and vice versa) the two effects tend to cancel out, with the gas mass being only slightly more influential than metallicity.The difference in the model predictions is even smaller between Fig. 3 upper and lower right panels, where the DtS overlap for models with different metallicities.The almost identical predicted dust masses for metal-rich and metal-poor model galaxies is not only a problem from the point of view of galaxy metal enrichment at high-redshift, but also for the understanding of the relative contributions between "dust-cycle" processes in the buildup of the dust mass.To show an example, in Fig. 4 we plot the ratios between stardust production and grain growth timescales 5 (both dominating over destruction in this case) for the dust-rich models presented in Fig. 3 right panels.As expected, metal-poor models (solid lines) exhibit a dominance of stardust production over dust growth at all galactic evolutionary times, since the lack of metals in the ISM halts the metal accretion on dust seeds (e.g.Mattsson et al. 2012;Asano et al. 2013).On the other hand, metal-rich models (dashed lines) show that the dust growth rate can be comparable or larger than stardust production rate (see, e.g.Mancini et al. 2015;Graziani et al. 2020;Di Cesare et al. 2023).The trends also resemble the SFHs shown in Fig. 1, and therefore highlight the role of the SFR evolution together with the metallic content in shaping this ratio.
In conclusion, the current constraints we have at our disposal for high-redshift galaxies are not sufficient to give a compelling view on galactic dust evolution in the early Universe.Information on the galactic metallicity is fundamental to start overcoming such limitations.Metallicity itself allows to better constrain the metal enrichment and gas infall history, while comparison of the metallicity and dust constraints tells us how efficiently metals are being locked up in dust grains and can give an indication of what grain formation processes are needed to reproduce the observations.In this way, we 5 for dust timescale we mean the ratio between the total dust mass and the respective rate (of stardust production, dust growth or dust destruction).The lower the dust timescale, the larger is the importance of a dust process in regulating the dust budget.
can constrain the metal and dust evolutionary scenarios in more detail similar to studies in local up to intermediate-redshift ( < 4) galaxies (e.g.Rémy-Ruyer et al. 2015;De Looze et al. 2020).

DUST VS. METALLICITY: PREDICTIONS BY MODELS
In this Section, we investigate several metal and dust enrichment scenarios that will become possible to constrain with future spectroscopic analyses of JWST data of REBELS and other high-z galaxies.
In particular, we focus on the predictions of various dust scaling relations under the assumption of different dust and metallicity scenarios tested in this work.In doing this, we also provide a first important test for recently inferred, metallicity dependent  [C II] −  gas calibrations for high-redshift galaxies.All these findings are then framed in the context of the new mass-metallicity relations at  > 4 obtained by means of JWST large programs in order to provide the first suggestions on the most likely evolutionary pictures.

Dust scaling relations
In the following, we look at the three most common dust diagnostic diagrams used for the analysis of dust evolution at lower redshift, i.e., the dust-to-gas (DtG), dust-to-stellar (DtS) and dust-to-metal (DtM) mass ratios as a function of the metallicity (in terms of the oxygen abundance 12 + log(O/H)).The results of these predictions are compared with the dust, stellar, and gas masses presented in Section 2 (see Tab. 1) at different metallicities.
As done in Fig. 3, we bin the data according to the stellar mass of the galaxies as computed by Topping et al. (2022).It is worth mentioning that this choice of data binning leads to a smaller scatter than the average error for the studied quantities (see Fig. 3).In addition, here we group the galactic evolution models for individual REBELS sources accordingly, to ensure direct comparison between models and data of similar stellar mass and SFH.

Metal-poor scenario
We start by looking at the results from the models with the largest amount of gas accreting onto the galaxies, i.e.  gas,tot = 15  * , indicative of a metal-poor scenario.
In Fig. 5 we show the predictions of the chemical and dust evolution models in the case dust evolution parameters have intermediate values (intermediate-dust scenario, top panels) or dust production is particularly favoured (dust-rich scenario, bottom panels).The results of the models are compared with the DtG, DtS and DtM ratios computed from the values in Tab. 1 by assuming a metallicity typical of the output of the models, i.e. 1/10 Z ⊙ , and a metallicity of 1/3 Z ⊙ (both indicated by filled squares).Due to the uncertainties in the gas mass calibrations in EoR galaxies, we show two additional data points for each metallicity in the left and right panels of Fig 5 .Crosses should be seen as a sort of lower limit for the gas mass by considering only the [C ii]-H 2 calibration of Vizgan et al. (2022) for the total gas mass.The plus symbols instead indicate an intermediate scenario, where UV/[C ii] weighted gas masses, i.e. gas masses scaled to the UV-to-[C ii] effective radii ratio in REBELS (see Appendix C), are considered.
We note that the models that belong to a dust-poor setup are not shown in Fig. 5.This is due to the extremely low values in the dust ratios (see also Fig. 3).An analog situation is also found for chemical evolution models evolving to larger metallicities (see .Dust-to-gas (DtG), dust-to-stellar (DtS) and dust-to-metal (DtM) mass ratios predictions from galaxy evolution models assuming a metal-poor scenario.Top panels: predictions for an intermediate-dust scenario (see 3.2.1).Bottom panels: predictions for a dust-rich scenario (see 3.2.1).Light orange and dark red shaded regions show the theoretical model predictions for REBELS galaxies belonging to the low and high stellar mass bins, respectively.Squares with same color code represent low and high stellar mass REBELS galaxies data bins, respectively.Plus and crosses are the same as squares, but for  UV / [C II] weighted gas masses and H 2 gas masses only (see 5. 5.1.2).Therefore, we can exclude that models with such a slow dust buildup can be a reliable solution to explain REBELS targets.
By looking at Fig. 5 (upper panels), we note that the simulations for the intermediate-dust scenario are usually within the errorbars of the inferred DtG, DtS and DtM for REBELS targets.However, it is worth noting that the DtS predicted for the low  * bin is at the lower edge of the large data errorbar.In general, the predictions of the models are below the data bins, therefore suggesting a slightly faster dust buildup to match the observed DtS and the DtG and DtM using HI+H 2 gas estimates from Heintz et al. (2021) and Vizgan et al. (2022).
The models assuming a "dust-rich scenario" agree well with the DtG and DtM computed using the weighted UV/[C ii] gas masses (see Fig. 5 bottom panels).The evolutionary tracks for these models also agree well with the DtG and DtM vs. metallicity trends observed in Damped Lyman- systems (DLAs) at redshifts  ∼ 1 − 4 (Péroux & Howk 2020).By looking at the DtS, the models for the high-mass REBELS galaxies overestimate the quantities inferred for the correspondent bin.Such a discrepancy cannot be attributed to a significant scatter in the data, since the observed DtS have low rms (< 10 −3 ) relative to the mean observational uncertainties (see also Fig. 3).At the same time, the models for low-mass galaxies are in better agreement with the DtS data with respect to the ones shown in Fig. 5 upper panels.However, we still notice an underestimation of the dust mass relative to the stellar mass by the models for the low-mass REBELS targets.
The adoption of a top-heavy IMF (i.e.favouring the formation of massive stars) may help alleviating the discrepancies between the observed trends for the DtG and the DtS for low-mass REBELS targets.In fact, a top-heavy IMF increases the level of chemical enrichment for a given star formation history, due to the larger number of CC-SNe enriching the ISM with metals.The faster chemical enrichment also corresponds to a faster dust mass buildup and, as a consequence, to larger DtS ratios, but not to larger DtG ratios (see Fig. D1 in  Appendix D).This happens because the top-heavy IMF enhances dust production from SNe and favours dust growth due to the larger amount of metals available, but at the same time boosts the amount of gas available in the ISM due to the larger mass ejection rates from evolved stellar populations, and in particular from massive stars (see, e.g.Palla et al. 2020a).
Another explanation for the very high dust-to-stellar mass ratio for the low-mass subsample could also be driven by a significant underestimate of the stellar mass  * for these targets.In fact, Hubble Space Telescope (HST) counterparts of REBELS-14 and REBELS-39 (which are part of the low-mass subsample) show very strong emission by rest-frame optical nebular lines (Bowler et al. 2017), which can be explained only by the presence of a very recent, large burst of star formation.In turn, this can limit the ability of nonparametric SFHs to overcome the so-called outshining problem, i.e. the underestimation of the contribution of older stellar populations due to overwhelming light form recent bursts (e.g.Leja et al. 2019).The hypothesis of a stellar mass underestimation is also supported by the very large dynamical mass estimated from the [C ii] line width and size6 (see Schouws et al., in prep;Topping et al. 2022) of REBELS-39 (> 10 11 M ⊙ ).The derived stellar mass for this galaxy contributes only to ∼ 1% of the total dynamical mass within the [C ii]-emitting region, therefore leaving room for the presence of older stellar populations, significantly enhancing the total stellar mass.

Metal-rich scenario
Moving to the models suited to reproduce a metal-rich scenario for REBELS targets (see Fig. 6), we show the predictions of the models in the case dust evolution parameters have intermediate values (intermediate-dust scenario, top panels) or dust production is particularly favoured (dust-rich scenario, bottom panels).The predictions of the models are compared with the DtG, DtS and DtM computed by assuming metallicities of 1/3 Z ⊙ and Z ⊙ , with the latter being a typical output of the models.
The models assuming the intermediate-dust scenario do reproduce the DtG and DtM ratios estimated by considering HI+H 2 mass calibrations.However, these models underestimate the DtS ratio.The reproduction of the estimated DtG and the underprediction of the DtS by the models may be driven by an overestimation of the gas mass through the [C ii] calibrations: at variance with what happens for lowmass REBELS galaxies in the low-metallicity case, here the problem stands for both low-and high-mass targets, with the latter showing large stellar masses (log( * /M ⊙ ) ≳ 10) for galaxies at  ≳ 6, which are less likely to be underestimated, despite the uncertainties.As a further proof, the comparison with the local DustPedia sample (pale grey crosses, Casasola et al. 2017Casasola et al. , 2020)), redshift  ≃ 2 strongly star forming objects (green diamonds, Shapley et al. 2020) and the galaxy GNz7q at  ≃ 7.2 (Fujimoto et al. 2022) show comparable DtS ratios to the REBELS ones, but almost an order of magnitude larger DtG ratios.We will come back to this point in 5.2.
In Fig. 6 lower panel, we see instead what happens when assuming a dust-rich scenario for the models reaching ∼ solar metallicity.The galactic evolution models do reproduce the observed DtS, especially for the low stellar mass bin.However, as also seen for models in the upper panel, we do not find unambiguous agreement between data and models in the different dust scaling relations.Here, we do not find any agreement with the data for DtG and DtM, even considering only the H 2 estimated mass contributing to the galactic gas budget.The models show in fact a steep rise in the predicted DtG and DtM between log(O/H) + 12 = 8.2 dex and 8.5 dex that allows reproducing the constraints for the bulk of local galaxies,  ≃ 2 dusty starbursts and the galaxy GNz7q at  ≃ 7.2, but at the same time prevents any agreement with the data presented in Section 2, further suggesting an overestimation of the gas mass for the REBELS targets.The jump to large values in the DtG and DtM is due to the transition between two regimes.In the first, the dust mass buildup is primarily set by stardust production, while in the second one dust growth also plays a fundamental role, leading the dust amount to saturate since no gas-phase metals available to form dust are left (e.g.Hirashita 2015; Aoyama et al. 2017Aoyama et al. , 2018Aoyama et al. , 2020)).This transition corresponds to the so-called "critical metallicity" regime (Asano et al. 2013, see also Gioannini et al. 2017;Graziani et al. 2020;Di Cesare et al. 2023), which happens at metallicities ≳ 0.1 ⊙ , with the exact value depending on different galactic evolution parameters, e.g. the galactic timescale/efficiency of star formation and the typical dust growth timescale (see, e.g.Hirashita 2015;Aoyama et al. 2017;Hirashita & Murga 2020).It is worth noting that the transition between the regimes is also present when looking at models with mixed dust prescriptions (see Appendix D2): while for metal-poor models SN dust production is the main driver of the resulting dust mass, in the metal-rich case it is the balance between dust growth and destruction that is discriminant for the level of dust enrichment.

Discussion
By looking at Fig. 5 and 6, a striking feature in the dust evolution tracks is the similarity, at similar levels of efficiency of star formation and therefore metallicity, between models with rather different SFHs (see Fig. 1).Such an outcome highlights how the dust evolution modelled in high-redshift galaxies is mostly influenced by parameters setting the efficiency of the different dust processes, i.e. production, growth and destruction (see also Appendix D2).Popping & Péroux (2022) showed that the assumption of different dust parameters and/or parametrisation is the most critical ingredient influencing the predictions of dust evolution at different redshifts in the framework of several SAMs and hydro simulations (e.g.Popping et al. 2017;Li et al. 2019;Triani et al. 2020;Graziani et al. 2020), however not considering the effect of the different modeling of the rest of the baryonic physics in the simulations.Here, we further support this conclusion by adopting a simpler framework on a galactic scale (instead of cosmic), by showing the very weak dependence on the particular SFH shapes for the inferred dust scaling relations.
The comparison of our model tracks for different metal enrichment scenarios with the REBELS data also highlights other important features.On the one hand, our metal-poor models generally reproduce the constraints from REBELS galaxies without assuming extreme prescriptions for dust evolution (see Fig. 5), in agreement with Di Cesare et al. (2023); Schneider et al. (in prep.).Only for the REBELS galaxies with the lowest stellar masses we discuss the possibility of a top-heavy IMF to reconcile the DtG and DtS measurements.These results are partly in tension with what was found in Dayal et al. (2022), where some of the REBELS galaxies where only reproduced by a "maximal dust model" for which extreme dust recipes (e.g.extremely fast dust growth timescale, no dust destruction by SNe) were adopted.However, it is worth noting that the stellar masses adopted here (from Topping et al. 2022) are different from those adopted in Dayal et al. (2022) (from Stefanon et al., in prep.).In turn, the different stellar masses lead to different dust masses in the methodology employed by Sommovigo et al. (2022a,b).Therefore, it is likely that the different dataset of stellar masses might lie at the heart of this need for extreme dust prescriptions and therefore of the discrepancy with our results.
On the other hand, the models tracks mimicking a metal-rich scenario (see Fig. 6) are unable to simultaneously reproduce the data for the DtG, DtM and DtS ratios independently of the assumed dust evolution scenario.Regarding this latter point, by adopting the most recent dust mass estimates by Algera et al. (2023) for REBELS-25 (log( dust /M ⊙ ) = 8.13 +0.38  −0.27 ) and REBELS-38  (log( dust /M ⊙ ) = 7.90 +0.39 −0.35 )7 based on additional sampling of the FIR continuum (see Section 2), the model difficulties in simultaneously reproducing the DtG and DtS are partially alleviated.This is because of the lower  dust in Algera et al. (2023), leading to larger dust masses.The comparison between the models and this other dataset can be observed in Fig. 7, where the DtG and DtS derived using newer dust mass estimates are shown.A relatively good agreement between data and model DtGs for the metal-rich case can be recovered for REBELS-25 by assuming molecular gas mass alone, whereas by adding HI masses the models still largely overestimate the dust-togas ratios.However, the larger dust masses estimated by Algera et al. (2023) give a good agreement with the predicted dust-to-stellar mass evolution tracks for the dust-rich scenario.This same conclusion is reached when the comparison is made with dustygadget, supporting the fact that at least some of these galaxies may have lower dust temperature and higher dust masses compared to what is estimated on the basis of single band observations (Algera et al. 2023).
This agreement also reinforces our previous statement of no need for extreme assumptions for the buildup of dust in UV-bright, dusty galaxies as the ones in REBELS.The prescriptions and dust evolution parameters used to reproduce the dust content in these high-redshift galaxies do not seem to differ much from the framework used for low-redshift galaxies (e.g.Asano et al. 2013;Gioannini et al. 2017;Galliano et al. 2018;Calura et al. 2023).This can be clearly seen in Fig. 8, where we show the timescales of stardust production, dust grain growth and dust destruction in the three REBELS galaxies with different SFHs shown in Fig. 1, i.e.REBELS-14, 29 and 25, compared with dust timescales found in local studies (De Looze et al. 2020;Galliano et al. 2021).In particular, we show the evolution of the different dust timescales in the two antipodal scenarios shown in Figs. 5 and 6, i.e. the metal-poor, intermediate-dust (top panels) and the metal-rich, dust-rich one (bottom panels).Despite of the large differences, i.e. of the order of one magnitude, between the timescale tracks for dust grain growth and SN dust destruction in upper and lower panels, we observe that in both metal and dust scenarios all dust timescales are well above few tenths of Myr at all evolutionary times.In particular, a striking feature for all panels is that the timescales found for individual REBELS galaxies are always comparable to the ones found in local studies, either in the case we have a higher dust destruction per SN event (upper panels) or we favour stardust production and dust grain growth (lower panels).

A problem of gas mass at high metallicity?
The marginal global agreement between models and data within the metal-rich scenario raises the question about the origin of this discrepancy: even allowing for additional combinations in the dust evolution parameters (see Appendix D2), none of the models is able to simultaneously recover the different dust scaling relations obtained by using the data presented in Section 2.
To further explore this, we look at the gas fraction  gas, [C II] /( gas, [C II] +  * ) inferred from observations, and compare these with the predictions of our chemical evolution models (see Fig. 9 top and bottom panels for the metal-poor and metal-rich scenario, respectively).
The gas fractions obtained using [C ii] calibrated HI+H 2 masses yield very large values, i.e. of the order of 0.9.While this is consistent with model predictions at low metallicities (Fig. 9 upper panel), large gas fractions do not match with what is obtained by galaxy evolution models at high metallicities (Fig. 9 lower panel).The motivation behind this discrepancy can be searched in some basic consideration of the chemical enrichment process that is irrespective of the mechanism of the stellar mass buildup in galaxies, from gas accretion from the cosmic web to galaxy mergers.In fact, to allow for a significant ISM metal enrichment, a large number of subsequent stellar generations are needed.At the same time, the short cosmic timescales for these high-redshift galaxies (≲ 0.85 Gyr at  ≃ 6.5) only allow for intermediate-massive stars ( ≳ 3M ⊙ ) to contribute to the chemical enrichment process.Due to the negative power-law nature of the stellar IMF at such masses, the large majority of the gas mass accreting the galactic system remains trapped in long-lived stars, therefore lowering significantly the available gas fraction.This leads to the suggestion that the gas masses estimated through [C ii] may significantly overestimate the gas content at high metallicities and, as a consequence, be at the basis of the inconsistency between models and data in the different dust scaling relations.As a proof of this, Fig. 9 shows that if we only assume the [C ii]-calibrated H 2 mass to contribute to the galactic gas budget we better match typical values observed in other intermediate to high-redshift sources with high metallicities (Shapley et al. 2020;Fujimoto et al. 2022).Moreover, the modeled gas mass fractions reproduce the values in the REBELS high-mass bin when considering H 2 masses only.Considering the acceptable match between data and models for the DtG and DtS in Fig. 7 when adopting only H 2 calibration for the gas mass, this clearly highlights the need for lower gas masses derived through [C ii] calibration at high-metallicities.
Therefore, calibrations from additional FIR lines (e.g.[OI]63m, see Wilson et al. 2023) may be needed to depict a clearer view on the atomic (but also total) gas contribution in these galactic sources, which is also crucial to understand their dust buildup.

Findings in the context of high-redhsift MZR relations
Several JWST programs have already started probing fundamental galactic scaling relations, such as the mass-metallicity relation (MZR) at redshifts  ≳ 4. In particular, Nakajima et al. (2023) studied the evolution of the MZR at  = 4 − 10 derived from ∼ 130 galaxies observed with JWST/NIRSpec taken from three major public spectroscopy programs (ERO, GLASS, and CEERS), using a similar stellar mass determination method to the one employed in this work (see also Harikane et al. 2023).They found that galaxies show a small evolution relative to the MZR at  ∼ 3 (Sanders et al. 2021), deriving this linear relation: log(O/H) + 12 = (8.24± 0.05) + (0.25 ± 0.03) log(  * 10 10 M ⊙ ), (11) in the stellar mass regime log( * /M ⊙ ) ≃ 7.5 − 9.5.Curti et al. (2023) found a shallower MZR, especially for redshift  > 6 galaxies.However, the sample collected by Curti et al. (2023) considered targets at the low-mass end of the MZR ( * /M ⊙ ∼ 10 6.5 − 10 8.5 , from the JADES program) that can significantly influence the inferred slope.Moreover, as also mentioned by the authors, the  > 6 redshift bin in Curti et al. (2023) is sparsely populated and probably subject to strong selection biases.
By extrapolating the relation by Nakajima et al. (2023) up to ∼ 10 10 M ⊙ , i.e. slightly above the upper mass limit of validity of Eq. ( 11), we find that our galaxies should have a metallicity between log(O/H) + 12 ≃ 8 dex (for the low-mass end of REBELS) and log(O/H) + 12 ≃ 8.35 dex (for the high-mass end).This means that low-mass REBELS galaxies should have a global metallicity of around 0.2Z ⊙ , i.e. between the low and intermediate metallicity data shown in Fig. 5. Therefore, low-mass REBELS galaxies should be in agreement with the metal-poor, dust-rich scenario outlined in 5.1.1,where most of the dust is the result of stardust production by SNe, with other dust processes as destruction and especially dust growth playing a minor role in shaping the dust budget (see also Appendix D2).Moreover, the comparison of the different dust scaling relations may suggest the presence of an IMF favoring the formation of massive stars or a slight revision of the computed stellar masses.
The most massive REBELS galaxies instead should have a metallicity of ∼ 0.5 Z ⊙ , that would be compatible with the model predictions presented in 5.1.2for high-metallicity galaxies.This would indicate either i) rapid dust buildup in case gas masses are overestimated by [C ii]- gas calibrations, or ii) slower dust formation in case stellar masses are underestimated.The former scenario seems favoured according to the set of stellar masses adopted, which shows larger values than those determined through other SED fitting methods (Topping et al. 2022).In addition, the larger dust masses by Algera et al. (2023) for two high-mass REBELS galaxies would also make the scenario i) more appealing.Therefore, different indications point towards a metal-rich, dust-rich scenario, where the large dust-to-gas, stellar and metal ratios are the result of a concurrent production of dust by stars and ISM growth processes.In such a case, the positive difference between the dust growth and the dust destruction rate is comparable to the stardust production rate (see also Fig. 8 to look at this in terms of dust timescales).
With this section, we have tried to place the dust evolution of REBELS galaxies in context of the first JWST results about the MZR relation for galaxies across the first Gyr of cosmic time.Drawing a detailed analysis about the relation between REBELS targets and the high-redshift galaxy scaling relations (which are prone to be revised in the next JWST cycles) is beyond the scope of this work.Conversely, we want to provide insights on the scenarios that JWST spectroscopy may open in the analysis of UV-bright, dusty galaxies in the EoR.In this way, we want to stress the urgent need for JWST followup for a statistically consistent number of FIR detected sources to characterize in much greater detail the dust evolution in the highredshift Universe.

CONCLUSIONS
In this paper, explored the possible scenarios of metal and dust evolution that apply to the ALMA REBELS (Bouwens et al. 2022) galaxy sample during the Epoch of Reionization (EoR).We took advantage of the sources with both [C ii]158m line and FIR continuum detections, which makes this sample of galaxies the largest with simultaneous stellar and dust mass determinations in the EoR (Sommovigo et al. 2022a;Inami et al. 2022), and extracted the gas masses from recent, metallicity dependent calibrations based on the [C ii] line luminosity (Heintz et al. 2021;Vizgan et al. 2022).We then modeled the chemical and dust evolution in these galaxies by means of detailed models adopting the non-parametric star formation histories (SFHs) derived by Topping et al. (2022) for each galaxy of the REBELS sample.In particular, we allow the gas flows and dust evolutionary history to vary in the model galaxies in order to simulate the different scenarios of metal and dust enrichment that may be found for such objects.
In this way, we wanted to provide clues on the dominant dust production mechanisms in EoR galaxies in different scenarios of metal enrichment, which can be probed by means of rest-frame optical ISM line spectroscopy from JWST.Moreover, the analysis performed throughout the paper also provides a test for [C ii]-gas mass calibrations (Zanella et al. 2018;Heintz et al. 2021;Vizgan et al. 2022), which play a crucial role in the study of galaxies at redshift  ≳ 6.
Our results can be summarized as follows: (i) at similar levels of efficiency in star formation and therefore metallicity, the predicted dust evolutionary tracks for different SFHs are very similar when assuming the same parametrizations for the different dust processes, i.e. stardust production, dust growth and dust destruction.This highlights that the modelling of dust evolution in high-redshift galaxies is primarily influenced by the assumptions and parameters setting the efficiency of the different dust processes in the ISM, in agreement with Popping & Péroux (2022); (ii) as expected from previous dust modelling at  ≳ 5 (Dayal et al. 2022;Di Cesare et al. 2023), models showing a slow dust buildup (dust-poor scenario) are far from reproducing the dust constraints from REBELS galaxies and therefore should be excluded when studying "normal" star-forming galaxies in the EoR.At the same time, there is no need for extreme assumptions about dust formation and evolution for these sources, with dust timescales comparable to what claimed by several studies focused on local galaxies (e.g.De Looze et al. 2020;Galliano et al. 2021).However, current observational constraints, i.e. the galactic dust and stellar masses, are still insufficient to discriminate between different galactic evolution scenarios that in turn influence the relative importance of the different dust processes in reproducing the observed dust masses.Information on ISM metallicity will therefore have a crucial importance to start placing tighter constraints on the contribution of dust evolution processes; (iii) in case JWST observations will reveal low metal enrichment ( ≃ 0.1  ⊙ ), we demonstrate that dust evolution is mainly driven by stardust production, whose rate dominates over ISM dust processes (growth, destruction).Either models where dust parameters are set to have a fast dust buildup (dust-rich scenario) or a quieter evolution (intermediate-dust scenario) are compatible with the constraints from dust-to-gas (DtG), dust-to-stellar (DtS) and dust-to-metal (DtM) mass ratios.However, models for REBELS targets with the lowest stellar masses ( * ∼ 10 9 M ⊙ ) generally underpredict the observed DtS.A possible way to reconcile models with observations for these sources is to adopt a top-heavy IMF (i.e.favouring the production of massive stars, see e.g.Palla et al. 2020a).However, indications from HST (Bowler et al. 2017) and dynamical mass estimates from the [C ii] line (Topping et al. 2022) may also indicate an underestimation of the stellar mass for these targets; (iv) in the case of sources evolving up to solar metallicity, the level of dust enrichment in galaxies is mostly driven by the balance between ISM dust growth and destruction processes.The comparison with REBELS galaxy measurements reveal that different galaxy models struggle to reproduce simultaneously the estimated DtG/DtM and DtS.In particular, models with a intermediate-dust scenario do reproduce the estimated DtG/DtM ratios while they underestimate the observed DtS.At the same time, models with a fast dust buildup do reproduce the observed DtS, while they significantly overestimate the DtG/DtM.This may be due to an overestimation of the gas mass by means of [C ii] line calibrations (Heintz et al. 2021;Vizgan et al. 2022) at high metallicities: gas fractions using these calibrations yield very large values which are not compatible with the presence of a chemically evolved galaxy basing on general consideration of galactic chemical enrichment.To alleviate the disagreement, lower gas masses relative to the reference dataset should be considered; (v) looking at the first JWST probes of the mass-metallicity relation (MZR) at high-redshift (e.g.Nakajima et al. 2023), REBELS galaxies should exhibit metallicities between ∼ 0.2 and ∼ 0.5 ⊙ .The galaxies with the lowest stellar masses should be representative of the metal-poor case mentioned in point iii) of this Section.For these galaxies, our modeling suggests a dust-rich scenario.Galaxies in the upper mass end should instead represent a case compatible to the metal-rich one.For such galaxies, either a dust-rich scenario in case gas masses are overestimated or an intermediate-dust scenario in case stellar masses are underestimated may be suggested, with the former more likely according to point (iv) of this Section.
Present and upcoming JWST/NIRSpec programs at intermediateto high-resolution (R ∼ 1000 − 2700) are allowing to resolve several nebular emission lines in galaxy spectra in  > 6 galaxies, such as [OII]3727, 3729 and [OIII]4960, 5008 doublets, [NII]6585 and H lines (for  < 6.7 galaxies) H, H and possibly even the auroral [OIII]4363 transition.This opens up the possiblity to determine the level of chemical enrichment in REBELS sources as we get access to the most common strong line metallicity indicators (e.g.O3, O32, R2, R23), even when direct-  measurements will not be possible.Considering uncertainties of the order of 0.2-0.3dex within log(O/H)+12 determinations (e.g.Nakajima et al. 2023;Bunker et al. 2023;Curti et al. 2023), such metallicity measurements will be sufficient to start discriminating between the different scenarios illustrated in this paper.
Therefore, the results presented here may be used as a viable tool to interpret and constrain the metal and dust evolution in EoR galaxies.At the same time, these results highlight even more the great potentials of the ALMA-JWST synergy to unveil some of the most debated issues in contemporary astronomy.
for a forsterite compound (Mg 2 SiO 4 ) we need 2 O atoms per Mg one (with the latter fixing maximum number of molecules that can form).For every Mg atom, more than 12 O atoms are available: therefore, the number of O atoms trapped in forsterite are around 15% of the total.However, the 20% of Si remains spare.The latter can form ferrosilite (FeSiO 3 ) or fayalite (Fe 2 SiO 4 ) molecules.To maximise the O depletion, let us assume that fayalite is formed.In this case around 10% of oxygen atoms remain trapped in this molecule.Summing it with the fraction of O in forsterite, we get a O depletion level around 25%.
C dust is instead mostly found in the ISM in the form of graphite and amorphous carbon.However, a significant fraction of C is locked in the very stable CO molecule (around 40% locally, Lodders 2010; Agúndez & Wakelam 2013).If we consider a lower limit of 20% of C atoms locked CO, then the maximum level of C depletion can arrive to 80%.
Up to now we have worked using number ratios between elements.However, the DtM ratio is a mass ratio: therefore, we need to convert the number ratios into mass ratios.By doing and this and then by summing up all the elemental contributions to ISM dust, we obtain that around 50% of all ISM metals can be depleted into dust.
Our calculations approximate dust as formed by a limited number of compounds.Other dust molecules (e.g.silicon carbide, SiC) are also present, however being much less abundant.At the same time, it was shown that interstellar ices can accrete silicate grains and form mantles around them.Nonetheless, several works indicate that these mantles are unlikely to survive to the harsh radiation field from new stellar generations, especially at high redshift (Ferrara et al. 2016;Gall & Hjorth 2018).This indicates very short lifetimes and therefore a very small contribution to the total dust mass budget in galactic evolution.For these reasons, it is safe to assume for each model timestep a slightly larger limit for the total dust mass, which translate into a maximum DtM of 0.6.It is also worth noting that such limit was also found in ISM dust depletion studies (e.g.Konstantopoulou et al. 2022) in several environments, from the MW galaxy up to  ≳ 2 Damped Lyman- systems.
On another note, the interstellar abundance pattern considered for the above calculations is the solar one.Nonetheless the values displayed can be easy generalised in case of a subsolar chemical enrichment, as the one tested in this paper.In fact, while the magnesium and carbon to oxygen ratios are similar to those of the solar system, Si and mainly Fe are underproduced.This means less seeds for oxygen atoms depleted into dust and therefore slightly smaller maximal DtM.
Table B1.Prescriptions adopted throughout this work for the model galaxies representing different metallicity scenarios.The acronyms in the stellar yields column stand for: Limongi & Chieffi (2018) (LC18), Ventura et al. (2013Ventura et al. ( , 2018Ventura et al. ( , 2020Ventura et al. ( , 2021) )    Effect of SN dust production prescriptions on the dust-to-gas (DtG) and dust-to-stellar (DtS) mass ratios predictions from galaxy evolution models for REBELS-14 assuming a metal-poor scenario (upper panels) and REBELS-25 assuming a metal-rich scenario (lower panels).Thick dashed and dash-dotted lines show results for the dust-rich and intermediate-dust scenario (see 3.2.1)analysed in Section 5.The thin dashed and dash-dotted lines represent a dust-rich scenario with smaller SN dust production and an intermediate-dust scenario with larger SN dust production, respectively.Data legend is as in Fig. 5.

Figure D3.
Effect of dust growth prescriptions on the dust-to-gas (DtG) and dust-to-stellar (DtS) mass ratios predictions from galaxy evolution models for REBELS-14 assuming a metal-poor scenario (upper panels) and REBELS-25 assuming a metal-rich scenario (lower panels).Thick dashed and dash-dotted lines show results for the dust-rich and intermediate-dust scenario (see 3.2.1)analysed in Section 5.The thin dashed and dash-dotted lines represent a dust-rich scenario with larger dust growth and an intermediate-dust scenario with smaller dust growth, respectively.Data legend is as in Fig. 5.

Figure D4.
Effect of SN dust destruction prescriptions on the dust-to-gas (DtG) and dust-to-stellar (DtS) mass ratios predictions from galaxy evolution models for REBELS-14 assuming a metal-poor scenario (upper panels) and REBELS-25 assuming a metal-rich scenario (lower panels).Thick dashed and dash-dotted lines show results for the dust-rich and intermediate-dust scenario (see 3.2.1)analysed in Section 5.The thin dashed and dash-dotted lines represent a dust-rich scenario with larger SN dust destruction and an intermediate-dust scenario with smaller SN dust destruction, respectively.Data legend is as in Fig. 5.

Figure 1 .
Figure 1.Non-parametric star formation histories (SFHs) for REBELS galaxies from Topping et al. (2022).Left, central and right panels show the SFH of REBELS-14, REBELS-29, and REBELS-25.The shaded areas represent the 1- confidence interval for the SED-modelled SFHs.The orange lines are the smoothed SFHs used as inputs for our galactic chemical evolution models.

Figure 2 .
Figure2.Evolution of the metallicity (both in terms of the oxygen abundance log(O/H)+12 and solar scaled metallicity Z/Z ⊙ ) in modelled REBELS galaxies assuming different gas infall masses ( gas,tot = ×3, ×10, ×15) and therefore metallicity evolution scenarios (see also Tab.B1).Left, central and right panels show the evolution for REBELS-14, REBELS-29 and REBELS-25.The thin lines in the left and central panels represent models in which multiple gas accretion episodes are allowed, as described in 3.2.

Figure 3 .
Figure 3. Evolution of the dust-to-stellar (DtS) mass ratio as a function of galactic stellar mass  * predicted by models for REBELS-14, REBELS-29 and REBELS-25 for different metallicity and dust evolutionary scenarios.Top panels show the results for galaxy evolution models assuming a metal-poor scenario, while lower panels for models assuming a metal-rich scenario (see also Tab.B1).Left, central and right panels indicate different dust evolutionary scenarios, i.e. dust-poor, intermediate-dust and dust-rich, respectively (see also Tab.B2).The color code for the galaxies with different star formation histories is the same as in Fig. 2. Small red squares are REBELS data for individual galaxies, whereas light orange and dark red large squares with errorbars represent low and high stellar mass REBELS galaxy bins with the average error within the bin, respectively.Additional data are from Casasola et al. (2020) (DustPedia galaxies), Shapley et al. (2020) (luminous SF galaxies), Watson et al. (2015); Hashimoto et al. (2019); Reuter et al. (2020); Bakx et al. (2021); Fujimoto et al. (2022) ( ∼ 7 galaxies).Average errorbars for these data sets are indicated in the top-left part of each panel.

Figure 4 .
Figure 4. Evolution of the ratio between stardust production and dust growth timescales for gas-rich (i.e. more metal-poor, solid lines) and gas-poor (i.e. more metal-rich, dashed lines) models in the dust-rich scenario (see 3.2.1).The color code for the galaxies with different star formation history is the same as in Fig.2.The thin, grey dashed line indicates a ratio equal to 1, corresponding to the equilibrium between the two dust processes.

Figure 5
Figure5.Dust-to-gas (DtG), dust-to-stellar (DtS) and dust-to-metal (DtM) mass ratios predictions from galaxy evolution models assuming a metal-poor scenario.Top panels: predictions for an intermediate-dust scenario (see 3.2.1).Bottom panels: predictions for a dust-rich scenario (see 3.2.1).Light orange and dark red shaded regions show the theoretical model predictions for REBELS galaxies belonging to the low and high stellar mass bins, respectively.Squares with same color code represent low and high stellar mass REBELS galaxies data bins, respectively.Plus and crosses are the same as squares, but for  UV / [C II] weighted gas masses and H 2 gas masses only (see 5.1.1).High stellar mass bin data are shifted by 0.05dex in metallicity for sake of clarity.Additional data are from Casasola et al. (2020) (DustPedia galaxies), Shapley et al. (2020) (luminous SF galaxies), Péroux & Howk (2020) (DLA systems), Heintz et al. (2023) (I Zw 18, S04590), Fujimoto et al. (2022) (GNz7q).
Figure5.Dust-to-gas (DtG), dust-to-stellar (DtS) and dust-to-metal (DtM) mass ratios predictions from galaxy evolution models assuming a metal-poor scenario.Top panels: predictions for an intermediate-dust scenario (see 3.2.1).Bottom panels: predictions for a dust-rich scenario (see 3.2.1).Light orange and dark red shaded regions show the theoretical model predictions for REBELS galaxies belonging to the low and high stellar mass bins, respectively.Squares with same color code represent low and high stellar mass REBELS galaxies data bins, respectively.Plus and crosses are the same as squares, but for  UV / [C II] weighted gas masses and H 2 gas masses only (see 5.1.1).High stellar mass bin data are shifted by 0.05dex in metallicity for sake of clarity.Additional data are from Casasola et al. (2020) (DustPedia galaxies), Shapley et al. (2020) (luminous SF galaxies), Péroux & Howk (2020) (DLA systems), Heintz et al. (2023) (I Zw 18, S04590), Fujimoto et al. (2022) (GNz7q).

Figure 6 .
Figure6.Dust-to-gas (DtG), dust-to-stellar (DtS) and dust-to-metal (DtM) ratios predictions from galaxy evolution models assuming a metal-rich scenario.Top panels: predictions for an intermediate-dust scenario (see 3.2.1).Bottom panels: predictions for a dust-rich scenario (see 3.2.1).Models and data legend are as in Fig.5.

Figure 7 .
Figure 7. Dust-to-gas (DtG) and dust-to-stellar (DtS) mass ratios predictions from galaxy evolution models for the high stellar mass subsample assuming a metal-rich, dust-rich scenario (dark red shaded regions).Purple and cyan symbols refer to the dust mass estimates by Algera et al. (2023) for REBELS-25 and REBELS-38, respectively.REBELS-38 data are shifted by 0.05dex in metallicity for the sake of clarity.Dark red symbols and additional data are as in Fig. 5.

Figure 8 .
Figure 8. Evolution of dust formation and destruction timescales predicted by models for REBELS-14, REBELS-29 and REBELS-25 for different metallicity and dust evolutionary scenarios.Top panels show the results for galaxy evolution models assuming a metal-poor, intermediate-dust scenario, while lower panels for models assuming a metal-rich, dust-rich scenario (see also Tabs.B1 and B2).Solid lines represent stardust production, dashed lines dust growth, and dash-dotted lines dust destruction.For the sake of comparison, we overlay the average dust formation and destruction timescales obtained for local galaxies by De Looze et al. (2020) and Galliano et al. (2021).

Figure 9 .
Figure 9. Gas fraction evolution with metallicity as predicted by galaxy evolution models compared with gas fraction derived through different [C ii] calibrations for the gas mass.Data legend is as in Figs. 5. Top panel: prediction for the metal-poor scenario.Bottom panel: prediction for the metal-rich scenario.For both panels, the gas fraction estimated for individual REBELS galaxies using their dynamical mass is shown in the right-side panels according to their stellar mass (see bottom colorbar).
Figure D1.Effects of IMF variation on the dust-to-gas (DtG) and dust-to-stellar (DtS) mass ratios predictions from galaxy evolution models for REBELS-14 (representative of the low-mass REBELS sample) assuming a metal-poor scenario.Solid lines show results for a Chabrier (2003) IMF, while dashed lines show results for a top-heavy IMF.Thick lines stand for a dust-rich scenario (see 3.2.1),while thin lines for an intermediate-dust scenario (see 3.2.1).Data legend is as in Fig. 5.

Figure D2.
Figure D2.Effect of SN dust production prescriptions on the dust-to-gas (DtG) and dust-to-stellar (DtS) mass ratios predictions from galaxy evolution models for REBELS-14 assuming a metal-poor scenario (upper panels) and REBELS-25 assuming a metal-rich scenario (lower panels).Thick dashed and dash-dotted lines show results for the dust-rich and intermediate-dust scenario (see 3.2.1)analysed in Section 5.The thin dashed and dash-dotted lines represent a dust-rich scenario with smaller SN dust production and an intermediate-dust scenario with larger SN dust production, respectively.Data legend is as in Fig.5.
.  gas, [C II] ranges are computed by summing the sum of HI and H 2 masses assuming metallicities between 0.1  ⊙ and  ⊙ .