# A3COSMOS: Measuring the cosmic dust-attenuated star formation rate density at $4 < z < 5$

Benjamin Magnelli, Sylvia Adscheid, Tsan-Ming Wang, Laure Ciesla, Emanuele Daddi, Ivan Delvecchio, David Elbaz, Yoshinobu Fudamoto, Shuma Fukushima, Maximilien Franco, et al.

► To cite this version:

Benjamin Magnelli, Sylvia Adscheid, Tsan-Ming Wang, Laure Ciesla, Emanuele Daddi, et al. A3COSMOS: Measuring the cosmic dust-attenuated star formation rate density at $4 < z < 5$. Astronomy & Astrophysics - A&A, 2024, 688, pp.A55. 10.1051/0004-6361/202450081 . cea-04666674

HAL Id: cea-04666674

https://cea.hal.science/cea-04666674v1

Submitted on 1 Aug 2024

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
# A³COSMOS: Measuring the cosmic dust-attenuated star formation rate density at 4 < z < 5

Benjamin Magnelli¹, Sylvia Adscheid², Tsan-Ming Wang², Laure Ciesla³, Emanuele Daddi¹,
Ivan Delvecchio⁴, David Elbaz¹, Yoshinobu Fudamoto⁵,⁶,⁷, Shuma Fukushima⁸, Maximilien Franco⁹,
Carlos Gómez-Guijarro¹, Carlotta Gruppioni¹⁰, Eric F. Jiménez-Andrade¹¹, Daizhong Liu¹²,¹³
Pascal Oesch¹⁴,¹⁵, Eva Schinnerer¹⁶, and Alberto Traina¹⁰,¹⁷

*(Affiliations can be found after the references)*

Received 22 March 2024 / Accepted 27 May 2024

## ABSTRACT

**Context.** In recent years, conflicting results have provided an uncertain view of the dust-attenuated star-forming properties of $z \ge 4$ galaxies.

**Aims.** To solve this, we need to accurately measure the mean dust-attenuated properties of star-forming galaxies (SFGs) at $4 < z < 5$ and therefore constrain the cosmic dust-attenuated star formation rate density (SFRD) of the Universe 1.3 Giga-years after the Big Bang.

**Methods.** We used the deepest optical-to-near-infrared data publicly available in the Cosmic Evolution Survey (COSMOS) field to build a mass-complete ($>10^{9.5} M_\odot$) sample of SFGs at $4 < z < 5$. Then, we measured their mean dust-attenuated properties (i.e., infrared luminosity, $\langle L_{IR} \rangle$; dust-attenuated star formation rate, $\langle SFR_{IR} \rangle$) by dividing our sample in three stellar mass ($M_*$) bins (i.e., $10^{9.5} < M_*/M_\odot < 10^{10}$, $10^{10} < M_*/M_\odot < 10^{10.5}$, and $10^{10.5} < M_*/M_\odot < 10^{11.5}$) and by stacking in the uv domain all archival Atacama Large Millimeter/submillimeter Array (ALMA) band 6 and 7 observations available for these galaxies. Then, we combined this information with their mean rest-frame ultraviolet (UV) emission measured from the COSMOS2020 catalog (i.e., UV luminosity, $\langle L_{UV} \rangle$; UV spectral slope, $\langle \beta_{UV} \rangle$; and unattenuated SFR, $\langle SFR_{UV} \rangle$), and constrained the IRX ($\equiv L_{IR}/L_{UV}$)-$\beta_{UV}$, IRX-$M_*$, and SFR-$M_*$ relations at $z \sim 4.5$. Finally, using these relations and the stellar mass function of SFGs at $z \sim 4.5$, we inferred the unattenuated and dust-attenuated SFRD at this epoch.

**Results.** SFGs follow an IRX-$\beta_{UV}$ relation that is consistent with that observed in local starbursts. Our measurements favors a steepening of the IRX-$M_*$ relation at $z \sim 4.5$, compared to the redshift-independent IRX-$M_*$ relation observed at $z \sim 1-3$. Our galaxies lie on a linear SFR-$M_*$ relation, whose normalization varies by 0.3 dex, when we exclude or include from our stacks the ALMA primary targets (i.e., sources within 3'' from the ALMA phase center). The cosmic SFRD($>M_*$) converges at $M_* \lesssim 10^9 M_\odot$, with SFGs at $10^8 < M_*/M_\odot < 10^9$ contributing already less than 15% of the SFRD from all SFGs with $M_* > 10^8 M_\odot$. The cosmic SFRD at $z \sim 4.5$ is dominated by SFGs with a stellar mass of $10^{9.5-10.5} M_\odot$. Finally, the fraction of the cosmic SFRD that is attenuated by dust, SFRD$_{IR}$(>M$_*$)/SFRD(>M$_*$), is $90 \pm 4\%$ for $M_* = 10^{10} M_\odot$, $68 \pm 10\%$ for $M_* = 10^{8.9} M_\odot$ (i.e., $0.03 \times M^*; M^*$ being the characteristic stellar mass of SFGs at this epoch) and this value converges to $60 \pm 10\%$ for $M_* = 10^8 M_\odot$.

**Conclusions.** A non-evolving IRX-$\beta_{UV}$ relation suggests that the grain properties (e.g., size distribution, composition) of dust in SFGs at $z \sim 4.5$ are similar to those in local starbursts. However, the mass and geometry of this dust result in lower attenuation in low-mass SFGs ($\lesssim 10^{10} M_\odot$) at $z \sim 4.5$ than at $z \lesssim 3$. Nevertheless, the fraction of the cosmic SFRD that is attenuated by dust remains significant ($\sim 68 \pm 10\%$) even at such an early cosmic epoch.

**Key words.** dust, extinction – galaxies: evolution – galaxies: high-redshift – infrared: galaxies

## 1. Introduction

One of the greatest achievements of the last 20 years in the field of galaxy evolution is the high precision measurement of the cosmic star formation rate density (SFRD) up to $z \sim 4$ (Madau & Dickinson 2014). This accomplishment was made possible by combining rest-frame infrared (IR) and ultraviolet (UV) observations probing, respectively, the dust-attenuated and unattenuated star formation rate (SFR) of star-forming galaxies (SFGs). At $z > 4$, however, estimates of the cosmic SFRD become highly uncertain because current measurements rely primarily on rest-frame UV observations to measure the SFRs of galaxies, accounting for dust attenuation using the local relations between the UV spectral slope ($\beta_{UV}$) and $L_{IR}/L_{UV}$ (≡IRX; where $L_{IR}$ and $L_{UV}$ are the IR and UV luminosities, respectively) or between the stellar mass ($M_*$) and IRX (Madau & Dickinson 2014). While such dust attenuation corrections are known to fail for a large fraction of massive galaxies at low redshifts ($z \lesssim 2.5$; Wuyts et al. 2011), initial Lyman break studies argued that the

low dust content of $z > 4$ galaxies made this method reliable in the early Universe (e.g., Bouwens et al. 2012). Observations from the Atacama Large Millimeter/submillimeter Array (ALMA) have since cast doubt on the ability of rest-frame UV data alone to provide an accurate view on the cosmic dust-attenuated SFRD, despite somewhat contradictory literature results. On the one hand, serendipitous detections in the ALMA large program ALPINE suggest that dust-attenuated star formation still dominates the cosmic SFRD up to $z \sim 6$ and that initial rest-frame UV measurements have underestimated its contribution in the early Universe (Gruppioni et al. 2020; Khusanova et al. 2021). On the other hand, ALMA continuum follow-ups of UV-selected galaxies at $z > 4$ have revealed a redshift evolution of the dust attenuation, with galaxies at $z > 4$ having even lower IR luminosities than expected from their UV spectral slope (Fudamoto et al. 2020a), although opposite results have recently appeared (e.g., Bowler et al. 2024). These conflicting results provide a highly uncertain view on the dust attenuation and dust-attenuated star formation properties of galaxies
at $z > 4$: the IRX-$M_*$ and IRX-$\beta_{UV}$ relations could be consistent with their low redshift incarnation or significantly steeper and flatter, respectively (Fudamoto et al. 2020a; Bowler et al. 2022, 2024); the SFR-$M_*$ relation is unknown within a factor ~4 (Khusanova et al. 2021; Popesso et al. 2023; Goovaerts et al. 2024), although it is commonly used to select “normal” SFGs at these epochs; and the cosmic SFRD is uncertain by a factor ~10 at $z \sim 4$ (Gruppioni et al. 2020; Zavala et al. 2021; Algera et al. 2023; Traina et al. 2024).

Current confusion on the dust attenuation properties of high-redshift galaxies stems from the past use of biased samples of $z > 4$ SFGs (UV-, IR-, or $z_{spec}$-selected), and from the fact that these past studies are also based on few objects, so that average measurements are dominated by the intrinsic scatter of the SFR-$M_*$ and IRX-$M_*$ relations. To constrain statistically these relations and the cosmic SFRD at $z > 4$, it is important to measure the mean IR and UV luminosities of a large sample of mass-selected SFGs at these redshifts. To this end, we are taking advantage of the deepest optical-to-near-IR data publicly available in the Cosmic Evolution Survey (COSMOS; i.e., COSMOS2020 catalog; Weaver et al. 2022) and we are capitalizing on the availability of a large number of ALMA archival data on this field, collected by the A$^3$COSMOS project (Liu et al. 2019a; Adscheid et al. 2024). We selected SFGs at $4 < z < 5$ from the COSMOS2020 catalog, a deep near-infrared izYJHK$_S$-selected catalog (with $5\sigma$ limits of 27 and 24.7 mag$_{AB}$ in the deepest *i* and shallowest *K$_S$* bands, respectively) that allows for detection of dusty galaxies that are absent from classical UV/optical-selected catalogs (Weaver et al. 2023). We applied an uv-domain stacking analysis to the archival ALMA data of this sample. Our stacking analysis is crucial as it provides measurements of the mean IR luminosities of these SFGs down to a stellar mass of ~$10^{9.5} M_\odot$, even though most of these galaxies are not detected individually in the ALMA archival data.

The structure of this paper is as follows. In Sect. 2, we present the A$^3$COSMOS database used in our analysis. In Sect. 3, we describe how we built our mass-complete sample of SFGs at $4 < z < 5$. In Sect. 4, we describe how we measured the mean IR and UV properties of our galaxies using a uv-domain stacking analysis of the A$^3$COSMOS database and the COSMOS2020 catalog, respectively. In Sects. 5.1–5.4, we present the main results of this study, that is, the IRX-$\beta_{UV}$ and IRX-$M_*$ relations, the SFR-$M_*$ relation, the cosmic dust-attenuated SFRD, and the “total” cosmic dust-attenuated SFRD at $z \sim 4.5$, respectively. Finally, in Sect. 6, we summarize and conclude.

We assume a flat $\Lambda$ cold dark matter cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.30$, and $\Omega_\Lambda = 0.70$. A Chabrier (2003) initial mass function (IMF) is assumed for all stellar masses and SFRs.

## 2. Data

We used all archival ALMA band 6 and band 7 observations publicly available in COSMOS (RA = $10^h00^m28.6^s$, Dec = $+02^\circ12'21.0''$) as of September 1, 2022. We have limited our analysis to these bands, as they are close to the peak of the dust emission of $z \sim 4.5$ SFGs and therefore provide good proxies of their infrared luminosities. In contrast, longer-wavelength bands probe the Rayleigh-Jeans emission of dust at these redshifts, which is a poor proxy for the infrared luminosity. Shorter-wavelength bands would in principle be better proxies of the infrared luminosity of SFGs at $z \sim 4.5$, but the number of such observations in the ALMA archive is so small that they would have no impact on our analysis. This database of all

archival ALMA band 6 and band 7 observations was assembled by the A$^3$COSMOS project, which aims to process all ALMA projects targeting the COSMOS field in a homogeneous way, and to provide calibrated visibilities, cleaned images, and value-added source catalogs (targeted and serendipitously detected) via a single access portal (Liu et al. 2019a; Adscheid et al. 2024). We limited our analysis to projects from cycles $\ge$3 because measurement sets from earlier cycles have different definitions of visibility weights, which makes their use in our uv-domain stacking analysis problematic (Wang et al. 2022). Our final database contains 87 ALMA projects (53 in band 6 and 34 in band 7), and has a total of 2100 images (equivalently ALMA pointings; 946 in band 6 and 1154 in band 7).

## 3. Sample

To build our mass-selected sample of SFGs, we used the deepest data publicly available in COSMOS (i.e., COSMOS2020 catalog; Weaver et al. 2022), as this wide, deep survey is ideal for studying representative samples of massive, high-redshift galaxies. Following Weaver et al. (2023), we used the FARMER COSMOS2020 catalog, which is the areal union of the deep near-infrared UltraVISTA DR4 imaging and the Subaru Suprime-Cam intermediate bands PDR2 that covers $1.279 \text{ deg}^2$ after removing contamination due to bright stars. Galaxies in the COSMOS2020 catalog were selected from the near-infrared izYJHK$_S$ co-added detection image, with $5\sigma$ limits of 27 and 24.7 mag$_{AB}$ (2'' aperture) in the deepest *i* and shallowest *K$_S$* bands, respectively. As in Weaver et al. (2023), we selected SFGs from this FARMER COSMOS2020 catalog using a standard NUV-$r-J$ criterion (Ilbert et al. 2013) and only kept SFGs with robust stellar mass estimates ($m_{IRAC3.6} < 26 \text{ AB}$ and $\chi^2_{\text{reduced}} < 10$) and accurate photometric redshifts ($\Delta z/(1+z) < 0.1$) at $4 < z < 5$. As detailed in Weaver et al. (2023), these photometric redshift and stellar mass estimates were obtained by fitting with LePhare the FARMER photometry in 30 bands from UV to near-infrared. We adopted the photometric redshift and stellar mass estimates from LP_zPDF and LP_MASS_MED, defined as the median of their likelihood distributions as they are generally less susceptible to template fitting systematics than those taken at the minimum $\chi^2_{\text{reduced}}$. From this initial sample, we further excluded AGNs, identified through a X-ray Chandra detection within $0.6^\circ$ (Weaver et al. 2022) and/or a spectral energy distribution (SED) best fit with an AGN template ($\chi^2_{\text{reduced}}(\text{AGN}) < 0.5 \times \chi^2_{\text{reduced}}(\text{Gal.})$; Weaver et al. 2022). These AGN exclusions concern less than 0.1% of the SFGs in our $z \sim 4.5$ parent sample, and has a negligible impact on our results. Finally, we kept only galaxies above the stellar mass completeness limit for SFGs in the COSMOS2020 catalog, that is, $M_* > 10^{9.5} M_\odot$ at $4 < z < 5$ (Weaver et al. 2023). This mass-complete sample, identical to that used in Weaver et al. (2023) to build the stellar mass function at $z \sim 4.5$, contains 5810 galaxies. We note that, while this sample may still miss part of the population of optically-dark galaxies that are heavily dust attenuated and have recently been discovered through Spitzer-IRAC and ALMA observations (e.g., Wang et al. 2019; Xiao et al. 2023), their contribution to the dust-attenuated SFRD remains modest as discussed in Sect. 5.4.

Finally, we cross-matched this mass-complete sample with the A$^3$COSMOS database, and excluded galaxies that were not covered by any ALMA observations (i.e., located at regions with a primary beam response of less than 0.2). This reduces our sample to 440 galaxies, with redshift, stellar mass, UV luminosity, and UV spectral slope distributions that are indistinguishable
<table><caption>**Table 1.** Measured dust-attenuated luminosities and derived properties of our six uv-domain stacks.</caption><thead><tr><th>Stellar mass<br>log<sub>10</sub>(M<sub>☉</sub>)</th><th>N</th><th>log<sub>10</sub>⟨M<sub>*</sub>⟩</th><th>⟨z⟩</th><th>log<sub>10</sub>⟨L<sup>rest</sup><sub>150</sub>⟩</th><th>S/N</th><th>log<sub>10</sub>⟨SFR<sub>IR</sub>⟩</th><th>log<sub>10</sub>⟨SFR<sub>UV</sub>⟩</th><th>⟨β<sub>UV</sub>⟩</th><th>log<sub>10</sub>⟨IRX⟩</th></tr></thead><tbody><tr><td colspan="11" style="text-align: center;">Including all SFGs with a A<sup>3</sup>COSMOS coverage</td></tr><tr><td>9.5–10</td><td>302</td><td>9.7<sup>+0.1</sup><sub>-0.1</sub></td><td>4.4<sup>+0.1</sup><sub>-0.1</sub></td><td>31.60<sup>+0.12</sup><sub>-0.16</sub></td><td>6</td><td>1.19<sup>+0.14</sup><sub>-0.17</sub></td><td>0.91<sup>+0.04</sup><sub>-0.05</sub></td><td>-1.53<sup>+0.10</sup><sub>-0.10</sub></td><td>0.45<sup>+0.14</sup><sub>-0.22</sub></td></tr><tr><td>10–10.5</td><td>96</td><td>10.3<sup>+0.1</sup><sub>-0.1</sub></td><td>4.5<sup>+0.1</sup><sub>-0.1</sub></td><td>31.98<sup>+0.07</sup><sub>-0.09</sub></td><td>8</td><td>1.56<sup>+0.09</sup><sub>-0.11</sub></td><td>0.75<sup>+0.13</sup><sub>-0.20</sub></td><td>-1.03<sup>+0.21</sup><sub>-0.21</sub></td><td>1.00<sup>+0.15</sup><sub>-0.23</sub></td></tr><tr><td>10.5–11.5</td><td>42</td><td>10.8<sup>+0.1</sup><sub>-0.1</sub></td><td>4.5<sup>+0.1</sup><sub>-0.1</sub></td><td>33.16<sup>+0.25</sup><sub>-0.68</sub></td><td>37</td><td>2.74<sup>+0.26</sup><sub>-0.68</sub></td><td>0.91<sup>+0.12</sup><sub>-0.17</sub></td><td>-0.31<sup>+0.62</sup><sub>-0.62</sub></td><td>2.01<sup>+0.27</sup><sub>-0.85</sub></td></tr><tr><td colspan="11" style="text-align: center;">Excluding the ALMA primary targets</td></tr><tr><td>9.5–10</td><td>271</td><td>9.7<sup>+0.1</sup><sub>-0.1</sub></td><td>4.4<sup>+0.1</sup><sub>-0.1</sub></td><td>31.28<sup>+0.33</sup><sub>-∞</sub></td><td>3</td><td>0.87<sup>+0.34</sup><sub>-∞</sub></td><td>0.70<sup>+0.06</sup><sub>-0.07</sub></td><td>-1.43<sup>+0.15</sup><sub>-0.15</sub></td><td>0.34<sup>+0.35</sup><sub>-∞</sub></td></tr><tr><td>10–10.5</td><td>79</td><td>10.3<sup>+0.1</sup><sub>-0.1</sub></td><td>4.5<sup>+0.1</sup><sub>-0.1</sub></td><td>31.86<sup>+0.12</sup><sub>-0.16</sub></td><td>5</td><td>1.44<sup>+0.14</sup><sub>-0.17</sub></td><td>0.44<sup>+0.13</sup><sub>-0.19</sub></td><td>-0.79<sup>+0.23</sup><sub>-0.23</sub></td><td>1.17<sup>+0.17</sup><sub>-0.28</sub></td></tr><tr><td>10.5–11.5</td><td>26</td><td>10.8<sup>+0.1</sup><sub>-0.1</sub></td><td>4.7<sup>+0.1</sup><sub>-0.1</sub></td><td>32.66<sup>+0.21</sup><sub>-0.41</sub></td><td>15</td><td>2.25<sup>+0.22</sup><sub>-0.41</sub></td><td>0.50<sup>+0.24</sup><sub>-0.58</sub></td><td>0.51<sup>+0.50</sup><sub>-0.50</sub></td><td>1.92<sup>+0.29</sup><sub>-1.43</sub></td></tr></tbody></table>

**Notes.** Errors on each quantity were obtained from a bootstrap analysis. We used the 50th percentile as the best estimate and the 16th and 84th percentiles of this bootstrap analysis as errors. For $\langle L_{150}^{\text{rest}} \rangle$, we also provide the signal-to-noise ratio (S/N) in the stacks. Errors on $\langle L_{150}^{\text{rest}} \rangle$ from the bootstrap analysis are much larger than these, as they take into account not only the photometric noise in the stacks but also the intrinsic dispersion in the luminosities of the stacked population.

(Kolmogorov-Smirnov test) from those of their parent sample of 5810 galaxies (see Appendix A). Of these 440 galaxies, 64 can be considered as the primary target of a given ALMA observation, as they are located within 3" of its phase center. In Appendix A, we found that these ALMA primary targets are biased toward bright, massive galaxies in a redshift range favorable to [CII] observations, but do not appear to be biased toward a particular UV spectral slope. These primary targets can thus potentially introduce complex and uncontrollable selection biases into our final ALMA-covered mass-complete galaxy sample. On the other hand, according to the analysis presented in Appendix A, excluding these ALMA primary targets provides us with a fairly representative sample of our parent sample of $z \sim 4.5$ SFGs, albeit perhaps slightly biased against bright galaxies. In Sect. 5, we therefore investigate the impact on our results of including or excluding these ALMA primary targets in our stacks. The number of galaxies in each of our stellar mass bins, and the number of galaxies that can be considered as the ALMA primary target, are given in Table 1.

# 4. Methods

## 4.1. Measuring $\langle L_{IR} \rangle$ by stacking ALMA data in the uv domain

Stacking the A$^3$COSMOS dataset is a challenge, as this database is heterogeneous in terms of observed frequencies, angular resolution, and integration time. As demonstrated in Wang et al. (2022), this challenge can, however, be overcome. The problem of frequency heterogeneity is solved by prior rescaling of each individual dataset to a common rest-frame frequency using submillimeter SED templates, while the problem of angular resolution heterogeneity is solved by performing our stacking analysis in the uv domain. Finally, the problem of integration time heterogeneity is simply solved by keeping track of the weight of each galaxies in the final stacks.

To perform the uv-domain stacking analysis on the A$^3$COSMOS dataset, we followed the same approach as Wang et al. (2022). Details regarding these steps can be found in their Sects. 3.1 and 3.2. First, we scaled the observed ALMA visibility amplitudes of each galaxy to its rest-frame luminosity at 150 µm (i.e., $L_{150}^{\text{rest}}$) by using the redshift of the source, the SED

template of Béthermin et al. (2017) and the CASA tasks gencal and applycal. This particular SED template was chosen to be consistent with the $L_{IR}$-to-$L_{150}^{\text{rest}}$ conversion used subsequently (see below) and because it provides a reasonable fit to the Herschel stacks of $4 < z < 5$ SFGs performed in Béthermin et al. (2020). The reference wavelength of 150 µm was chosen to limit the impact of the assumed SED template on our results, as for galaxies located at $z \sim 4.5$, it actually corresponds to the rest-frame wavelength probed by the ALMA band 7, which dominates in number the A$^3$COSMOS database. For each galaxy, we then shifted the phase center of its visibilities to its coordinate using a CASA-based package, the STACKER (Lindroos et al. 2015). The visibilities of the galaxies in each of our stellar mass bins were then stacked together using the CASA task concat. The cleaned image was generated from the stacked measurement set with the CASA task tclean, using natural weighting and cleaning the image down to $3\sigma$. Finally, we measured the stacked $L_{150}^{\text{rest}}$ (i.e., $\langle L_{150}^{\text{rest}} \rangle$) using the central pixel value, as at the angular resolution of our stacked measurement set (i.e., $\sim 0.5'5-1''$), our stacked population remains basically spatially unresolved. We verified, however, that using instead 2D Gaussian fits at the phase center of the stacked images gave consistent $\langle L_{150}^{\text{rest}} \rangle$ measurements.

To assess the uncertainty on $\langle L_{150}^{\text{rest}} \rangle$, we used a bootstrap analysis. For each stellar mass bin, we made 100 realizations of the stacking analysis, using for each realization a different sample, drawn from the original one, with the same number of sources but allowing for replacement (i.e., a galaxy can be picked several times). In what follows, we use the 50th percentile as the best estimate and the 16th and 84th percentiles of this bootstrap analysis as errors on $\langle L_{150}^{\text{rest}} \rangle$.

These $\langle L_{150}^{\text{rest}} \rangle$ measurements were converted into mean IR luminosities, $\langle L_{IR} \rangle$ [8–1000 µm], following Khusanova et al. (2021),

$$ \frac{\langle L_{IR} \rangle}{L_{\odot}} = (3.80 \pm 0.63) \times 10^{-21} \times \frac{\langle L_{150}^{\text{rest}} \rangle}{\text{erg s}^{-1} \text{Hz}^{-1}} \quad (1) $$

This conversion factor was calculated in Khusanova et al. (2021) by averaging the $L_{IR}$-to-$L_{150}^{\text{rest}}$ conversion factors of all SED templates providing a reasonable fit (i.e., $\chi^2_{\text{reduced}} < 1.5$) to the Herschel stacks of $4 < z < 5$ SFGs performed in Béthermin et al. (2020). It corresponds approximately to that of
<figure><img src="image_2.png" /><figcaption>Fig. 1. Cutouts of the stacked images. These cutouts correspond to our $10^{9.5} < M_*/M_\odot < 10^{10}$ (left), $10^{10} < M_*/M_\odot < 10^{10.5}$ (center), and $10^{10.5} < M_*/M_\odot < 10^{11.5}$ (right) stellar mass bins, including (top) or excluding (bottom) the ALMA primary targets. The number of stacked galaxies is given in each cutout. The red dashed ellipse shows the ALMA synthesized beam. Color scales are similar for a given stellar mass bin. Each cutout has a size of $4'' \times 4''$.</figcaption></figure>

a modified blackbody function with a dust temperature of 41 K (Khusanova et al. 2021). Uncertainties on the $L_{IR}$-to-$L_{150}^{rest}$ conversion factor are mitigated with respect to the assumed dust temperature, since at this rest-frame wavelength we are probing close to the dust SED peak. Errors in these $L_{IR}$-to-$L_{150}^{rest}$ conversion factors, as evaluated in Khusanova et al. (2021), are naturally propagated to all the measurements presented hereafter. We note that the SED template of Béthermin et al. (2017) used here to rescale each individual visibility amplitudes to a common rest-frame frequency is one of the templates providing the best $\chi^2_{reduced}$ from comparison to the Herschel stacks of Béthermin et al. (2020). Also, it yields a $L_{IR}$-to-$L_{150}^{rest}$ conversion factor that is within 5% to that of Eq. (1). Our rescaling to a common rest-frame frequency and our $L_{IR}$-to-$L_{150}^{rest}$ conversion are thus fully consistent.

Finally, these mean IR luminosities were converted into dust-attenuated SFRs following Madau & Dickinson (2014) for a Chabrier (2003) IMF,

$$ \frac{\langle \text{SFR}_{\text{IR}} \rangle}{M_{\odot} \text{yr}^{-1}} = 1 \times 10^{-10} \times \frac{\langle L_{\text{IR}} \rangle}{L_{\odot}}. \qquad (2) $$

The results of this uv-domain stacking analysis (i.e., $\langle L_{150}^{rest} \rangle$, $\langle L_{IR} \rangle$, and $\langle \text{SFR}_{\text{IR}} \rangle$) applied to our mass-complete sample of $4 < z < 5$ SFGs divided into three stellar mass bins (i.e., $10^{9.5} < M_*/M_\odot < 10^{10}$, $10^{10} < M_*/M_\odot < 10^{10.5}$, and $10^{10.5} < M_*/M_\odot < 10^{11.5}$), are given in Table 1 and shown in Fig. 1. This particular choice of stellar mass bin was made to maximize the signal-to-noise ratio of the stacks while providing a

good sampling of the stellar mass range probed by our sample. Figure 1 shows that indeed in each of our stellar mass bins we obtain robust detection in terms of signal-to-noise ratio (S/N). Naturally, errors on $\langle L_{150}^{rest} \rangle$ from our bootstrap analysis are much larger, as they take into account not only the photometric noise in the stacks but also the intrinsic dispersion in the luminosities of the stacked population.

We note that the dust-attenuated properties given in Table 1 are weighted values, as a galaxy with deeper ALMA coverage has a higher weight in our uv-domain stacking analysis. The weight of each galaxy is simply given by the number of visibilities to which it corresponds in our stacks (ALMA achieved greater sensitivity by increasing the number of observed visibilities, rather than by increasing the integration time per visibility). Therefore, in what follows, when we compare these dust-attenuated properties to other physical quantities (e.g., redshift, stellar mass, UV luminosity), these are also weighted using the number of visibilities of each galaxy.

## 4.2. Measuring $\langle L_{UV} \rangle$ and $\langle \beta_{UV} \rangle$ using the COSMOS2020 catalog

In our study, the UV luminosity refers to monochromatic luminosity at 1600 Å (i.e., $L_{UV} \equiv v_{1600}L_{v_{1600}}$), not corrected for dust attenuation. For each of our galaxies, $L_{UV}$ was taken from the COSMOS2020 catalog and corresponds to the UV luminosity of their best-fit SED found with LePhare assuming the Bruzual & Charlot (2003) model with a delayed star formation history (for more details see Weaver et al. 2022). The robustness
of these $L_{UV}$ measurements is ensured by the fact that all our galaxies have at least one (but more for most; see below) photometric detection (i.e., $S/N > 3$) in the rest-frame UV (i.e., $\lambda_{rest} < 3300$ Å) and that they have accurate SED fits with LePhare (i.e., $\chi^2_{reduced} < 10$; see Sect. 3). From these $L_{UV}$, we then measured $\langle L_{UV} \rangle$, accounting for the weight of each galaxy in our uv-domain stacking analysis (see Sect. 4.1). We assumed that the uncertainties on $\langle L_{UV} \rangle$ were dominated by the scatter within the population rather than by fitting uncertainties. Therefore, the uncertainties on $\langle L_{UV} \rangle$ were estimated by bootstrapping over the distributions of $L_{UV}$ in each stellar mass bin. These mean UV luminosities were then converted into unattenuated SFRs following Madau & Dickinson (2014) for a Chabrier (2003) IMF,

$$ \frac{\langle \text{SFR}_{\text{UV}} \rangle}{M_{\odot} \text{yr}^{-1}} = 1.5 \times 10^{-10} \times \frac{\langle L_{\text{UV}} \rangle}{L_{\odot}}. \quad (3) $$

Finally, from the best-fit SED of each of our galaxies in COSMOS2020 catalog, we also inferred their UV spectral slope (i.e., $\beta_{UV}$ with $f_{\lambda} \propto \lambda^{\beta_{UV}}$). To this end, we used their rest-frame luminosities at 1600 Å and 2300 Å (e.g., Calzetti et al. 1994). From these $\beta_{UV}$, we then measured $\langle \beta_{UV} \rangle$ and associated uncertainties using a bootstrap analysis and accounting for the weights of each galaxy in our uv-domain stacking analysis (see Sect. 4.1).

The mean UV luminosities $\langle L_{UV} \rangle$, UV spectral slopes $\langle \beta_{UV} \rangle$ and unattenuated SFRs $\langle \text{SFR}_{UV} \rangle$ of the SFGs in our three stellar mass bins are given in Table 1. In addition, this table also provides their mean stellar mass and redshift, as inferred using the COSMOS2020 catalog, and accounting for the weights of each galaxy in our uv-domain stacking analysis (see Sect. 4.1).

To test the robustness of our UV luminosity and UV spectral slope measurements and especially to ensure that they are not biased because they are based on SED fits, we turned to the observed photometry of our galaxies, after applying to the COSMOS2020 magnitudes the appropriate zero-point offsets and foreground extinction corrections (see Laigle et al. 2016; Weaver et al. 2022). To measure the UV spectral slope (i.e., $\beta_{UV}^{\text{phot.}}$), we required for each galaxy at least two detections (i.e., $S/N > 3$) in the rest-frame UV (i.e., $\lambda_{rest} < 3300$ Å) separated by 500 Å. For our two lower stellar mass bins, this requirement is met by 100% of our galaxies, while in our highest stellar mass bin it is met by only 80% of our galaxies, probably due to higher dust attenuation. We then calculated the UV luminosity (i.e., $L_{UV}^{\text{phot.}}$) of these galaxies using the photometry closest to the rest-frame 1600 Å and applying a k-correction based on the previously measured UV spectral slope. For all galaxies for which this analysis was possible (i.e., ~95%), we find very good agreements between $L_{UV}$ and $L_{UV}^{\text{phot.}}$. The median and 16th and 84th percentiles of the $L_{UV}^{\text{phot.}}$-to-$L_{UV}$ ratio is $1.03^{+0.11}_{-0.12}$, with no significant dependence on stellar mass. Our fiducial UV luminosities $L_{UV}$ can therefore be considered accurate to within ~10% and only slightly overestimated by ~3%. As these systematics and uncertainties are much smaller than those introduced by the dispersion in $L_{UV}$ within the galaxy population of each of our stellar mass bins, we decided, for the sake of completeness, to use $L_{UV}$ instead of $L_{UV}^{\text{phot.}}$ in the rest of our analysis. Similarly, we find good, although poorer, agreements between $\beta_{UV}^{\text{phot.}}$ and $\beta_{UV}$, with a median $\beta_{UV}^{\text{phot.}}$-to-$\beta_{UV}$ ratio and 16th and 84th percentiles of $0.74^{+0.35}_{-0.51}$. Despite these slight differences, for the sake of completeness and consistency, we decided to use $\beta_{UV}$ instead of $\beta_{UV}^{\text{phot.}}$ in the rest of our analysis. However, in Sect. 5.1, we illustrate the

impact of using $\beta_{UV}^{\text{phot.}}$ instead of $\beta_{UV}$ on the observed IRX-$\beta_{UV}$ relation.

## 5. Results

Using our uv-domain stacking analysis, we are able to measure for the first time the mean dust attenuation properties (i.e., $\langle L_{IR} \rangle$ and $\langle \text{SFR}_{IR} \rangle$) of a mass-complete sample of SFGs at $4 < z < 5$ divided in three stellar mass bins, that is, $10^{9.5} < M_{*}/M_{\odot} < 10^{10}$, $10^{10} < M_{*}/M_{\odot} < 10^{10.5}$, and $10^{10.5} < M_{*}/M_{\odot} < 10^{11.5}$ (see Table 1). Combining this information with the emission of these galaxies in the rest-frame UV (i.e., $\langle L_{UV} \rangle$, $\langle \beta_{UV} \rangle$, and $\langle \text{SFR}_{UV} \rangle$), we now study the IRX-$\beta_{UV}$, IRX-$M_{*}$, and SFR-$M_{*}$ relations at $z \sim 4.5$ and infer the total, dust-attenuated, and unattenuated cosmic SFRD at this epoch.

### 5.1. The IRX-$\beta_{UV}$ and IRX-$M_{*}$ relations

The IR excess (IRX $\equiv L_{IR}/L_{UV}$) and its relation with $\beta_{UV}$ and $M_{*}$ is often used to account for the dust-attenuated fraction of the SFR seen in the UV. Indeed, in the absence of IR photometry, the wealth of optical-to-near-IR data available for high-redshift galaxies allows us in most cases to accurately measure $L_{UV}$, $\beta_{UV}$, and/or $M_{*}$, which, combined with the (locally) calibrated IRX-$\beta_{UV}$ and IRX-$M_{*}$ relations, can be used to estimate $L_{IR}$ (and consequently SFR$_{IR}$).

The IRX-$\beta_{UV}$ relation is the most fundamental of these two relations. Indeed, assuming a universal, dust unattenuated UV spectral slope $\beta_0$ for SFGs, the frequency dependence of the attenuation curve directly determines the position of these galaxies in the IRX-$\beta_{UV}$ plane (see, e.g., Meurer et al. 1999; Fudamoto et al. 2020a),

$$ \log_{10} \text{IRX} \equiv \log_{10} \left( \frac{L_{\text{IR}}}{L_{\text{UV}}} \right) = 1.5 \times \left( 10^{0.4 \times \frac{dA_{1600}}{d\beta_{\text{UV}}} \times (\beta_{\text{UV}} - \beta_0)} - 1 \right), \quad (4) $$

where $dA_{1600}/d\beta_{UV}$ depends on the frequency-dependence of the attenuation curve (and thus mostly on the dust grain size distribution and composition), and where the factor 1.5 is the bolometric correction of $L_{UV}$ and can be obtained by comparing Eq. (3) and Eq. (2), as $L_{IR}$ is approaching a bolometric luminosity (i.e., for dust $L_{IR}$ [8–1000 µm] $\sim L_{\text{bol}}[0-\infty]$).

The distribution in the IRX-$\beta_{UV}$ plane of the mean $\langle \text{IRX} \rangle$ and $\langle \beta_{UV} \rangle$ measurements for our mass-complete sample of SFGs at $4 < z < 5$ divided in three stellar mass bins is shown in the left panel of Fig. 2. This distribution is overall very consistent with the local reference for starbursts from Meurer et al. (1999, i.e., with $\beta_0 = -2.22$ and $dA_{1600}/d\beta_{UV} = 1.99$; see also Calzetti et al. 2000) and its updated version from Overzier et al. (2011, i.e., with $\beta_0 = -1.96$ and $dA_{1600}/d\beta_{UV} = 1.96$). On the contrary, our measurements disfavor a flatter IRX-$\beta_{UV}$ relation, such as that observed in the Small Magellanic Cloud (SMC; here parametrized with $\beta_0 = -2.22$ and $dA_{1600}/d\beta_{UV} = 1.10$, following Fudamoto et al. 2020a). Considering $\langle L_{UV}^{\text{phot.}} \rangle$ and $\langle \beta_{UV}^{\text{phot.}} \rangle$ instead of $\langle L_{UV} \rangle$ and $\langle \beta_{UV} \rangle$, and excluding instead of including the ALMA primary targets in our stacks, does not qualitatively change our results. These various measurements are indeed consistent with each other, except for our lowest stellar mass bin where the photometric estimates slightly shift toward redder UV spectral slopes. Due to these uncertainties, we are unfortunately unable to impose strict constraints on the value of $\beta_0$.

Overall, our results suggest that the grain size distribution and composition of the dust in $z \sim 4.5$ SFGs are very similar to those of local starbursts, extending to a higher redshift the
<figure><img src="image_3.png" /><figcaption>Fig. 2. Dust attenuation properties of SFGs at $z \sim 4.5$. Left: IRX-$\beta_{UV}$ relation derived at $4 < z < 5$ by stacking all SFGs with a $A^3\text{COSMOS}$ coverage (dark-green hexagons) or by stacking all SFGs except the ALMA primary targets (lime-green circles). Hexagons and circles with faded colors correspond to our $\langle I_{\text{UV}}^{\text{phot.}} \rangle$ and $\langle \beta_{\text{UV}}^{\text{phot.}} \rangle$ estimates, i.e., inferred from a COSMOS2020 photometry as opposed to their best-fit SED (see Sect. 4.2). For each stellar mass bin, the fraction of galaxies for which we were able to make such photometric measurement (see text for details) is given in the bottom part of the panel. In parenthesis, we provide the same number but after excluding the ALMA primary targets. The gray shaded squares present measurements for UV-selected SFGs at $z \sim 4-8$ from Bowler et al. (2024). The black solid, dash-dotted, dashed, and dotted lines show the IRX-$\beta_{UV}$ relation observed in local starbursts by Meurer et al. (1999) and by Calzetti et al. (2000), an updated version for local galaxies (Overzier et al. 2011), and an SMC-like dust attenuation relation, respectively. The orange dash-dot-dotted line shows the relation observed in a UV-selected sample of SFGs at $z \sim 4.5$ (Fudamoto et al. 2020a). Right: IRX-$M_*$ relation at $4 < z < 5$. Symbols are the same as in the left panel. Dark-green and lime-green solid lines are linear fits to these data points, while the associated shaded regions show the $1\sigma$ uncertainties on these fits (i.e., the 16th to 84th ranges). The pink dotted line shows the redshift-independent IRX-$M_*$ relation derived by Pannella et al. (2015) for $z \le 4$ SFGs. The black dash-dotted line shows the IRX-$M_*$ relation at $z \sim 4$ inferred by Koprowski et al. (2018). The blue dashed line shows the IRX-$M_*$ relation at $z \sim 3$ and $z \sim 4.5$ found in Fudamoto et al. (2020a,b).</figcaption></figure>

results for $z \sim 3$ SFGs obtained by Fudamoto et al. (2020b). The inconsistency with the rather flat IRX-$\beta_{UV}$ relation found at $z \sim 4.5$ in Fudamoto et al. (2020a) could lie on the fact that this latter study is based on a UV-selected sample of SFGs that missed a fraction of highly dust-attenuated galaxies. However, more recently, by re-analysing the ALPINE measurements used in Fudamoto et al. (2020a) and combining them with $z \sim 7$ measurements from REBEL, Bowler et al. (2024) find that using a consistent methodology, UV-selected SFGs at $z \sim 4-8$ appear to follow the same IRX-$\beta_{UV}$ as local starbursts, in perfect agreement with our findings.

The IRX-$M_*$ is less fundamental than the IRX-$\beta_{UV}$ relation, as it involves two quantities that are only indirectly related, the stellar mass being a proxy for the integrated past star formation and hence dust production. Therefore, unlike the IRX-$\beta_{UV}$ relation that mainly depends on the grain size distribution and composition of the dust, the IRX-$M_*$ relation also depends on its mass per unit of SFR and its geometry relative to the newly-formed stars. Strikingly, the IRX-$M_*$ relation is found to evolve very weakly from $z \sim 0$ to $z \sim 3$ (Pannella et al. 2015), despite the significant evolution in this redshift range in the star formation efficiency (SFR/$M_{\text{gas}}$, e.g., Liu et al. 2019b; Wang et al. 2022), metallicity (e.g., Bellstedt et al. 2021), and size of SFGs (e.g., van der Wel et al. 2014) at a given stellar mass (see also discussion in Shapley et al. 2022).

The distribution in the IRX-$M_*$ plane of the $\langle \text{IRX} \rangle$ and $\langle M_* \rangle$ measurements for our mass-complete sample of SFGs at $4 < z < 5$ is shown on the right panel of Fig. 2. Overall, this distribution (including or excluding the ALMA primary targets in our stacking analysis) is consistent with a steepening of the IRX-$M_*$ relation, our measurements being only marginally consistent with the $z = 1-4$ relation of Pannella et al. (2015) at low mass. These findings reinforce those of Fudamoto et al. (2020b), who

also find a Meurer-like IRX-$\beta_{UV}$ relation but a steeper IRX-$M_*$ relation in a mass-complete sample of SFGs at $z \sim 3$. These results are also very consistent with the $z \sim 4-8$ measurements of Bowler et al. (2024), obtained by combining the UV-selected sample of SFGs from the ALPINE and REBELS surveys. Therefore, although the dust in $z \sim 4.5$ SFGs appears to have the same composition as in local starbursts, its mass and geometry result in lower attenuation in $\le 10^{10} M_{\odot}$ SFGs at $z \sim 4.5$ than at $z \le 3$ (e.g., $-0.5$ dex in IRX at $M_* \sim 10^{9.75} M_{\odot}$). We note that at these low stellar masses, Shapley et al. (2023) do not find significantly lower dust attenuation as measured by the JWST using the Balmer decrement of SFGs, perhaps suggesting a difference between continuum and nebular reddening.

A linear fit to our data points yields the following IRX-$M_*$,

$$ \log_{10} \text{IRX} = 1.1^{+0.4}_{-0.4} \times \log_{10} \left( \frac{M_*}{10^{10} M_{\odot}} \right) + 0.67^{+0.13}_{-0.14}, \quad (5) $$

when including the ALMA primary targets in our stacks, or,

$$ \log_{10} \text{IRX} = 1.8^{+1.0}_{-0.7} \times \log_{10} \left( \frac{M_*}{10^{10} M_{\odot}} \right) + 0.45^{+0.30}_{-0.42}, \quad (6) $$

when excluding the ALMA primary targets. These two fits are very consistent overall, but the latter gives a less constrained and slightly steeper IRX-$M_*$ relation, due to the fact that by excluding the ALMA primary targets, our lowest stellar mass bin becomes an upper limit. In the rest of our analysis, we consider the whole set of IRX-$M_*$ relations found in these two fits (shaded regions in Fig. 2), and use as the fiducial IRX-$M_*$ relation at $z \sim 4.5$ the linear combination in the log-space of these two relations,

$$ \log_{10} \text{IRX} = 1.45 \times \log_{10} \left( \frac{M_*}{10^{10} M_{\odot}} \right) + 0.56. \quad (7) $$
## 5.2. The SFR–M* relation

In the SFR–M* plane, SFGs lie on a tight locus, known as the main sequence (MS) of SFGs and characterized by an almost linear increase in SFR with stellar mass (i.e., SFR ∝ M*; e.g., Noeske et al. 2007; Elbaz et al. 2007; Schreiber et al. 2015; Leslie et al. 2020; Daddi et al. 2022; Popesso et al. 2023; Goovaerts et al. 2024; Koprowski et al. 2024, and references therein). The small scatter of the MS suggests that secular evolution is the dominant mode of growth in SFGs, where gas inflow, outflow, and star formation are in equilibrium (e.g., Bouché et al. 2010; Davé et al. 2012; Lilly et al. 2013; Peng & Maiolino 2014; Rathaus & Sternberg 2016; Magnelli et al. 2020). Over the last decade, the MS has been a cornerstone of galaxy evolution studies, enabling the selection of “normal” SFGs at a given (z, M*) and the study of galaxy properties as a function of their position relative to the MS. While well established up to z ~ 3, the MS (i.e., its normalization, slope, and dispersion) is not well constrained at higher redshifts, even though it is commonly used to select normal SFGs at these epochs. Indeed, current mass-complete samples of SFGs, including those from JWST, lack systematic detections in the IR and thus rely on highly uncertain dust-attenuated SFRs, while even state-of-the-art stacking analysis with Herschel are limited to the most massive end of the SFG population (i.e., ≥10¹¹ M☉; e.g., Schreiber et al. 2015; Koprowski et al. 2024). These latter stacking estimates are also hampered by clustering signal, which has to be removed beforehand due to the coarse angular resolution of Herschel-SPIRE (>>14′′). Our ALMA stacks are not affected by such clustering bias thanks to their high angular resolution.

The distribution in the SFR–M* plane of the ⟨SFR⟩ and ⟨M*⟩ measurements for our mass-complete sample of SFGs at 4 < z < 5 is shown in Fig. 3. The measurements including the ALMA primary targets lie almost perfectly on the MS of Schreiber et al. (2015, i.e., MS<sub>All,Fit</sub> ≡ MS<sub>S15</sub>). This agreement is somewhat surprising as at these redshifts, Schreiber et al. (2015) are only able to probe the most massive SFGs (i.e., ≥10¹¹ M☉) with their Herschel stacks. In contrast to the previous section, the exclusion of the ALMA primary targets in our stacks has a significant impact on the inferred MS relation. Indeed, although these two measurements formally agree within their uncertainties, those obtained by excluding the ALMA primary targets give systematically lower SFRs; this tendency increasing at higher stellar mass where the weight of the ALMA primary targets in our stacks increases (see numbers at the bottom of Fig. 3). By fitting these latter measurements with the MS of Schreiber et al. (2015) leaving its normalization as a free parameter, we find MS<sub>Offcenter,Fit</sub> ≡ MS<sub>S15</sub> – 0.3 dex. Naturally, this difference stems from the fact that a fraction of the principal investigators of ALMA has targeted high-redshift starbursts (e.g., submillimeter-selected galaxies) or at least galaxies located on the upper envelop of the MS dispersion. Consequently, the SFR–M* relation obtained here by including the ALMA primary targets is probably slightly biased toward high values, while the SFR–M* relation obtained by excluding the ALMA primary targets is probably biased toward low values. In the rest of our analysis, we take into account this uncertainty on the normalization of the MS (i.e., ±0.15 dex) and use as our fiducial MS the linear combination in the log-space of these two relations, that is, MS<sub>fid</sub> ≡ MS<sub>S15</sub> – 0.15 dex.

## 5.3. The cosmic SFRD

Combining the results of the previous sections with the stellar mass function (SMF) of SFGs at z ~ 4.5, we can calcu-

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

**Fig. 3.** SFR–M* relation at 4 < z < 5. Dark-gray hexagons show the dust-attenuated SFRs (i.e., SFR<sub>IR</sub>) inferred by stacking all SFGs with a $A^3\text{COSMOS}$ coverage, while dark-green hexagons show their total SFRs, i.e., adding the contribution of their unattenuated SFRs measured in the UV (i.e., SFR<sub>UV</sub>; symbols are slightly offset along the stellar mass axis for clarity). Light-gray and lime-green circles display the same quantities but for our stacking analysis which excludes the ALMA primary targets (symbols are slightly offset along the stellar mass axis for clarity). The dark-green solid and red dotted lines represent the MS and its ±0.3 dex dispersion from Schreiber et al. (2015), whose normalization happens to be perfectly consistent with our dark-green hexagons, i.e., MS<sub>All,Fit</sub> ≡ MS<sub>S15</sub>. The lime-green line is obtained by renormalizing the MS of Schreiber et al. (2015) to fit our lime-green circles, i.e., MS<sub>Offcenter,Fit</sub> ≡ MS<sub>S15</sub> – 0.3 dex. For each stellar mass bin, we give the fraction of ALMA primary targets and in parenthesis their weights in our stacks. Finally, pink squares and diamonds show the total and dust-attenuated SFRs found in Khusanova et al. (2021) by stacking in three stellar mass bins, ALMA observations of all (detected and undetected) z ~ 4.5 SFGs from the ALMA ALPINE survey.

late the dust-attenuated and unattenuated cosmic SFRD at this epoch and the evolution of these quantities with stellar mass. In the top left panel of Fig. 4, we display the SMF of SFGs used in our analysis and inferred in Weaver et al. (2023). The choice of this particular SMF was relatively straightforward as it is arguably the best estimate available to date and, more importantly, it was measured from the same parent sample as that used here. As shown in Fig. 4 and also noted in Weaver et al. (2023), at the massive end, their data points deviate, however, from a canonical Schechter function (Schechter 1976) and thus from their best fit. After an in-depth analysis of this population of massive SFGs, Weaver et al. (2023) conclude that their luminosity, redshift, mass, and number density are not compatible with AGN-dominated galaxies, but rather with the population of optically-dark galaxies recently discovered thanks to Spitzer-IRAC and ALMA observations (e.g., Wang et al. 2019; Xiao et al. 2023). They argue that the $K_s$ ~ 24.7 AB depth of the UltraVISTA DR4 used to build the izYJHK<sub>s</sub> detection image on which the COSMOS2020 catalog was selected could indeed be sufficient to reach out into this population missed by previous optical or even H-band selections. To account analytically for this massive population, we fit simultaneously all the SMF data points (3.5 < z < 4.5 and 4.5 < z < 5.5) of Weaver et al. (2023) with the combination of a Schechter function dominating at low stellar masses and a power-law dominating at high stellar masses. The result of this fit and its associated uncertainties is represented by the solid line and the gray shaded region in the top-left panel of Fig. 4. In the rest of the analysis, we use this fit
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Evolution with stellar mass of the cosmic SFRD at $4 < z < 5$. **Top left:** fiducial SMF used in our calculations is represented by the gray line and shaded region. This is the combination of a Schechter and power-law function fit to the SMF of Weaver et al. (2023, gray and black circles). A Schechter function fit to these data points is shown by the gray dotted line. **Bottom left:** cosmic SFRD (thick black line) and dust-attenuated SFRD (red line) above a given stellar mass, as inferred using the fiducial SMF. The dotted gray and red lines present the same quantities but using a simple Schechter SMF in our calculations. **Top right:** fraction of the $M_{*} > 10^{8} M_{\odot}$ SFRD (gray histogram) and dust-attenuated SFRD (red histogram) that is attributed to a given stellar mass bin. **Bottom right:** fraction of the cosmic SFRD above a given stellar mass that is dust-attenuated (red region) and unattenuated (blue region). The dotted gray line present the same quantity but using a simple Schechter SMF in our calculations. In all panels, the vertical line represents the $0.03 \times M^{*}$ lower integration limits commonly used to infer the “total” cosmic SFRD (e.g., Madau & Dickinson 2014). All quantities in the bottom-left, top-right, and bottom-right panels were calculated using our fiducial SMF, IRX-$M_{*}$, MS relations. The propagation of the $1-\sigma$ uncertainties on the SMF and the IRX-$M_{*}$ relation are represented by shaded regions. The uncertainties on the MS would move the black and red lines of the bottom-left panel up and down by 0.15 dex, but would have no effect on the other quantities displayed in this figure.</figcaption></figure>

as our fiducial SMF, but also discuss the impact of using instead the Schechter fit of Weaver et al. (2023) on our calculations.

First, we calculated the cosmic SFRD above a given stellar mass (i.e., SFRD($>M_{*}$)) and the contribution of each stellar mass bin to SFRD($>10^{8} M_{\odot}$) by multiplying and then integrating our fiducial SMF with our fiducial MS (bottom-left and top-right panels of Fig. 4, respectively). Then, using our fiducial IRX-$M_{*}$ relation, we calculated the cosmic dust-attenuated SFRD above a given stellar mass (i.e., SFRD$_{IR}$($>M_{*}$)) and the contribution of each stellar mass bin to SFRD$_{IR}$($>10^{8} M_{\odot}$; bottom-left and top-right panels of Fig. 4, respectively). Finally, we calculated the fraction of the cosmic SFRD above a given stellar mass that is dust attenuated ($f_{\text{attenuated}}$; bottom-right panel of Fig. 4). In these calculations, the complex propagations of the uncertainties on the SMF and the IRX-$M_{*}$ relation were inferred using Monte Carlo realizations and are represented by shaded regions. In contrast, using MS$_{All,Fit}$ or MS$_{Offcenter,Fit}$ instead of

our fiducial MS simply move the black and red lines of the bottom-left panel of Fig. 4 (i.e., SFRD($>M_{*}$) and SFRD$_{IR}$($>M_{*}$)) up or down by 0.15 dex, respectively, but would have no effect on the other quantities displayed in this figure. Finally, we note that our calculations does not account for the contribution of star-bursts (i.e., SFGs located $>0.6$ dex above the MS) to the cosmic SFRD, since by construction we have assumed that at a given mass, all SFGs are MS galaxies. However, this rare population of galaxies is known to have a modest ~10% contribution to the SFRD (Rodighiero et al. 2011; Schreiber et al. 2015, i.e., at most +0.05 dex to our estimates) and part of this contribution is actually taken into account by our MS$_{All,Fit}$ relation that is supposably biased towards this population (see Sect. 5.2).

Three main conclusions can be drawn from Fig. 4. Firstly, the combination of the MS and SMF slopes implies that the cosmic SFRD converges at low stellar masses: the number of SFGs at low stellar mass does not increase sufficiently to counterbalance
the decrease in their SFRs. This is highlighted by the fact that, based on our extrapolations, SFGs at $10^8 < M_*/M_\odot < 10^9$ already contribute less than 15% of the SFRD($>10^8 M_\odot$). The steep slope of the IRX-$M_*$ relation makes this result even more pronounced for SFRD<sub>IR</sub>($>M_*$), with SFGs at $10^8 < M_*/M_\odot < 10^9$ contributing less than 5% of the SFRD<sub>IR</sub>($>10^8 M_\odot$). This implies that the study of SFGs much less massive than $\le 10^9 M_\odot$ is almost irrelevant for our understanding of the cosmic SFRD at $z \sim 4.5$. Secondly, the cosmic SFRD is mainly dominated by SFGs with a stellar mass of $10^{9.5-10.5} M_\odot$: galaxies that are more massive and therefore more star-forming, are too rare to contribute significantly to the SFRD; while if less massive galaxies are numerous, they do not form enough stars to contribute significantly to the SFRD. This characteristic stellar mass of the SFG population is consistent with that observed at lower redshifts ($\sim 10^{10.3} M_\odot$; e.g., Karim et al. 2011; Leslie et al. 2020). Thirdly, the fraction of the cosmic SFRD that is attenuated by dust remains significant even at such an early cosmic epoch, with $f_{\text{attenuated}} \equiv \text{SFRD}_{\text{IR}}(>M_*)/\text{SFRD}(>M_*)$ converging to a value of $\sim 60 \pm 10\%$ for $M_* = 10^8 M_\odot$.

Using a Schechter SMF instead of our fiducial SMF would not qualitatively change any of these results (see dotted lines in Fig. 4). Indeed, while a massive population is responsible for a large deviation of our fiducial SMF from a Schechter function, its number density remains too low to have a significant impact on the global star formation activity in the Universe. In fact, excluding this population from our calculation shifts the characteristic mass of the SFG population to lower values by only ~0.3 dex and decrease the cosmic SFRD by ~0.1 dex, this latter value being in perfect agreement with the finding of, for example, Wang et al. (2019) for optically-dark galaxies. As these variations are negligible compared to those introduced by uncertainties on the MS and IRX-$M_*$ relations, in what follows we restrict our calculation to our fiducial SMF (i.e., accounting for this massive population of SFGs).

Finally, as detailed in Appendix B, we also tested the robustness of our calculations against the most recent JWST-derived SMF from Weibel et al. (2024). Because the SMF of Weibel et al. (2024) differs from that of Weaver et al. (2023) only at $M_* \le 10^9 M_\odot$ where it has a slightly steeper slope ($-1.79$ versus $-1.56$), our results remain qualitatively unchanged using this alternative SMF: the convergence of the cosmic SFRD occurs but at lower stellar mass ($10^7 M_\odot$ vs. $10^8 M_\odot$); at the convergence, the fraction that is dust attenuated is lower but still significant (43% vs. 60%); the contribution of very low stellar mass galaxies increases but remains sub-dominant (35% vs. 15% for $M_* < 10^9 M_\odot$); while the “total” cosmic SFRD (i.e. at $10^{8.9} M_\odot$; see Sect. 5.4) and the fraction of it that is attenuated by dust remains mostly unchanged from those deduced using the COSMOS2020-derived SMF of Weaver et al. (2023). Because of these very good agreements over the stellar mass range for which we are able to measure the dust-attenuated properties of $z \sim 4.5$ SFGs, and because using the COSMOS2020-derived SMF ensure consistency with our stack samples, we have decided to retain the COSMOS2020-derived SMF from Weaver et al. (2023) as our fiducial SMF. We defer to future work the possibility of extending the constraints on the dust-attenuation of SFGs to lower stellar masses using JWST-based catalogs.

## 5.4. The “total” cosmic SFRD

In the literature, the IR and UV luminosity functions are commonly integrated down to $0.03 \times L^*$ to infer the “total” cosmic

SFRD (e.g., Madau & Dickinson 2014), where $L^*$ is the characteristic luminosity of the Schechter function. The slope of our fiducial MS being equal to one (Schreiber et al. 2015), integrating the SMF down to $0.03 \times M^*$ (i.e., $10^{8.9} M_\odot$; vertical lines in Fig. 4) should provide consistent results. In Fig. 5, we compare the total cosmic SFRD and dust-attenuated SFRD at $z \sim 4.5$ calculated here by integrating the SMF down to this limit with the literature results. The symbols correspond to the calculation using our fiducial SMF, MS and IRX-$M_*$ relations, while the range on the cosmic SFRD corresponds to the $\pm 0.15$ dex uncertainties on the MS relation (i.e., excluding or not the ALMA primary targets) and the range on the cosmic dust-attenuated SFRD takes into account the uncertainties on both the MS relation (i.e., excluding or not the ALMA primary targets) and the IRX-$M_*$ relations. Hereafter these total cosmic SFRD and total cosmic dust-attenuated SFRD are referred as SFRD<sub>UV+IR</sub> and SFRD<sub>IR</sub>, respectively. We find $\text{SFRD}_{\text{UV+IR}} = 0.012-0.028 M_\odot \text{yr}^{-1} \text{Mpc}^{-3}$ and $\text{SFRD}_{\text{IR}} = 0.007-0.021 M_\odot \text{yr}^{-1} \text{Mpc}^{-3}$.

Our estimate of SFRD<sub>UV+IR</sub> is consistent with Madau & Dickinson (2014), in particular the upper value of our range, which corresponds to $MS_{\text{All,Fit}} \equiv MS_{\text{S15}}$. Our estimate of SFRD<sub>IR</sub> is also in agreement with previous measurements based on Herschel (Liu et al. 2018), SCUBA (Maniyar et al. 2018) or ALMA (e.g., Zavala et al. 2021; Traina et al. 2024).

In contrast, our estimate of SFRD<sub>IR</sub> disagrees with the recent ALMA measurements by Gruppioni et al. (2020) and Khusanova et al. (2021), both of which found large SFRD<sub>IR</sub> values, about two times and four times larger than the SFRD<sub>UV+IR</sub> and SFRD<sub>IR</sub> inferred in Madau & Dickinson (2014), respectively. This disagreement (although our estimates are formally consistent within the uncertainties) could naturally come from our sides, owing to the fact that our method is based on a near-infrared-selected sample (i.e., COSMOS2020) that may miss a population of heavily dust-attenuated SFGs only recoverable using an ALMA-selected sample such as that used in, e.g., Gruppioni et al. (2020). The existence of a population of optically-dark galaxies (e.g., Wang et al. 2019; Xiao et al. 2023) supports this view. However, the impact of this population, while undeniable, remains modest with the most recent estimates ranging from SFRD<sub>IR</sub> ~0.002 $M_\odot \text{yr}^{-1} \text{Mpc}^{-3}$ to ~0.012 $M_\odot \text{yr}^{-1} \text{Mpc}^{-3}$ (e.g. Xiao et al. 2023; Barrufet et al. 2023; Williams et al. 2024). In addition, their impact on our measurements would be further mitigated because at the $K_S \sim 24.7$ AB depth of the UltraVISTA DR4 used to build the izYJHK<sub>s</sub> detection image, the COSMOS2020 catalog should detect part of this population (Weaver et al. 2023). Therefore, while the inclusion of optically-dark galaxies could lead to an increase of our SFRD<sub>UV+IR</sub> and SFRD<sub>IR</sub> measurements by at most ~0.012 $M_\odot \text{yr}^{-1} \text{Mpc}^{-3}$ (+0.25 dex compared to our fiducial SFRD<sub>IR</sub>), such increase is well within the uncertainties associated with our measurements and not enough to explain the disagreements (~+0.7 dex compared to our fiducial SFRD<sub>IR</sub> value) with Gruppioni et al. (2020) and Khusanova et al. (2021). Future ALMA stacks of JWST-selected SFGs at $z \sim 4.5$ will enable us to assess the exact contribution of this population of galaxies to the cosmic SFRD.

Instead, we argue that the large SFRD<sub>IR</sub> values found in Gruppioni et al. (2020) and Khusanova et al. (2021) could be due to low number statistic and extrapolation of the IRX-$M_*$ relation, respectively. On the one hand, using the same approach as Gruppioni et al. (2020) but using a much larger number of ALMA observations drawn from the $A^3\text{COSMOS}$ database, Traina et al. (2024) measure a SFRD<sub>IR</sub> that is more consistent with our measurements. This suggests that a larger statistic
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Redshift evolution of the “total” cosmic SFRD. The cosmic SFRD (UV+IR) and dust-attenuated SFRD (IR) derived in this work are shown by black and orange circles, respectively (offset slightly along the redshift axis for clarity). To be consistent with Madau & Dickinson (2014), these measurements were obtained by integrating the SMF down to $0.03 \times M^*$ (i.e., $10^{8.9} M_{\odot}$). Orange open symbols are a collection of dust-attenuated SFRDs (i.e., SFRD<sub>IR</sub>) from the literature (Sanders et al. 2003; Takeuchi et al. 2003; Schiminovich et al. 2005; Dahlen et al. 2007; Reddy & Steidel 2009; Magnelli et al. 2011, 2013; Cucciati et al. 2012; Bouwens et al. 2012; Gruppioni et al. 2013; Finkelstein et al. 2015; Bouwens et al. 2015; Liu et al. 2018; Khusanova et al. 2020, 2021; Zavala et al. 2021; Algera et al. 2023; Traina et al. 2024). Light blue open symbols are a collection of unattenuated SFRDs (i.e., SFRD<sub>UV,uncorr</sub>) from the literature, while dark blue open symbols correspond to the same measurements after correction for dust attenuation (i.e., SFRD<sub>UV,corr</sub>). The dark-gray line is the redshift evolution of the cosmic SFRD inferred in Madau & Dickinson (2014). Orange dotted and dashed lines show the redshift evolution of the dust-attenuated SFRD inferred in Koprowski et al. (2017) and Maniyar et al. (2018), respectively.</figcaption></figure>

and the combination of Herschel and ALMA in Traina et al. (2024) have lowered the inferred IR luminosity function and, subsequently, SFRD<sub>IR</sub> at $z \sim 4.5$. On the other hand, for $\ge 10^{8.35} M_{\odot}$ where they have direct constraints on the $L_{IR}-M_{*}$ relation, Khusanova et al. (2021) inferred a SFRD<sub>IR</sub> of $0.023 M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$, very consistent with our measurements. This suggests that the disagreement with Khusanova et al. (2021) stems from their extrapolation to lower stellar mass, as their total cosmic SFRD<sub>IR</sub> is for SFGs with $>10^6 M_{\odot}$. The reason for these very low integration limits and for the increase of SFRD<sub>IR</sub> from $10^{8.35} M_{\odot}$ to $10^6 M_{\odot}$ when our measurements have already converged at these masses (see Fig. 4) is threefold: Khusanova et al. (2021) used a MS relation with a sub-unity slope of ~0.85, which implies that the $0.03 \times L^*$ integration limit is converted to a very low-mass limit of $10^6 M_{\odot}$; Next, Khusanova et al. (2021) used a $L_{IR}-M_{*}$ relation with a sub-unity slope of ~0.85, which implies an unrealistic mass-independent value for IRX; Finally, all these elements combined imply that as stellar mass decreases, the SFR of their

SFGs decreases more slowly than in our case (flatter MS relation) and more importantly the dust-attenuated fraction of these SFRs remains significant. In brief, an flat IRX-$M_{*}$ relation combined with a shallow MS relation result in unrealistic extrapolations to lower stellar masses and an increase of SFRD<sub>IR</sub>($>M_{*}$) from $0.023 M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ for $M_{*} = 10^{8.35} M_{\odot}$ to $0.085 M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ for $M_{*} = 10^6 M_{\odot}$.

Finally, in Fig. 6, we present the cosmic dust attenuation in magnitudes, $A_{UV} \equiv 2.5 \times \log_{10}(1 + \text{SFRD}_{\text{IR}}/\text{SFRD}_{\text{UV}})$ (see, e.g., Madau & Dickinson 2014) and compare it to various measurements from the literature. Our estimate that about ~68% of the total cosmic SFRD is attenuated by dust at $M_{*} \ge 10^{8.9} M_{\odot}$ (and thus $A_{UV} = 1.14$) is in line with previous results from the literature, which range from ~35% (Khusanova et al. 2020), ~45% (Zavala et al. 2021), to ~70% (Finkelstein et al. 2015). The fraction of the total cosmic SFRD that is attenuated by dust remains thus significant at this early cosmic epoch, even though this fraction has clearly decreased from $z \sim 1-2$, where it peaks, to $z \sim 4.5$.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Cosmic dust attenuation in magnitudes ($A_{UV} \equiv 2.5 \times \log_{10}(1 + \text{SFRD}_{\text{IR}}/\text{SFRD}_{\text{UV}})$) as a function of redshift. The cosmic dust attenuation inferred in this work is shown by an orange circle. Blue and orange symbols are a collection of cosmic dust attenuation measurements from the literature.</figcaption></figure>

## 6. Summary

Using the wealth of multiwavelength observations publicly available over the COSMOS field, we measured for the first time the mean UV and IR emission properties of a mass-complete sample of SFGs at $4 < z < 5$ divided in three stellar mass bins, that is, $10^{9.5} < M_*/M_{\odot} < 10^{10}$, $10^{10} < M_*/M_{\odot} < 10^{10.5}$, and $10^{10.5} < M_*/M_{\odot} < 10^{11.5}$. We built this mass-complete sample of SFGs and measured their rest-frame UV properties (i.e., $\langle L_{UV} \rangle$, $\langle \beta_{UV} \rangle$, and $\langle \text{SFR}_{UV} \rangle$) using the deepest optical-to-near-IR data publicly available in COSMOS (i.e., the COSMOS2020 catalog; Weaver et al. 2022). Their mean IR properties (i.e., $\langle L_{IR} \rangle$ and $\langle \text{SFR}_{IR} \rangle$) were then measured by stacking in the uv domain all archival ALMA band 6 and 7 observations publicly available for these galaxies (A$^3$COSMOS; Liu et al. 2019a; Adscheid et al. 2024). With this unique approach, we find the following:

1. The relation between the IR excess ($\text{IRX} \equiv L_{\text{IR}}/L_{\text{UV}}$) and the UV spectral slope ($\beta_{UV}$) in our mass-complete sample of SFGs at $4 < z < 5$ is overall very consistent with the local IRX-$\beta_{UV}$ relation from Meurer et al. (1999) and Calzetti et al. (2000), and disfavors a flatter IRX-$\beta_{UV}$ relation, such as that observed in the SMC. This result, consistent with the most recent $z \sim 4-8$ measurements from Bowler et al. (2024), suggests that the grain size distribution and composition of the dust in SFGs with $M_* > 10^{9.5} M_{\odot}$ and at $z \sim 4.5$ are very similar to those of local starbursts.

2. Our measurements favor a slight steepening of the IRX-$M_*$ relation at $z \sim 4.5$, when compared to the redshift-independent IRX-$M_*$ relation observed at $z \sim 1-3$ (e.g., Pannella et al. 2015). Thus, although the dust in $z \sim 4.5$ SFGs appears to have the same composition than in local starbursts, its mass and geometry result in lower attenuation in $\le 10^{10} M_{\odot}$ SFGs at $z \sim 4.5$ than at $z \le 3$ (e.g., $-0.5$ dex in IRX at $M_* \sim 10^{9.75} M_{\odot}$).

3. In the SFR-$M_*$ plane, our galaxies lie almost perfectly on the MS of Schreiber et al. (2015), while they lie 0.3 dex below this relation when excluding from our stacks the ALMA primary targets (i.e., sources within 3'' from the ALMA phase center). ALMA primary targets are probably slightly biased toward galaxies located on the upper envelope of the MS dispersion. We set our fiducial MS to $MS_{S15} - 0.15$ dex.

4. The combination of the MS and SMF slopes implies that the SFRD($M_*$) converges at $M_* \lesssim 10^9 M_{\odot}$, as the number of SFGs at lower stellar mass does not increase sufficiently to counterbalance the decrease in their SFRs. For example, SFGs with $10^8 < M_*/M_{\odot} < 10^9$ contribute already less than 15% of the SFRD($>10^8 M_{\odot}$), and less than 5% of the SFRD$_{IR}$($>10^8 M_{\odot}$).

5. The “total” cosmic SFRD inferred here at $z \sim 4.5$ for SFGs with $M_* > 0.03 \times M^*$ ($M^*$ being the characteristic stellar mass of SFGs at this epoch) is consistent with Madau & Dickinson (2014) and is dominated by SFGs with a stellar mass of $10^{9.5-10.5} M_{\odot}$. At $z \sim 4.5$, SFRD$_{UV+IR} = 0.012-0.028 M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ and SFRD$_{IR} = 0.007-0.021 M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$. The population of optical-dark galaxies potentially missed by our study could add at most $0.012 M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ to these values (Xiao et al. 2023; Barrufet et al. 2023; Williams et al. 2024). Future ALMA stacks of JWST-selected SFGs at $z \sim 4.5$ will assess the exact contribution of this population of galaxies to the cosmic SFRD.

6. The fraction of the cosmic SFRD that is attenuated by dust remains significant even at such an early cosmic epoch, with $f_{\text{attenuated}} \equiv \text{SFRD}_{\text{IR}}(> M_*)/\text{SFRD}(> M_*)$ having a value of $68 \pm 10\%$ for $M_* \ge 10^{8.9} M_{\odot}$ (i.e., $0.03 \times M^*$) and converging to a value of $60 \pm 10\%$ for $M_* \ge 10^8 M_{\odot}$.

**Acknowledgements.** We would like to thank the referee for their comments that have helped to improve our paper. BM, DE, and CGG acknowledges support from CNES. S.A. gratefully acknowledges the Collaborative Research Center 1601 (SFB 1601 sub-project C2) funded by the Deutsche Forschungsgemeinschaft (DFG. German Research Foundation) - 500700252. E.F.-J.A. acknowledges support from UNAM-PAPIIT project IA102023, and from CONAHCYT Ciencia de Frontera project ID: CF-2023-I-506. MF acknowledges support from NSF grant AST-2009577 and NASA JWST GO program 1727. CG acknowledges the support from grant PRIN MIUR 2017-20173ML3WW_001 “Opening the ALMA window on the cosmic evolution of gas, stars, and supermassive black hole”, and co-funding by the European Union – NextGenerationEU within PRIN 2022 project n.20229YBSAN “Globular clusters in cosmological simulations and in lensed fields: from their birth to the present epoch”. BM acknowledges the following open source software used in the analysis: Astropy (Astropy Collaboration 2022), photutils (Bradley et al. 2022), and NumPy (Harris et al. 2020).

## 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

Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJ, 946, L12

Barrufet, L., Oesch, P. A., Weibel, A., et al. 2023, MNRAS, 522, 449

Bellstedt, S., Robotham, A. S. G., Driver, S. P., et al. 2021, MNRAS, 503, 3309

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

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

Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJ, 754, 83

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

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

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, https://zenodo.org/records/6825092

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

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

Chabrier, G. 2003, PASP, 115, 763

Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, A&A, 539, A31

Daddi, E., Delvecchio, I., Dimauro, P., et al. 2022, A&A, 661, L7

Dahlen, T., Mobasher, B., Dickinson, M., et al. 2007, ApJ, 654, 172

Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98

Eisenstein, D.J., Willott, C., Alberts, S., et al. 2023, ApJS, submitted [arXiv:2306.02465]
Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33

Finkelstein, S. L., Ryan, R. E., Jr, Papovich, C., et al. 2015, ApJ, 810, 71

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

Fudamoto, Y., Oesch, P. A., Magnelli, B., et al. 2020b, MNRAS, 491, 4724

Goovaerts, I., Pello, R., Burgarella, D., et al. 2024, A&A, 683, A184

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8

Harris, C. R., Millman, K. J., Walt, S. J. V. D., et al. 2020, Nature, 585, 357

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Karim, A., Schinnerer, E., Martínez-Sansigre, A., et al. 2011, ApJ, 730, 61

Khusanova, Y., Le Fèvre, O., Cassata, P., et al. 2020, A&A, 634, A97

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

Koprowski, M. P., Dunlop, J. S., Michałowski, M. J., et al. 2017, MNRAS, 471, 4155

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

Koprowski, M. P., Wijesekera, J. V., Dunlop, J. S., et al. 2024, A&A, submitted [arXiv:2403.06575]

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58

Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119

Lindroos, L., Knudsen, K. K., Vlemmings, W., Conway, J., & Martí-Vidal, I. 2015, MNRAS, 446, 3502

Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172

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

Liu, D., Schinnerer, E., Groves, B., et al. 2019b, ApJ, 887, 235

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35

Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132

Magnelli, B., Boogaard, L., Decarli, R., et al. 2020, ApJ, 892, 66

Maniyar, A. S., Béthermin, M., & Lagache, G. 2018, A&A, 614, A39

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

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43

Overzier, R. A., Heckman, T. M., Wang, J., et al. 2011, ApJ, 726, L7

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

Peng, Y.-J., & Maiolino, R. 2014, MNRAS, 443, 3643

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Rathaus, B., & Sternberg, A. 2016, MNRAS, 458, 3168

Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778

Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJ, 739, L40

Sanders, D. B., Mazzarella, J. M., Kim, D. C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607

Schechter, P. 1976, ApJ, 203, 297

Schiminovich, D., Ilbert, O., Arnouts, S., et al. 2005, ApJ, 619, L47

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Shapley, A. E., Sanders, R. L., Salim, S., et al. 2022, ApJ, 926, 145

Shapley, A. E., Sanders, R. L., Reddy, N. A., Topping, M. W., & Brammer, G. B. 2023, ApJ, 954, 157

Takeuchi, T. T., Yoshikawa, K., & Ishii, T. T. 2003, ApJ, 587, L89

Traina, A., Gruppioni, C., Delvecchio, I., et al. 2024, A&A, 681, A118

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28

Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Nature, 572, 211

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

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11

Weaver, J. R., Davidzon, I., Toft, S., et al. 2023, A&A, 677, A184

Weibel, A., Oesch, P. A., Barrufet, L., et al. 2024, MNRAS, submitted [arXiv:2403.08872]

Williams, C. C., Alberts, S., Ji, Z., et al. 2024, ApJ, 968, 34

Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, ApJ, 742, 96

Xiao, M. Y., Elbaz, D., Gómez-Guijarro, C., et al. 2023, A&A, 672, A18

Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, ApJ, 909, 165

<footnote n="1">¹ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France
e-mail: benjamin.magnelli@cea.fr</footnote>

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

<footnote n="3">³ Aix-Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="4">⁴ INAF – Osservatorio Astronomico di Brera 28, 20121, Milano, Italy and Via Bianchi 46, 23807 Merate, Italy</footnote>

<footnote n="5">⁵ Waseda Research Institute for Science and Engineering, Faculty of Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan</footnote>

<footnote n="6">⁶ National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo, Japan</footnote>

<footnote n="7">⁷ Center for Frontier Science, Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan</footnote>

<footnote n="8">⁸ Department of Pure and Applied Physics, Graduate School of Advanced Science and Engineering, Faculty of Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan</footnote>

<footnote n="9">⁹ The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA</footnote>

<footnote n="10">¹⁰ Istituto Nazionale di Astrofisica (INAF) – Osservatorio di Astrofisica e Scienza dello Spazio (OAS), Via Gobetti 101, 40129 Bologna, Italy</footnote>

<footnote n="11">¹¹ Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Antigua Carretera a Pátzcuaro #8701, Ex-Hda. San José de la Huerta, Morelia, Michoacán 58089, Mexico</footnote>

<footnote n="12">¹² Max-Planck-Institut für extraterrestrische Physik (MPE), Giessen-bachstrasse 1, 85748 Garching, Germany</footnote>

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

<footnote n="14">¹⁴ Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland</footnote>

<footnote n="15">¹⁵ Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen, Jagtvej 128, København N 2200, Denmark</footnote>

<footnote n="16">¹⁶ Max Planck Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="17">¹⁷ Dipartimento di Fisica e Astronomia (DIFA), Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy</footnote>
# Appendix A: Representativeness of the stacked samples

As our stacked samples include only a small fraction of the sources in our parent sample of $z \sim 4.5$ SFGs (those with an ALMA coverage), their representativeness may not be optimal. To check this, we compared the distributions of key physical properties (redshift, stellar mass, UV luminosity, and UV spectral slope) in our parent and stacked samples, i.e., all galaxies with an ALMA coverage ("All") and galaxies with an ALMA coverage but excluding the ALMA primary targets ("Off"; Fig. A.1). To highlight the potential bias of the ALMA primary target subsample, we also examined its distributions, although this subsample has not been stacked separately in our analysis.

For all these physical parameters, our stacked samples ("All" and "Off") have very similar distributions to our parent sample. A two-sample Kolmogorov-Smirnov (KS) analysis indicates that our "All" and "Off" samples are both consistent with being randomly drawn from the same distribution as our parent sample ($p$-value > 0.05). Moreover, the similarity between our parent and the "Off" samples is consistently better than between our parent and the "All" samples. This is naturally explained by the fact that the "All" sample contains the small subsample of ALMA primary targets (64 out of 440), which can, with high confidence, be ruled out as being randomly drawn from our parent sample (p-value < 0.05). As expected, these ALMA primary targets are biased toward bright, massive galaxies, in a redshift range favorable for [CII] observations. However, this subsample does not appear to be biased toward galaxies with a particular UV spectral slope. This last finding probably explains why there is no significant difference between the $IRX-\beta_{UV}$ and $IRX-M_*$ relations inferred by including or excluding these ALMA primary targets from our stacks.

In short, this analysis reinforces our choice to perform our stacking analysis by including and excluding these ALMA primary targets, as these latter are slightly biased towards bright, massive galaxies. Excluding these ALMA primary targets gives us a fairly representative sample of our parent sample of $z \sim 4.5$ SFGs.

<figure><img src="image_8.png" /><figcaption>**Fig. A.1.** Representativeness of the stacked samples. Normalized distribution of some key physical properties, i.e., redshift (top left), stellar mass (top right), UV luminosity (bottom left), and UV spectral slope (bottom left) in our parent sample (green; "Par."), our stacked samples, including (red; "All") and excluding the ALMA primary targets ("Off"), and the sub-sample of only the ALMA primary targets (yellow; "Targ."). In each panel, we provide the *p*-value from the two-sample Kolmogorov-Smirnov (KS) analysis performed between these different samples and our parent sample.</figcaption></figure>
## Appendix B: An alternative fiducial SMF

Here we test the robustness of our calculations against the most recent SMF of Weibel et al. (2024), derived by combining the public JWST/NIRCam imaging programs from CEERS (Bagley et al. 2023), PRIMER (PI: J. Dunlop), and JADES (Eisenstein et al. 2023), covering a total area of ~500 arcmin². Using this alternative SMF, we repeated the analysis performed in Sect. 5.3, that is, combining it with our IRX-$M_*$ and MS relations to calculate the dust-attenuated and unattenuated cosmic SFRD at $z \sim 4.5$ and the evolution of these quantities with stellar mass. The result of this analysis is shown in Fig. B.1.

In the common stellar mass range, there is very good agreement in the SMFs of Weaver et al. (2023) and Weibel et al. (2024). In particular, as in Weaver et al. (2023), the JWST-derived SMF deviates at the massive end from the canonical Schechter function and requires a combination of a Schechter function and a power-law function to be fit. Moreover, such fit yields to a characteristic mass, $M^*$, of $10^{10.4} M_\odot$, which is in perfect agreement with that deduced from the Weaver et al. (2023) data points. In fact, the only significant difference between the SMF of Weibel et al. (2024) and that of Weaver et al. (2023) appears at stellar masses not probed by the COSMOS2020 catalog (i.e., at $\lesssim 10^9 M_\odot$), where the JWST-

derived SMF has a slightly steeper faint-end slope, -1.79 versus -1.56.

The very good agreement between the SMF of Weibel et al. (2024) and that of Weaver et al. (2023), especially down to $0.03 \times M^*$ where the “total” cosmic SFRD is commonly defined, implies that our results remain qualitatively unchanged using this alternative SMF: the convergence of the cosmic SFRD occurs but at lower stellar mass ($10^7 M_\odot$ vs. $10^8 M_\odot$); at the convergence, the fraction that is dust attenuated is lower but still significant (43% vs. 60%); the contribution of very low stellar mass galaxies increases but remains sub-dominant (35% vs. 15% for $M_* < 10^9 M_\odot$); while the “total” cosmic SFRD (i.e. at $M_* = 10^{8.9} M_\odot$) and the fraction of it that is attenuated by dust remains mostly unchanged from those deduced using the COSMOS2020-derived SMF of Weaver et al. (2023).

Because this alternative SMF does not alter our results over the range of stellar mass for which we are able to measure the dust-attenuated properties of $z \sim 4.5$ SFGs (i.e., $\gtrsim 10^9 M_\odot$), we have decided to retain as our fiducial SMF the COSMOS2020-derived SMF from Weaver et al. (2023). This ensures also consistency between the SMF and our stack samples. We defer to future work the possibility of extending the constraints on the dust-attenuation of SFGs to lower stellar masses using JWST-based catalogs.

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

**Fig. B.1.** Evolution with stellar mass of the cosmic SFRD at $4 < z < 5$. This evolution was inferred using the JWST-derived SMF of Weibel et al. (2024, red and orange circles) instead of the COSMOS2020-derived SMF of Weaver et al. (2023). The symbols are the same as in Fig. 4. However, to allow a direct comparison with the results displayed in Fig. 4, in the top left panel we show with black and gray circles the SMF of Weaver et al. (2023), in the bottom left panel we show with black and red circles the cosmic SFRD and dust-attenuated SFRD at $0.03 \times M^*$ (i.e., “total” cosmic SFRDs) measured using the SMF of Weaver et al. (2023), and in the bottom right panel we show with the black line the fraction of the SFRD that is attenuated measured using the SMF of Weaver et al. (2023). Also, compared to Fig. 4, the x-axis has been extended down to $10^7 M_\odot$ and the top right panel is now the fraction of the $10^7 M_\odot$ SFRD.