# Dust emission from the bulk of galaxies in the Epoch of Reionisation

L. Ciesla, S. Adscheid, B. Magnelli, M. Boquien, N. Laporte, M. Béthermin,

C. Carvajal, E. Schinnerer, D. Liu

► To cite this version:

L. Ciesla, S. Adscheid, B. Magnelli, M. Boquien, N. Laporte, et al. Dust emission from the bulk of galaxies in the Epoch of Reionisation. Astronomy & Astrophysics - A&A, 2025, 700, pp.A277. 10.1051/0004-6361/202453278 . insu-05263011

HAL Id: insu-05263011

https://insu.hal.science/insu-05263011v1

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dust emission from the bulk of galaxies in the Epoch of Reionisation

L. Ciesla<fnref n="1" />,*, S. Adscheid<fnref n="2" />, B. Magnelli<fnref n="3" />, M. Boquien<fnref n="4" />, N. Laporte<fnref n="1" />, M. Béthermin<fnref n="5" />, C. Carvajal<fnref n="1" />,
E. Schinnerer<fnref n="6" />, and D. Liu<fnref n="7" />,<fnref n="8" />

Received 3 December 2024 / Accepted 14 July 2025

## ABSTRACT

The excess of UV-bright galaxies observed at $z > 10$ has been one of the major surprises of the early JWST observations. Several explanations have been proposed to understand the mild change in space density of the UV-bright galaxies at these high redshifts, among them an evolution of dust attenuation properties in galaxies. However, our view of dust in primordial galaxies is limited towards a few tens of $z \sim 7$ galaxies, pre-selected from UV-optical observations, and are thus not necessarily representative of the bulk of the sources at these redshifts. In this work, we aim to constrain the dust properties of galaxies at $6 < z < 12$ by making the most of the A$^3$COSMOS database in the JADES/GOODS-South field. We stacked ALMA band 6 and 7 observations of 4464 JADES galaxies covered by the A$^3$COSMOS database and used the measurements as constraints to perform UV-to-FIR SED modelling. We obtained tentative signals for the brightest UV galaxies ($M_{UV} < -19$ mag) as well as for the most massive ones ($\log M_*/M_{\odot} > 9$) at $6 < z < 7$, and upper limits for fainter ($M_{UV} > -19$ mag), lower-mass sources ($\log M_*/M_{\odot} < 9$), or ones at higher redshift ($z > 7$). Fitting these $6 < z < 7$ galaxies with ALMA constraints results in lower star formation rates (-0.4 dex) and FUV attenuation (-0.5 mag) for galaxies with $\log M_*/M_{\odot} > 8$, compared to the fit without FIR. We extend the $L_{IR}$ vs $M_{UV}$ relation down to $M_{UV} = -19$ mag and show a tentative breakdown of the relation at fainter UV magnitudes. The positions of the JADES $z \sim 6.5$ sample on the infrared excess (IRX) versus $\beta$ and IRX versus $M_*$ diagrams are consistent with the ones of the ALPINE ($z \sim 5.5$) and REBELS ($z \sim 6.5$) samples, suggesting that the dust composition and content of our mass-selected sample are similar to the ones of these UV-selected galaxies. Extending our analysis of the infrared properties to $z > 7$ galaxies, we find a non-evolution of $\beta$ in the $M_{UV}$ range probed by our sample (-17.24$^{+0.54}_{-0.62}$), and highlight the fact that samples from the literature are not representative of the bulk of galaxy populations at $z > 6$. We confirm a linear relation between $A_V$ and $sFIR^{-1}$ with a flatter slope than previously reported due to the use of ALMA constraints. Our results suggest that rapid and significant dust production has already happened at $z \sim 7$.

**Key words.** galaxies: evolution – galaxies: fundamental parameters – galaxies: high-redshift

## 1. Introduction

In the interstellar medium (ISM) of galaxies, dust plays a crucial role, catalysing the transformation of atomic hydrogen into molecular hydrogen, the raw material for star formation (Wolfire et al. 1995). Dust also enables gas to cool and condense into new stars by absorbing ultraviolet (UV) radiation from nearby young stars (Draine 1978; Dwek 1986; Hollenbach & Tielens 1997). This absorbed energy is then re-emitted in the infrared (IR) range, where dust grain thermal emission dominates the spectral energy distribution (SED) of galaxies, covering wavelengths from approximately 8 to 1000 µm. Because of its important role in the ISM and its tight link with the other components of galaxies, the study of dust emission is paramount for a better understanding of all of the processes at play, especially galaxies’ star formation activity.

At $z > 4$, most of the constraints on star formation in galaxies are obtained via UV emission, which is highly sensitive to dust attenuation and needs to be corrected to obtain a complete view (e.g. Calzetti et al. 2000; Salim & Narayanan 2020). Understanding how dust attenuation varies with redshift and other galaxy properties is thus essential for accurately determining star formation rates (SFRs), and thereby constructing a precise account of galaxy evolution throughout cosmic history. At $z > 4$, two programmes targeted UV-selected galaxies to be observed with the Atacama Large Millimeter/Submillimeter Array (ALMA). ALPINE (ALMA Large Program to INvestigate CII at Early times; Le Fèvre et al. 2020; Béthermin et al. 2020; Faisst et al. 2020) provided 23 far-infrared continuum-detected galaxies at $z \sim 5$. At a slightly higher redshift, the REBELS survey (Reionization Era Bright Emission Line Survey, Bouwens et al. 2022) provided dust continuum ALMA observations of 40 UV-selected bright galaxies at $6.5 < z < 7.7$. At $z > 7$, a few galaxies have been detected

<footnote n="1">¹ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="2">² Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany</footnote>

<footnote n="3">³ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, 06000 Nice, France</footnote>

<footnote n="5">⁵ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote n="6">⁶ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="7">⁷ Max-Planck-Institut für Extraterrestische Physik (MPE), Giessenbachstr., 85748 Garching, Germany</footnote>

<footnote n="8">⁸ Purple Mountain Observatory, Chinese Academy of Sciences, 10 Yuanhua Road, Nanjing 210023, PR China</footnote>

<footnote>* Corresponding author: laure.ciesla@lam.fr</footnote>
in dust emission (e.g. Watson et al. 2015; Laporte et al. 2017, 2019; Bowler et al. 2018; Tamura et al. 2019; Bakx et al. 2020; Sugahara et al. 2021; Endsley et al. 2022; Witstok et al. 2022; Schouws et al. 2022; Algera et al. 2023), providing hints that significant dust reservoirs are already in place at early times. Although they provide an indispensable and detailed view of dust emission at $z \ge 4$, these UV-selected samples are limited in terms of statistics and might only represent the tip of the infrared luminosity distribution at these redshifts. To overcome this, with a different approach, Magnelli et al. (2024) used the A³COSMOS ALMA database (Liu et al. 2019; Adscheid et al. 2024), gathering a wealth of ALMA archival observations to statistically probe dust properties of a mass selected sample at $4 < z < 5$.

The infrared excess (IRX$\equiv L_{IR}/L_{UV}$) versus $\beta$ relation is widely used as a dust attenuation diagnostic (Calzetti et al. 1994, 2000; Meurer et al. 1999). While this relation provides accurate estimates of the attenuation for local starburst galaxies (Meurer et al. 1999), the position of normal star-forming galaxies relative to it is still debated (e.g. Capak et al. 2015; Barisic et al. 2017; Hashimoto et al. 2019; Bakx et al. 2020; Bowler et al. 2022, 2024; Schouws et al. 2022; Boquien et al. 2022; Magnelli et al. 2024). The position of galaxies with respect to the IRX-$\beta$ relation depends on the geometry of dust and stars in the galaxy (Popping et al. 2017; Narayanan et al. 2018; Ferrara et al. 2022; Pallottini et al. 2022; Vijayan et al. 2022) and an effect due to star formation history (SFH) is also expected, since the UV $\beta$ slope is a combination of dust attenuation and the intrinsic UV slope, $\beta_0$ (Boquien et al. 2012). Magnelli et al. (2024) used averaged measurements of ALMA galaxies using A³COSMOS to investigate the dust attenuation properties of star-forming galaxies at $4 < z < 5$. They find that the average properties of their sample are compatible with the relations of Meurer et al. (1999) and Calzetti et al. (2000) for local starbursts. Similar conclusions were reached by Bowler et al. (2024) on stacked properties of the REBELS galaxies at $z \sim 7$. Both results point towards non-evolving grain properties (e.g. size distribution, composition). However, using the ALPINE sample of galaxies at $4.4 < z < 5.5$, Boquien et al. (2022) found a wider variety of attenuation curves among these sources from relations steeper than the Small Magellanic Cloud (SMC) extinction curve to relations shallower than for local starbursts. In their fitting process, the slope of the attenuation curve, $\delta$, was let free, which could explain the larger variety of curves that they find.

The relation between IRX and the stellar mass is another diagnostic of the dust attenuation properties, with the IRX quantifying the dust obscured fraction, while the stellar mass is a proxy of the past star formation activity, and hence of dust production. Recently, there have been claims of a steepening of the relation between IRX and stellar mass with increasing redshift (Fudamoto et al. 2020; Magnelli et al. 2024); that is, a decreasing obscured fraction with redshift. This would imply an evolution of dust geometry relative to new-formed stars. Using the REBELS sample of galaxies at $z \sim 7$, Bowler et al. (2024) also found a decrease in the obscured fraction with redshift but associated with a IRX-stellar mass relation with a relatively flat slope, in contradiction with the results of Fudamoto et al. (2020) and Magnelli et al. (2024), adding to the debate.

The goal of this work is to extend these studies in terms of statistics but also redshift with a statistical analysis using the A³COSMOS (A³GOODSS) archival data, pushing it to its limits by stacking galaxies at $z = 6$ and beyond, observed as part of the JADES survey (Rieke et al. 2023).

Indeed, dust at higher redshifts ($z \ge 9-10$) is also central in the discussion to understand the origin of the excess of UV-bright galaxies observed at $z > 10$ (e.g. Castellano et al. 2022; Naidu et al. 2022; Donnan et al. 2023; Harikane et al. 2024; Arrabal Haro et al. 2023; Curtis-Lake et al. 2023; Robertson et al. 2023; Roberts-Borsani et al. 2024). The evolution of attenuation properties is mentioned as a possible explanation to the abundance of such ‘blue monsters’ (e.g. Ferrara et al. 2023). Constraints on the infrared emission of early galaxies could provide benchmark properties for these studies.

The paper is organised as follows. In Sect. 2 we describe the JADES sample used as a prior for the ALMA stacking and in Sect. 2.2 the SED fitting procedure to obtain the UV magnitudes and stellar mass from which we bin them. The ALMA stacking of the sources using A³COSMOS is detailed in Sect. 3. We compare the properties of the galaxies using and not using the ALMA constraint in the SED modelling in Sect. 4 and study the dust properties of the sample in Sect. 5. Constraints on galaxies at $z > 7$ are provided in Sect. 6. Throughout the paper, we use a Salpeter (1955) initial mass function and WMAP7 cosmology (Komatsu et al. 2011).

# 2. The sample

## 2.1. Sample selection

We used the second release of the JWST Advanced Deep Extra-galactic Survey (JADES, Rieke et al. 2023; Bunker et al. 2024; Eisenstein et al. 2023a; Hainline et al. 2024) of the Hubble Ultra Deep Field covering 68 sq.arcmin presented in Eisenstein et al. (2023b). Although the main part of the study focuses on galaxies with $6 < z < 7$, we selected galaxies with $6 < z < 12$ in the JADES catalogue. We refer to Rieke et al. (2023), Eisenstein et al. (2023b), and Robertson et al. (2024) for more details on the data reduction process, photometry procedures, and redshift estimates. This catalogue combines nine broad and five medium JWST/NIRCam bands from the JADES and JEMS (Williams et al. 2023) surveys, respectively, and existing HST imaging, for a total of 23 photometric bands. The reddest NIR-Cam filter, where the galaxies were selected (F444W), allows us to probe the 0.63 µm and 0.34 µm rest frame at $z = 6$ and $z = 12$, respectively. In this work we have used the photometry catalogue made from KRON apertures. In the broad NIRCam bands, the $5\sigma$ flux depths are between 1.9 and 3.0 nJy, while they range from 1.6 to 3.3 nJy in the medium bands. Photometric redshifts available in the JADES catalogue were computed with the code EAZY (Brammer et al. 2008). The average offset between photometric redshifts and a compilation of spectroscopic redshifts is 0.05 with a $\Delta z/(1+z)$ of 0.024. We used spectroscopic redshifts when available, although this concerns only 0.6% of the sample at $6 < z < 7$, and 0% at $z > 7$. We also checked for spectroscopic data through the Dawn JWST Archive<fnref n="1" /> (Heintz et al. 2025) and found 37 NIRSpec spectra. Out of these, nine are galaxies with $z_{spec} < 5.5$, which we excluded from our analysis.

As in Ciesla et al. (2024), we imposed that each galaxy should be detected with a signal-to-noise (S/N) higher than 3 in at least three JWST filters. This criterion results in a subsample of 8810 galaxies with redshifts between 6 and 12, 6000 of them between 6 and 7. Given the uncertainty on photometric redshift, several studies impose a criterion on the quality of the photometric redshift, based on its PDF (e.g. Cole et al. 2025;

<footnote n="1">¹ https://dawn-cph.github.io/dja/</footnote>
Finkelstein et al. 2023). However, in the following, we consider the full galaxy sample. Indeed, after redoing the analysis of this work imposing the same photo-z quality criteria as Cole et al. (2025), we find no strong impact on our results, a conclusion also reached by Cole et al. (2025).

A particular focus on the $6 < z < 7$ bin is made in the paper, driven by the tentative signal obtained from stacking the ALMA observations that is more constraining at $6 < z < 7$ than at higher redshifts. We extend our analysis at higher redshifts in Sect. 6.

## 2.2. SED modelling with CIGALE and physical properties

We used the SED modelling code CIGALE<fnref n="2" /> (Code Investigating GALaxy Emission, Boquien et al. 2019) on the HST+JWST set of data, probing the UV to NIR rest-frame wavelength range. The fitting set-up of CIGALE was the same as the one tested and used in Ciesla et al. (2024), which adopts a non-parametric SFH specifically tailored for early galaxies. The non-parametric approach usually assumes a prior to handle the relation between two consecutive SFR bins of the SFH (continuity, bursty continuity, etc.). In this work, we used the non-parametric SFH from Ciesla et al. (2024), who used no prior and showed that it provides more flexibility to model the bursty SFH expected in early galaxies. We note that this SFH library results in a wide range of $\beta_0$, the intrinsic UV slope, ranging from -2.78 to 2.65 with a median value of -2.20. We used the stellar population models of Bruzual & Charlot (2003), and a Calzetti et al. (2000) dust attenuation law that is found to be suitable for galaxies at $z > 6$ (Bowler et al. 2022; Markov et al. 2025).

CIGALE assumes an energy balance: the energy absorbed in the UV-optical is re-emitted in the IR. Some studies have shown a spatial offset between the stellar (UV) and dust (IR) emission in some high-redshift sub-millimetre galaxies (e.g. Hodge et al. 2016; Elbaz et al. 2018) resulting in discussions on the validity of UV-to-IR SED modelling based on energy balance. Using the disturbed morphology of M 82 as a test-bed, Seillé et al. (2022) show that the estimates of the SFR and stellar mass for the whole galaxy are found to be consistent with the sum of the same parameters obtained for individual regions. The same conclusion is reached by Li et al. (2024) at $4 < z < 6$ from a resolved analysis of ALMA-CRISTAL galaxies, suggesting that there is no bias due to outshining by young stars on the derived global properties of these sources. Furthermore, using FIRE zoom-in simulations of dusty and high-redshift galaxies, Haskell et al. (2023) demonstrate that when the quality of the fit is acceptable, the fidelity of SED modelling estimates using energy balance is independent from the degree of UV/FIR offset, with a similar quality as was previously reported for local galaxies. Therefore, we are confident in using the ALMA constraint and performed the SED modelling from UV to FIR using the energy balance.

To model the IR emission, we used the dust library of Schreiber et al. (2018), which has two free parameters. The first is the fraction of dust mass locked in the form of polycyclic aromatic hydrocarbons (PAHs), which can vary from 0 to 100%. The second is the dust temperature, which varies from 15 to 60 K. The data used in this work does not allow us to constrain the PAH fraction. Therefore we chose to fix this input parameter to 3%, which corresponds to the estimate of Schreiber et al. (2018) at $z \sim 3$, the highest redshift constraint provided in their study. We note that Shivaei et al. (2024) found a constant value of 3.4% for the PAH fraction at metallicities of ~0.5 $Z_\odot$, which decreases to values lower than 1% at ~0.3 $Z_\odot$. However, due to

the lack of observational constraints, the PAH fraction does not have a significant impact on the results and we kept a fixed value of 3%. Regarding the dust temperature, following the results of Sommovigo et al. (2022) from REBELS galaxies, we set it to 47 K. In Sect. 5.3 we test the impact of this dust temperature assumption on our results.

To estimate the physical properties, the models were fitted to the observations and the likelihoods were computed from the $\chi^2$. The galaxy properties were taken as the likelihood-weighted means and the uncertainties are the likelihood-weighted standard deviations. One important aspect to note is that upper limits were taken into account in the fitting procedure (see Sect. 4.3 of Boquien et al. 2019, for the presentation of the mathematical implementation). In this work, we take $3\sigma$ observational uncertainties as the upper limits.

The parameters used for the fit are provided in Table 1. In their analysis, Ciesla et al. (2024) show that the physical properties, such as stellar mass and SFR, of galaxies with $6 < z < 12$ are well constrained with the HST+JWST set of filters of the JADES survey. As we are using the same combination of filters, we rely on their results. We first applied this CIGALE set-up to the JADES sample using only HST and JWST observations. The resulting normalised distributions of UV magnitude, stellar mass, SFR, and UV slope ($\beta$) are shown in Fig. 1 for the $6 < z < 7$ galaxies. Our sample allows us to probe relatively faint UV magnitudes ($-17.11_{0.63}^{0.55}$ mag), a median stellar mass of $\log M_\star/M_\odot = 7.82_{0.22}^{0.27}$, a median SFR of $\log SFR/M_\odot \text{yr}^{-1} = -0.05_{0.21}^{0.30}$, and relatively blue UV $\beta$ slopes ($-2.24_{0.14}^{0.18}$). The sample that we use for the ALMA stacking was built from this HST+JWST CIGALE run. We selected galaxies with a reduced $\chi^2$ no larger than 4, removing 26 galaxies for which many sources have SEDs that are similar to brown dwarves' SEDs. The second run of CIGALE using the ALMA constraints is presented in Sect. 4.

## 3. ALMA stacking of JADES galaxies

We used the A³GOODSS ALMA archival database, part of the A³COSMOS database (version 20220606; Liu et al. 2019; Adscheid et al. 2024), to perform a stacking analysis on our galaxy sample, choosing ALMA images from bands 6 (~243 GHz) and 7 (~324 GHz). Out of the 5974 galaxies in our $6 < z < 7$ sample, 4464 are covered by one or multiple ALMA pointings (i.e. within an area of a primary beam attenuation > 0.5). Regarding astrometry, the most critical optical versus ALMA astrometric errors are <0.2" (e.g. Franco et al. 2018). JWST astrometry are tied to Gaia and should therefore be better. As the ALMA observations we are using have a point spread function full width at half maximum greater than 0.2", we do not consider astrometry to be a possible issue. This reduces our sample to 4464 galaxies, with redshift, stellar mass, UV luminosity, and UV spectral slope distributions that are indistinguishable (Kolmogorov-Smirnov test) from the ones of their parent sample of 5974 galaxies.

Because of this large number of galaxies, it was impractical to stack them in the uv domain in order to deal with the heterogeneity in angular resolution (ranging from ~0.1" to ~2") of these multiple ALMA pointings (see, e.g. Wang et al. 2022, 2024; Magnelli et al. 2024). To solve this problem, we instead used the fact that at ~1 mm, galaxies are compact in size (~0.1"; e.g. Gómez-Guijarro et al. 2022) and can therefore be considered spatially unresolved at resolutions $\ge 0.2''$, so that $S_{\text{peak}} = S_{\text{total}}$. Therefore, we excluded images with synthesised beam

<footnote n="2">² http://cigale.lam.fr</footnote>
<table><caption>**Table 1.** CIGALE input parameters used to fit the JADES sample.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Definition</th></tr></thead><tbody><tr><td colspan="3">Non-parametric SFH – SFHNLEVELS_FLAT</td></tr><tr><td>age (Myr)</td><td>200, 300, 500, 750</td><td>maximum possible age</td></tr><tr><td>sfr_max</td><td>500 M<sub>☉</sub> yr<sup>-1</sup></td><td>max. value of the SFR</td></tr><tr><td>N<sub>bins</sub></td><td>10</td><td># of time bins of the SFH</td></tr><tr><td>N<sub>SFH</sub></td><td>500</td><td># of SFH per age value</td></tr><tr><td>Last bin size</td><td>10 Myr</td><td></td></tr><tr><td colspan="3">Stellar population – BC03</td></tr><tr><td>metallicity</td><td>0.004, 0.02</td><td></td></tr><tr><td colspan="3">Emission lines – NEBULAR</td></tr><tr><td>logU</td><td>-3, -2</td><td>ionisation parameter</td></tr><tr><td>f<sub>esc</sub></td><td>0, 0.1</td><td>Lyman continuum photons escape fraction</td></tr><tr><td>f<sub>dust</sub></td><td>0, 0.1</td><td>Lyman continuum photons fraction absorbed by dust</td></tr><tr><td colspan="3">Dust attenuation – DUSTATT_MODIFIED_STARBURST</td></tr><tr><td>E(B-V) lines</td><td>[0; 0.5]</td><td>8 values linearly sampled</td></tr><tr><td>E(B-V)s factor</td><td>0.5</td><td>colour excess ratio between continuum & nebular</td></tr><tr><td colspan="3">Dust emission – SCHREIBER2016</td></tr><tr><td>f<sub>PAH</sub></td><td>3%</td><td>fraction of PAH</td></tr><tr><td>T<sub>dust</sub></td><td>47 K</td><td>dust temperature</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Normalised distributions of the physical properties of the 6 < z < 7 JADES selected sample obtained from UV-to-NIR rest-frame SED modelling (without IR). From left to right, and top to bottom: UV magnitude, stellar mass, SFR, and UV slope, β.</figcaption></figure>

sizes <0.2" from our analysis and 'stacked' only primary beam-corrected flux densities measured from the central pixel at each source position. However, we checked that our results and conclusions are not affected by this cut-off when applying a larger one (<0.4"). The measured fluxes were re-normalised from their observed wavelength (within the ALMA bands 6 or 7) to a common 1.1 mm observed wavelength using the star-forming galaxy SED template of Béthermin et al. (2012). The final stacked flux

<table><caption>**Table 2.** Results of the ALMA stacking analysis at 6 < z < 7.</caption><thead><tr><th></th><th>N<sub>galaxies</sub></th><th>N<sub>pointings</sub></th><th>⟨S<sub>v</sub>⟩<br>[μJy]</th></tr></thead><tbody><tr><td>-22 &lt; M<sub>UV</sub> &lt; -20</td><td>45</td><td>198</td><td>5.4 ± 3.7</td></tr><tr><td>-20 &lt; M<sub>UV</sub> &lt; -19</td><td>154</td><td>515</td><td>5.2 ± 3.0</td></tr><tr><td>-19 &lt; M<sub>UV</sub> &lt; -18</td><td>629</td><td>2056</td><td>&lt;4.8</td></tr><tr><td>-18 &lt; M<sub>UV</sub> &lt; -17</td><td>1627</td><td>5943</td><td>&lt;2.5</td></tr><tr><td>8 &lt; logM<sub>★</sub> &lt; 9</td><td>1352</td><td>5341</td><td>&lt;2.5</td></tr><tr><td>9 &lt; logM<sub>★</sub> &lt; 10</td><td>59</td><td>226</td><td>8.7 ± 3.8</td></tr><tr><td>10 &lt; logM<sub>★</sub> &lt; 11</td><td>4</td><td>14</td><td>38.5 ± 29.0</td></tr></tbody></table>

density is the noise-weighted mean of all individual flux densities, S<sub>i</sub>:

$$S_{1.1 \text{ mm,stack}} = \frac{\sum S_{1.1,i} \cdot \sigma_i^{-2}}{\sum \sigma_i^{-2}},$$

with σ<sub>i</sub> the pixel RMS noise at the position of each galaxy.

We first binned galaxies by UV magnitude in four bins: -22 < M<sub>UV</sub> < -20 mag, -20 < M<sub>UV</sub> < -19 mag, -19 < M<sub>UV</sub> < -18 mag, and -18 < M<sub>UV</sub> < -17 mag. Selecting galaxies with M<sub>UV</sub> < -17 mag yields a total of 2455 galaxies. We provide in Table 2 the resulting ALMA measurements for each bin as well as the corresponding number of galaxies and pointings.

We show these measurements in the left panel of Fig. 2 along with the predicted ALMA fluxes obtained from UV-NIR only SED fitting. The ALMA flux of the brightest sources (M<sub>UV</sub> < -20 mag) is slightly overestimated in the absence of IR constraints although still compatible within the errors. At fainter magnitudes, the upper-limit predictions are in agreement with the ALMA measurements. We note that if we assume the ALMA fluxes predicted from the UV-NIR SED modelling, a S/N = 3 detection from the ALMA stacking would require approximately 10 660 and 52 000 sources, in the -19 < M<sub>UV</sub> < -18 mag and
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Left panel: Flux densities obtained by stacking ALMA observations in four UV magnitude bins (black symbols). Circles with downward arrows are $3\sigma$ upper limits. The contours show the spread of predicted fluxes without the ALMA constraint. The grey circles with black borders are the median values of these predicted fluxes in each UV magnitude bin. The blue squares are the stack values of the REBELS sample (Bowler et al. 2024). Right panel: Same as left panel but using three stellar mass bins. The white circle is the expected flux density assuming a simple scaling to the stellar mass of the $9 < \log M_* < 10$ bin.</figcaption></figure>

$-18 < M_{UV} < -17$ mag bins, respectively, or deeper ALMA observations of the 629 and 1627 sources observed in these bins, respectively.

We also stacked the ALMA observations of the JADES/A$^3$COSMOS galaxies by bins of stellar mass ($8 < \log M_*/M_{\odot} < 9$, $9 < \log M_*/M_{\odot} < 10$, and $10 < \log M_*/M_{\odot} < 11$), and show the results in the right panel of Fig. 2. Selecting galaxies with $\log M_*/M_{\odot} > 8$ yields a total of 1415 galaxies. The predictions from UV-NIR SED fitting are consistent with the constraints obtained from the ALMA stacking. Despite being the most populated bin (1352 galaxies), the ALMA measurement of the lowest stellar mass bin ($8 < \log M_*/M_{\odot} < 9$) results in a $3\sigma$ upper limit of 2.5 µJy. More sources, by a factor of 5, would be needed according to SED modelling predictions to reach a S/N of 3 on this measurement. These galaxies being less massive, it is expected that they contain less dust than more massive galaxies. We cannot constrain the dust mass of the galaxies of our sample since we only have one ALMA data point probing the 137–157 µm rest-frame emission. To understand if this non-detection is compatible with the fact that lower-mass galaxies should have a lower dust content, we took the stellar mass as a proxy for galaxy SED normalisation and scaled down the ALMA flux density obtained from galaxies in the intermediate mass bin ($9 < \log M_*/M_{\odot} < 10$). We obtained a prediction of S$_{1.1\,\text{mm}$ of $3.2^{+0.6}_{-0.5}$ µJy, which is in tension with our limit, suggesting that the dust content per unit of stellar mass of these low-mass systems is indeed lower than that observed in more massive objects.

**4. Impact of the ALMA constraint on physical properties of $6 < z < 7$ galaxies**

In this section, we investigate the impact of the ALMA constraints on the inference of the physical properties of the JADES sub-sample. We ran CIGALE a second time using the ALMA stacking measurements as constraints for the IR part of the SED. We attributed to each galaxy the ALMA flux density obtained from stacking according to the derived UV magnitude of the

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Example of SED for a source at $z = 6.78$ in the $9 < \log(M_{\star}^{\text{no-IR}}/M_{\odot}) < 10$ stellar mass bin. The open black squares are the fluxes from HST and JWST and the blue line is the best fit model obtained fitting these data points. The filled black square indicates the ALMA measurement from stacks obtained in the same stellar mass bin, and the red line indicates the best fit model obtained when adding this IR constraint to the SED fitting. Parameters obtained from the two fits are indicated using the same colours. As an indication, the best SEDs obtained from the rest of the galaxies in the same stellar mass bin are shown in light grey.</figcaption></figure>

first CIGALE run. For example, each galaxy with a UV magnitude within $-22 < M_{UV} < -20$ was attributed a S$_{1.1\,\text{mm}$ of $5.4 \pm 3.7$ µJy.

In Fig. 3, we show an example of a $z = 6.78$ galaxy with $M_{UV} = -20.5$ mag and $\log M_{\star}^{\text{no-IR}} = 9.3$. We show on the figure the resulting best fits with and without the ALMA constraint in red and blue, respectively. Clearly, without the ALMA constraint, the IR luminosity is over-predicted by the UV-NIR SED modelling, resulting in a significant $A_V$ of 1.1 mag. Using ALMA stacked measurements as a constraint, the $A_V$ decreases to 0.3 mag. This difference in $A_V$ results in a difference in age,
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison between the physical properties obtained when fitting the UV-to-NIR rest-frame observations only (HST+JWST) and the ones obtained when adding the ALMA constraints. The solid black line is the one-to-one relationship. Results for the stellar mass, SFR, FUV attenuation, and $\beta$ are shown from left to right. The black circles with error bars indicate the median error in each panel. The inset panels show the same figure for the most massive ($\log M_{\star}^{\text{no-IR}}/M_{\odot} > 8$) sources in the sample.</figcaption></figure>

<table><caption>**Table 3.** Stacked physical properties at $6 < z < 7$ derived from the SED modelling using the ALMA constraint.</caption><thead><tr><th></th><th>⟨redshift⟩</th><th>log⟨M<sub>★</sub>⟩<br>[M<sub>☉</sub>]</th><th>⟨SFR⟩<br>[M<sub>☉</sub> yr<sup>-1</sup>]</th><th>⟨β⟩</th><th>⟨M<sub>UV</sub>⟩<br>[mag]</th><th>⟨IRX⟩</th><th>log ⟨L<sub>IR</sub>⟩<br>[L<sub>☉</sub>]</th></tr></thead><tbody><tr><td>-22 &lt; M<sub>UV</sub> &lt; -20</td><td>6.18<sup>+0.51</sup><sub>-0.06</sub></td><td>9.43<sup>+0.0</sup><sub>-0.14</sub></td><td>14.6<sup>+1.3</sup><sub>-6.7</sub></td><td>-2.05<sup>+0.01</sup><sub>-0.03</sub></td><td>-20.47<sup>+0.38</sup><sub>-0.07</sub></td><td>0.05<sup>+0.12</sup><sub>-0.09</sub></td><td>10.53<sup>+0.06</sup><sub>-0.06</sub></td></tr><tr><td>-20 &lt; M<sub>UV</sub> &lt; -19</td><td>6.64<sup>+0.09</sup><sub>-0.42</sub></td><td>8.85<sup>+0.12</sup><sub>-0.51</sub></td><td>12.6<sup>+6.8</sup><sub>-4.6</sub></td><td>-2.13<sup>+0.05</sup><sub>-0.18</sub></td><td>-19.44<sup>+0.19</sup><sub>-0.25</sub></td><td>0.30<sup>+0.16</sup><sub>-0.15</sub></td><td>10.42<sup>+0.25</sup><sub>-0.18</sub></td></tr><tr><td>-19 &lt; M<sub>UV</sub> &lt; -18</td><td>6.68<sup>+0.20</sup><sub>-0.38</sub></td><td>8.23<sup>+0.31</sup><sub>-0.07</sub></td><td>2.0<sup>+0.7</sup><sub>-1.1</sub></td><td>-2.34<sup>+0.07</sup><sub>-0.07</sub></td><td>-18.40<sup>+0.25</sup><sub>-0.25</sub></td><td>&lt;-0.09</td><td>&lt;9.59</td></tr><tr><td>-18 &lt; M<sub>UV</sub> &lt; -17</td><td>6.62<sup>+0.12</sup><sub>-0.36</sub></td><td>7.99<sup>+0.21</sup><sub>-0.22</sub></td><td>0.7<sup>+0.3</sup><sub>-0.2</sub></td><td>-2.28<sup>+0.10</sup><sub>-0.12</sub></td><td>-17.38<sup>+0.19</sup><sub>-0.28</sub></td><td>&lt;-0.10</td><td>&lt;9.17</td></tr><tr><td>8 &lt; log M<sub>★</sub> &lt; 9</td><td>6.65<sup>+0.14</sup><sub>-0.36</sub></td><td>8.27<sup>+0.15</sup><sub>-0.12</sub></td><td>0.4<sup>+0.4</sup><sub>-0.2</sub></td><td>-2.16<sup>+0.22</sup><sub>-0.14</sub></td><td>-17.47<sup>+0.59</sup><sub>-0.69</sub></td><td>&lt;0.11</td><td>&lt;9.20</td></tr><tr><td>9 &lt; log M<sub>★</sub> &lt; 10</td><td>6.19<sup>+0.53</sup><sub>-0.07</sub></td><td>9.27<sup>+0.14</sup><sub>-0.19</sub></td><td>8.8<sup>+8.3</sup><sub>-7.2</sub></td><td>-2.03<sup>+0.08</sup><sub>-0.13</sub></td><td>-20.06<sup>+0.43</sup><sub>-0.48</sub></td><td>0.13<sup>+0.05</sup><sub>-1.82</sub></td><td>10.59<sup>+0.12</sup><sub>-2.30</sub></td></tr><tr><td>10 &lt; log M<sub>★</sub> &lt; 11</td><td>6.31<sup>+0.00</sup><sub>-0.00</sub></td><td>10.86<sup>+0.00</sup><sub>-0.44</sub></td><td>9.2<sup>+0.0</sup><sub>-0.0</sub></td><td>-0.42<sup>+0.00</sup><sub>-0.67</sub></td><td>-18.37<sup>+0.00</sup><sub>-0.96</sub></td><td>1.28<sup>+0.00</sup><sub>-0.24</sub></td><td>10.96<sup>+0.15</sup><sub>-0.00</sub></td></tr></tbody></table>

with an older age obtained from the run using ALMA of 235 Myr compared to 105 Myr when ALMA is not used.

To understand the impact of the IR constraint on the derived physical properties of the galaxies of the whole sample, we compare in Fig. 4 the stellar masses, SFRs, FUV attenuations, and UV continuum $\beta$ slopes obtained with and without the ALMA constraint.

The stellar mass is weakly affected by the addition of the ALMA constraints. However, for galaxies with $8 < \log(M_{\star}^{\text{no-IR}}/M_{\odot}) < 10$ the use of ALMA leads to slightly higher stellar masses, although it is well within the typical error on $M_{\star}$ which is 0.24 dex. Regarding the SFR, the estimates with and without IR measurements are on average in agreement although there are sources, the most massive ones ($\log(M_{\star}^{\text{no-IR}}/M_{\odot}) > 8$, for which the SFR is clearly lower using IR measurements, with a median difference of -0.37 dex. For the FUV attenuation, the average trend shows a lower $A_{\text{FUV}}$ when ALMA is used but only by 5% on the average value, which is negligible compared to the typical errors on this measurement. However, as for the SFR, we find that for the most massive sources, $A_{\text{FUV}}$ is also lower when using IR constraints, with a median value of -0.49 mag. Finally, the UV slope, $\beta$, tends to be bluer when using ALMA constraints. The sources that deviate the most from the one-to-one relationship for each parameter are galaxies with $\log(M_{\star}^{\text{no-IR}}/M_{\odot}) > 8$. For these sources, without IR constraints, the code is using stronger attenuation to fit the observations, while the fit with IR tends towards lower attenuation but older ages instead. This comparison suggests biases when estimating the physical properties of $z > 6$ galaxies without FIR constraints.

## 5. Dust attenuation properties

In this section, we study the dust attenuation properties of $6 < z < 7$ galaxies, taking into account the ALMA constraint obtained from stacking. The physical properties of this sample are derived from their distribution (median, 25th, and 75th percentiles), using the same weights as were used to estimate the ALMA fluxes, and are provided in Table 3. For bins containing both galaxies with high S/N values and galaxies with upper limits on these physical properties, we derived upper limits of the median value using a parametric approach; that is, by fitting a normal distribution to the data and determining the $z$ score at a 84% confidence level. The upper limit of the bin wavs then derived by adding to the mean of the distribution the $z$ score times the standard deviation of the distribution.

Several pre-JWST and JWST studies have reported a relation between the UV $\beta$ slope and UV magnitude such as UV slopes shift towards bluer $\beta$ at low UV luminosities (e.g. Finkelstein et al. 2012; Bouwens et al. 2014; Bhatawdekar & Conselice 2021; Topping et al. 2024). Although still debated at $z > 8$, this relation would indicate a trend between dust attenuation, which is the dominant component affecting the UV slope, and UV luminosity that may hint at a relationship between luminosity and metals in early galaxies (Cullen et al. 2023; Topping et al. 2024). In Fig. 5, we compare the UV properties ($\beta$ and $M_{\text{UV}}$) of our sample with sources from the literature at the same redshift.

Combining our measurement and the stack estimates from REBELS (Bowler et al. 2024; Austin et al. 2024; Roberts-Borsani et al. 2024), we obtain the following
<figure><img src="image_6.png" /><figcaption>Fig. 5. β as a function of UV magnitude. Big circles with black edges are weighted medians in four UV magnitude bins. Samples from the literature are indicated by symbols with different shades of blue (Bowler et al. 2024; Cullen et al. 2023; Austin et al. 2023; Roberts-Borsani et al. 2024; Mitsuhashi et al. 2024), white edges showing individual galaxies and black edges stacked or median values. The relation from Topping et al. (2024) is shown with the solid line.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Infrared luminosity as a function of UV magnitude. The contours show the position of the whole sample with $L_{IR}$ obtained from the fit using the ALMA constraint. Big circles with black borders are the weighted median values in four UV magnitude bins. Stacked or median values for samples from the literature are shown in various shades of blue symbols (Bowler et al. 2024; Khusanova et al. 2021).</figcaption></figure>

relation from a linear regression fit:

$$ \beta = -0.08 \pm 0.03 \times M_{UV} - 3.70 \pm 0.11. \quad (1) $$

The resulting slope is in agreement with the relations from Topping et al. (2024), who found a $d\beta/dM_{UV}$ slope of $-0.11 \pm 0.02$ and $-0.12 \pm 0.02$ at $z \sim 5.9$ and $z \sim 7.3$, respectively, as well as the slopes obtained at $z \sim 5$ by Rogers et al. (2014) and Bouwens et al. (2014) within errors ($-0.12 \pm 0.02$ and $-0.14 \pm 0.02$, respectively).

The ALMA stacking provides a constraint on the IR luminosity of the galaxies in our sample. In Fig. 6, we show the $L_{IR}$ as a function of the UV magnitude. The JADES sample combined with A$^3$COSMOS observations allows us to probe fainter magnitudes and extend the $L_{IR}$ versus $M_{UV}$ relation down to $M_{UV} \sim -17$ mag. Our sample has $\log L_{IR}/L_{\odot}$ spanning from 7.4 to 12.3 with a median value of 9.3. For comparison, we added the stacked measurements of the REBELS sample as well as the ones from the ALPINE sample at $z \sim 5.5$ from Khusanova et al. (2021). Both these samples probe higher UV magnitude ranges and higher IR luminosities. A fit to the REBELS+ALPINE stacked signals combined with our two highest UV magnitude bins yields the following relation:

$$ \log (L_{IR}/L_{\odot}) = -0.31 \pm 0.05 \times M_{UV} + 4.31 \pm 0.42. \quad (2) $$

Extended to fainter UV magnitude, this relation predicts higher $L_{IR}$ at $M_{UV} > -19$ than our estimates using ALMA. At $M_{UV} = -18.4$ mag and $-17.4$ mag, the relation derived from higher UV magnitudes overestimates the $L_{IR}$ by 0.45 dex and 0.56 dex, respectively. This is a possible hint at a breakdown in the $L_{IR}$ versus $M_{UV}$ relation at low UV magnitudes that could be due to lower dust content. For the two UV-faintest bins, which have a median $M_{UV}$ of $-18.4$ and $-17.4$, the average stellar mass is $\log M_{\star}/M_{\odot} = 8.2$ and 8.0, respectively. These two UV bins also have low dust attenuation: 0.31 and 0.33, respectively. Therefore, these low UV magnitudes are dominated by low-mass, low-attenuation systems, leading to a turnover of $L_{IR}$, which has already been observed at lower redshifts ($z \sim 4$) by Bernhard et al. (2014).

## 5.1. IRX versus β

We place our measurements, obtained from ALMA stacking in UV magnitude bins, on the IRX-β diagram in Fig. 7. Our averaged data points are close to the stacked measurements from REBELS, especially the brightest UV magnitude bin ($-22 < M_{UV} < -20$), which corresponds to the UV magnitude range probed by the REBELS galaxies. They disfavour the expected relation from a SMC extinction curve but are compatible with a Calzetti et al. (2000) relations used in Bowler et al. (2024) with $\beta_0 = -2.3$. Although our SFH library has a wide range of $\beta_0$, galaxies of our sample favour blue $\beta_0$ with a median value of $-2.5$, which is consistent with the assumed $\beta_0$ of Reddy et al. (2018). Therefore, including the effect of dust, our sample has relatively blue $\beta$ with a median value of $-2.25$.

The galaxies in our sample are also compatible with both the starburst and SMC relations of Reddy et al. (2018), who used them with a different $\beta_0$ than Bowler et al. (2024) ($\beta_0 = -2.6$). These relations were found by Boquien et al. (2022) to be suitable for the ALPINE sample. In our process, we tested letting the $\delta$ slope of the attenuation curve free, as in Boquien et al. (2022). However, the improvement of the fit, using a $\Delta$BIC (Bayesian information criterion) analysis, proved to be statistically significant for only 3.7% of the sample. Therefore, we kept $\delta$ fixed to 0 (no modification of the attenuation curve). Although they seem to favour the Reddy et al. (2018) relations, most estimates at low IRX were obtained from the ALMA upper limits on the IR and are thus upper limits on the IRX, as is indicated by the two faintest UV magnitude bins. However, at higher IRX, our estimates are compatible with the relations of both Calzetti et al. (2000) with $\beta_0 = -2.3$ and Reddy et al. (2018) (starbursts and SMC, using $\beta_0 = -2.6$). These findings are consistent with the results that grain size and dust composition at very high redshift are similar to local starbursts and SMC, indicating no fundamental evolution. Our results confirm that at $z > 6$, there is no clear sign of a breakdown of the Calzetti et al. (2000) relation.

## 5.2. IRX versus stellar mass

We place our galaxies on the IRX-$M_{\star}$ plane in Fig. 8. Specifically in this figure, we use the results of the SED fitting using
<figure><img src="image_8.png" /></figure>

**Fig. 7.** IRX as a function of $\beta$. Coloured circles with black borders are weighted medians derived in four UV magnitude bins from our sample. The blue squares are the stacked values of Bowler et al. (2024) for the REBELS sample. The light blue downward triangles are the ALPINE $z > 5$ from Boquien et al. (2022), while the dark blue triangles are ALPINE stacks at $z > 5$ values from Fudamoto et al. (2020). The solid thick black line is the relation assuming a Calzetti et al. (2000) attenuation law, while the dotted line is the one assuming a SMC-like extinction curve. The thin black lines are the corresponding relations from Reddy et al. (2018).

the ALMA constraints obtained from stacking in bins of stellar mass ($8 < \log M_{\star}/M_{\odot} < 9$, $9 < \log M_{\star}/M_{\odot} < 10$, and $10 < \log M_{\star}/M_{\odot} < 11$). Along with our data points, we add the relations from Pannella et al. (2015) at $z \sim 1-3$, Magnelli et al. (2024) at $4 < z < 5$, Fudamoto et al. (2020) at $z \sim 5.5$, and Bowler et al. (2024) at $z \sim 7$ (see also Álvarez-Márquez et al. 2019; Koprowski et al. 2018; Heinis et al. 2014; Bouwens et al. 2016).

Our low-mass bin upper limit on IRX does not provide any useful constraint on the position of these galaxies on the IRX-$M_{\star}$ plane. The intermediate-mass median value is compatible with the relations of Bowler et al. (2024) and Magnelli et al. (2024), and marginally with the relation of Fudamoto et al. (2020). There seems to be a better agreement with the relation of Bowler et al. (2024), despite the large scatter indicated by the error bars (25th and 75th percentile). The high-mass bin data point is in agreement with the relations of both Magnelli et al. (2024) and Fudamoto et al. (2020), and marginally with the Bowler et al. (2024) relation, but here again the large scatter prevents us from discriminating a relation. For comparison, we show the brightest UV magnitude bin, $-22 < M_{UV} < -20$ (red star), for consistency with the REBELS sample and find that it falls on the relation from Bowler et al. (2024), and is thus in agreement with their result. Finally, we show the relation obtained by Pannella et al. (2015) at lower redshift ($z < 4$), indicating that there is a clear decrease in IRX at fixed stellar mass with increasing redshift.

### 5.3. Impact of dust temperature assumption

At $6 < z < 7$, ALMA 1.1 mm probes the 137–157 µm rest frame, which provides a good constraint on the peak of the IR SED. However, with only one flux measurement in the IR, the $L_{IR}$ strongly depends on the assumption made about the IR SED shape. In the dust emission model of Schreiber et al. (2018), this shape is handled through two free parameters: the PAH fraction and the dust temperature. The PAH fraction has little

<figure><img src="image_9.png" /></figure>

**Fig. 8.** IRX as a function of stellar mass. Circles with black edges are weighted medians colour-coded according to the stellar mass bin used for stacking. The red star is the position of the $-22 < M_{UV} < -20$ bin for direct comparison with the stacked values of Bowler et al. (2024) for the REBELS sample (blue squares). The light blue downward triangles are the ALPINE $z > 5$ stacks values from Boquien et al. (2022). The solid light purple line is the relation from Magnelli et al. (2024), the dashed line is from Fudamoto et al. (2020), the dash-dotted line is the relation from Bowler et al. (2024), and the dotted light purple line the relation from Pannella et al. (2015) at lower redshifts. The darker the color of these relations, the higher the redshifts of the galaxies used to derive them.

effect on the $L_{IR}$, but the dust temperature, which we originally fixed at 47 K Sommovigo et al. (2022), is a key assumption that we can test in this section. Indeed, given the large masses and more evolved ISM properties found for REBELS galaxies (e.g. Algera et al. 2024; Rowland et al. 2025), dust temperature of REBELS galaxies are not necessarily appropriate for the broader JADES sample that contains lower-mass and metal-poor galaxies. Indeed, observations suggest that such galaxies could have higher dust temperatures, as high as 80 K (e.g. Laporte et al. 2017; Bakx et al. 2020). (Algera et al. 2024) also suggest that the low dust temperatures in REBELS galaxies could result from inefficient heating in a massive, dust-rich ISM with less bursty star formation, a scenario that may not apply to the lower-mass JADES galaxies.

We thus ran CIGALE several times, each time fixing the dust temperature to a different value. We varied $T_{dust}$ from the highest value allowed by Schreiber et al. (2018) library, 60 K, down to a minimum value of 30 K, which corresponds to the lowest temperature of REBELS galaxies found by Sommovigo et al. (2022) (REBELS32, 39 K minus the lower error of 9 K). We note that Mitsuhashi et al. (2024) found temperatures as low as $26_{-6.5}^{+14.5}$ K in the SERENADE sample, which is typical of the dust temperature at $z < 1$ (e.g. Schreiber et al. 2018), and recently Algera et al. (2024) obtained a dust temperature of $32_{-6}^{+9}$ K for REBELS25, a galaxy at redshift 7.31. Our lowest assumed temperature of 30 K is thus also consistent with the temperatures derived from the SERENADE sample and REBELS25.

We looked at the distribution of $L_{IR}$ using either the UV magnitude bins for the ALMA stacking (Fig. 9) or stellar mass bins (Fig. 10). We show the results of the SED modelling run for each dust temperature assumption (30, 33, 36, 39, 42, 45, 47, and 60 K). As was expected, the effect of decreasing dust temperature assumption is to shift the $L_{IR}$ distributions towards lower values. However, temperatures of ~36–47 K result in similar $L_{IR}$. For galaxies with $\log M_{\star} < 9$ and $M_{UV} < -19$, assuming a
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** $L_{IR}$ distributions of galaxies in the four UV magnitude bins considered in this work, assuming different dust temperatures.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** $L_{IR}$ distributions of galaxies in the three stellar mass bins considered in this work, assuming different dust temperatures.</figcaption></figure>

temperature of 60 K pushes the estimate of the $L_{IR}$ to high values ($\log L_{IR} = 11$.). Despite these uncertainties on the IR SED, the distributions shown in Figs. 9 and 10 provide average constraints on the range of $L_{IR}$ of the bulk of galaxies at $6 < z < 7$.

The impact of $T_{dust}$ on the different physical properties studied in Sect. 5 is presented in Fig. 11. This figure shows on the left panel the IRX-$M_*$ plane using ALMA observations stacked in stellar mass bins, assuming the same set of dust temperatures than in Fig. 9. As was expected, the IRX parameter is affected. Very small variations are also found in $M_*$ due to the fact that CIGALE performs SED modelling using an energy balance over the whole wavelength range, another indication of the error on the measurement that we can have on this parameter. The intermediate-mass bin galaxies are still compatible with the Bowler et al. (2024) relation as long as the temperature is lower than 47 K. A high dust temperature would bring the observations closer to the Pannella et al. (2015) relation at lower redshifts. A

strong effect is, however, observed for the most massive bin, but this bin contains only four galaxies and is thus more subject to statistical variations.

In the right panel of Fig. 11, we show the impact of the dust temperature assumption on the IRX-$\beta$ diagram. The conclusions of Sect. 5 are still valid. Our data points are compatible with a Calzetti et al. (2000) relation using $\beta_0 = -2.3$ or the Reddy et al. (2018) relation for starbursts (SB) and the SMC with $\beta_0 = -2.6$. However, for high temperatures, the observations would prefer the SB relation of Reddy et al. (2018).

## 6. Infrared constraints at $z > 7$

We applied the ALMA stacking analysis described in Sect. 3 to JADES galaxies at $z > 7$. The measurements result only in upper limits; however, they provide first constraints on the IR emission of the bulk of the population at these redshifts. We provide these upper limits in Table 4 in UV magnitude and stellar mass bins. We then perform the same SED modelling procedure on the $z > 7$ sources as for the $6 < z < 7$ redshift bin, assuming the same dust temperature of $T_{dust} = 47$ K.

### 6.1. Evolution of $\beta$

We examine $\beta$ as a function of redshift in Fig. 12. Although the galaxies in our sample are slightly redder than the averaged measurements of Roberts-Borsani et al. (2024) at the same redshifts, they do not show extreme values and lie within the range spanned by individual galaxies in the literature. As is discussed in Franco et al. (2024) (see also Topping et al. 2024), there is a caveat in using SED modelling to measure $\beta$ as models are limited to a given $\beta$ value towards the blue side of the distribution. In our case this value is defined by the lowest intrinsic value allowed by our model: $\beta_0 = -2.77$. On the other hand, despite being more flexible, a direct measurement of $\beta$ from the data can also introduce systematics due to the different rest-frame wavelength range used to measure $\beta$.

To understand the extent to which our results could be affected by the lower limit on $\beta_0$ induced by SED modelling, we show in the inset panel of Fig. 12 the distributions of $\beta_0$ allowed by our models (dark blue), $\beta_0$ obtained from the SED fitting run (cyan), and the $\beta$ for our galaxies (purple). For the total sample, the three redshift bins, only 15% of the galaxies select the lowest possible value of $\beta_0$. A Kolmogorov-Smirnov test indicates that these galaxies have similar redshift and $A_{FUV}$ distributions to the total sample, but different $M_*$ and age distributions, resulting in lower values for both properties. The fact that galaxies with the lowest $\beta_0$ are not at a particular redshift means that these galaxies would not necessarily drive a trend that we would be missing between $\beta$ and redshift. Furthermore, although we would expect to miss galaxies with the bluest $\beta$ in the highest redshift bin from the expected trend, only 5% of these $z = 9-12$ galaxies have the lowest $\beta_0$ value compared to 18 and 13% at $z = 6-7$ and $z = 7-9$, respectively. Therefore, despite the fact that we may miss a fraction of galaxies with extreme blue colours, our results seem to favour a relatively flat relation between $\beta$ and redshift.

The $\beta$ distribution of our sample points to the absence of a trend between $\beta$ and redshift. A trend would be expected with galaxies having bluer colours with increasing redshift, since they become progressively metal- and dust-poor (e.g. Topping et al. 2024; Cullen et al. 2023). Recently, Cullen et al. (2023) found extremely blue galaxies at $z \sim 12$ with a median $\beta$ of $-2.5$, while Topping et al. (2024) obtained a median $\beta$ of $-2.48_{-0.12}^{+0.14}$ for a sample drawn from JADES. At the same redshift, we find a
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Left panel: IRX as a function of stellar mass (symbols and lines following Fig. 8, except the upper limits that are shown here with coloured downward triangles). Filled coloured circles are the stacked measurements obtained for each temperature assumption. Right panel: IRX-$\beta$ diagram (symbols and lines following Fig. 7). Galaxies are colour-coded as in the left panel.</figcaption></figure>

<table><caption>**Table 4.** Results of the ALMA stacking analysis at $7 < z < 9$ and $9 < z < 12$.</caption><thead><tr><th></th><th>N<sub>galaxies</sub></th><th>N<sub>pointings</sub></th><th>⟨S<sub>v</sub>⟩<br>[μJy]</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center;">7 &lt; z &lt; 9</td></tr><tr><td>-22 &lt; M<sub>UV</sub> &lt; -20</td><td>17</td><td>55</td><td>&lt;26.4</td></tr><tr><td>-20 &lt; M<sub>UV</sub> &lt; -19</td><td>90</td><td>267</td><td>&lt;13.6</td></tr><tr><td>-19 &lt; M<sub>UV</sub> &lt; -18</td><td>302</td><td>1058</td><td>&lt;6.2</td></tr><tr><td>-18 &lt; M<sub>UV</sub> &lt; -17</td><td>893</td><td>3337</td><td>&lt;3.2</td></tr><tr><td>8 &lt; log M<sub>★</sub> &lt; 9</td><td>863</td><td>3527</td><td>&lt;2.9</td></tr><tr><td>9 &lt; log M<sub>★</sub> &lt; 10</td><td>35</td><td>130</td><td>&lt;14.8</td></tr><tr><td>10 &lt; log M<sub>★</sub> &lt; 11</td><td>4</td><td>18</td><td>&lt;72.0</td></tr><tr><td colspan="4" style="text-align: center;">9 &lt; z &lt; 12</td></tr><tr><td>-22 &lt; M<sub>UV</sub> &lt; -20</td><td>7</td><td>18</td><td>&lt;62.3</td></tr><tr><td>-20 &lt; M<sub>UV</sub> &lt; -19</td><td>40</td><td>170</td><td>&lt;11.0</td></tr><tr><td>-19 &lt; M<sub>UV</sub> &lt; -18</td><td>159</td><td>479</td><td>&lt;9.4</td></tr><tr><td>-18 &lt; M<sub>UV</sub> &lt; -17</td><td>403</td><td>1447</td><td>&lt;4.5</td></tr><tr><td>8 &lt; log M<sub>★</sub> &lt; 9</td><td>452</td><td>1734</td><td>&lt;3.8</td></tr><tr><td>9 &lt; log M<sub>★</sub> &lt; 10</td><td>33</td><td>127</td><td>&lt;14.4</td></tr><tr><td>10 &lt; log M<sub>★</sub> &lt; 11</td><td>2</td><td>4</td><td>&lt;188.0</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** $\beta$ as a function of redshift. Grey contours show the full JADES/A$^3$COSMOS sample, while purple circles show the median position in several redshift bins. Different samples from the literature are shown in shades of blue (Laporte et al. 2017, 2019; Tamura et al. 2019; Inami et al. 2022; Mitsuhashi et al. 2024; Schouws et al. 2025; Topping et al. 2024). The inset panel shows the distribution of the intrinsic UV slope, $\beta_0$, allowed by the models (dark blue), as well as $\beta_0$ and $\beta$ obtained from the SED fitting.</figcaption></figure>

median value of $-2.14 \pm 0.21$, which is slightly redder. Although both the studies of Cullen et al. (2023) and Topping et al. (2024) found a decrease in $\beta$ with increasing redshift, the latter mention that the rate of this evolution is slowest at low UV luminosities with minor differences between galaxies with $M_{UV} = -18$ at $z = 5-9$ and $z > 9$. The bulk of galaxies in our sample are UV-faint (see Fig. 6) and the absence of strong evolution that we observe for the whole sample is consistent with the finding from Topping et al. (2024). Interestingly, we note that $\beta$ is slightly redder in the highest redshift bin ($z > 9.5$), possibly confirming the results of Saxena et al. (2024) who found such a red $\beta$ using a spectroscopic sample at the same redshift. They interpret this redder value as resulting from a rapid build-up of dust content in the very early Universe or a significant contribution

from the nebular continuum emission to the observed UV spectra, depending on gas temperatures and densities.

We now try to put constraints on the evolution of $L_{IR}$ with redshift. At these high redshifts, the temperature of the cosmic microwave background (CMB) increases rapidly, reaching 19 K, 22 K, 27 K, and 30 K at $z = 6, 7, 9$, and 10, respectively (da Cunha et al. 2013). If not taken into account in the SED fitting procedure, the impact of the CMB results in an under-estimation of the $L_{IR}$ that depends on the dust temperature and redshift. da Cunha et al. (2013) demonstrated that for $T_{dust} = 40$ K the $L_{IR}$ is under-estimated by $-0.018, -0.019, -0.045$, and $-0.070$ dex at $z = 6, 7, 9$, and 10, respectively. We thus corrected our $L_{IR}$ measurement accordingly in Fig. 13, where we show them as a function of redshift. We note that at $6 < z < 7$ the correction is negligible and was not taken into account in Fig. 6.
<figure><img src="image_15.png" /><figcaption>Fig. 13. IR luminosity as a function of redshift for our sample (black contours, downward arrows indicate upper limits) and different samples from the literature (shades of blue, Laporte et al. 2017, 2019; Tamura et al. 2019; Inami et al. 2022; Mitsuhashi et al. 2024; Schouws et al. 2025). Symbols with downward arrows are upper limits. Circles with black borders are our fiducial UV magnitude bins.</figcaption></figure>

For comparison, we added to Fig. 13 individual measurements from the literature: the REBELS sample (Inami et al. 2022), the SERENADE sample (Mitsuhashi et al. 2024), and the individual sources of Laporte et al. (2017, 2019), Tamura et al. (2019), and Schouws et al. (2025). Fig. 13 highlights the fact that the bulk of the population is still out of reach and that the samples on which models are basing their assumptions have extreme IR properties.

## 6.2. Dust formation in early galaxies

Combining V-band attenuation and sSFRs measured from JWST/NIRSpec spectra of 631 galaxies at $3 < z < 14$, Langeroodi et al. (2024) found that dust must have appeared rapidly and that attenuation follows star formation on timescales shorter than ~30 Myr, favouring the scenario in which dust is produced by supernovae (SNe). They found a linear correlation between $A_V$ and the inverse sSFR (sSFR$^{-1}$, a proxy for the age of the galaxies) for galaxies with redshifts from 3 to 14, and $\beta$ from -3 to 0. Through a model, they show that the slope of this relation, $\gamma$, depends on several dust properties such as the efficiency of SN in producing dust, the fraction of dust destroyed by SN reverse and forward shocks, the fraction of dust removed from galaxy by feedback-driven outflows, dust geometry, and dust grain properties. They obtained $\gamma = 0.67 \pm 0.2 \text{ cm}^2 \text{ g}^{-1}$, 30 times lower than the value expected if the entire SN dust yield were preserved in the ISM and had Milky-Way-like grain properties. Here, we use our photometry-based estimates of $A_V$, $M_\star$, and SFR to compare the results obtained with the two different methods.

In Fig. 14, we show the V-band attenuation as a function of the inverse sSFR for the whole sample ($6 < z < 12$) binned by the UV slope, $\beta$, similarly to Langeroodi et al. (2024). Consistently with the results of the previous sections, we confirm a rapid dust production in early galaxies with non-negligible dust attenuation at sSFR$^{-1} > 30$ Myr. Despite using different methods (spectroscopy vs photometry-based SED modelling), our data points also confirm the linear relation between $A_V$ and sSFR$^{-1}$

found by Langeroodi et al. (2024), implying that dust production is driven by star formation activity and stellar mass build-up.

However, the intercept and slope of the relation obtained by our sample are lower and flatter, respectively:

$$ \log A_V = (0.61 \pm 0.05) \times \log \text{sSFR}^{-1} - (0.30 \pm 0.03). \quad (3) $$

First we checked if the flatter slope that we observe is due to the SED modelling by performing a mock analysis (e.g. Buat et al. 2015; Ciesla et al. 2017; Boquien et al. 2019). Basically, we ran CIGALE once on the JADES galaxies. We took the best fit SED, convolved it in the JADES+ALMA set of filters, and added noise to the photometry. We again ran CIGALE on the mock photometry. In this way we can compare estimated parameters to the true ones used to build the mock photometry. From this test, we find that sSFR$^{-1}$ is recovered within the error bars with a systematic overestimate (60%) and that $A_V$ is systematically underestimated by -0.4 mag. However, applying corrections derived from this mock analysis does not yield a slope consistent with the observations from Langeroodi et al. (2024) ($0.58 \pm 0.12$). Although we looked at the impact of dust temperature assumption in Sect. 5.3 and concluded that it does not affect strongly our conclusions for the $6 < z < 7$ redshift bin, studies predict that dust temperature keeps increasing with redshift (e.g. Sommovigo et al. 2022; Mitsuhashi et al. 2024). The temperature of 47 K that we assumed in this work is consistent with the results and predictions of Sommovigo et al. (2022) and Mitsuhashi et al. (2024). However, it can reach up to 60-65 K at $z > 7$. We thus tested the impact of increasing the assumption on $T_{\text{dust}}$ for galaxies in the $9 < z < 12$ to 60 K (the maximum value available in Schreiber et al. (2018) library) and looked at the impact of it on the slope of the $A_V$-sSFR$^{-1}$ relation. We obtain a slope of $0.63 \pm 0.05$ and conclude that the assumption on $T_{\text{dust}}$ is not causing the difference observed between our study and the results of Langeroodi et al. (2024).

We then tested if the use of ALMA constraints is causing the discrepancy. Indeed, the third panel Fig. 4 shows that $A_{\text{FUV}}$ is lower when ALMA information is used in the SED modelling. However, using $A_V$, $M_\star$, and SFR obtained from the UV-to-NIR SED modelling run (Fig. 4, right panel), our binned points are in perfect agreement with Langeroodi et al. (2024). We obtain a slope of $0.91 \pm 0.17$ if we discard the reddest $\beta$ bin, which seems to be discrepant and has much lower statistics than the other bins. Therefore, we conclude that the addition of ALMA constraints is at the origin of the different slopes between the spectroscopic sample of Langeroodi et al. (2024) and our photometric one. Our flatter slope increases even more the discrepancy with $\gamma$ obtained in Milky-Way conditions, but also points to a relation between sSFR$^{-1}$ and the SFR surface density ($\Sigma_{\text{sSFR}}$) for which we do not have access for our sample.

## 7. Conclusions

To place constraints on dust emission in galaxies in the Epoch of Reionisation, we have used the $A^3\text{COSMOS}$ ($A^3\text{GOODSS}$) database to stack all JADES galaxies at $6 < z < 12$ with available ALMA coverage in bands 6 and 7. Focusing on the $6 < z < 7$ redshift bin, our results are summarised as follows:

- ALMA stacks result in upper limits for galaxies with $M_{\text{UV}} > 19$ and $\log M_\star < 9$, and tentative detections at brighter UV magnitudes and higher stellar masses at $6 < z < 7$.

- Upper limits for UV-faint galaxies are in agreement with predicted ALMA fluxes from UV-to-NIR SED fitting, while ALMA fluxes are lower than predicted for UV-bright $M_{\text{UV}} < 19$.
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Left: $A_V$ attenuation as a function of $sSFR^{-1}$. Galaxies with $6 < z < 12$ are binned with the UV $\beta$ slope (circles). The solid black line is the linear fit of these binned points. Predictions of the linear $A_V-sSFR^{-1}$ relation at a fixed SFR surface density from Langeroodi et al. (2024) are indicated with dashed coloured lines. The stacked points obtained by Langeroodi et al. (2024) are shown as cyan hexagons. Right: Same as in the left panel but using the $A_V$ and $sSFR^{-1}$ obtained with the UV-NIR SED fit; that is, without ALMA constraints.</figcaption></figure>

–20 galaxies. For galaxies binned with stellar mass, the predicted fluxes without ALMA constraints and the ALMA fluxes are in agreement. However, for galaxies with $\log M_{\star} < 9$, there is a tension between the ALMA-derived upper limit ($3\sigma$) and the flux expected from scaling down the ALMA flux of the upper stellar mass bin. This is a hint at a possible lower dust-to-stellar-mass ratio at these low stellar masses.

- The addition of the ALMA constraints for the SED modelling results in a lower SFR (0.4 dex), lower $A_{FUV}$ (0.5 dex) for galaxies with $\log M_{\star} > 8$, as well as slightly higher stellar masses and lower $\beta$ for the same sources, implying an older stellar population (25% older on average).

- Using infrared constraints for the SED modelling, our $\beta$ versus $M_{UV}$ is consistent with measurements from the literature, and in perfect agreement with the relation from Topping et al. (2024).

- Our measurements extend the $L_{IR} - M_{UV}$ relation down to $M_{UV} \sim -19$ and hint at a possible breakdown of the relation at fainter UV magnitudes.

- When the ALMA constraints are used in the fitting, galaxies are consistent with the literature (REBELS and ALPINE) when placed on the IRX against $\beta$, indicating a similar dust composition and content as in these UV-selected samples. Our stack measurements computed in bins of stellar masses do not provide strong constraints on the IRX-$M_{\star}$ relation. However, the position of our UV-bright bin ($M_{UV} < -20$) is in good agreement with the relation from Bowler et al. (2024). Nevertheless, we find that a given stellar mass, the IRX decreases with increasing redshift.

We tested the impact of dust temperature assumption on our results and found that temperature assumptions between 36 K and 47 K yield similar $L_{IR}$. However, pushing to a higher dust temperature (60 K) results in a higher $L_{IR}$, as was expected. The conclusions obtained with a dust temperature of 47 K are still valid when the assumed dust temperature is lower. In the IRX-$\beta$ diagram, our points are compatible with a Calzetti et al. (2000) relation using $\beta_0 = -2.3$ or the Reddy et al. (2018) relation for starburst and SMC ($\beta_0 = -2.6$). However, a high dust temperature (60 K) would bring our observations closer to the Reddy et al. (2018) relation for starbursts.

Combining the whole sample ($6 < z < 12$), we obtain the following results:

- We did not find a clear evolution of $\beta$ with redshift for the galaxies in our sample and showed that this absence of variation is not due to the lower limit in $\beta_0$ imposed by our SED modelling procedure as only 5% of galaxies at $z > 9$ have the lowest value, and 15% of the whole sample at $z > 6$. Considering that the median $M_{UV}$ of our sample is $-17.24^{+0.54}_{-0.62}^*$, this non-evolution is consistent with the results of Topping et al. (2024) at a similar UV magnitude range.

- We provide upper limits on the IR luminosity of NIRCam-selected galaxies at $z > 7$, and show that the detected samples from the literature (REBELS, SERENADE, and individual detections) are not representative of the bulk of galaxies at these redshifts.

- We confirm a linear relation between the V-band attenuation and the inverse sSFR, showing a non-negligible attenuation at $sSFR^{-1} > 30$ Myr and a flatter slope than the one obtained from Langeroodi et al. (2024) from a spectroscopic sample, due to the use of ALMA constraints in our physical property measurements.

Our results focused on $6 < z < 7$ galaxies indicate that the dust properties at the Epoch of Reionisation are mostly similar to the ones in lower-redshift galaxies, implying that significant dust production must have already happened at $z \sim 6.5$. However, the breakdown observed in the $L_{IR}$ versus $M_{UV}$ relation at $M_{UV} \ge -19$ is a possible hint of a different regime in the UV-faint galaxies. These galaxies have, on average, low masses ($\log M_{\star} \sim 8$), resulting in shallower gravitational potential, allowing metals to escape, and resulting in less dust. Or the breakdown could be also due to the lower SFR measured in these sources ($\sim 0.5-2 \text{ M}_\odot \text{ yr}^{-1}$) resulting in a lower dust production. Nevertheless, the rapid build-up concluded from the $6 < z < 7$ galaxies is also confirmed by the non-negligible amount of attenuation observed at $sSFR^{-1} > 30$ Myr for galaxies with $6 < z < 12$.

*Acknowledgements.* LC warmly thanks A. Ferrara, J. Lewis, M. Trebitsch, and D. Elbaz for insightful discussions that improved greatly this work. The authors thank the JADES team for the huge work and effort put in the preparation, observation, and production of the data used in this work. This work received support from the French government under the France 2030 investment plan, as part of the Initiative d'Excellence d'Aix-Marseille Université – A*MIDEX AMX-22-RE-AB-101. This project has received financial support from CNRS and CNES through the MITI interdisciplinary programmes. SA gratefully acknowledges the Collaborative Research Center 1601 (SFB 1601 sub-project C2) funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) – 500700252. LC and MB acknowledges the funding of the French Agence Nationale de la Recherche for the project iMAGE (grant ANR-22-CE31-0007). NIRCam was built by a team at the University of Arizona (UofA) and Lockheed Martin's Advanced Technology Center, led by Prof. Marcia Rieke at UoA. ES acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 694343). This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. MB gratefully acknowledges support from the ANID BASAL project FB210003 and from the FONDECYT regular grant 1211000. This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence of Université Côte d'Azur under reference number ANR-15-IDEX-01.

## References

Adscheid, S., Magnelli, B., Liu, D., et al. 2024, A&A, 685, A1

Algera, H. S. B., Inami, H., Oesch, P. A., et al. 2023, MNRAS, 518, 6142

Algera, H. S. B., Inami, H., De Looze, I., et al. 2024, MNRAS, 533, 3098

Álvarez-Márquez, J., Colina, L., Marques-Chaves, R., et al. 2019, A&A, 629, A9

Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, ApJ, 951, L22

Austin, D., Adams, N., Conselice, C. J., et al. 2023, ApJ, 952, L7

Austin, D., Conselice, C. J., Adams, N. J., et al. 2024, arXiv e-prints [arXiv:2404.10751]

Bakx, T. J. L. C., Tamura, Y., Hashimoto, T., et al. 2020, MNRAS, 493, 4294

Barisic, I., Faisst, A. L., Capak, P. L., et al. 2017, ApJ, 845, A1

Bernhard, E., Béthermin, M., Sargent, M., et al. 2014, MNRAS, 442, 509

Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23

Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2

Bhatawdekar, R., & Conselice, C. J. 2021, ApJ, 909, 144

Boquien, M., Buat, V., Boselli, A., et al. 2012, A&A, 539, A145

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Boquien, M., Buat, V., Burgarella, D., et al. 2022, A&A, 663, A50

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2014, ApJ, 793, 115

Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, ApJ, 833, 72

Bouwens, R. J., Smit, R., Schouws, S., et al. 2022, ApJ, 931, 160

Bowler, R. A. A., Bourne, N., Dunlop, J. S., McLure, R. J., & McLeod, D. J. 2018, MNRAS, 481, 1631

Bowler, R. A. A., Cullen, F., McLure, R. J., Dunlop, J. S., & Avison, A. 2022, MNRAS, 510, 5088

Bowler, R. A. A., Inami, H., Sommovigo, L., et al. 2024, MNRAS, 527, 5808

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Buat, V., Oi, N., Heinis, S., et al. 2015, A&A, 577, A141

Bunker, A. J., Cameron, A. J., Curtis-Lake, E., et al. 2024, A&A, 690, A288

Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455

Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJ, 938, L15

Ciesla, L., Elbaz, D., & Fensch, J. 2017, A&A, 608, A41

Ciesla, L., Elbaz, D., Ilbert, O., et al. 2024, A&A, 686, A128

Cole, J. W., Papovich, C., Finkelstein, S. L., et al. 2025, ApJ, 979, 193

Cullen, F., McLure, R. J., McLeod, D. J., et al. 2023, MNRAS, 520, 14

Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nature Astronomy, 7, 622

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Draine, B. T. 1978, ApJS, 36, 595

Dwek, E. 1986, ApJ, 302, 363

Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023a, arXiv e-prints [arXiv:2306.02465]

Eisenstein, D. J., Johnson, B. D., Robertson, B., et al. 2023b, arXiv e-prints [arXiv:2310.12340]

Elbaz, D., Leiton, R., Nagar, N., et al. 2018, A&A, 616, A110

Endsley, R., Stark, D. P., Fan, X., et al. 2022, MNRAS, 512, 4248

Faisst, A. L., Schaerer, D., Lemaux, B. C., et al. 2020, ApJS, 247, 61

Ferrara, A., Sommovigo, L., Dayal, P., et al. 2022, MNRAS, 512, 58

Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986

Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, ApJ, 756, 164

Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13

Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152

Franco, M., Akins, H. B., Casey, C. M., et al. 2024, ApJ, 973, 23

Fudamoto, Y., Oesch, P. A., Faisst, A., et al. 2020, A&A, 643, A4

Gómez-Guijarro, C., Elbaz, D., Xiao, M., et al. 2022, A&A, 658, A43

Hainline, K. N., Johnson, B. D., Robertson, B., et al. 2024, ApJ, 964, 71

Harikane, Y., Nakajima, K., Ouchi, M., et al. 2024, ApJ, 960, 56

Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019, PASJ, 71, 71

Haskell, P., Smith, D. J. B., Cochrane, R. K., Hayward, C. C., & Anglés-Alcázar, D. 2023, MNRAS, 525, 1535

Heinis, S., Buat, V., Béthermin, M., et al. 2014, MNRAS, 437, 1268

Heintz, K. E., Brammer, G. B., Watson, D., et al. 2025, A&A, 693, A60

Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103

Hollenbach, D. J., & Tielens, A. G. G. M. 1997, ARA&A, 35, 29

Inami, H., Algera, H. S. B., Schouws, S., et al. 2022, MNRAS, 515, 3126

Khusanova, Y., Bethermin, M., Le Fèvre, O., et al. 2021, A&A, 649, A152

Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18

Koprowski, M. P., Coppin, K. E. K., Geach, J. E., et al. 2018, MNRAS, 479, 4355

Langeroodi, D., Hjorth, J., Ferrara, A., & Gall, C. 2024, arXiv e-prints [arXiv:2410.14671]

Laporte, N., Ellis, R. S., Boone, F., et al. 2017, ApJ, 837, L21

Laporte, N., Katz, H., Ellis, R. S., et al. 2019, MNRAS, 487, L81

Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1

Li, J., Da Cunha, E., González-López, J., et al. 2024, ApJ, 976, 70

Liu, D., Lang, P., Magnelli, B., et al. 2019, ApJS, 244, 40

Magnelli, B., Adscheid, S., Wang, T.-M., et al. 2024, A&A, 688, A55

Markov, V., Gallerani, S., Ferrara, A., et al. 2025, Nat. Astron., 9, 458

Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJS, 521, 64

Mitsuhashi, I., Harikane, Y., Bauer, F. E., et al. 2024, ApJ, 971, 161

Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJ, 940, L14

Narayanan, D., Conroy, C., Davé, R., Johnson, B. D., & Popping, G. 2018, ApJ, 869, 70

Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621

Pannella, M., Elbaz, D., Daddi, E., et al. 2015, ApJ, 807, 141

Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152

Reddy, N. A., Shapley, A. E., Sanders, R. L., et al. 2018, ApJ, 869, 92

Rieke, G. M. J., Robertson, B., Tacchella, S., et al. 2023, ApJS, 269, 16

Roberts-Borsani, G., Treu, T., Shapley, A., et al. 2024, ApJ, 976, 193

Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, Nat. Astron., 7, 611

Robertson, B., Johnson, B. D., Tacchella, S., et al. 2024, ApJ, 970, 31

Rogers, A. B., McLure, R. J., Dunlop, J. S., et al. 2014, MNRAS, 440, 3714

Rowland, L. E., Stefanon, M., Bouwens, R., et al. 2025, arXiv e-prints [arXiv:2501.10559]

Salim, S., & Narayanan, D. 2020, ARA&A, 58, 529

Salpeter, E. E. 1955, ApJ, 121, 161

Saxena, A., Cameron, A. J., Katz, H., et al. 2024, arXiv e-prints [arXiv:2411.14532]

Schouws, S., Stefanon, M., Bouwens, R., et al. 2022, ApJ, 928, 31

Schouws, S., Bouwens, R. J., Ormerod, K., et al. 2025, ApJS, 988, 19

Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30

Seillé, L. M., Buat, V., Haddad, W., et al. 2022, A&A, 665, A137

Shivaei, I., Alberts, S., Florian, M., et al. 2024, A&A, 690, A89

Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2022, MNRAS, 513, 3122

Sugahara, Y., Inoue, A. K., Hashimoto, T., et al. 2021, ApJ, 923, 5

Tamura, Y., Mawatari, K., Hashimoto, T., et al. 2019, ApJ, 874, 27

Topping, M. W., Stark, D. P., Endsley, R., et al. 2024, MNRAS, 529, 4087

Vijayan, A. P., Wilkins, S. M., Lovell, C. C., et al. 2022, MNRAS, 511, 4999

Wang, T.-M., Magnelli, B., Schinnerer, E., et al. 2022, A&A, 660, A142

Wang, T.-M., Magnelli, B., Schinnerer, E., et al. 2024, A&A, 681, A110

Watson, D., Christensen, L., Knudsen, K. K., et al. 2015, Nature, 519, 327

Williams, C. C., Tacchella, S., Maseda, M. V., et al. 2023, ApJS, 268, 64

Witstok, J., Smit, R., Maiolino, R., et al. 2022, MNRAS, 515, 1751

Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152