Observational and theoretical constraints on the formation and early evolution of the ﬁrst dust grains in galaxies at 5 < z < 10

Context. The ﬁrst generation of stars were born a few hundred million years after the big bang. These stars synthesise elements heavier than H and He, which are later expelled into the interstellar medium, initiating the rise of metals. Within this enriched medium, the ﬁrst dust grains were formed. This event is cosmologically crucial for molecule formation, as dust plays a major role by cooling low-metallicity star-forming clouds, which can fragment to create lower mass stars. Collecting information on these ﬁrst dust grains is di ﬃ cult because of the negative alliance of large distances and low dust masses. Aims. We aim to combine the observational information from galaxies at redshifts 5 (cid:46) z (cid:46) 10 to constrain their dust emission and theoretically understand the ﬁrst evolutionary phases of the dust cycle. Methods. Spectral energy distributions (SEDs) are ﬁtted with CIGALE and the physical parameters and their evolution are modelled. From this SED ﬁtting, we built a dust-emission template for this population of galaxies in the reionisation epoch. Results. Our new models explain why some early galaxies are observed and others are not. We follow in time the formation of the ﬁrst grains by supernovae later destroyed by other supernova blasts and expelled in the circumgalactic and intergalactic media. Conclusions. We ﬁnd evidence for the ﬁrst dust grains formed in the universe. But above all, this work underlines the need to collect more data and to develop new facilities to further constrain the dust cycle in galaxies in the reionisation epoch.


Introduction
Understanding the characteristics of the dust cycle is still an issue. This statement is even truer for Lyman break galaxies at all redshifts (LBGs, e.g. Burgarella et al. 2007) and at redshifts z > 5 (Bouwens et al. 2011;Finkelstein et al. 2010) at the epoch of reionization (EoR). This is due to the severe limits due to submillimetre (submm) and millimetre (mm) observations but also because we have a poor idea on the chemical conditions, the dust grain characteristics and even the stellar populations.
Fundamentally, as in the local universe, dust grains in highredshift (Hi-z) galaxies absorb far ultraviolet (FUV) photons from young and massive stars. By doing so, these grains warm up to dust temperatures of a few tens of degrees and emit farinfrared (FIR) photons.
This paper re-evaluates the roles of the relevant physical processes in the formation and early evolution of dust grains in these LBGs: formation in supernovae (SNe) explosion and outflows from lower mass stars (e.g. Todini & Ferrara 2001;Matsuura et al. 2015;Ventura et al. 2012;Nanni et al. 2013Nanni et al. , 2014Marassi et al. 2019), destruction by SNe shocks, growth by accretion (e.g. Hirashita & Omukai 2009;Dwek & Cherchneff 2011;Asano et al. 2013), and galactic outflows (e.g. Jones et al. 2018;Ohyama et al. 2019). Within the context of this paper, two main questions can be asked: What is the dust cycle in high-redshift objects?. Two types of stellar sources form and eject dust grains in the interstellar medium (ISM) of galaxies: stars with an initial mass in the 1-8 M range during the asymptotic giant branch (AGB) phase, and stars with an initial mass in the 8-40 M range during the core-collapse SNe phase. Dust grains also form seeds that grow in the ISM (e.g. Dwek 1998;Draine 2009;Jones & Nuth 2011;Asano et al. 2013). The AGB dominates dust production in galaxies old enough to allow low-mass stars to evolve to the AGB phase (Valiante et al. 2009). The AGB's contribution to the global dust production in metal-poor environments reaches at most 30% after about 2 Gyr (Dell'Agli et al. 2019). This is confirmed in other works (e.g. Dwek & Cherchneff 2011;). However, SNe blasts produce a shockwave able to destroy dust grains in the ISM (e.g. Slavin et al. 2015;Matsuura et al. 2019;Nozawa et al. 2003;Dwek & Cherchneff 2011).
There are numerous works (e.g. review by Rupke 2018) showing that galactic winds are ubiquitous at low redshifts. These multi-phase winds carry neutral and ionised gas and dust grains outside the galaxies. For instance, in M 82, dust is found well beyond the radius where gas is thought to be outflowing from the galaxy. This suggests that dust escapes from the gravitational potential of galaxies into the intergalactic space (e.g. Engelbracht et al. 2006;Yoshida et al. 2019). Dust is also observed in the outflows from other local star-forming galaxies and strong evidence for dust in the outflowing neutral clouds is also found at z = 3 (Shapley et al. 2003). Very recently, in a very detailed analysis, Jones et al. (2018) found that their measurements favour a picture in which the majority of heavy elements are ejected in a predominantly low-ionisation outflow, which regulates galactic chemical evolution.
The origin and the evolution of dust at high redshift has become a prominent issue since the early works of Watson et al. (2015). At z > 5, large amounts of dust (∼10 8 M ) are found in galaxies (e.g. Riechers et al. 2013;Mortlock et al. 2011;Venemans et al. 2012;Michałowski 2015) or quasars (e.g. Bertoldi et al. 2003;Priddey et al. 2003;Robson et al. 2004;Beelen et al. 2006;Michałowski et al. 2010;Cullen et al. 2017;Leśniewska & Michałowski 2019). Recently, we started to collect new and interesting data on the galaxy host of γ-ray bursts that will be very useful for the study of the properties of dust and the comparison of dust extinction and dust attenuation laws in these objects (e.g. Stratta et al. 2011;Hjorth et al. 2013;Bolmer et al. 2018).
How can we explain the rest-frame UV colour properties and detectability of Hi-z LBGs? Most Hi-z LBGs observed in FIR are not detected (e.g. Capak et al. 2015;Bouwens et al. 2016;Hirashita et al. 2017). Explanations are proposed to understand the non-detections. Ferrara et al. (2017) suggested that these objects can contain a large molecular gas mass. While dust in the diffuse ISM attains relatively high temperatures (T ≈ 70−80 K), dust embedded in dense gas remains cold (T ≈ 30−40 K). Faisst et al. (2017) suggested that radiation pressure causes a spatial offset between dust clouds and young star-forming regions within the lifetime of O/B stars. These offsets modify the radiation balance and create viewing-angle effects that can change UV colours at fixed IRX (IRX = log (L dust / L FUV )).
We use a seven-year WMAP7 cosmology (Komatsu et al. 2011) as adopted by Astropy, of Ω M = 0.273, Ω Λ = 0.727 and H 0 = 70.4. We assume a Chabrier initial mass function (IMF, Chabrier 2003) for the SED fitting because top-heavy IMF (m −α with α = 1.0, 1.35, 1.5, Dabringhausen et al. 2010;Gall et al. 2011) single stellar population (SSPs) models are not available in CIGALE. The Chabrier IMF for the SED fitting is fixed in the one assumed by Bruzual & Charlot (2003), that is, the lower and upper mass cut-offs m Low = 0.1 M and m Up = 100 M . Chabrier and top-heavy IMFs are assumed for the dust modelling lower and upper mass cut-offs m Low = 1.3 M and m Up = 40 M (with Kobayashi yields). In agreement with Sect. 5, with the exception of the pop III stars used between 10 and 30 M .

Building spectral energy distributions
The Hi-z LBG sample is inhomogeneous and collected from several papers as listed in Table 1. It corresponds to an FUV selection. The criteria are as follows: we selected all LBGs to have spectroscopic redshifts in the 5 z 10 range for which a mm detection was attempted (successful or not) in the rest-frame FIR (Table 1). Practically speaking, this means that we focused on LBGs observed in ALMA's bands 6 and 7. This redshift range (at EoR) provides constraints in rest-frame FIR, meaning, observed in the mm range. As an additional criterion, we only use LBGs with at least five data points in the rest-frame UV and optical.
We collected photometric data to build SEDs and derive stellar masses (M star ), dust masses (M dust ), star formation rates (SFRs), FUV and luminosities (L FUV and L dust ), the age of the dominant stellar population (Age main ), and the UV slope β (Calzetti et al. 1994). The set of data is consistently fitted with the same code, CIGALE (Burgarella et al. 2005;Noll et al. 2009;Boquien et al. 2019), and using the same dust models (Draine & Li 2007;Draine et al. 2014). The CIGALE modules and input parameters used for all the fits are described in Table 2. To model the dust emission, we also tested a modified, single dust temperature grey body, plus a mid-IR power law that approximates the hot-dust emission from the AGN (Active Galactic Nuclei) heating or clumpy hot star bursting regions (Casey 2012). The impact is negligible, and we base this paper on the former to constrain our models. The dust models are calibrated with a sub-sample of dwarf (compact with sizes ranging from 0.08 to 1.2 , Madden et al. 2013) low-metallicity galaxies (Low-zZ) with 12 + log(O/H) ≤ 8.4 in the Dwarf Galaxy Survey (DGS, Madden et al. 2013;Rémy-Ruyer et al. 2015). The DGS was originally selected from the Herschel Dwarf Galaxy Survey. It is a compilation of a SPIRE SAG2 guaranteed time (GT) key program and SPIRE GT2 observations, using the PACS and SPIRE instruments onboard the Herschel Space Observatory to obtain 55-500 µm photometry and spectroscopy of 48 dwarf galaxies in our local Universe, chosen to cover a broad range of physical conditions. We add data from the NASA Extragalactic Database (NED) in the UV (GALEX) and in the optical (B and R bands). When selecting the FUV, NUV, B, and R data in NED, we kept the data that provide total fluxes for our galaxies, and, when several total fluxes were available, we favoured new or newly reprocessed ones. If these criteria are not respected, no data were added for the corresponding objects, which in the end meant that we did not keep them in the Low-zZ sample.
The IR information on these Hi-z LBGs is limited. Our approach relies on the assumption that the physical conditions for dust grains are shared by these LBGs, which means that these LBGs present the same FIR SED shape within uncertainties. Here, we used luminosity-weighted SEDs normalised at λ = 200 µm, where the dust is optically thin (e.g. Casey 2012). Our LBG sample presents different FIR-to-FUV ratios, and the dust attenuations are different. We leverage the redshift range to combine and build the IR SEDs of our ALMA-detected LBGs, and thus gain access to better spectral information for these objects (Fig. 1). However, we have no observational information on the dust emission, at wavelengths shorter than the peak, that is, near-and mid-IR. From the assumption that these Hi-z LBGs are metal-poor (Castellano et al. 2014;Yuan et al. 2019), we use the minimum value for q PAH (fraction of the dust mass in the PAH) because q PAH is correlated with the metallicity (e.g. Ciesla et al. 2014). From our SEDs, γ, the fraction of the dust heated by starlight above the lower cutoff U min (the minimum of the distribution of starlight intensity relative to the local interstellar radiation field), in Draine and Li's models cannot be constrained. We assume that the best value for the DGS sample is valid for the Hi-z LBGs (but the impact of this parameter on our result is limited to ±8%, estimated using the range found by Álvarez-Márquez et al. (2019) on a sample of 22000 LBGs at z ∼ 3.0). In addition to the dust continuum information, we used the emission in the rest-frame FUV and optical range to model the stellar populations. Whenever available, [CII]158 µm and [OIII]88.3 µm were used to constrain the SED fitting. The nebular module in CIGALE uses nebular templates based on Inoue (2011), which were generated using CLOUDY 08.00 (Ferland et al. 1998). Each line has a Gaussian shape with A32, page 2 of 16 Notes. The Hi-z LBG sample contains 18 objects, while the Low-zZ sample contains 31 objects for which the final fits are adequate, i.e. χ 2 reduced ≤ 5. 9 Hi-z LBGs are detected in ALMA in Band 6 or Band 7, they provide 15 detected measurements, but only 11 measurements with S /N submm > 3 are used to build the IR template. a user-defined line width to take rotation into account, which can be especially important for narrowband filters and high-redshift objects due to the apparent line broadening with redshift in the observed frame. This normalised nebular emission-line spectrum is rescaled to the appropriate flux level by multiplying by the ionising photon luminosity that was computed along with the composite stellar population.

Methodology to derive the physical parameters
The SED analysis is performed in five phases.
Phase 0. A necessary preliminary phase was to perform a lot of initial fits on individual Low-zZ and Hi-z LBGs to understand the limits of the physical parameters to be explored (star formation history (SFH), dust attenuation, AGN fraction, etc.) in the CIGALE fits. Because of the exploration nature of these fits, these initial fits are not reported in Table 2. In order to perform this exploration of the parameter space, we first used different types of SFHs available in CIGALE: a periodic SFH, that is the same star formation event regularly occurring during the galaxy lifetime and a single or double exponential SFH. None of them provided an improvement of the fits. We also attempted to modify the metallicity of the SSPs but we did not manage to get any strong contraint. Finally, we tested various dust emission models within the set available into CIGALE. But the limited number of FIR data points for individual galaxies prevented us to carry on this study and led us to the template approach described below.
Phase 1. CIGALE fits individual Hi-z LBGs with a large distribution of priors for the dust parameters and SFH (delayed 1 , delayed + burst, and constant SFHs are assumed) with average reduced χ 2 ν < 5. On the other hand, constant SFHs provide results very close to the delayed ones. This was expected, given the ages of the stellar populations and the age of the universe. Fits with a delayed + burst SFH are not as good. As none of the stellar models available in CIGALE feature a topheavy IMF, we used a Chabrier IMF. Using a top-heavy IMF in the SED fitting instead of a Chabrier one could change the results of the SED fitting in the UV-optical range, but less so in the IR range. Future development of CIGALE will address this issue. Similarly, stellar models that include binaries (e.g. BPASS Stanway & Eldridge 2018) would change the UV-optical emission of these galaxies. We will address this possibility in a future paper. These initial fits provide an estimated rest-frame flux 200 µm to normalise the SEDs (with FIR detections) at λ = 200 µm where the dust is optically thin (e.g. Casey 2012). Detected LBGs at 5 z 10 made it possible to sample the spectrum of dust emission (Fig. 1). Since we do not have data beyond the rest-frame wavelength λ = 200 µm, normalizing at any longer wavelength would not be meaningful. Although only the Rayleigh-Jeans range is covered, dust temperatures in the range 40 K to 70 K fit the data better.
Phase 2. Local Low-zZ galaxies are often said to be analogues to early galaxies Capak et al. 2015;Hou et al. 2019). We compiled a sample of Low-zZ galaxies from Rémy-Ruyer et al. (2013, 2015. We kept DGS galaxies with 12 + log(O/H) < 8.4. With the same assumptions on the priors as for the LBGs, these galaxies were fitted with CIGALE. We checked that the main dust physical parameters found for the Low-zZ sample were in agreement with the ones from the Hi-z LBG sample (which confirms the partial analogy), and we leveraged the former fits to set the value of the γ parameter in the models from Draine & Li (2007) and Draine et al. (2014). We also checked that the value assumed for q PAH is consistent with our assumption of low metallicity (Fig. 2).
Phase 3. The previous phases suggested that a common template could be used to fit most of the objects in our samples. We interpret this possible generalization of the template as a similarity in the physical conditions in the Hi-z LBG sample. Note however, that there does not seem to be an universality as some objects present different different dust emission (for instance SBS0335-052) and it would interesting to explore the origin of the difference). In this phase, q PAH and γ are fixed to the value given in Table 2, and we fitted the IR-combined SEDs of the Hi-z LBGs (Fig. 3). Phase 3 provides physical parameters for the average combined SEDs of the Hi-z LBGs used to build the Hi-z IR template. The fit with a modified black body gives: β submm = 1.51 (slope on rest-frame sub-mm side of the dust emission) and T dust = 56.6 K (Fig. 4).
Phase 4. We re-fitted all the Hi-z LBGs over the FUV-to-FIR range with this Hi-z IR template (Tables 3 and 4). As a safety check, we also fit the Low-zZ galaxies with the same IR template. The quality of the fit is shown in the Annex.
In conclusion, after studying the individual dust emission of the Hi-z LBG sample, we found some similarities with the local Low-zZ sample. Furthermore, we find that most of the Hi-z LBG sample show a similar dust emission. We have characterised this dust emission using a modified blackbody and a single model from the library of Draine & Li (2007), Draine et al. (2014). We use the latter in the rest of the paper as it presents a more physical understanding of the physical conditions in the ISM of galaxies in the early universe. The above methodology was successful in fitting the majority of the Hi-z and Low-zZ samples.

Results
From Table 3 and 4, we compute the specific SFR (sSFR = SFR/M star ), the specific dust mass (sM dust = M dust /M star ), and IRX = L dust /L FUV . In Figs. 5 and 6, we define the dust-formation rate diagram (DFRD): sM dust vs. sSFR. The x-axis provides information on the stellar population and the y-axis on dust A32, page 4 of 16 Fig. 1. Dust emission of the detected Hi-z LBGs with S /N submm > 3 (11 data points) are normalised at λ rest−frame = 200 µm and combined to form a single Hi-z LBG template. For each detected LBG in our sample, we use the number of S /N submm > 3 ALMA detections, which are validated in agreement with our criteria, as described in Table 1. The redshift range (5 z 10) makes it possible to cover a wide wavelength range and explore dust emission. We over-plot colour-coded curves for SEDs based on modified blackbodies in the range 30 ≤ T dust [K] ≤ 85. grains. For some of the galaxies, the fit is not good, for example, SBS 035-052 (see Sect. 4). The two possible origins for these bad fits are the quality of the data, and/or the fact that our template does not match the dust properties of these objects. These galaxies are not used in our analysis. The physical parameters derived for the Hi-z LBG and the Low-zZ DGS samples by CIGALE are listed in Tables 3, 4, and C.1.
By definition, CIGALE checks that the observed SED is consistent with the assumed SFH, and we tried several SFHs before selecting the delayed one. A first result from the SED fitting, both for the Hiz-LBGs and the Low-zZ DGS objects is that both τ main and the range of ages are quite short (Tables 3, 4 and C.1) with ages in the range of 100-500 Myrs for the Hiz-LBGs, and 35-1300 Myrs for the Low-zZ sample. From this result, two important points should be noted: (i) we know that some young, and even very young, stellar populations are present in some of these objects. The most extreme one is the galaxy at z = 9.1 (Hashimoto et al. 2018) and observational evidences ([OIII]88.3 µm line, for instance) suggest the presence of a very young stellar population (a few Myrs). However, this paper deals with dust characteristics and any dust produced/destroyed over such short timescales do not impact on the dust budget. This is why we assume a delayed (no final bursts) star formation history in the SED fitting; (ii) this age range is in agreement with the properties of this sample as they are dwarf, low-metallicity, star-forming galaxies (Maiolino et al. 2008;Madden et al. 2013;Cormier et al. 2019). This is the reason why they are often cited as being high-redshift LBG analogues.   Figure 5 is colour-coded by age/τ main of the stellar population, assuming a delayed SFH. Both panels of Fig. 5 exhibit an age-evolutionary sequence with younger objects at higher sSFR. To the right-hand side of both panels, sM dust presents (or reaches) an apparent maximum at sSFR ∼10 −8 yr −1 with a decline at lower sSFRs. Figure 6 (top and central panels) is colour-coded by IRX. Galaxies have low IRX to the right-hand side, meaning at large sSFRs, but also to the left-hand side, meaning at low sSFRs. In between, objects have larger IRX values.
Combining the information from Figs. 5 and 6, a qualitative scenario seems to emerge. We first focus on the Hi-z LBG sample, which is the core objective of this work: Hi-z LBGs with

Notes.
We note that the parameters for which the uncertainty (e.g. param_err) is ≥ parameter (i.e. param) should be considered as upper limits. (a) Estimated from the relation described in the text.
However, this method being very indirect (see text), these values should only be considered as indicative. low IRX, low sM dust , and low ages are found at high sSFR. We also find Hi-z LBGs with low IRX and even lower sM dust , but large ages are found at low sSFR. In between, we see an apparent maximum observed in IRX and sM dust for the present Hi-z LBG sample at sSFR ∼10 −8 yr −1 that needs to be confirmed by future observational data. Both IRX and sM dust decline, and age increases to a locus at sSFR∼10 −10 -10 −9 yr −1 , where we find the LBGs with upper limits like those discussed in Capak et al. (2015), Ferrara et al. (2017) and Faisst et al. (2017). This seems to define an evolutionary sequence from the right of the DFRD, where the first dust grains form to the left and where they are removed/destroyed. In Fig. 7, we confirm the trends observed in Figs. 5 and 6 by plotting sM dust vs. age main . The decline of sM dust is interpreted as an efficient dust removal (destruction or outflow) with a timescale of ∼500 Myrs. A word of caution, though: if the dust grains are carried away into the wider ISM, they could cool down or warm up from the injection of kinetic energy, and could therefore be more difficult to detect given our wavelength coverage (Fig. 1). We can try to quantify the removal trends: the distribution of Bayesian values t/τ found by the CIGALE SED fitting are in the range 0.5 ≤ age/τ main ≤ 7.5. The galaxies with the maximum sM dust are in the range 0.4-2.9, corresponding to an integrated SFR = 12% and 80%, respectively. The maximum SFR occurs at age/τ main = 1.0, where about 26% of the stars have formed. This means (predictably) that sM dust correlates with SFR and therefore with the dust production by stars and not the ISM.
If we now focus on the Low-zZ sample, our interpretation is that separate parallel sM dust vs. sSFR sequences could be stratified from bottom to top in addition to from right to left. Accounting for the horizontal age sequence for these low-zZ galaxies from the lower panel of Fig. 5, it could be possible that these local galaxies undergo an increase of sM dust as they get older. However, this interpretation is not valid for the Hi-z sample, where we do not see low-IRX objects at inter-  Fig. 6. Hi-z LBGs (filled dots, top) and Low-zZ galaxies (filled boxes, middle) dust formation rate diagram, sM dust vs. sSFR, colour-coded in IRX. Like in Fig. 5, colour-coding the data points helps towards interpreting the sequence as an evolutionary sequence with low-IRX objects at high sSFR and again at low sSFR. In the middle range, IRX is higher. Low-zZ galaxies (filled dots, bottom) dust formation rate diagram, sMdust vs. sSFR, colour-coded in metallicity. mediate sSFR, and there are no high-IRX objects on the lefthand side of the plot, at low sSFR, while we still observe the evolutionary sequence from right to left. Under this hypothesis, we can suggest that the Low-zZ and Hi-z LBG samples are not perfect analogues. It is possible that we observe the effect of another parameter, the metallicity, for the Low-zZ galaxies (Fig. 6, lower panel), because these galaxies formed at much later cosmic time than the Hi-z LBGs. We already knew that the metallicity and IRX of Low-zZ galaxies are correlated (e.g. Heckman et al. 1998;Boissier et al. 2004). In Fig. 8, the linear fit we estimate for the Low-zZ sample is close to, but different, from the starburst galaxy relation found by Heckman et al. (1998). We also added two data points corresponding to galaxies at z = 3.3 by Pannella et al. (2015). Using the updated derived relation between IRX and 12+log 10 (O/H) for the Low-zZ sample (log 10 (IRX) = 1.856 × (12 + log 10 (O/H)), we can estimate the metallicity of the Hi-z LBG sample. The maximum of the normalised distribution of the latter shown in the lower panel of Fig. 8 Fig. 7. sM dust vs. Age main diagrams colour-coded in sSFR for Hi-z LBGs. Low-sSFR objects are the oldest objects and low-sM dust objects are the youngest ones. The decline of sM dust is interpreted as an efficient dust removal (destruction or outflow) with a time-scale of ∼500 Myrs. This diagram makes it possible to obtain an estimate of the dust removal grain in these Hi-z galaxies. Symbols coded as in Fig. 1. From left to right, the vertical dotted lines show 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 90%, and 99% of the integrated SFR, assuming a delayed SFH. highest metallicity of these Low-zZ galaxies, as is suggested by several papers for the DGS sample (e.g. Asano et al. 2013). Asano et al. (2013) define a metallicity as the critical metallicity Z cr , the switching metallicity point at which the increasing rate of dust mass by grain growth exceeds the dust production rate by stars. This Z cr ∼ 0.1−0.5, which approximately corresponds  Fig. 9. Top: the evolution of Hi-z LBGs in the IRX vs.β diagram is an evolutionary sequence, with younger galaxies being bluer and older galaxies redder. The continuous line is the Calzetti dust attenuation law, and the dotted line is the SMC one. The trend observed here explains the location of galaxies not detected by Capak et al. (2015). Bottom: the same objects are colour-coded with sM dust . Hi-z LBGs with low dust masses are found below Calzetti and SMC dust attenuation laws.
to 12+log 10 (O/H)) = 8.0-8.6, in agreement with the Fig. 5 in Ma et al. (2016). In other words, the Low-zZ sample might not be as clean and homogeneous as the Hi-z LBG one, and they are not complete analogues to galaxies in the early universe.
In the upper panel of Fig. 9, we see that Hi-z LBGs with young stellar ages are found to the left of, and above, Meurer's law (the local starburst law linking IRX and β, Meurer et al. 1999) while older and lower sM dust LBGs are found below this law, in the region identified from ALMA-undetected LBGs by Capak et al. (2015). Therefore, the origin of the position of these Hi-z LBGs is consistent with having both low-IRX values and red UV slopes β.

Modelling the dust evolutionary sequence
It is beyond the scope of this paper to describe the models in detail. They will be fully explained in an associated paper (Nanni et al. 2020). In a nutshell, we performed the calculations for the metal evolution using the OMEGA code (Onezone Model for the Evolution of GAlaxies, Côté et al. 2017). We assumed the metal yields for Type II SNe from Kobayashi et al. (2006) computed up to 40 M , and from the FRUITY database for low-mass stars with 1.3 < M/M < 6 evolving through the thermally pulsing AGB phase and developing stellar winds (Cristallo et al. 2011(Cristallo et al. , 2015Piersanti et al. 2013). The yields for population III stars are from Heger & Woosley (2010) and are limited to the mass range 10 < M/M < 30 in OMEGA.
Dust removal from the galaxy through galactic outflow follows a rate proportional to the SFR through the "mass-loading factor": ML × SFR. This assumes that the galactic outflow is generated by the feedback of stars on the gas in the ISM (e.g. Murray et al. 2005). The SFH is a delayed one with τ assumed to be 83 Myrs (in agreement with the average value found from the SED fitting). Two kinds of IMFs are tested: top-heavy IMFs: m −α with α = 1.0, 1.35, 1.5 defined above in this paper and a Chabrier IMF.
We assumed that a fraction of silicates (olivine and pyroxene), iron and carbon grains ejected in the ISM are condensed. In different works (Ventura et al. 2012;Nanni et al. 2013Nanni et al. , 2014, it has been shown that for low-mass stars, the condensation fraction ( f i ) of silicate and carbon dust can be up to 50-60% or slightly more during the super wind phase, when most of the mass is lost by the stars through stellar winds.
In a recent work in which dust dust condensation in SNe remnants is modelled (Marassi et al. 2019), a fair estimate of the mass of dust over the mass of metals is found to be 0.3 < f i < 0.6. In this calculation, however, the effect of the reverse shock on dust destruction is not taken into account. There is a tension on the amount of dust produced by SNe in models and what is actually estimated from observations. Bocchio et al. (2016) argued that only a small fraction (a few per cent) of the dust produced initially condensed in SNe remnants would survive the reverse shock, while Gall et al. (2014) found evidence that large dust grains (>1 µm) in SNe remnants might survive the passage of the reverse shock. Matsuura et al. (2019) also found that dust might be reformed after the passage of a shock.
In this work, we define the condensation fraction at each time-step as the ratio between the number of atoms locked into dust grains over the maximum that can condense according to its stoichiometric formula and to the chemical composition of the gas. On the basis of theoretical calculations, we assume, for TP-AGB stars: f sil = 0.6 for silicates, f iron = 0.01 for solid iron, and f car = 0.5 for carbon dust grains. For SNe, we consider two cases: a very high value of the condensation fraction: f sil = 1.0 for silicates, f iron = 1.0 and f car = 0.5 for carbon dust grains. Note that we only considered pyroxene in silicates (MgSiO3). We also add an additional test case with half of the previously assumed condensation fraction and IMF with α = 1. These assumptions correspond to two very favourable scenarios for dust condensation in SNe remnants.
For each dust species, the destruction rate by SNe is computed as: SNe destr = M dust,i /τ destr , where M dust,i is the dust mass of each dust species, τ destr is the destruction timescale evaluated as: τ destr = M gas /(τ destr R SNe M swept ), where M gas is the gas mass in the galaxy that changes as a function of the outflow and star formation, and destr is the efficiency of the destruction. Here, destr = 0.1 (Hirashita & Aoyama 2019) R SNe is the SNe rate and M swept is the mass of gas swept by the SNe blast wave.
We explore values for M swept = 1000 and 6800 M . The latter case represents the typical value assumed in the literature (e.g. Dwek et al. 2007), while the former assumption considers that much lower values for M swept are also possible (Nozawa et al. 2006). The average τ SFH = 83 Myrs found by fitting individual LBGs is used for the SFH in our modelling. The calculations were performed by normalising the total stellar mass formed in the galaxy after 13 Gyr to 1 M . The best chemical evolution models selected by computing the χ 2 with respect to sM dust , sSFR, and age are presented in Fig. 10. The blue models show that we can explain the structure of the DFRD diagram. However, even though the models are consistent with the data if we account for the uncertainties, we cannot reach the top data points. This may or may not mean that we are still missing some important physics about the dust cycle in these early galaxies. The fact that the models are low in the DFRD could be relevant, given the uncertainties in the theoretical metal yields presently available in the literature. Further observational constraints and deeper theoretical studies need to be devoted to the understanding of the dust cycle in the first galaxies. Fig. 10. Left: dust formation rate diagram: sM dust vs. sSFR in the Hi-z LBGs (red data points). Only models with a delayed star formation history are plotted. Consistently, a delayed star formation history only was used in the SED fitting with CIGALE to estimate sM dust and sSFR. The models shown do not have any grain growth in the ISM. We select the best models (shown in blue) by computing the χ 2 of the models with respect to sM dust , sSFR, and age, and finally use the combined total χ 2 as quality of the fit to select the best models. We also plot the other models in light grey to plot which parameter space the models cover. Right: specific dust mass vs. age of the dominant stellar population. Models with a larger condensation fraction provide better fits, on average. The Hi-z LBGs are young galaxies. The steep decline of sM dust for Hi-z LBGs with respect to Low-zZ galaxies suggest a fast removal (destruction + outflow) of grains in about 0.5 Gyr for Hi-z LBGs that can explain the various detectability of these Hi-z LBGs. Colours have the same meaning as in the left panel.
We adopted a different initial total baryonic mass of the galaxy between M gas = 20 M and 100 M . The larger the mass of the gas, the lower the dust destruction efficiency of SNe. Different values of the mass loading factor (ML) have also been tested, ranging from 0.5 to 0.98 × M gas . For the best models shown in blue in Fig. 10, we computed the mean and standard deviation values: M gas = 57.9 ± 26.7, α = −1.14 ± 0.18, ML = 38.3± 18.2, and M swept = 3951±2900. In the dust removal budget, outflows are dominant ( 70-80%) over SN destruction ( 20-30%).

Conclusions
The major result from this work is a consistent model that explains the (sub-)mm detections (or non-detections) of Hi-z LBGs in the EoR. The non-detections are explained by dust destruction from shocks produced by SNe plus removal in the circum-and intergalactic media by outflows (Fig. 10). Such large-scale outflows of interstellar material are a generic feature of LBGs (Shapley et al. 2003;Pettini et al. 2002). The models presented here do not include grain growth in the ISM of these galaxies since the efficiency of such a process can be highly inhibited at the very low temperature characterising the ISM (Ferrara et al. 2016;Ceccarelli et al. 2018). Our modelling suggests that a Chabrier IMF cannot produce enough dust mass to reach the upper part of the sequence. For these LBGs, a topheavy IMF is preferred.
The possible rise, suggested from data, and decline of sM dust with sSFRs would be the formation and removal (dust destruction and outflow) of the first dust grains in the universe. If a coeval and massive burst of star formation occurred in the early universe (population III), we could predict an earlier bump at higher redshift in the evolution of the cosmic average dust attenuation (Burgarella et al. 2013). However, this bump is likely to be blurred by a non-coeval galaxy formation.
There are still large uncertainties in our modelling due to a lack of knowledge regarding the physical conditions in the early universe. It is crucial to further constrain the models by securing new observational data, and to invest in work on dust modelling to fully explain the dust cycle in Hi-z LBGs (and in other objects).      To build these plots, we use the best-fit modelled SED for each of the Hi-z LBGs. Then, we add the observational noise to the observed data points and we refit the noised best models. Finally, we compare the input and the output parameters. These plots in this figure shows that CIGALE is able to re-derive the input parameters quite nicely and also provide an estimate of the degeneracy between these parameters.