# HSC-CLAUDS survey: The star formation rate functions since z ~ 2 and comparison with hydrodynamical simulations

V. Picouet, S. Arnouts, E. Le Floc'H, T. Moutard, K. Kraljic, O. Ilbert, M.

Sawicki, G. Desprez, C. Laigle, D. Schiminovich, et al.

► To cite this version:

V. Picouet, S. Arnouts, E. Le Floc'H, T. Moutard, K. Kraljic, et al. HSC-CLAUDS survey: The star formation rate functions since z ~ 2 and comparison with hydrodynamical simulations. Astronomy and Astrophysics - A&A, 2023, 675, 10.1051/0004-6361/202245756. insu-04169720

HAL Id: insu-04169720

https://insu.hal.science/insu-04169720

Submitted on 13 Aug 2023

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
# HSC-CLAUDS survey: The star formation rate functions since $z \sim 2$ and comparison with hydrodynamical simulations

V. Picouet<fnref n="1" />,<fnref n="2" />, S. Arnouts<fnref n="1" />, E. Le Floc'h<fnref n="3" />, T. Moutard<fnref n="1" />, K. Kraljic<fnref n="1" />,<fnref n="4" />, O. Ilbert<fnref n="1" />, M. Sawicki<fnref n="5" />, G. Desprez<fnref n="5" />, C. Laigle<fnref n="6" />, D. Schiminovich<fnref n="2" />, S. de la Torre<fnref n="1" />, S. Gwyn<fnref n="7" />, H. J. McCracken<fnref n="6" />, Y. Dubois<fnref n="6" />, R. Davé<fnref n="8" />,<fnref n="9" />,<fnref n="10" />, S. Toft<fnref n="11" />,<fnref n="12" />, J. R. Weaver<fnref n="13" />, M. Shuntov<fnref n="11" />, and O. B. Kauffmann<fnref n="1" />

Received 21 December 2022 / Accepted 5 April 2023

## ABSTRACT

**Context.** Star formation rate functions (SFRFs) give an instantaneous view of the distribution of star formation rates (SFRs) in galaxies at different epochs. They are a complementary and more stringent test for models than the galaxy stellar mass function, which gives an integrated view of the past star formation activity. However, the exploration of SFRFs has been limited thus far due to difficulties in assessing the SFR from observed quantities and probing the SFRF over a wide range of SFRs.

**Aims.** We overcome these limitations thanks to an original method that predicts the infrared luminosity from the rest-frame UV/optical color of a galaxy and then its SFR over a wide range of stellar masses and redshifts. We applied this technique to the deep imaging survey HSC-CLAUDS combined with near-infrared and UV photometry. We provide the first SFR functions with reliable measurements in the high- and low-SFR regimes up to $z = 2$ and compare our results with previous observations and four state-of-the-art hydrodynamical simulations.

**Methods.** The SFR estimates are based on the calibration of the infrared excess (IRX = $L_{IR}/L_{UV}$) in the NUVrK color-color diagram. We improved upon the original calibration in the COSMOS field by incorporating Herschel photometry, which allowed us to extend the analysis to higher redshifts and to galaxies with lower stellar masses using stacking techniques. Our *NrK* method leads to an accuracy of individual SFR estimates of $\sigma \sim 0.25$ dex. We show that it reproduces the evolution of the main sequence up to $z = 2$ and the behavior of the attenuation (or $\langle \text{IRX} \rangle$) with stellar mass. In addition to the known lack of evolution of this relation up to $z = 2$ for galaxies with $M_* \le 10^{10.3} M_\odot$, we observe a plateau in $\langle \text{IRX} \rangle$ at higher stellar masses that depends on redshift.

**Results.** We measure the SFR functions and cosmic SFR density up to $z = 2$ for a mass-selected star-forming galaxy sample (with a mass limit of $M_* \ge 2.10^9 M_\odot$ at $z = 2$). The SFR functions cover a wide range of SFRs ($0.01 \le \text{SFR} \le 1000 M_\odot \text{ yr}^{-1}$), providing good constraints on their shapes. They are well fitted by a Schechter function after accounting for the Eddington bias. The high-SFR tails match the far-infrared observations well, and show a strong redshift evolution of the Schechter parameter, SFR*, as $\log_{10}(\text{SFR}^*) = 0.58z + 0.76$. The slope of the SFR functions, $\alpha$, shows almost no evolution up to $z = 1.5-2$ with $\alpha = -1.3 \pm 0.1$. We compare the SFR functions with predictions from four state-of-the-art hydrodynamical simulations. Significant differences are observed between them, and none of the simulations are able to reproduce the observed SFRFs over the whole redshift and SFR range. We find that only one simulation is able to predict the fraction of highly star-forming galaxies at high $z$, $1 \le z \le 2$. This highlights the benefits of using SFRFs as a constraint that can be reproduced by simulations; however, despite efforts to incorporate more physically motivated prescriptions for star-formation and feedback processes, its use remains challenging.

**Key words.** galaxies: evolution – galaxies: star formation – galaxies: statistics – surveys – ultraviolet: galaxies – infrared: galaxies

## 1. Introduction

Spectroscopic and multiwavelength imaging surveys have provided insights into galaxy properties and their evolution across

cosmic time. A major result was the determination of the history of the cosmic star formation rate density (SFRD) thanks to analyses of galaxy luminosity functions in different wavelengths (from the UV/optical to the far-infrared and radio). After

<footnote n="1">¹ Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France
e-mail: vincent@picouet.fr</footnote>

<footnote n="2">² Department of Astronomy, Columbia University, 550 W. 120th Street, New York, NY 10027, USA</footnote>

<footnote n="3">³ Laboratoire AIM, CEA/DSM-CNRS-Université Paris Diderot, IRFU/Service d'Astrophysique, Bât. 709, CEA-Saclay, 91191 Gif-sur-Yvette Cedex, France</footnote>

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

<footnote n="5">⁵ Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary's University, 923 Robie Street, Halifax, Nova Scotia B3H 3C3, Canada</footnote>

<footnote n="6">⁶ Institut d'Astrophysique de Paris, UMR 7095, CNRS, UPMC Univ. Paris VI, 98 bis boulevard Arago, 75014 Paris, France</footnote>

<footnote n="7">⁷ Herzberg Astronomy and Astrophysics, National Research Council of Canada, 5071 West Saanich Rd., Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="8">⁸ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="9">⁹ University of the Western Cape, Bellville, Cape Town 7535, South Africa</footnote>

<footnote n="10">¹⁰ South African Astronomical Observatories, Observatory, Cape Town 7925, South Africa</footnote>

<footnote n="11">¹¹ Cosmic Dawn Center (DAWN), Copenhagen, Denmark</footnote>

<footnote n="12">¹² Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark</footnote>

<footnote n="13">¹³ Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA</footnote>
an increase from early time up to $z \sim 3-4$ (Smit et al. 2012; Mashian et al. 2016), the SFRD reaches a maximum at cosmic noon, $z \sim 1.5-2.5$ (Gruppioni et al. 2013), followed by a decline of an order of magnitude until today (Schiminovich et al. 2005; Karim et al. 2011; Madau & Dickinson 2014) despite a considerable amount of neutral and atomic gas available (Péroux & Howk 2020). This is corroborated by analyses of the integrated quantity and the galaxy stellar mass function (GSMF), which shows that half of the stellar mass density has already been assembled since $z \sim 2$ (Arnouts et al. 2007), with remarkably little evolution of the GSMF of the star-forming population since then (Moutard et al. 2016a). Star-forming galaxies (SFGs) gradually transition toward quiescent systems and have contributed to the buildup and evolution of the passive GSMF up to now (Ilbert et al. 2013; Davidzon et al. 2017).

Characterizing the main mechanisms involved in the evolution of the SFG population and understanding what triggers their star formation and what ultimately causes their migrations into passives, is thus a major challenge. The SFGs appear to lie on a tight sequence that links their stellar mass to their star formation rate (SFR), the so-called main sequence (MS; Noeske et al. 2007; Salim et al. 2007; Speagle et al. 2014; Whitaker et al. 2012). The tightness of the MS suggests that the SFGs grow through a secular evolution in an equilibrium between gas accretion, star formation, and outflows (Bouché et al. 2010; Davé et al. 2012; Lilly et al. 2013). The slope of the MS also suggests that lower-mass SFGs tend to be more efficient at forming stars with a higher specific star formation rate (sSFR; the SFR per stellar mass) than more massive ones. The depletion timescale based on the molecular gas content reveals that galaxies rapidly consume their gas, $t_{dep} = M_{mol}/\text{SFR} \sim 1-2$ Gyr (Bigiel et al. 2008; Tacconi et al. 2020), which must then be replenished for galaxies to stay on the MS. In the Λ cold dark matter (CDM) framework, it has long been claimed that galaxies can be fueled by cold gas thanks to cold mode accretion from cosmic web filaments, without the gas being gravitationally shock-heated (Kereš et al. 2005; Dekel et al. 2009), which also contributes to the acquisition of their angular momentum (e.g., Pichon et al. 2011). This cold accretion mode from the intergalactic medium is expected to be ubiquitous in the early Universe and in low-mass galaxies, while the hot accretion mode may dominate at lower redshifts for high-mass systems (Van de Voort et al. 2011; Snedden et al. 2016). Furthermore, CO observations show that SFGs appear to contain three to ten times more gas at redshift $z = 1-2$ than their local counterparts (Tacconi et al. 2010; Daddi et al. 2010). These gas-rich systems show more disturbed thick disks, with more dispersion-dominated kinematics. This suggests that the disks are less settled than their local counterparts (Kassin 2010; Kassin et al. 2012) and prone to violent disk instabilities that trigger enhanced SFRs (Cacciato et al. 2012), with SFGs moving up and down the MS following their gas accretion episodes and disk perturbations. Smooth and continuous gas flows from the cosmic web appear to be the key ingredient for conveying large quantities of cold gas at high redshifts and triggering intense star formation episodes.

Outflows, on the other hand, are the other components that alter the evolution of SFGs. Galactic winds from massive stars and supernova (SN) explosions can expel a fraction of the gas outside the disk, which reduces the star formation efficiency and contributes to the metal enrichment of the interstellar medium (ISM) and the circumgalactic medium (Davé et al. 2011a; Hopkins et al. 2014; Fontanot et al. 2017). While these processes may be efficient for low-potential-well

systems, they may not be sufficient for massive ones. At high masses, active galactic nucleus (AGN) feedback appears more effective at producing high-velocity winds, ejecting a large fraction of the gas, and preventing its cooling on a short timescale, thereby halting the star formation activity of the host galaxy (Hopkins & Beacom 2006; Cattaneo et al. 2009).

All these inflow and outflow processes that govern the evolution of galaxies are incorporated into the most recent semi-analytical models (SAMs) via empirical recipes and into numerical simulations as sub-grid physics (e.g., Somerville & Davé 2015; Davé et al. 2011b; Dubois et al. 2016). Their comparison with observations is crucial to constraining the influence of feedback processes. Such comparisons show that stellar winds on the low-mass end and AGN feedback on the massive end can explain the shape of the observed GSMFs and their deviation from the theoretical dark matter (DM) halo mass functions (e.g., Silk & Mamon 2012), and they are necessary for reproducing the stellar-to-halo ratio (Shuntov et al. 2022).

The star formation rate function (SFRF) is another independent constraint. In contrast to the GSMF, which provides an integrated view of the past star formation activity, the SFRF gives an instantaneous view of the distribution of the in situ SFR and its relative evolution with cosmic time. It provides insights into the importance of stellar winds, as some observations suggest a correlation between outflow velocities and SFRs (Heckman et al. 2015). By implementing four different stellar wind recipes in hydrodynamical simulations, Davé et al. (2011b) show that models can reproduce the faint end of the SFRF, but they all fail to suppress high SFRs at low redshifts ($z \le 2$). Katsianis et al. (2017a), using EAGLE simulations with different AGN and SN stellar wind implementations, show that SN winds play an essential role in reproducing the SFRFs at high redshifts and that AGN feedback becomes prominent at low redshifts. However, some discrepancies arise depending on the assumed SFRF measurements, based on UV, $H_\alpha$, or far-infrared (FIR) estimators. The simulations are in better agreement with the high-end SFR functions from UV/$H_\alpha$ and underestimate the number of high-SFR systems observed with the FIR SFRFs.

Accurately measuring SFR functions is a difficult task. It relies on observational SFR estimates, which are timescale dependent, subject to different dust attenuation effects, and sensitive to different selection effects. This can affect the shape of the SFRFs. All the known tracers (from the far-UV to the radio) have different benefits and drawbacks.

As the emission of galaxies in the rest-frame UV is dominated by young, short-lived ($t \sim 10^8$ yr), massive stars (Kennicutt 1998), UV represents a direct tracer of the SFR (Bouwens et al. 2009; Schiminovich et al. 2005). It is easily accessible over the entire history of the Universe, but UV light is efficiently absorbed and scattered by dust grains, which heat up and re-emit the absorbed energy at FIR wavelengths. A correlation between the infrared excess (IRX; $IRX = L_{IR}/L_{UV}$), a measurement of the UV attenuation ($A_{UV}$), and the slope of the UV continuum ($\beta$ slope) has been observed for starburst galaxies (Meurer et al. 1999; Calzetti et al. 2000). This dust correction is abundantly used to derive the SFR of high redshift galaxies, as the UV slope is the only accessible quantity (Smit et al. 2012; Katsianis et al. 2017b). However, a large scatter in the $IRX-\beta$ relation is observed for the SFG population, spanning a range between a Calzetti- and a Small Magellanic Could-like attenuation law (Seibert et al. 2005; Salim et al. 2007), which depends on galaxies' ages and metallicities (Boquien et al. 2009; Shivaei et al. 2020).

The total infrared luminosity ($L_{IR}$) is produced by the dust continuum emission and is a direct probe of the SFR
(Kennicutt 1998). It is defined as the integrated luminosity between 8 and 1000 µm and can be assessed either by combining multiband photometry or via monochromatic wavelengths, where tight correlations are observed between monochromatic and total luminosities (Bavouzet et al. 2008; Goto et al. 2011). However, short wavelength monochromatic luminosities (≤30 µm) can be impacted by the presence of AGNs and the heating of the dust by old stellar populations for evolved galaxies (Cortese et al. 2008). A constraint from the Rayleigh-Jeans part of the FIR spectral energy distribution (SED) is required to minimize their impacts on SFR estimates. Finally, FIR observations suffer from limited instrumental sensitivity and angular resolution, restricting detections to luminous distant infrared galaxies. While sensitivity can be partly compensated for by different stacking techniques (Heinis et al. 2013), the resolution leads to important confusion issues (Bethermin et al. 2012).

The $H_\alpha$ luminosity ($L_{H_\alpha}$) is produced from the gas ionized by short-lived massive stars ($t \sim 10^7$ yr) and is thus an excellent tracer of the instantaneous SFR (Kennicutt 1998). It is difficult, however, to observe at high redshifts as the line is redshifted into the near-infrared (NIR) domain. Alternatively, narrowband imaging surveys can be efficiently used to detect $H_\alpha$ line emitters from their color excess (Ly et al. 2011; Sobral et al. 2013). However, the derived luminosity is subject to several uncertainties: the contribution of the adjacent [NII] line, which is sensitive to the galaxy metallicities; and dust attenuation effects, which are based on an empirical relation with a large scatter (Ly et al. 2011; Sobral et al. 2013). Furthermore, the density of sources is sensitive to mismatched line contamination.

The SFR functions derived from the luminosity functions of the above tracers have recently been compiled over a wide redshift range (Katsianis et al. 2017a,b), each with specific caveats for converting their luminosities into dust-free SFRs. Below $z \sim 2$, all the SFRFs derived from UV luminosities show a shortage of high-SFR sources (SFR > $100 M_\odot \text{yr}^{-1}$), while FIR selection reveals sources up to SFR = $1000 M_\odot \text{yr}^{-1}$. In the low-SFR regime, the slope of the SFRFs from $H_\alpha$ and UV luminosities shows a large range of values, $-1.4 \le \alpha \le -1.8$, while the FIR observations do not have access to this regime.

In this work we make use of a new approach, combining UV and FIR observations to derive the dust correction to be applied to the UV luminosities of a mass-selected sample of SFGs. We then derive the SFRFs up to $z \sim 2$, based on the contribution of galaxies with $M_\star \ge 10^9 M_\odot$.

Arnouts et al. (2013, hereafter A13) found that the IRX shows a remarkable behavior in the rest-frame color-color diagram (NUV $-$ r) versus (r $-$ K). They identified a single vector by combining the two colors, $NrK$, which captures the behavior of the IRX over a large dynamical range with a small dispersion, $\sigma$(IRX) $\sim 0.2$ dex, and no mass dependence. By using a two-component dust model (i.e., birth clouds and the diffuse ISM; Charlot & Fall 2000) and a full distribution of galaxy inclinations (Tuffs et al. 2004; Chevallard et al. 2013), this model can reproduce the IRX distribution in the color-color diagram, confirming that it encodes information about energy transfer between starlight and dust. This method, by combining the UV and infrared luminosities, provides a direct measurement of the energy budget that can be used to predict the infrared luminosity with a simple optical estimator and assess the total SFR (Bell et al. 2005), without any assumption on the shape of the attenuation law.

While in A13 the IRX measurement was based on the $L_{IR}$ derived from the Spitzer 24 µm observations, in this work we make use of the COSMOS2020 catalog (Weaver et al. 2021),

which gathers Spitzer/MIPS (24 µm) and Herschel/PACS and SPIRE (100, 160, 250, 350, and 500 µm) bands (Jin et al. 2018), allowing us to extend the calibration up to $z \sim 2$. In addition, we use the stacking technique to extend the calibration into the low stellar mass regime. We then apply the new COSMOS2020 IRX versus $NrK$ calibration to the HSC-CLAUDS deep survey with very deep U ($u \sim 27$) and optical imaging ($i \sim 27$; Sawicki et al. 2019), where robust photometric redshifts have been estimated (Desprez et al. 2023). Thanks to its depth, the CLAUDS-HSC survey is ideally suited to measure the unobscured UV luminosity functions up to $z = 2$, providing the best constraints on both ends of the UV luminosity function (Moutard et al. 2020). In this work we restrict the analysis to the regions covered by deep NIR imaging ($\sim 5.5$ deg$^2$ in the COSMOS-E and XMM-LSS fields).

The paper is organized as follows. In Sect. 2 we describe both the COSMOS2020 data set used to derive the IRX calibration and the HSC-CLAUDS survey, to which the calibration is applied to measure the SFRFs. In Sect. 3 we describe the estimates of the physical parameters, and in Sect. 4 we perform the calibration of the IRX from direct FIR observations and extend it to low-mass systems using the FIR stacking technique. Section 5 presents the SFR functions derived between $0 \le z \le 2$ and a comparison with previous results from the literature, as well as a comparison with four hydrodynamical simulations, TNG100 from the ILLUSTRISTNG project (Pillepich et al. 2018; Nelson et al. 2019), EAGLE (Crain et al. 2015; Schaye et al. 2015; McAlpine et al. 2016), HORIZON-AGN (Dubois et al. 2014), and SIMBA (Davé et al. 2019). Finally, Sect. 6 presents the cosmic SFRD for different stellar mass and SFR regimes, and we conclude in Sect. 7. Appendix A describes the stacking analysis with the Spitzer-24 µm and Herschel-SPIRE data, and Appendix B gives a more detailed description of the four hydrodynamical simulations used in this work. Throughout this paper, we use a Chabrier (2003) initial mass function, all magnitudes are in the AB system (Oke 1974), and we adopt a flat $\Lambda$CDM cosmology with $\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$ and the Hubble constant $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2. Data

In this work we used the latest version of the Cosmic Evolution Survey (COSMOS) catalog provided by Weaver et al. (2021, COSMOS2020), combined with the super-deblended FIR photometry of the Spitzer and Herschel data from Jin et al. (2018). We performed the calibration of the IRX as a function of $NrK$ vector and stellar mass parameters to assess the SFR of individual galaxies, first based on direct measurements with detected FIR sources then using stacking techniques to extend the calibration to lower stellar masses. This calibration was then applied to the sources detected in the HSC-CLAUDS-NIR catalog (Desprez et al. 2023, hereafter D23) to derive the SFR functions up to $z \sim 2$.

### 2.1. COSMOS and far-infrared catalogs

#### 2.1.1. COSMOS2020 Catalog

COSMOS is a major extragalactic field with a large multi-wavelength photometric coverage over a 2 deg$^2$ field. It collects ground-based optical observations with intermediate and broadband filters, NIR photometry for a total of 31 filter bands as described in Laigle et al. (2016, COSMOS2015). The main improvement of the latest COSMOS2020 catalog with respect to COSMOS2015 is the depth of broadband imaging.
It combines the deeper CFHT *U*-band imaging from CLAUDS survey (Sawicki et al. 2019), the ultra-deep data from Public Data Release 2 (PDR2) of the HSC Subaru Strategic Program (HSC-SSP; Aihara et al. 2019) and the latest release of UltraVISTA (DR4; McCracken et al. 2012) as well as all the mid-infrared (MIR) imaging available with the *Spitzer*/IRAC channels.

COSMOS2020 provides a set of four photometric redshift catalogs. They rely on two different photometric extractions, based on SEtractor software (Bertin & Arnouts 1996) and the Farmer (a package running the Tractor code on multiwavelength images, Weaver et al. 2022, Sect. 3.2) and the two photometric redshift codes LePhare (Arnouts et al. 2002; Ilbert et al. 2006) and EAZY (Brammer et al. 2008). For more consistency with our catalog, in the following, we only refer to the SEtractor-LePhare catalog for both the calibration of the *NrK* versus IRX relations and the comparisons of the photometric redshifts and physical parameters derived with our HSC-CLAUDS catalog.

### 2.1.2. Far-infrared catalog

We use the “super-deblended” FIR to millimeter photometric catalog from Jin et al. (2018). They use the *Spitzer*/MIPS 24 µm images from the COSMOS-*Spitzer* survey (Le Floc’h et al. 2009), the *Herschel*/PACS (100 and 160 µm) images from PEP survey (Lutz et al. 2011) and *Herschel*/SPIRE (250, 350, and 500 µm) from the HerMES survey (Oliver et al. 2012). Point spread function (PSF) prior-fitting multiband photometry was performed by adopting the prior positions of the 24 µm, radio, and *K*<sub>s</sub>-band mass-selected sources (Liu et al. 2018). The typical fluxes at *S*/*N* = 5 correspond to 50 µJy in *Spitzer*/24 µm, 8.3, 22.9 mJy in PACS-100 and 160 µm, 7.6, 11.0, 13.2 Jy in SPIRE-250, 350, and 500 µm.

We restricted the FIR population to the 24 µm sources with a signal-to-noise *S*/*N* ≥ 5 and rejected a few anomalous sources with *m*<sub>AB</sub>(3.6 µm) > 22.5, leading to a catalog of ~26 000 galaxies. Our sample is driven by the 24 µm, and only the most luminous sources are detected with *Herschel*. When adopting a *S*/*N* = 3 for *Herschel*, ~12%, 10%, 34%, 8%, and 4% are detected respectively in PACS-100 µm, PACS-160 µm, and SPIRE-250, 350, and 500 µm.

The FIR sources are matched with the COSMOS2020 catalog, adopting a small positional uncertainty of 0.25 arcsec (as all the high S/N sources of interest are attached to a *K*-band counterpart in the COSMOS2015 catalog). Limiting the redshift range in between 0 ≤ *z* ≤ 2, the final FIR catalog contains ~23 200 sources.

## 2.2. The HSC-CLAUDS survey

### 2.2.1. Observations

HSC-CLAUDS combines the PDR2 of the deep and ultra-deep layers of the HSC-SSP (Aihara et al. 2019) with the deep *U*-band observations carried out with the MegaCam instrument at CFHT (CLAUDS; Sawicki et al. 2019). This data set is a unique combination of depth and area, reaching 26–27th magnitude over a total area of ~20 deg², split into four separate regions (E-COSMOS, XMM-LSS, ELAIS-N1, and DEEP2-3). In the E-COSMOS and XMM-LSS fields, HSC-CLAUDS is combined with the publicly available VIRCAM NIR observations from the VIDEO and UltraVISTA surveys (Jarvis et al. 2013; McCracken et al. 2018, respectively). The different bands are shown in Fig. 1.

Two catalogs have been produced for the source extraction and flux measurements, which are described in the companion paper by D23. They discuss the processing steps to combine all the above data sets into the HSC grid and the addition of the external bands into the dedicated photometric HSC pipeline (Bosch et al. 2018). In addition to the catalog produced with the HSC pipeline, a second catalog based on SEtractor software (Bertin & Arnouts 1996) is produced using a multiband $\chi^2$ image in dual mode. In the SEtractor catalog we also include the far-UV (FUV; $\lambda$ ~ 1500 Å) and near-UV (NUV; $\lambda$ ~ 2300 Å) imaging from the GALEX satellite (Martin 2005). The deblending of the UV photometry was performed with the EMPHOT code (Conseil et al. 2011) by using optical *u*-band detections as prior down to *u* ~ 25 (e.g., Zamojski et al. 2007; Moutard et al. 2016b).

D23 gives all pieces of information regarding magnitude, color, depth, and photometric redshift measurements, and detailed comparisons between the two catalogs and external data sets (COSMOS2020, CANDELS, and spectroscopic samples) are discussed. The two catalogs are publicly available<fnref n="1" />. We restricted our analysis to the E-COSMOS ($\alpha$ ~ 150°, $\delta$ ~ 2.5°) and XMM-LSS ($\alpha$ ~ 35.5°, $\delta$ ~ -5°) fields, where the deep NIR photometry is available. Figure 2 shows the layouts of the CLAUDS (with *u* and *u*<sup>*</sup> filters), HSC-SSP, and VIRCAM observations in these two regions. The 5$\sigma$ depths in 2 arcsec apertures for the deep and ultra-deep components as well as their respective area are given in Table 1. According to those expected depths, we restrict our analysis to galaxies with $i_{AB} \le 26.5$ and $K_{AB} \le 24.5$, which also corresponds to the limits where the number counts in the two fields are consistent (D23). The large area covered benefits the present analysis by reducing the impact of cosmic variance in the estimate of the SFR functions and cosmic SFRD.

### 2.2.2. Photometric redshifts

The photometric redshifts are derived with LePhare code using all available passbands. They are described in D23 and compared to the extensive spectroscopic redshift samples available in the two fields. A good agreement is observed at magnitudes brighter than $i_{AB} \le 25$, with a scatter $\sigma \le 0.03$ ($\sigma$ being defined as $\sigma = 1.48 \times \text{Median}(|z_p - z_s|/(1 + z_s))$ and a small fraction of outliers $\eta \le 7\%$ ($\eta$ being defined as the fraction of galaxies with $\Delta z > 0.15 \times (1 + z)$).

At fainter magnitude, $i_{AB} \ge 25.5$, the scatter remains small but the outlier fraction doubles with the majority being catastrophic redshifts between low and high redshifts.

## 3. Physical parameter estimates

### 3.1. Stellar masses and luminosities

To derive the physical parameters and rest-frame luminosities we use LePhare code (Arnouts et al. 2002; Ilbert et al. 2006) with the stellar population synthesis model from Bruzual & Charlot (2003, hereafter BC03), following a similar procedure as Ilbert et al. (2015). We fix the redshift to its spectroscopic value if available or to the photometric redshift value based on the median of the marginalized probability distribution function. Our BC03 library includes six exponentially declining star formation histories, following $\tau^{-1}e^{-l/\tau}$ with $\tau$ varying between 0.1 and 30 Gyr, and two delayed star formation histories; two stellar

<footnote n="1">¹ https://www.clauds.net/</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Transmission curves of the photometric bands used to derive galaxy properties (left: UV to NIR) and infrared luminosities (right: Spitzer/MIPS + Herschel/PACS-SPIRE) for the calibration of the IRX. Transmission curves are arbitrarily scaled to one. Two spectra of SFGs, without dust (blue) and with dust (red), are shown at $z = 0, 0.5, 1, 1.5, 2$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Deep (solid lines) and ultra-deep (dashed lines) footprints of the HSC-CLAUDS and VIRCAM observations, as indicated in the inset and overlaid on the background detection images. Starting in 2015, the *u*\* filter (blue lines) was replaced by a new filter (*u*̄, with a slightly bluer effective wavelength; light blue line).</figcaption></figure>

metallicities ($Z/Z_⊙$ = 0.4, 1); two extinction laws with a maximum dust reddening of $E(B-V) = 0.7$. We also imposed the prior $E(B-V) < 0.15$ if age/$\tau > 4$ (a low extinction is imposed for galaxies that have a low SFR) and include the contribution of emission lines using an empirical relation between the UV light and the emission line fluxes (Ilbert et al. 2009).

The physical parameters are derived by computing the median of the marginalized likelihood for each parameter and the errors corresponding to the 68% confidence level. To derive the rest-frame luminosities (or absolute magnitudes), we adopted the same approach as Ilbert et al. (2005), using the photometry in the nearest rest-frame broadband filter to minimize the dependence on the k-correction. We note that in this work, the SFR measurement does not rely on the SED fitting ingredients, such as the dust attenuation law or the reddening excess. They are only used to best fit the observed multiband photometry and

derive the luminosities in different passbands. The SFR is then estimated from a combination of the NUV, *r*, $K_s$ luminosities, and redshift (see Sect. 4.1).

As we are only interested in SFGs in this analysis, we use the rest-frame color-color diagram ($NUV_{abs} - r_{abs}$) versus ($r_{abs} - K_{abs}$) (Fig. 3) to classify quiescent and SFGs (Arnouts et al. 2013), following the redshift dependence proposed by Moutard et al. (2016a; see their Fig. 8 and Sect. 5.1). All the galaxies not classified as quiescent are classified as star-forming.

In Fig. 4 we compare the stellar mass and the *NrK* vector (defined in the next section) derived with our HSC-CLAUDS photometry and with the COSMOS2020 catalog, as a function of redshift and stellar mass. For this comparison, we restrict the sample to the SFGs with photo-z within $|z_{COSMOS} - z_{CLAUDS}|/ (1 + z̄) \le 0.1$. This cut will always be used when comparing the physical parameters. Our stellar masses (top panels) appear
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Behavior of the IRX in the NUVrK color-color diagram. Left: Mean IRX (⟨IRX⟩ color-coded in logarithmic scale) in four redshift bins. The thin parallel lines show the modeled evolution of the ⟨IRX⟩ stripes with the norm of the NrK vector perpendicular to them. In each panel, the dotted line indicates the region of passive galaxies (Moutard et al. 2016a) that are not included in this analysis. Right: Same as the left panel but for the dispersion around the mean (σ(IRX)), color-coded in a logarithmic scale.</figcaption></figure>

<table><caption>**Table 1.** Summary of UV-optical-infrared observations with the depth (corresponding to the magnitude at 5σ in 2" aperture) and area for the deep (DD) and ultra-deep (UDD) regions.</caption><thead><tr><th>Telescope<br>(Survey)</th><th>Filter</th><th>Central<br>λ [Å]</th><th>Depth<br>DD(UDD)</th><th>Surface<br>deg<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="2">GALEX<br>(GALEX-DIS)</td><td>FUV</td><td>1526</td><td>25.5</td><td>5.3</td></tr><tr><td>NUV</td><td>2307</td><td>25.5</td><td>4.6</td></tr><tr><td rowspan="2">CFHT<br>(CLAUDS)</td><td><i>u</i></td><td>3709</td><td>27.2</td><td>1.4</td></tr><tr><td><i>u</i>*</td><td>3858</td><td>26.0 (27.5)</td><td>5.6 (1.7)</td></tr><tr><td rowspan="5">SUBARU<br>(HSC-SSP)</td><td><i>g</i></td><td>4847</td><td>27.0 (27.7)</td><td>5.6 (2.8)</td></tr><tr><td><i>r</i></td><td>6219</td><td>26.6 (27.5)</td><td>5.6 (2.8)</td></tr><tr><td><i>i</i></td><td>7699</td><td>26.5 (27.2)</td><td>5.6 (2.8)</td></tr><tr><td><i>z</i></td><td>8894</td><td>26.0 (26.5)</td><td>5.6 (2.8)</td></tr><tr><td><i>y</i></td><td>9761</td><td>25.0 (26.0)</td><td>5.6 (2.8)</td></tr><tr><td>VISTA</td><td><i>Y</i></td><td>10216</td><td>26.0/25.5</td><td>4.2/1.4</td></tr><tr><td rowspan="3">VIDEO/U-VISTA</td><td><i>J</i></td><td>12525</td><td>25.7/25.3</td><td>2.9/1.4</td></tr><tr><td><i>H</i></td><td>16466</td><td>25.2/25.0</td><td>4.2/1.4</td></tr><tr><td><i>K</i><sub>s</sub></td><td>21557</td><td>24.8/24.8</td><td>4.2/1.4</td></tr></tbody></table>

**Notes.** The two regions are entirely covered by NIR data, except for the *J* band in one pointing.

in good agreement with COSMOS2020 estimates with a relative difference lower than 25% (which means less than ~0.1 dex) and no significant trend is observed with redshift and stellar mass. Since the *K*-band luminosity relies on some extrapolation of the SEDs at high redshift, in the bottom panels we compare our *NrK* vectors with the COSMOS2020 estimates, which are better-constrained thanks to the use of the IRAC MIR photometry. The *NrK* relative differences show no bias with redshift and stellar mass and typical variation of less than ~25% except in the highest redshift bins and lowest stellar mass bins where it increases up to 40%.

The stellar mass completeness (smallest mass at which most of the objects would still be observable) of our SFG sample

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Relative differences of the stellar masses and the *NrK* vectors (see Sect. 4) between HSC-CLAUDS and COSMOS2020 as a function of redshift and stellar mass.</figcaption></figure>

($i_{AB} < 26.5$ and $K_{AB} < 24.5$) is empirically computed following the commonly used method (Ilbert et al. 2013; Weaver et al. 2021) developed by Pozzetti et al. (2010) based on the mass-to-light ratio $M/L$ dependence. At each redshift, we define a minimum mass $M_{min}$ above which the stellar mass function of each subpopulation is complete. To do so, we rescale each galaxy's stellar mass computed via template-fitting to the stellar mass limit ($M_{lim}$) it would have at its redshift if its $K_{AB}$ apparent magnitude was equal to the limiting magnitude of the survey ($K_{lim} = 24.5$): $log_{10}(M_{lim}) = log_{10}(M_{med}) + 0.4(K - K_{lim})$. At
<figure><img src="image_6.png" /><figcaption>Fig. 5. Difference between the infrared luminosity, $L_{IR}$, estimated from the 24 µm photometry alone and by adding the Herschel photometry up to $z \sim 2$.</figcaption></figure>

each redshift, to derive a representative limit, we select the 20% faintest galaxies of our sample. The stellar mass completeness can then be determined in a given redshift bin from the distribution of the rescaled masses of this faint subsample: the 95th percentile of the distribution defines a mass at which most of the objects would still be observable.

## 3.2. Far-infrared luminosities and far-infrared sample properties

The infrared luminosity ($L_{IR}$) is defined as the luminosity integrated from 8 to 1000 µm:

$$L_{IR} = \int_{\lambda=8 \mu\text{m}}^{1000 \mu\text{m}} L(\lambda)d\lambda. \qquad (1)$$

It is derived by using the code LePhare (Arnouts et al. 2002; Ilbert et al. 2006) combined with the FIR SED templates of Dale & Helou (2002). In A13, the $L_{IR}$ was derived by extrapolating the 24 µm observed flux density with the FIR templates in such a way that the templates follow the locally observed dust temperature-luminosity relationship (Chary & Elbaz 2001; Goto et al. 2010). By using the Herschel observations, Elbaz et al. (2011) confirmed the tight correlation between the infrared luminosity extrapolated from 24 µm photometry ($L_{IR}^{24\mu\text{m}}$) and the infrared luminosity ($L_{IR}$) measured with the PACS and SPIRE photometry. However, they noticed a bias at high luminosity ($L_{IR}^{24\mu\text{m}} > 10^{12} L_{\odot}$), where the $L_{IR}^{24\mu\text{m}}$ systematically overestimates the $L_{IR}$. In Fig. 5 we compare the infrared luminosity measured with the six passbands (24 µm to 500 µm) with the one derived from the 24 µm alone. All the galaxies in this sample have at least one Herschel flux with $S/N > 3$ and the scaling and SED shape are left free in the fitting procedure for the $L_{IR}$ estimate. We confirm the excellent agreement for the SFGs with low infrared luminosity, with a small scatter ($\sigma \sim 0.07$ dex) in between the two estimates, while at high luminosity we confirm the overestimation of $L_{IR}^{24\mu\text{m}}$ compared to $L_{IR}$. We therefore included the Herschel photometry in the analysis to improve the estimate at high infrared luminosity; it does not impact the $L_{IR}$ estimate for the bulk of the 24 µm sample.

The properties of the FIR population are shown in Fig. 6. In the bottom panel, we show the SFR distribution as a function of photometric redshift. The SFR is defined as the sum of the UV and FIR contribution, as in A13:

$$\text{SFR}(M_{\odot} \text{yr}^{-1}) = 8.6 \times 10^{-11} \times (L_{\text{IR}} + 2.3 \times L_{\text{NUV}}), \qquad (2)$$

<figure><img src="image_7.png" /><figcaption>Fig. 6. Stellar mass distributions (top panel) and SFR distributions (bottom panel) as a function of redshift for the 24 µm sample. The solid blue line represents the 50% mass completeness limit of the FIR sample (see text), while the orange line represents the stellar mass completeness of the HSC-CLAUDS K-selected sample ($K \le 24.5$; see text). The horizontal gray lines in the bottom panel indicate the SFR threshold that corresponds to the luminous (LIRG, $L_{IR} \ge 10^{11} L_{\odot}$) and ultra-luminous (ULIRG, $L_{IR} \ge 10^{12} L_{\odot}$) infrared galaxies.</figcaption></figure>

where $L_{NUV}$ is the monochromatic NUV luminosity: $L_{NUV}/L_{\odot} = \nu L_{\nu}(2300 \AA)$. As already shown by Le Floc'h et al. (2005), below $z = 0.5$, the population is composed of moderately SFGs with SFR $\le 10 M_{\odot} \text{yr}^{-1}$.

The fraction of luminous infrared galaxies (LIRGs) gradually increases from $z = 0.5$ to $z = 1$ and dominates at $z \ge 1$. At all redshifts, the fraction of ultra-luminous infrared galaxies (ULIRGs) is negligible.

In the top panel, we show the stellar mass distribution as a function of photometric redshift. The FIR sample is dominated by galaxies with $M_{\star} \ge 5 \times 10^9 M_{\odot}$. To characterize how representative the 24 µm sample is with respect to the entire star-forming population at a given mass and redshift, we define a 50% completeness stellar-mass limit with the ratio $\int_{M_{\text{lim}}}^{\infty} \Phi_{\text{SF}}^{\text{FIR}}(M_{\star}, z)dM / \int_{M_{\text{lim}}}^{\infty} \Phi_{\text{SF}}^{\text{all}}(M_{\star}, z)dM_{\star} \sim 0.5$; where $\Phi_{\text{SF}}^{\text{FIR}}(M_{\star})$ and $\Phi_{\text{SF}}^{\text{all}}(M_{\star})$ are the $V_{\text{max}}$ weighted comoving volume densities of FIR ($f_{24\mu\text{m}} \ge 50 \mu\text{Jy}$) and $K_{\text{s}}$-selected ($K \le 24.5$) samples of SFGs, respectively. Above this limit (shown as a solid blue line), we consider the physical properties of the FIR population to be representative of the whole star-forming sample. This corresponds to a stellar mass of $M_{\star} \sim 1-2 \times 10^{10} M_{\odot}$ at $z \ge 0.5$. The lower completeness near $z \sim 1.4$ reflects the dip in between the polycyclic aromatic hydrocarbon features passing in the 24 µm passband (Fig. 1).

## 4. The infrared excess in the NUVrK diagram

For the majority of galaxies, the summation of the infrared and UV luminosities is a reliable indicator of the bolometric luminosity coming from young stars and therefore a good proxy for SFR (Buat et al. 2002). As the UV luminosity is of easy reach up to high redshift, measuring the IRX, defined as $IRX = L_{IR}/L_{NUV}$, allows for assessing the obscured star-formation contribution and offers an interesting alternative to SED fitting
<figure><img src="image_8.png" /><figcaption>Fig. 7. Evolution of the IRX with NrK (upper plots) and redshift (bottom plots) in different redshift and NrK bins. The back dots represent the IRX computed on the 24 µm-selected sample, and the dotted line corresponds to the minimum least square fitting at the center of each bin. The back lines around each dotted line represent the calibration for the edge of each bin.</figcaption></figure>

derived SFR, which depends strongly on the adopted attenuation laws<fnref n="2" />.

In contrast, the IRX is weakly dependent on the age of the stellar population, dust geometry, and nature of the extinction law (Witt & Gordon 2000). We note that our definition of the IRX differs from the literature, which usually adopts the FUV luminosity. Indeed, the NUV luminosity can be impacted by the dust attenuation bump at 2175 Å, which falls in the blue side of the NUV passband. But we decided to use the NUV luminosity for practical reasons as it proves to be more reliable thanks to the deep GALEX NUV and CFHT u-band observations in the whole redshift range considered in this work. As shown by Hao et al. (2011), this choice does not impact the reliability of the SFR estimates.

Based on a 24 µm-selected sample, A13 have measured a tight correlation between the NrK vector and the IRX values up to redshift $z = 1.3$, with almost no dependence on the stellar mass. In this section we revisit this analysis and extend it to higher redshift, $z \sim 2$, by using the FIR photometry (including MIPS/Spitzer and PACS and SPIRE/Herschel) and the latest COSMOS2020 catalog (Weaver et al. 2021). We also extend to lower-mass populations, based on an FIR stacking technique.

## 4.1. IRX calibration with detected FIR sources

Figure 3 shows the volume-weighted mean IRX (⟨IRX⟩) in the NUVrK diagram. In all the redshift bins, we observe an increase of the ⟨IRX⟩ by 1.5–3 dex from the bottom left corner to the upper right one with a small scatter around the mean. We also observe constant IRX stripes, which allows us to describe the variation in IRX by a single vector perpendicular to those stripes:

$$NRK = \sin(\phi) \times (\text{NUV} - r) + \cos(\phi) \times (r - K). \quad (3)$$

Since the dispersion $\sigma[\text{IRX}(\phi)]$ reaches a minimum when the vector $\text{NRK}(\phi)$ is perpendicular to the stripes, we minimize $\sigma[\text{IRX}(\phi)]$ to find the best angle. We obtain an angle of $\phi = 21^\circ$, consistent with A13 ($\phi = 18^\circ$), but we note that a difference up to 10 degrees does not impact the final calibration. The right side of the Fig. 3 shows that the scatter at all redshift stays below 0.2 dex over the whole color diagram.

With the definition of the NrK, we can now derive the relationship between ⟨IRX⟩ and NrK. In Fig. 7 we show the mean ⟨IRX⟩ and the associated scatter for the FIR-detected sources as a function of NrK and redshift. In all panels, a tight linear correlation is observed with a small scatter ($\sigma \le 0.3$ dex) compared to the dynamical range covered by ⟨IRX⟩ (~2 dex). The entire calibration does not depend on the stellar mass, and while the slope of the relation does not change with redshift, the normalization increases with redshift. At a fixed redshift, the populations in different stellar mass bins follow the same relation but occupy

<footnote n="2">² This is not the case for SED fitting codes preserving the energy budget, which gives an SFR consistent with the addition of infrared and UV contribution, and also corrects for the contribution of the infrared emission due to dust heating by old stars and not directly connected to star formation.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Comparison of the infrared luminosity, $L_{IR}$, estimated from the COSMOS2020 FIR sample with that estimated with the IRX-NrK calibration ($L_{IR}^{NrK}$) as a function of redshift and infrared luminosity.</figcaption></figure>

different regions as their NrK distributions shift to higher values at higher stellar mass.

We first adopted the same parametrization for the 〈IRX〉 versus NrK relation as A13, with two separated quantities 〈IRX〉 = $f(z) + \alpha \cdot NrK$, where $f(z)$ is a third-order polynomial function describing the redshift evolution, and $\alpha$ is a constant for the NRK dependence. We performed a linear least-square fit to derive the four free parameters and their standard deviation uncertainty. For the redshift evolution $f(z) = a_0 + a_1 \cdot z + a_2 \cdot z^2 + a_3 \cdot z^3$ we derived $a_0 = -0.73 \pm 0.02$; $a_1 = 2.02 \pm 0.06$; $a_2 = -1.33 \pm 0.06$; $a_3 = 0.33 \pm 0.02$; for the NrK term, $\alpha = 0.63 \pm 0.003$.

The result of this calibration is added as dashed black lines in Fig. 7, and the predicted $L_{IR}$ ($L_{IR}^{NrK}$) that comes from it are compared to the reference $L_{IR}$ in Fig. 8. The predicted and reference $L_{IR}$ are in good agreement with almost no bias with redshift and a global dispersion $\sigma \sim 0.24$ dex. While this scatter is slightly larger than in the analysis by A13, it covers a larger redshift range by extending the analysis up to $z = 2$. Finally, the comparison as a function of infrared luminosity shows almost no bias except in the extreme regimes (for ULIRGs and low infrared luminosity galaxies), where our calibration slightly underestimates or overestimates the luminosity by less than a factor of 2.

## 4.2. IRX calibration based on FIR stacking

Surveys conducted in the thermal infrared regime are generally not sensitive enough for detecting low-mass sources individually. To explore the validity of the IRX-NrK relationship in the low-to-intermediate stellar-mass regime we resort to the FIR stacking technique (Bethermin et al. 2010). However, while the stacking technique is often applied to stellar mass-selected samples (Whitaker et al. 2014), at a fixed stellar mass the 〈IRX〉 spans a large dynamical range, which can be better constrained by dividing the sample as a function of the unobscured UV luminosity (Heinis et al. 2014) where the most extinguished galaxies are expected to be the less UV luminous in a fixed stellar mass bin, or alternatively (and hereafter), as a function of NrK values

as it best follows the evolution of 〈IRX〉 in order to minimize the dispersion around the mean. We thus split the sample into NrK, stellar mass, and redshift bins and co-add the images from the Spitzer-24 µm and SPIRE-Herschel FIR channels for all the selected sources to get their average FIR emissions. Some stacks are shown in Fig. A.1 and the details about the procedure are given in Appendix A.

The results of the stacking procedure are shown in Fig. 9 with the mean infrared luminosity derived from the mean 24 µm flux (blue diamonds) or by combining the 24, 250, 350, and 500 µm fluxes (red circles). Despite all the potential, hard-to-control biases that such a technique can introduce, as discussed in the appendix, we observe an overall good agreement between the stacking results and the individual calibration (black circles). We obtain the same trend with the stacking technique, namely the tight correlation between 〈IRX〉 and NrK, with a similar slope and small scatter around the mean values. However, noticeable changes in the normalization in the different bins are observed. In the regime where the FIR-detected sources are the dominant population (panels below the thick black line), the stacking is in agreement within the respective uncertainties. For the panels above the black line, the individual detections lead to a systematically higher 〈IRX〉 at a fixed NrK. By nature, the FIR detections biased the samples toward the most obscured ones and did not reflect the behavior of the whole population in a given stellar mass and NrK bin.

We also observe that at high redshift and stellar mass bins (bottom-right panels), the $L_{IR}$ measured from the 24 µm mean flux alone tends to be slightly higher than the one estimated with 24+SPIRE mean fluxes. This is consistent with the expected overestimation of the $L_{IR}$ at high luminosity discussed in Sect. 3.2.

To take the stacking results into account in the following analyses, we considered two cases. In the first case, we disregarded the stacking and kept the original IRX-NrK calibration as derived with the FIR-detected sources. In the second case, we used a mixed calibration. For galaxies below the FIR stellar mass completeness, we applied a linear correction with stellar mass and redshift (shown as red lines in Fig. 9) and parameterized as

$$ \text{IRX}_{\text{stack}} = \text{IRX}_{\text{NrK}} - (0.34z - 0.10M_{\star} + 0.69), \quad \text{if } M_{\star} < M_{\text{lim}}(z). \qquad (4) $$

In each ($M_{\star}, z$) panel of Fig. 9, we report the corresponding shifts between the original and stacking calibrations. While the former will lead to an overestimate of the predicted SFR, the latter will lead to lower SFR values but maybe yet more appropriate at low mass or high redshift for the whole population.

## 4.3. Dust attenuation evolution with stellar mass and redshift

The IRX, or equivalently the dust attenuation, has been found to be a strong function of stellar mass at low (Garn & Best 2010) and high (Pannella et al. 2009; Reddy et al. 2010; Whitaker et al. 2012; Heinis et al. 2014; Shapley et al. 2022) redshift, with the most massive galaxies being more highly obscured. Since our NrK calibration does not rely explicitly on the stellar mass, in Fig. 10, we show the behavior of the mean IRX as a function of stellar mass in four different redshift bins. The mean IRX values are measured by using the FIR-detected sources (solid dots), the stacked sample (empty diamonds) both from the COSMOS2020 sources and by using the NrK prediction with or without the stacking correction (solid and dashed lines, respectively) for the HSC-CLAUDS sample. In each redshift panel, the stellar mass distributions are shown in the bottom part. The
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Mean IRX as a function of NrK for the detected FIR sources (black circles and gray shaded histograms) and for the stacked populations (open blue diamonds with 24 µm alone and red circles with 24 µm+SPIRE images) in different redshift and mass bins. The IRX-NrK calibration based on the detected sources is shown as solid black lines, while the dashed red lines show the calibration based on the stacking results, which are considered only in the panels where the FIR-detected sources are not representative of the whole population (limits are established with the blue line of Fig. 6, and corresponding panels are colored in pale red) and delineated by the thick solid black line. The shift between the two calibrations is reported in each panel.</figcaption></figure>

bottom right panel combines the 〈IRX〉 versus stellar mass relationship for the four redshift bins.

The mean IRX derived from the different approaches is consistent with each other. Despite the absence of calibration with stellar mass, the NrK-based estimates reproduce the same behavior as a function of stellar mass and redshift. It should be noted that in Fig. 10 we do not show results from the literature as some assumptions regarding the dust attenuation law (which vary across the NUVrK plane and the specific SFR of galaxies, Arnouts et al. 2013) is required to convert A(Hα) (Garn & Best 2010; Shapley et al. 2022) or LIR/LFUV (Heinis et al. 2014; Whitaker et al. 2014) into our LIR/LNUV measurements, making the correction uncertain.

Two main features are observed. First, at low stellar mass ($M_{\star} \le 10^{10.3} M_{\odot}$), almost no evolution of the IRX versus stellar mass relationship is observed between $z=0$ and $z=2$ (in agreement with Whitaker et al. 2014; Shapley et al. 2022). Considering the well-established evolution of the atomic/molecular gas, dust content, and metallicity at fixed stellar mass with redshift, the lack of evolution of the attenuation with redshift is rather

puzzling (Bogdanoska & Burgarella 2020) and could underlie some changes in the dust properties such as grain size distribution and composition (e.g., Shapley et al. 2022).

Second, at high stellar mass ($M_{\star} \ge 10^{10.3} M_{\odot}$), a 〈IRX〉 plateau is observed despite the large scatter in the measurements. It appears at different stellar mass with redshift, near $M_{\star} \sim 10^{10.7} M_{\odot}$ at $z \sim 2$ and $M_{\star} \sim 10^{10.3} M_{\odot}$ at $z = 0.25$. This result is in qualitative agreement with the trend reported by Whitaker et al. (2014, their Fig. 5) based on the stacking of Spitzer-24 µm data.

The origin of the plateau and its redshift dependence is also unclear. It may be related to a change in the balance between the formation and destruction of dust grains in massive galaxies. SNe and AGB stars are efficient sources to produce and inject dust grains in the ISM (Gehrz 1989) while accretion growth in dense molecular clouds (Draine et al. 2009) contributes to increasing the dust mass in galaxies. On the other hand, hot gas ISM and SN shock waves can destroy dust grains (or modify dust grain size) by shattering or sputtering processes (e.g., Inoue 2011).
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Evolution of the mean IRX with stellar mass in four redshift bins based on the measurements from the COSMOS2020 FIR-detected sources (empty circles), FIR stacked sources (filled diamonds), and the HSC-CLAUDS sample with or without the IRX stacking correction (solid and dashed lines, respectively). In the bottom right panel, we show the behaviors for the COSMOS2020 FIR stacking in the four redshift bins.</figcaption></figure>

Alternatively, high redshift galaxies are more gas-rich, with more turbulent disks than their low-$z$ counterparts (Kassin 2010; Förster Schreiber & Wuyts 2020), and star formation is potentially more concentrated in long-live giant clumps as suggested by hydrodynamical simulations (Fensch & Bournaud 2021). At a fixed stellar mass and lower redshift, galaxies have a lower SFR and are less clumpy. This contributes to reducing the dust production by SNe and growth by accretion while keeping efficient dust destruction in the hot gas ISM by sputtering and SN shock waves, as well as the presence of X-ray feedback from AGNs (Choi et al. 2012; Davé et al. 2019), which may lead to a reduced dust mass and thus produces a lower mean dust attenuation. While it is beyond the scope of this paper to propose a coherent model to explain the behavior of the IRX versus stellar mass relation, it will be of interest to know if the origin of the plateau is related to the survival/destruction of dust and/or to the disk properties with cosmic time in massive galaxies.

## 4.4. Comparison of individual SFR estimates with literature values

The two previous 〈IRX〉–NrK calibrations have been applied to the HSC-CLAUDS sample. The total SFR is derived by computing respectively $L_{\text{IR}}^{\text{NrK}}$ ($L_{\text{IR}}^{\text{NrK}} = \langle\text{IRX}\rangle \cdot L_{\text{NUV}}$) and $L_{\text{IR}}^{\text{NrK,stack}}$ and by applying Eq. (2).

In Fig. 11 we compare the SFRs for the HSC-CLAUDS samples based on the NrK methods with the FIR-based SFR

measured with the COSMO2020 data set at high masses ($M_{\star} \ge 10^{10} M_{\odot}$), while at low mass we compare with the SFRs measured in the CANDELS fields (Barro et al. 2019), for the two fields included in the HSC-CLAUDS regions, to extend our comparison to low-mass galaxies. For the latter, we use their SFR-ladder, estimated from a combination of three tracers. They use FIR and MIR photometry when available to estimate the $L_{\text{IR}}$ and the $\beta$-slope when infrared data are not available. For extreme galaxies detected by both MIPS and Herschel, total SFRs are calculated with dust emission models fitting the infrared data points and adding the unobscured star formation. For sources having only 24 um flux measurements (undetected at longer wavelength), the obscured SFR is computed following Wuyts et al. (2008). For the remaining galaxies undetected by MIR and FIR surveys, the SFRs are estimated using the IRX-$\beta_{\text{UV}}$ relations (85% of the ~10K matched sources). In all the mass bins the different SFR estimates are consistent with each other within their relative uncertainties. At low mass, $M_{\star} \le 10^{10} M_{\odot}$, our results are in overall good agreement with CANDELS SFRs, with a small scatter ($\sigma < 0.2$) and bias. Adopting the calibration based on the stacking results slightly reduces the bias, suggesting that this calibration may be more appropriate in the low-mass regime. At high masses, with the COSMOS2020 FIR sample, we still have an overall good agreement. However, we observe a wavy shape of the mean difference with redshift, which contributes to enlarging the scatter reported in each panel, up to 0.3 dex. We checked that this was not due to the estimates of
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Comparison of the *NrK*-based SFR of HSC-CLAUDS with that of CANDELS based on the SFR-ladder estimates at low masses (top panels) and with COSMOS2020 based on the FIR SFR estimates at high stellar masses (bottom panels). The *NrK*-based SFRs are derived from the calibration without (black circles) or with (red crosses) corrections based on the staking results.</figcaption></figure>

the *NrK* vector between the HSC-CLAUDS and COSMOS2020 data set. It is rather due to the COSMOS2020 calibration, where our parametric form with redshift does not capture this residual variation at all masses<fnref n="3" />, which does not exceed ±0.1dex. Since it affects all masses it does not modify the SFR functions but we account for this redshift residual in the cosmic SFRD measurements (Sect. 6) by adding them in quadrature in the error bars. In the highest stellar mass bin, $M_* > 10^{11} M_\odot$, the global shift is due to the fitting procedure, which is slightly above the measurements in Fig. 9. We checked that changing the normalization in the highest mass range does not impact the bright end of the SFR functions and then the results discussed in Sect. 5.

In Fig. 12 we show the behavior of SFR-stellar mass relation as a function of redshift. This correlation is known to be tight and to evolve with redshift (Elbaz et al. 2007; Noeske et al. 2007; Speagle et al. 2014). The mean SFRs per stellar mass bins are estimated for the two 〈IRX〉 − *NrK* calibrations discussed above. We compare our estimates with the local measurements from GALEX-SDSS survey (Salim et al. 2007), the Herschel Reference Survey (Ciesla et al. 2016), deep Herschel observations (Rodighiero et al. 2010; Magnelli et al. 2014; Heinis et al. 2014; Schreiber et al. 2015) and radio 3 GHz stacking (Leslie et al. 2020). Our local, $z \sim 0$, estimate is in excellent agreement with GALEX (Salim et al. 2007) and HRS Ciesla et al. (2016) surveys.

At higher redshift, the normalization of the MS agrees with previous studies up to $z \sim 2$ and over the whole stellar mass range. As previously noticed (Salim et al. 2007; Ilbert et al. 2015; Magnelli et al. 2014; Leslie et al. 2020; Schreiber et al. 2015), a single slope does not provide the best fit as a flattening is observed around $M_* \sim 10^{10} M_\odot$. Indeed, the double slope fit of the MS proposed by Magnelli et al. (2014) or Schreiber et al. (2015) provides a good fit for our SFR-$M_*$ measurements.

# 5. Total SFR functions

In this section we present the construction of the total SFR functions and compare their evolutions to observations and simulations. We use the results from the previous section to correct for dust attenuation.

## 5.1. From the unobscured to total SFR functions

Figure 13 illustrates the impact of our dust correction on the SFR function. In the top panel, we show the unobscured SFR function (based on the observed NUV luminosity function) for different stellar mass bins (colored lines) and for the whole sample (black solid line) in the redshift range $0.6 \le z \le 1$. One noticeable feature is that the unobscured SFR distributions for the three high-mass bins saturate in SFR around a similar value ($L_{NUV} \sim 10^{9.7} L_\odot$). This saturation effect was already reported by Martin et al. (2005) in their analysis of the bi-variate FIR and UV luminosity function in the local Universe. They showed that at low luminosity, the FIR and FUV luminosities track each other. At high luminosity, a saturation effect happens around $L_{UV} \sim 10^{10} L_\odot$, while the $L_{FIR}$ continues to increase by a factor of $\sim 100$. The wide range of dust attenuation at a fixed UV luminosity prevents the derivation of a statistically meaningful dust correction that could be used to translate the unobscured SFR functions into dust-free SFR functions.

By splitting it into stellar mass bins, Heinis et al. (2014) showed that in addition to the trend between dust attenuation (or IRX) and stellar mass, a dependence with $L_{UV}$ is observed, so that fainter $L_{FUV}$ have higher 〈IRX〉 than brighter ones at fixed stellar mass, providing a better perspective to derive the SFR function by estimating the 〈IRX〉 in bins of $L_{UV}$ and stellar mass (see also Bourne et al. 2017). The IRX-*NrK* calibration incorporates those dependences at once. The dust-free SFR function is shown in Fig. 13 bottom panel. While the main effect is a clear separation with stellar mass, the residual correction at fixed stellar mass brings the asymmetric shapes into nearly Gaussian SFR functions per stellar mass bin at the origin of the MS and

<footnote n="3">³ Note that the comparison of our SFR-*NrK* with the SFR from CANDELS at low stellar mass is essentially based on SFRs derived with the *β* slope in UV and not FIR measurement.</footnote>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Mean NrK-based SFR versus stellar mass relation as a function of redshift. The two SFRs derived from the IRX-NrK relations with (red dots) and without (black dots) the stacking results are shown. The completeness limits in stellar mass are shown as vertical orange lines. Comparison with the literature is also shown, with the symbols specified in the inset. The MS at z ~ 0 from Rodighiero et al. (2010) is also reported in each panel.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Impact of dust correction on the shape of the SFR functions. The unobscured (upper panel) and dust corrected (lower panel) SFR functions are shown for different stellar mass bins (in color) and the whole sample (solid and dotted black lines) in the range 0.6 ≤ z ≤ 1.0. One noticeable feature of the dust correction is the ability to dissociate the SFR functions for each stellar mass bin with a larger correction at high masses and to modify the SFR functions per stellar mass bin from an asymmetric to a nearly Gaussian distribution.</figcaption></figure>

its moderate scatter. Adopting a dust correction based only on 〈IRX〉 - M* relationship would translate the unobscured SFR per stellar mass bins but would fail to produce the right shape of the SFR function. Finally, the dust-free SFR function leads to a flatter faint-end slope resulting from a larger correction for the higher stellar mass bin and a more concentrated SFR distribution in a fixed stellar mass bin.

## 5.2. Measurement of the SFR functions

After applying the above dust correction to each individual galaxy, we can derive the SFR functions in different redshift bins up to z = 2. We restrict our sample to SFGs with $i < 26.5$, $K < 24.5$.

The SFR functions are then measured using a $V_{max}$ estimator (Felten 1976; Ilbert et al. 2004) combining the i- and Ks-band selections. Each galaxy is then weighted by $1/V_{max}(z = \min(z_i, z_k))$, the maximum comoving volume within which the galaxy could have been observed in both filters given the limiting magnitudes of our sample and its best-fit template.

To estimate the SFR function uncertainties, we take into account the contribution from the Poissonian errors ($\sigma_{poisson}$) and the cosmic variance due to large-scale density fluctuations ($\sigma_{cv}$):

$$ \sigma_{\text{tot}}^2 = \sigma_{\text{poisson}}^2 + \sigma_{\text{cv}}^2 \quad (5) $$

$$ = \langle N \rangle + \langle N \rangle^2 \times \frac{1}{V^2} \int_V dV_1 dV_2 \xi_{\text{DM}} b^2(m, z) \quad (6) $$

$$ = \langle N \rangle + \langle N \rangle^2 \times b^2(m, \bar{z}) \times \sigma_{\text{DM}}^2(\bar{z}). \quad (7) $$

The cosmic variance term, $\sigma_{cv}$, is estimated by using the Moster et al. (2011) cookbook. The DM cosmic variance, $\sigma_{DM}(\bar{z})$, is given by their Eq. (10) and depends on the redshift and the considered surface. To scale this term to the surface of
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** SFR functions per redshift bin from $z = 0.05$ to $z = 2.05$. The black- and red-filled circles correspond to the HSC-CLAUDS data with and without the stacking correction. The short-dashed black and long-dashed red lines correspond to the best-fitted Schechter functions assuming a fixed slope parameter ($\alpha = -1.3$). The gray area corresponds at the bright end to the Eddington correction and at the faint end to the slope uncertainty ($\Delta\alpha = 0.1$). The orange areas represent the SFR functions based on the SFR derived from the COSMOS2020 FIR data set. In the top panels, the SFR range is translated by 1 dex to make the low-SFR regime of the SFRFs visible.</figcaption></figure>

our survey, we linearly extrapolated their estimate predicted for 2 deg$^2$ to 5.5 deg$^2$. The galaxy bias, $b(m, \bar{z})$, is parameterized in their Eq. (13) for different stellar mass bins that we converted into SFR bins assuming the MS relationship.

Figure 14 shows the SFR functions in ten redshift bins between $0.05 \le z \le 2.05$ and derived with the IRX-NrK

calibration without (black dots) and with (red dots) the stacking correction.

The orange shaded area shows the Vmax weighted SFR functions derived for the COSMOS2020 FIR-selected catalog ($f_{24\mu m} \ge 50\mu Jy$), with the total SFR derived according to Eq. (2).
<table><caption>**Table 2.** Parameter values of the SFR function Schechter fit for both individual FIR calibration and the verified calibration.</caption><thead><tr><th>z<sub>min</sub></th><th>z<sub>max</sub></th><th>Φ<sup>*</sup><sub>No-Stacking</sub><br>(10<sup>-3</sup> Mpc<sup>-3</sup>)</th><th>SFR<sup>*</sup><sub>No-Stacking</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>Φ<sup>*</sup><sub>Stacking</sub><br>(10<sup>-3</sup> Mpc<sup>-3</sup>)</th><th>SFR<sup>*</sup><sub>Stacking</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>Φ<sup>*</sup><sub>DP</sub><br>(10<sup>-3</sup> Mpc<sup>-3</sup>)</th><th>SFR<sup>*</sup><sub>DP</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>σ<sub>DP</sub></th></tr></thead><tbody><tr><td>0.05</td><td>0.25</td><td>3.98 ± 0.20</td><td>2.32 ± 0.18</td><td>3.98 ± 0.20</td><td>2.32 ± 0.18</td><td>12.32 ± 1.60</td><td>0.73 ± 0.25</td><td>0.46 ± 0.06</td></tr><tr><td>0.25</td><td>0.45</td><td>2.78 ± 0.16</td><td>5.90 ± 0.42</td><td>2.78 ± 0.16</td><td>5.90 ± 0.42</td><td>5.67 ± 0.24</td><td>5.97 ± 0.71</td><td>0.27 ± 0.02</td></tr><tr><td>0.45</td><td>0.65</td><td>2.36 ± 0.10</td><td>9.69 ± 0.42</td><td>2.27 ± 0.09</td><td>9.81 ± 0.38</td><td>4.90 ± 0.23</td><td>7.13 ± 0.81</td><td>0.33 ± 0.02</td></tr><tr><td>0.65</td><td>0.85</td><td>2.72 ± 0.13</td><td>13.48 ± 0.65</td><td>2.77 ± 0.12</td><td>13.31 ± 0.58</td><td>6.14 ± 0.43</td><td>9.16 ± 1.51</td><td>0.34 ± 0.03</td></tr><tr><td>0.85</td><td>1.05</td><td>2.14 ± 0.10</td><td>22.97 ± 1.11</td><td>2.35 ± 0.09</td><td>22.09 ± 0.86</td><td>5.75 ± 0.40</td><td>11.55 ± 1.91</td><td>0.39 ± 0.03</td></tr><tr><td>1.05</td><td>1.25</td><td>2.15 ± 0.13</td><td>23.14 ± 1.31</td><td>2.59 ± 0.10</td><td>21.72 ± 0.80</td><td>6.91 ± 0.50</td><td>8.61 ± 1.28</td><td>0.43 ± 0.02</td></tr><tr><td>1.25</td><td>1.45</td><td>1.38 ± 0.08</td><td>34.84 ± 1.75</td><td>1.82 ± 0.10</td><td>32.04 ± 1.63</td><td>3.73 ± 0.22</td><td>21.57 ± 2.70</td><td>0.35 ± 0.02</td></tr><tr><td>1.45</td><td>1.65</td><td>1.67 ± 0.13</td><td>36.72 ± 2.41</td><td>2.43 ± 0.12</td><td>33.80 ± 1.44</td><td>6.41 ± 0.52</td><td>13.23 ± 2.02</td><td>0.44 ± 0.02</td></tr><tr><td>1.65</td><td>1.85</td><td>1.06 ± 0.12</td><td>55.83 ± 4.74</td><td>1.64 ± 0.07</td><td>56.02 ± 1.86</td><td>2.90 ± 0.22</td><td>42.86 ± 5.90</td><td>0.34 ± 0.02</td></tr><tr><td>1.85</td><td>2.05</td><td>0.69 ± 0.07</td><td>71.34 ± 6.88</td><td>0.95 ± 0.03</td><td>83.77 ± 2.34</td><td>2.15 ± 0.14</td><td>47.90 ± 6.62</td><td>0.37 ± 0.02</td></tr></tbody></table>

Notes. The fit parameters of the double power law (noted DP) are also given for the original IRX calibration.

As can be seen, the high-SFR regimes of the HSC-CLAUDS SFR functions are in excellent agreement with the FIR-selected SFR functions. The IRX-NrK calibration does not over/under-predict the density of galaxies with high SFRs, validating our dust correction procedure. On the other hand, since our sample is optical, it allows us to probe the slope of the SFRFs in the low-SFR regime beyond what is currently reached by FIR observations.

## 5.3. Parametric fits of the SFR functions

We fit the SFR functions with a three-parameter Schechter function (Schechter 1976):

$$ \phi(\text{SFR})\text{dSFR} = \phi_{\text{SFR}}^{*} e^{-\frac{\text{SFR}}{\text{SFR}^{*}}} \left(\frac{\text{SFR}}{\text{SFR}^{*}}\right)^{\alpha} \frac{\text{dSFR}}{\text{SFR}^{*}}, \quad (8) $$

where $\alpha$ is the faint-end slope of the power-law regime and $\phi_{\text{SFR}}^{*}$ and $\text{SFR}^{*}$ the characteristic density and SFR separating the power-law and exponential regimes.

Before performing the fit, we account for the Eddington bias (Eddington 1913), due to the uncertainty in the SFR estimate and the exponential cutoff of the SFR function. As a result, more galaxies are shifted toward high SFRs than the reverse, producing a shallower decline in the high-SFR regime than the intrinsic one. To correct for it we adopt the same procedure as Ilbert et al. (2013). We first estimate the uncertainty according to Fig. 11. We adopt an uncertainty of $\sigma = 0.15$ at $z \le 1$ and $\sigma = 0.2$ at $z \ge 1$, assuming that the two SFR estimators used in the comparison have similar errors (i.e., by dividing the observed error by $\sqrt{2}$). We then convolve the Schechter function (Eq. (8)) by the SFR uncertainty:

$$ \phi(\text{SFR})^{\text{conv}} = \phi(\text{SFR}) \star G(\sigma = 0.15, 0.2). \quad (9) $$

A least square fit procedure is performed between $\phi(\text{SFR})^{\text{conv}}$ and the $V_{\text{max}}$ weighted SFR functions up to the completeness limit derived by converting our stellar mass limit into an SFR limit by assuming the MS SFR - $M_{\star}$ relations. The SFR limits are illustrated by a change in symbol size in the different panels.

To estimate the slope of the SFRFs, we fit the lowest redshift SFRF where we get the best constraint. We derive $\alpha = -1.3 \pm 0.1$. Since there is no evidence for the evolution of $\alpha$ up to redshift $z = 2$ in our data within the fit uncertainties (which is also consistent with Ilbert et al. 2015; Mancuso et al. 2015), we simply assume a fixed slope at all redshift. The slope uncertainty

$\sigma(\alpha) = 0.1$ is then propagated into the cosmic SFRD measurements (next section).

The best-fit parameters of the SFRFs are reported in Table 2 (before convolution by the SFR uncertainty) and the final SFRFs for the two calibrations are shown in Fig. 14 (black and red lines). The gray shaded area reflects the slope uncertainties ($\delta\alpha = \pm 0.1$) at the faint end and the Eddington correction at the bright end.

Finally, we also fit the SFR functions by a double power law with the same fixed faint-end slope:

$$ \Phi(\text{SFR}) = \Phi^{*} \left( \frac{\text{SFR}}{\text{SFR}^{*}} \right)^{(1-\alpha)} \exp \left( -\frac{\log_{10}^2 \left( 1 + \frac{\text{SFR}}{\text{SFR}^{*}} \right)}{2\sigma^2} \right). \quad (10) $$

The double-power-law fits give very similar results to the Schechter functions down to $\Phi(\text{SFR}) \sim 10^{-5}$. The parameters are reported in Table 2 but are not shown in Fig. 14 for clarity. As seen in Table 2, the characteristic SFR ($\text{SFR}_{\star}$) exhibits a monotonic decline from ~60 to ~2 $M_{\odot}$ yr$^{-1}$ from $z = 2$ to $z = 0$ and a smooth, though a noisier, increase of the normalization $\Phi_{\star}$ from ~0.001 to ~0.004 Mpc$^{-3}$.

## 5.4. Comparison with observations

In Fig. 14 we compare our results with previous SFR functions from the literature based either on the FIR+UV luminosities or the dust-corrected $H_{\alpha}$ emission lines, and we apply a correction factor for studies using a different initial mass function than Chabrier$^4$.

Bothwell et al. (2011) used the IRAS Faint Source Catalogue and the GALEX All-Sky Imaging Survey (AIS) to perform a combined weighted analysis to derive the SFR function in the local Universe. They also include deep *Spitzer* and GALEX imaging of galaxies in the Local Volume Legacy (LVL) survey ($\le 11$ Mpc) to constrain the faint end of the local SFR function down to SFR < 0.01 $M_{\odot}$ yr$^{-1}$.

Gruppioni et al. (2013, 2015) used the PEP and HerMES surveys of the *Herschel* mission, covering the passbands at 70, 100, and 160 µm (PACS) and 250, 350, and 500 µm (SPIRE), in the COSMOS and GOODS-South fields to measure the infrared luminosity functions up to redshift $z = 4$ with a flux-limited sample at 160 µm. In Gruppioni et al. (2015) they also perform a SED fitting from the UV to the submillimeter to subtract possible

<footnote>$^4$ Most of them were already compiled in Katsianis et al. (2017b).</footnote>
contributions from AGNs to the infrared luminosity and derive the SFRs by summing up the UV and infrared luminosities.

Reddy et al. (2008) used a sample of spectroscopically confirmed Lyman-break galaxies and *Spitzer* MIPS 24 µm observations to derive the SFR (UV+infrared) functions in between $1.9 \le z \le 2.7$ after correction for incompleteness effects. Marchetti et al. (2016) used Herschel observations to infer the infrared luminosity down to $L_{IR} = 10^9 L_{\odot}$ at $z < 0.2$. Similarly, Wang et al. (2016) derive a luminosity function at 250 µm up to $z = 0.5$. Their faint-end slope is consistent with Marchetti et al. (2016). In the following, we only refer to the Marchetti et al. (2016) value.

Ilbert et al. (2015) use a 24 µm-selected sample in the COSMOS and GOODS surveys up to $z = 1.4$. The SFR is estimated by combining the UV+infrared luminosities, and SFR functions are derived by summing up their sSFR functions split per stellar mass bins.

Ly et al. (2011) used emission line galaxies from narrowband imaging at 1.18 µm from the New H$\alpha$ Survey corresponding to H$\alpha$ at $z = 0.85$. They correct for incompleteness and [N II] flux contamination. The SFR is derived by applying a luminosity-dependent dust correction following Hopkins et al. (2001).

Sobral et al. (2013) used four narrowband imaging observations in the UDS and COSMOS fields to select H$\alpha$ emitters at $z = 0.40, 0.84, 1.47$ and $2.23$. The H$\alpha$ luminosity functions were then corrected for incompleteness, [NII] contamination, and dust extinction assuming an average attenuation of A(H$\alpha$) = 1 mag.

Parsa et al. (2016) used the deep fields (HUDF, CANDELS, and UltraVista-COSMOS) to measure the UV luminosity functions at $1.5 \le z \le 2.5$. To convert it into an SFR function, Katsianis et al. (2017b) adopt a luminosity-dependent evolution of the $\beta$-slope as proposed by Smit et al. (2012). These data are the only ones based on a UV selection in this compilation as the authors claimed that this data set provides the best constraint on the slope of the UV luminosity function at this redshift.

As a sanity check, we first compared our SFRFs with the COSMOS2020 ones derived with the same *NrK* method and adopted their photometric redshifts and luminosity estimates. The COSMOS2020-*NrK* SFRFs are shown as black dotted lines. At all redshifts, they are in excellent agreement with our HSC-CLAUDS SFRFs and within our uncertainties.

The COSMOS2020-(UV+FIR) SFRFs derived with the 24 µm flux-limited sample are shown as shaded orange histograms. The two SFRFs are also in good agreement up to the completeness limit. Even though the *NrK* method relies on the COSMOS2020 FIR data for the calibration, when applied to the entire HSC-CLAUDS sample, our *NrK* method reproduces well the SFRF at high SFRs at all redshifts. On the other hand, the UV-optical selection of the HSC-CLAUDS sample spans a wider range of SFRs, extending down to at least a factor of 10 in the low-SFR regime, allowing us to explore the faint end slope.

At low redshift, $0 \le z \le 0.5$, our *NrK* SFR function is in good agreement with the FIR+UV SFRF obtained in the local volume by Bothwell et al. (2011, blue line) and exhibits a comparable faint-end slope ($\alpha = -1.41$ with the Vmax estimator). It is also in good agreement with Marchetti et al. (2016, thin red stars) after converting their FIR luminosity function into a SFR function. While we observe a good match at the bright end, their faint-end slope is flatter, as expected since it neglects the contribution of faint UV sources. The 160 µm-selected SFRFs (red light and dark stars, respectively; Gruppioni et al. 2013, 2015) do not probe the faint end<fnref n="5" /> but their normalizations appear consistent

with us around SFR ~1–3 $M_{\odot}$ yr$^{-1}$. However, they overpredict the high-SFR end with respect to all the other FIR+UV and *NrK* measurements. This excess could be due to the FIR photometric extraction, where we adopt the super-deblended FIR photometry in the COSMOS field (Jin et al. 2018, see their Sect. 2.1.2), and/or to the photometric redshift estimates. We note that the Gruppioni et al. (2015) SFRF leads to a high SFRD (see Sect. 6) due also to a steep faint-end slope. This is not the case for the Gruppioni et al. (2013) SFRF, based only on the FIR SFRF with a slope $\alpha = -1.2$. The SFRF from H$\alpha$ by Sobral et al. (2013, dark green triangles) at $z \sim 0.4$ shows a steeper slope and a deficit at high SFRs, which could be attributed to a unique and averaged dust correction factor (A(H$\alpha$) = 1 mag) applied (see below).

At intermediate redshifts, $0.5 \le z \le 1.5$, all the FIR+UV SFRFs are consistent with each other as do our *NrK* SFRFs. To extend the SFRFs in the low-SFR regime, Ilbert et al. (2015, open blue circles) have included the contribution of low-mass galaxies by assuming that the shape of the sSFR function at low masses is the same as for their lowest stellar mass measurement ($\log_{10}(M_{\star}) = 9.5–10$) and by normalizing the sSFR with the density of the star-forming GSMFs at the appropriated redshift. This leads to a slope of the SFRF in excellent agreement with our estimate ($\alpha \sim -1.3$) with no sign of evolution up to $z \sim 1.5$. The H$\alpha$ SFRF from Ly et al. (2011, light green triangles) at $z \sim 0.85$ (shown in $z = 0.8$ and $z = 1.0$ panels) shows a good agreement with the other measurements. They reproduce the high-SFR distribution and the faint end slope remarkably well compared to the one from Sobral et al. (2013, dark green triangles). This comes from the different dust correction treatments. Ly et al. (2011) adopted a luminosity-SFR-dependent correction as proposed by Hopkins et al. (2001) with large or small correction for respectively high or low luminosity (similar to what we observe with stellar mass in Fig. 10), changing significantly the shape of the original H$\alpha$ luminosity function.

At high redshifts, $1.5 \le z \le 2$, our SFRFs are still in good agreement with Gruppioni et al. (2013, 2015) at the bright end, while H$\alpha$ SFRF slightly underestimate the high star-forming population. At $z \sim 2$, the SFRFs start to differ in different regimes. The UV-selected sample (Parsa et al. 2016, dark blue stars) shows a significant shortage of high SFRs. This shortage is most likely a consequence of the uncertainty in the dust correction, especially for the most luminous as discussed at the beginning of Sect. 5. Adopting a $\beta$-slope varying only with luminosity cannot properly capture the wide scatter of dust attenuation at high UV luminosities (Martin 2005) and then cannot properly reproduce the high end of the SFR function. On the other hand, the UV sample explores the low-SFR regime. They derive a faint end slope consistent with our value despite a higher density normalization. The Lyman break galaxies (LBG) selected sample (Reddy et al. 2008, light blue stars), with SFR derived from UV+FIR, also shows a higher normalization of the SFRF. In contrast to the UV and LBG samples, our HSC-CLAUDS sample has a stellar mass limit above $\log_{10}(M_{\star}/M_{\odot}) = 9.5$ at this redshift. We thus can miss the potential contribution of lower-mass galaxies in the SFRF around SFR~10 $M_{\odot}$ yr$^{-1}$ and produce a flattening of the faint-end slope. We also note that the LBG and UV SFR functions are derived in a much wider redshift bin ($\Delta z \sim 1$), which can impact the comparison, and a nontrivial correction for incompleteness is required for the LBG sample to assess the global SFR function.

Finally, despite the very deep data set used in this work, our optical/NIR-selected sample can potentially miss heavily

<footnote n="5">⁵ Due to the lower sensitivity at 160 µm, their flux-limited sample introduces a brighter SFR cutoff than the 24 µm-selected samples.</footnote>
obscured galaxies such as submillimeter galaxies (Chapman et al. 2005) or dark-HST galaxies (Wang et al. 2019). This population can contribute to the cosmic SFRD at high redshifts, but its comoving density is expected to be low in our redshift range of interest ($z \le 2$, Chapman et al. 2005). It could help us better match our SFR functions above SFR $\ge 100 M_{\odot}$ yr$^{-1}$ with those of Gruppioni et al. (2015). Considering our bright-end SFR functions as lower limits, it is an even more stringent test for the comparison with simulations discussed in the next section.

In conclusion, the *NrK* method presented in this work allows, for the first time, the SFR functions to be measured for a wide range of SFRs. The derived SFR functions can reproduce the number density of high-SFR galaxies as observed with FIR samples as well as the slope at low SFRs, which is found to be relatively shallow ($\alpha \sim -1.3$), with no evolution at least up to $z \sim 1.6$ and potentially $z \sim 2$, according to the UV-selected sample. This method overcomes the current limitations of the other approaches (i.e., at the faint end for the FIR samples due to instrumental sensitivity and the dust treatment especially at the high-SFR end for the $H_{\alpha}$- and UV-selected samples).

## 5.5. Comparison with simulations

Star formation rate functions give an instantaneous view of the distribution of the in situ star formation at different epochs. It is a more stringent test for the models than the GSMF since the latter captures an integrated view of the past star formation activity.

In terms of SAMs, Gruppioni et al. (2015) already made a comparison of their observed SFR functions with several SAMs and found an overall good agreement with the bright end of the SFRFs up to $z = 2$, while the models fail to reproduce high star-forming systems at $z > 2$. For hydrodynamical simulation, Katsianis et al. (2017b) used EAGLE simulation and observed a deficit of high star-forming simulated galaxies at $z < 2$.

We aim here to make a broader comparison with several state-of-the-art cosmological hydrodynamical simulations. In this section we confront our SFR functions with four hydrodynamical simulations: SIMBA (Davé et al. 2019), HORIZON-AGN (Dubois et al. 2014), EAGLE (Crain et al. 2015; Schaye et al. 2015; McAlpine et al. 2016), and TNG100 from the ILLUSTRISTNG project (Pillepich et al. 2018; Nelson et al. 2019).

### 5.5.1. Main ingredients in the simulations

All these simulations incorporate different prescriptions to form stars, treat the stellar and black hole (BH) feedback, and adopt different observables at $z = 0$ to fine-tune the sub-grid physics models. The simulations are described in more detail in Appendix B and their main features are summarized in Table 3. Here we highlight some of the main differences that can have an impact on the SFR functions, which is the main topic of this paper.

SIMBA is the only one to model on the fly the formation, growth, and destruction of dust, and to introduce X-ray heating from BHs in addition to the regular AGN feedback (see Appendix B). The importance of X-ray heating has been explored in zoom simulations by Choi et al. (2012), showing that it can potentially drive the quenching of massive galaxies. Indeed, as shown in Davé et al. (2019, 2020), while the X-ray feedback has a minimal effect on the galaxy mass function, it represents an important additional energy input to fully quench

massive galaxies. This leads to a quenched galaxy population (Rodríguez Montero et al. 2019) with reduced central molecular gas (Appleby et al. 2020), and the galaxy–BH co-evolution in good agreement with observations (Thomas et al. 2019).

While all of these simulations successfully reproduce the bulk massive galaxy properties, this is achieved with very different AGN feedback models. SIMBA employs kinetic two-mode feedback: The radiative mode, or the so-called quasar mode, is designed to model outflows of cold molecular and warm ionized gas. The jet mode, often referred to as the radio mode, drives high-velocity collimated jets of hot gas in a direction given by the angular momentum of the inner disk.

The energy used to drive the AGN feedback and that serves to quench galaxies originates from the accretion energy onto BHs. SIMBA is unique in using a two-mode accretion model consisting of the torque-limited accretion from cold gas and the more standardly used Bondi accretion from hot gas. Similarly, both HORIZON-AGN and TNG100 employ Bondi accretion with a two-mode feedback model; however, they adopt spherical thermal feedback at high Eddington growth rates. In HORIZON-AGN, the radio mode deposits kinetic energy into a bipolar jet along the accreted gas angular momentum (similar to SIMBA), while in TNG100 it is in a random direction. In contrast, EAGLE does not make a distinction between the two modes, it assumes Bondi accretion model onto BHs and thermal energy injection follows the quasar-mode feedback scheme.

Finally, while all simulations rely in one way or another on the Schmidt relation to form stars (i.e., SFR proportional to the gas density over the dynamical or free-fall timescale), the models differ in detail. While HORIZON-AGN and TNG100 use the density of atomic hydrogen to trigger star formation, SIMBA relies on the molecular hydrogen, with comparable density thresholds and efficiencies, and EAGLE adopts a metallicity-dependent density threshold for star formation.

### 5.5.2. Descriptions of the snapshots

To perform the comparison with the observed SFRFs, we analyzed the snapshots at five different redshifts, $z = 0, 0.5, 1, 1.5, 2$. We restricted the analysis to galaxies identified by the galaxy finders used for different simulations, and therefore composed of at least 30 star particles (50 for HORIZON-AGN). This corresponds to a stellar-mass limit of $\log_{10}(M_{\star}/M_{\odot}) = 8.7, 8.2, 7.7$, and 7.6 for SIMBA, HORIZON-AGN, EAGLE, and TNG100, respectively. The SFR is typically estimated based on the number of stars formed over a certain period of time (frequently 100 Myr; Dubois et al. 2014).

We explicitly verified for the HORIZON-AGN and TNG100 simulations that using shorter timescales (e.g., 10 Myr or 50 Myr) does not alter our conclusions. For SIMBA, SFR is computed from the gas particles, corresponding to the SFR computed from the stellar particles averaged over a timescale between 50 and 100 Myr (Davé et al. 2019). The choice of 100 Myr roughly corresponds to a minimum measurable SFR ($M_{\star,init}/100$ Myr) of $0.18 M_{\odot}$ yr$^{-1}$ for SIMBA and $0.02 M_{\odot}$ yr$^{-1}$ for the other simulations. In EAGLE, stellar mass and SFR are computed using particles within a fixed spherical aperture with a radius of 30 proper kpc (see, e.g., Crain et al. 2015). Given that the majority of star formation occurs in the central 30 kpc, even for massive galaxies, this aperture constraint has only a minimal effect on the SFRs (see, e.g., Furlong et al. 2015). The effect of the aperture on the GSMF is negligible for galaxies with $M_{\star} < 10^{11} M_{\odot}$; however, for more massive galaxies the aperture reduces the stellar masses (see, e.g., Schaye et al. 2015). Finally,
<table><caption>**Table 3.** Summary of the different simulations used in this work: SIMBA, HORIZON-AGN, EAGLE, and TNG100.</caption><thead><tr><th>Parameters</th><th>SIMBA</th><th>HORIZON-AGN</th><th>EAGLE</th><th>TNG100</th></tr></thead><tbody><tr><td>Ω<sub>m</sub> − Ω<sub>b</sub></td><td>0.3–0.048</td><td>0.272–0.045</td><td>0.307–0.048</td><td>0.309–0.048</td></tr><tr><td>H<sub>0</sub></td><td>68 km s<sup>−1</sup> Mpc<sup>−1</sup></td><td>70.4 km s<sup>−1</sup> Mpc<sup>−1</sup></td><td>67.77 km s<sup>−1</sup> Mpc<sup>−1</sup></td><td>67.74 km s<sup>−1</sup> Mpc<sup>−1</sup></td></tr><tr><td>σ<sub>8</sub> − n<sub>s</sub></td><td>0.82–0.97</td><td>0.81–0.97</td><td>0.83–0.96</td><td>0.82–0.97</td></tr><tr><td>Code</td><td>GIZMO</td><td>RAMSES</td><td>GADGET-3</td><td>AREPO</td></tr><tr><td>Box size</td><td>100 h<sup>−1</sup> Mpc</td><td>100 h<sup>−1</sup> Mpc</td><td>67.77 h<sup>−1</sup> Mpc</td><td>75 h<sup>−1</sup> Mpc</td></tr><tr><td>Particles</td><td>1024<sup>3</sup></td><td>1024<sup>3</sup></td><td>1504<sup>3</sup></td><td>1820<sup>3</sup></td></tr><tr><td>Star particle mass<sup>(a)</sup></td><td>1.82 × 10<sup>7</sup> M<sub>☉</sub></td><td>2 × 10<sup>6</sup> M<sub>☉</sub></td><td>1.81 × 10<sup>6</sup> M<sub>☉</sub></td><td>1.4 × 10<sup>6</sup> M<sub>☉</sub></td></tr><tr><td>Gas particle mass</td><td>1.82 × 10<sup>7</sup> M<sub>☉</sub></td><td>–</td><td>1.81 × 10<sup>6</sup> M<sub>☉</sub></td><td>1.4 × 10<sup>6</sup> M<sub>☉</sub></td></tr><tr><td>DM particle mass</td><td>9.6 × 10<sup>7</sup> M<sub>☉</sub></td><td>8 × 10<sup>7</sup> M<sub>☉</sub></td><td>9.7 × 10<sup>6</sup> M<sub>☉</sub></td><td>7.5 × 10<sup>6</sup> M<sub>☉</sub></td></tr><tr><td>Gravitational<br>softening/resolution</td><td>0.5 h<sup>−1</sup>kpc down to 1<br>proper kpc</td><td>0.8 proper h<sup>−1</sup>kpc</td><td>1.8 h<sup>−1</sup>kpc down to 0.7<br>proper kpc</td><td>0.5 h<sup>−1</sup>kpc</td></tr><tr><td>Galaxy finder (n<sub>p</sub>)<sup>(b)</sup></td><td>6D FoF (30)</td><td>ADAPTAHOP (50)</td><td>FoF &amp; SUBFIND (30)</td><td>FoF &amp; SUBFIND (30)</td></tr><tr><td>Mass completeness<sup>(γ)</sup></td><td>log<sub>10</sub>(M<sub>comp</sub>/M<sub>☉</sub>) = 8.7</td><td>log<sub>10</sub>(M<sub>comp</sub>/M<sub>☉</sub>) = 8.2</td><td>log<sub>10</sub>(M<sub>comp</sub>/M<sub>☉</sub>) = 7.7</td><td>log<sub>10</sub>(M<sub>comp</sub>/M<sub>☉</sub>) = 7.6</td></tr><tr><td>SNe feedback</td><td>Winds-kinetic</td><td>Winds-kinetic</td><td>Thermal</td><td>Winds-kinetic</td></tr><tr><td>AGN feedback</td><td>Radio mode-kinetic &amp;<br>X-ray heating, quasar<br>mode-kinetic</td><td>Radio mode-kinetic and<br>quasar mode-thermal</td><td>Thermal</td><td>Radio mode-kinetic,<br>quasar mode-thermal</td></tr><tr><td>Feedback calibration</td><td>z = 0 galaxy stellar mass<br>function</td><td>z = 0 black hole-galaxy<br>scaling relation</td><td>z = 0 GSMF, galaxy<br>sizes, M<sub>⋆</sub>−M<sub>BH</sub> relation</td><td>SFRD evolution,<br>GSMF, SHMR at z = 0</td></tr><tr><td>Star formation</td><td>Proceeds in n<sub>H₂</sub> ><br>0.13 cm<sup>−3</sup> regions,<br>following a Schmidt rela-<br>tion with 2% efficiency</td><td>Proceeds in n<sub>H</sub> ><br>0.1 cm<sup>−3</sup> regions, follow-<br>ing a Schmidt relation<br>with 2% efficiency</td><td>Metalllicity dependant<br>density threshold repro-<br>ducing KS relation</td><td>Proceeds in n<sub>H</sub> ><br>0.1 cm<sup>−3</sup> regions, fol-<br>lowing a Schmidt rela-<br>tion with 2% efficiency</td></tr><tr><td>SFR estimator</td><td>~50–100 Myr</td><td>100 Myr</td><td>30 kpc aperture</td><td>100 Myr<sup>(δ)</sup></td></tr></tbody></table>

**Notes.** The three blocks respectively recap the cosmological parameters, general simulation parameters, and baryonic physics, including models for feedback and star formation. (a) Initial mass particle. (b) Minimal number of particles for a galaxy to be found. (γ) Galaxy stellar mass completeness based on the initial stellar mass particle and the minimum number of particles used by the galaxy finder. (δ) Donnari et al. (2019); Pillepich et al. (2019).

we note that the SFR considered in each simulation is not corrected for the mass loss due to winds and SNe. However, we checked that taking this correction into account has a minimal impact on SFRF.

Before performing our comparisons, we exclude the quiescent population by using a unique redshift-dependent criterion for all four simulations as defined in Davé et al. (2019), and based on the sSFR: $log_{10}(sSFR/yr^{-1}) = -10.8 + 0.3z$. This is equivalent to removing galaxies with $\Delta$SFR ~ 1 dex below the MS. We verified that this cut is consistent with the selection of SFGs based on the NUVrK color-color diagram. We note that, because the impact of passive galaxies is much less significant on the SFR functions than on stellar mass functions, the exclusion of simulated passive galaxies is often not applied when performing the comparison of the SFR functions with observations (Katsianis et al. 2017a).

### 5.5.3. Comparison of the simulated SFRFs with observations

In Fig. 15 we compare the SFR functions of SFGs obtained for the four simulations with the Schechter fits of our HSC-CLAUDS sample in five redshift bins. The differences between the simulations are quite noticeable, especially at increasing redshift. To quantify those differences, we fit each SFRF by a Schechter function, leaving the slope as a free parameter. We perform the fit down to the SFR limit (indicated as thick marks in each panel) as implied by the mass limit and MS relation of each simulation (see Fig. 15, bottom row). The data points and Schechter fits of the simulations are shown with dots and

solid lines, respectively. The parameters of the Schechter fits are shown in Fig. 16, along with those derived from the observations. We note that for TNG100 and SIMBA simulations at $z = 0$ and $z = 0.5$, we adopt a double-power-law fitting function to better reproduce their bright ends, even if the parameters derived with a Schechter function are used for consistency in Fig. 16.

Figures 15 and 16 reveal three trends regarding the general and/or specific behaviors of the simulations with respect to the observations. First, the SFRFs from TNG100 and EAGLE simulations show very similar behaviors at all $z$ and are in between the most extreme behaviors of the HORIZON-AGN and SIMBA simulations.

Second, while all simulations reproduce the high-SFR side at low $z$ ($z \le 0.5$), SIMBA is the only one able to reproduce its evolution with redshift up to $z \sim 2$, as can be seen in the evolution of SFR\* in Fig. 16. This is consistent with the analysis from Lovell et al. (2021), which shows that SIMBA is the only simulation able to reproduce the submillimeter galaxy number counts at $z \sim 2$. In contrast, HORIZON-AGN shows a very mild evolution of SFR\* with an increase by less than a factor of 2, while the observations suggest a factor of ~15 between $z = 0$ and $z = 2$. TNG100 and EAGLE are in between and show a shortage of high-SFR galaxies, as reflected by a low SFR\*, at $z \ge 1$.

Finally, the comoving density of galaxies with low to intermediate SFRs varies significantly between simulations. At $z \sim 2$, they all agree with a slope $-1.5 \le \alpha \le -1.3$, consistent with the observations. At later times, while the slope in HORIZON-AGN
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** SFR functions, stellar mass functions, and MS predicted by simulations: Top row: Comparison of the SFR functions for the SFGs (gray-shaded area based on the faint end slope uncertainty) with those from the four simulations in five redshift bins: TNG100 (blue), HORIZON-AGN (red), SIMBA (green), and EAGLE (orange). On the top axis, the vertical marks correspond to the SFR limits used to fit the observations (black) and the simulations (colored marks) as described in Sects. 5.3 and 5.5. The data points from the simulations are fitted by a Schechter or double-power-law function (thick lines; see text). Second row from top: GSMF for the star-forming population of the four simulations (colored lines) with a unique star-forming criterion from Davé et al. (2019, see text). The gray-shaded area is based on the published star-forming GSMFs in the COSMOS field (see text). Second row from bottom: Same as above but for the whole (star-forming + quiescent) population. The gray area shows the observed total GSMF from the COSMOS field (Weaver et al. 2022) except at $z = 0$, where we use Bernardi et al. (2013) and Li & White (2009). Bottom row: Evolution of the star-forming MS for the COSMOS2020 data set with the SFR based on the *NrK* method (gray shaded area; see text) and the simulations (colored dots).</figcaption></figure>

shows almost no evolution, in EAGLE and TNG100 it gradually flattens, but yet within the range of the observations (except at $z \sim 0$, where $(\alpha + 1) > 0$ for TNG100). In contrast, the evolution is much more pronounced for SIMBA, with $(\alpha + 1) = 0$ already at $z = 1$ and ends up with a positive slope ($(\alpha + 1) \gg 0$), that is to say, with a shortage of galaxies with low SFRs at $z < 1$. We note that the lower mass resolution of SIMBA should not be responsible for this effect at least at $z \ge 0.5$.

In summary, none of the simulations analyzed here is capable of reproducing the observed SFR functions of SFGs at all redshifts from $z = 2$ down to $z = 0$. Three of them, HORIZON-AGN, TNG100, and EAGLE, fail to reproduce the density of highly SFGs at $z \ge 0.5$, but reproduce reasonably well the low to inter-

mediate SFR regime. On the contrary, SIMBA nicely reproduces the high-SFR regime but fails in the low-SFR regime below $z = 1.5$. This highlights that it remains challenging to achieve the balance between the star-forming “MS” galaxies and the population of quenched galaxies over cosmic time despite recent attempts at using a more physically motivated prescription for star formation and various feedback processes. To clarify the differences in the SFR functions, in Fig. 15 we also include additional information with the comparisons of the GSMFs for the SFGs (second row from top), the whole population (second row from bottom) as well as the evolution of the MS (bottom row).

**The global GSMF.** To ease the comparison with the simulations in Fig. 15 (2nd row from bottom), the observed total
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Redshift evolution of the Schechter parameters – $\alpha$ (top panel), SFR\* (middle panel), and $\Phi^*$ (bottom panel) – of the SFR functions for the observations (black dots) and the four hydrodynamical simulations (color dots as indicated in the top panel). The fits of the redshift evolution for each parameter from the observations are shown as dotted lines.</figcaption></figure>

GSMFs are represented as a gray-shaded region. At $z \sim 0$, the gray area encloses the GSMF from Li & White (2009; bottom limit) and Bernardi et al. (2013; top limit), while at higher redshifts it is based on the total GSMFs in the COSMOS survey (Weaver et al. 2022). The same color code is used for the simulations. The simulations are often tuned to reproduce the observed GSMFs in particular at $z \sim 0$. This is indeed the case of the simulations analyzed here. The only exception is the HORIZON-AGN simulation, which is not calibrated in the local Universe apart from the choice of BH feedback parameters that reproduce the local BH mass versus stellar velocity dispersion relation. This possibly explains why the HORIZON-AGN overpredicts at all epochs the GSMFs density (Kaviraj et al. 2017), and particularly below $10^{11} M_{\odot}$. At $z \ge 0.5$, EAGLE underpredicts the density of high stellar mass galaxies (see also Furlong et al. 2015) and slightly overpredicts (similarly to ILLUSTRISTNG) the density of low stellar mass galaxies. At all stellar masses, SIMBA agrees reasonably well with the observations.

The slight over-prediction of massive galaxies ($M \ge 10^{11} M_{\odot}$) by HORIZON-AGN and ILLUSTRISTNG simulation below $z = 0.5$ may be due to an inefficient AGN feedback that does not quench star formation enough in these massive systems, while the over-prediction of galaxies at the low-mass end and $z \ge 0.5$ for all simulations but SIMBA can instead be attributed to an insufficient SN feedback model.

**The star-forming GSMF.** An observable that provides us with better insight into the ability of simulations to properly model different galaxy populations is the GSMF for SFGs. As it has not been specifically calibrated against observations in any of the simulations analyzed here, it may also better discriminate between them. The GSMF for SFGs is shown in Fig. 15 (2nd from top). Simulations exhibit a different level of agreement with the observed GSMFs in a redshift-dependent way.

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Cosmic SFRD. The blue dots correspond to the averaged SFRD resulting from the two IRX calibrations, with the upper and lower limits of the blue-filled area corresponding to the calibration with and without stacking, respectively. The error bars on the blue dots are computed using the faint-end slope uncertainty of the SFR functions, $\alpha = -1.3 \pm 0.1$. In the top panel, SFRDs from other observational studies are added with colors depending on their estimator: UV (blueish), infrared (reddish), and radio (black). In the bottom panel, SFRDs measured from the four hydrodynamical simulations are added. In both subplots, we also add the contributions to the cosmic SFRD of different stellar-mass and SFR regimes (color-shaded areas).</figcaption></figure>

**HORIZON-AGN** has quite similar behavior at all redshifts. While it manages to reproduce the very massive end, it significantly overestimates the density of SFGs below $M_{\star} = 10^{11} M_{\odot}$. At $z = 0$, SIMBA, EAGLE, and TNG100 reproduce the low-mass end of the star-forming GSMF reasonably well; however, they underpredict the density of massive SFGs. At $z > 0.5$, and as the redshift increases, the trend of each simulation accentuates. EAGLE and TNG100 slightly overproduce the low-to intermediate-mass population and EAGLE underestimates the massive end. HORIZON-AGN continues to overestimate the GSMF in all but the most massive galaxies. Interestingly, SIMBA reproduces reasonably well all the mass regimes.

In summary, HORIZON-AGN overpredicts the number of low- to intermediate-mass SFGs at all redshifts, similar to the global population. This is consistent with a possible interpretation that SN feedback quenching is too inefficient for intermediate- to low-mass galaxies. EAGLE and TNG100 overpredict the low-mass end for the star-forming and the global population at $z \ge 0.5$, and EAGLE is the only one simulation to under-predict the number of SFGs at $z > 1.5$. Overall SIMBA appears to be the only simulation able to reproduce relatively well the GSMFs for the global and star-forming populations.

**The star-forming main sequence.** Additional key observable for galaxy formation models is the relation between the SFR and stellar mass of galaxies. The cosmological simulations have been known to under-predict the amplitude of the SFR–$M_{\star}$ relation at high redshifts ($z \sim 2$) by a factor of a few, while the agreement improves at low redshift ($z \le 0.5$). This can be seen in Fig. 15 (bottom row), where the MS for all the simulations is compared to the observations based on the COSMOS2020 data set. The SFR is derived from the *NrK* method and the width of the gray-shaded area corresponds to the observed scatter.
At $z = 0$, all the simulations are close to the observations except for HORIZON-AGN and EAGLE, which slightly underestimate the SFR at low masses, $M_{\star} \sim 10^{9.5} M_{\odot}$. At increasing redshift, the three simulations, HORIZON-AGN, EAGLE, and ILLUSTRISTNG, gradually deviate from the observed MS. The departure is more pronounced at low masses, where it reaches up to a factor of 10 underestimation of the SFR at $z = 2$ for galaxies with $M_{\star} \sim 10^{9} M_{\odot}$ and a factor of ~3 for galaxies with $M_{\star} \ge 10^{10} M_{\odot}$. This trend is much less pronounced with SIMBA, which remains consistent with the observed MS at high masses up to $z = 2$ and deviates from the MS mainly for the low-mass regime, $M_{\star} \le 10^{10} M_{\odot}$.

While a natural reason for this offset on the side of modeling could be a too strong stellar feedback in low-mass galaxies, and in the case of HORIZON-AGN, EAGLE and TNG100, possibly also AGN feedback at the high-mass end at high redshift ($z \ge 1$), this seems unlikely considering the excess of simulated galaxies at intermediate- to low-mass in the GSMFs. This is probably not the whole story as simply boosting the SFR would for example bring HORIZON-AGN in even bigger disagreement with the observed cosmic SFRD at all redshifts, and it would induce a disagreement for other simulations at $z \le 0.5$, where they match the average SFRD derived in this work (see Fig. 17, bottom panel). Figure 15 also suggests that systematic offset does not explain the discrepancy in the SFR-$M_{\star}$ relation as it would simply shift the SFR functions toward higher values, inducing a further disagreement with observed SFRF at all SFRs at $z = 0$ for all simulations, at the bright end for SIMBA at all redshifts, and at the faint end and at intermediate SFR for all simulations at $z \ge 0.5$. Alternatively, as suggested by Furlong et al. (2015), a potential solution to low SFRs is an insufficient burstiness of star formation. Making star formation burstier could result in a higher SFR over shorter periods compared to the current models without significantly modifying the stellar mass of galaxies. It remains to be seen whether this solution could at the same time solve the offset of the SFR-$M_{\star}$ relation, and the discrepancies at the bright and low end of SFR functions while keeping GSMFs in relatively good agreement with observations.

Another interesting aspect is the scatter of the star-forming MS. This scatter is significantly broader for SIMBA ($\sigma \sim 0.4$) compared to the other simulations ($\sigma \sim 0.3$). This trend tends to broaden the SFR functions per stellar mass bin and could be in part responsible for the underestimation of the low-to-intermediate SFR end and the flattening of the faint end slope in SIMBA's SFRF.

In conclusion, this set of simulations, although not exhaustive, reflects today's state-of-the-art galaxy formation models and shows the complexity of simulating distributions of galaxies' SFR that match observations. Indeed, compared to stellar mass, SFR is a more instantaneous parameter and is subject to much more stochasticity. It is therefore not surprising that simulations do not yet fully agree with observations. Further SFR function predictions from simulations will provide an efficient testing ground, in addition to more standard used stellar mass functions and star-forming MS, for converging toward more realistic feedback mechanisms implementation in simulation.

## 6. Cosmic SFRD

The cosmic SFRD, $\psi$, is derived by integrating the SFR Schechter functions as

$$ \psi(z) = \int_{0}^{\infty} \text{SFR } \Phi(\text{SFR}, z) \, \text{dSFR.} \qquad (11) $$

The redshift evolution of $\psi(z)$ is shown in Fig. 17. The light blue shaded region encloses the integration of the SFR functions based on the IRX calibration with and without stacking technique (lower and upper bounds, respectively). The open circles are the average SFRD resulting from the two calibrations, while the error bars reflect the impact of a changing slope around $\alpha = -1.3 (\Delta\alpha = \pm 0.1)$ and redshift residuals in the SFR calibration mentioned in Sect. 4.4.

Our measurements are compared with the SFRD from the literature obtained with different SFR estimators – radio-based (Malefahlo et al. 2022; Karim et al. 2011; Leslie et al. 2020), UV-based (Schiminovich et al. 2005; Dahlen et al. 2007; Cucciati et al. 2012), and FIR-based (Sanders et al. 2003; Magnelli et al. 2011, 2013) – and the compilation fit from Madau & Dickinson (2014, thick black line). There is an overall good agreement with a gradual decline of the cosmic SFRD since $z \sim 2$, the cosmic noon, up to the present day by a factor of 10. While Madau & Dickinson (2014) predicts a decline scaling as $\psi(z) \propto (1+z)^{2.7}$ since $z = 1.5 - 2$, our observations suggest a slightly steeper rate with $\psi(z) \propto (1+z)^{2.8}$ when using the same parametric form as their Eq. (15).

The redshift evolution of $\psi(z)$ might be better captured by the parametric form proposed by Katsianis et al. (2021), which closely mimics the evolution of the gas reservoir model (Bouché et al. 2010) described with two parameters (see Eq. (16) in Katsianis et al. 2021). The gray-filled area shows this alternative parametric form, with the upper bound corresponding to their parameters while the lower bound is an adaptation to fit our observations. This suggests a faster evolution in between $0 \le z \le 1$ followed by a plateau between $1 \le z \le 2$ rather than a peak at $z \sim 2$, in better agreement with our measurements and recent FIR estimates (Gruppioni et al. 2013, 2020; Katsianis et al. 2021).

We also show the contribution to the SFRD of different populations split into SFR (top panel) and stellar mass (bottom panel) bins. The star-forming populations contributing the most to the SFRD evolve with redshift. At low-$z$, $z \le 0.5$, the main contributors are galaxies with low star formation activity ($0.3 \le \text{SFR} \le 3 M_{\odot} \text{yr}^{-1}$) while moderate SFGs ($3 \le \text{SFR} \le 30 M_{\odot} \text{yr}^{-1}$) take over at higher redshift. The contribution of the high star-forming population ($\text{SFR} \ge 30 M_{\odot} \text{yr}^{-1}$) shows a sharp increase from $0.5 \le z \le 2$ contributing equally at $z \sim 2$ than the moderate star-forming population, in qualitative agreement with previous FIR studies (Le Floc'h et al. 2005). When split into stellar mass bins, the main contribution to the SFRD comes from intermediate-mass galaxies ($9.5 \le \log(M_{\star}/M_{\odot}) \le 10.5$), namely, the population below the knee of the GSMF, at all redshift. All the SFRDs from the three stellar mass bins show a global decline with cosmic time but at different paces. By Adopting the SFRD parametrization from Madau & Dickinson (2014), we find an evolution rate, $\psi(z) \propto (1+z)^{\alpha}$, with $\alpha = 2.2$ for the lowest-mass population and $\alpha = 3.2$ for the most massive one. These trends are consistent with the downsizing picture (Cowie et al. 1996; Juneau et al. 2005) in which most massive DM halos have a rapid accretion rate at an early time followed by a quenching phase after their star formation onset. Low-mass systems accrete at a later time but are also sensitive to the global decline of the accretion rate at all mass at a later time, preventing them from becoming more dominant at low $z$ (e.g., Bouché et al. 2010). This picture is consistent with Gruppioni et al. (2013), based on the infrared luminosity functions, but slightly differ in the relative contribution to the total SFRD. When adopting their stellar mass bins, we both find that the high-mass bin $11 \le \log(M_{\star}/M_{\odot}) \le 12$ never exceeds 15%
of the total SFRD. While we find that the two low-mass bins ($\text{Log}(M_*/M_\odot) \le 10$ and $10 \le \text{Log}(M_*/M_\odot) \le 11$) contribute equally at ~50–40%, Gruppioni et al. (2013) find a higher contribution from the intermediate-mass bin. But we note that the parameters of the infrared luminosity function in their lowest-mass bin are poorly constrained at increasing redshift.

In the bottom panel of Fig. 17, we show the SFRD of the simulations by integrating their respective Schechter functions discussed in Sect. 5.5. The evolution of $\psi(z)$ for HORIZON-AGN is in excellent agreement with the observations above $z = 0.5$, despite noticeable differences in the shape of the SFR functions with the observations. The lack of evolution in the SFR* parameter is compensated for by a high normalization. All the other simulations are slightly below the observations. For the TNG100 and EAGLE simulations, this is essentially due to the scarcity of high SFGs at high redshift, while for the SIMBA simulation, which is the only one to reproduce the density of high SFGs, it is caused by the lack of intermediate star-forming systems. This illustrates that integrated quantities such as the cosmic SFRD (or stellar mass density) are a necessary test to be passed by the simulations but not sufficient to understand the origin of any discrepancies without the full characterization of the SFR (or stellar mass) functions.

## 7. Conclusion

In this paper we investigate the evolution of the SFR functions up to $z \sim 2$, based on stellar-mass-selected samples from deep HSC-CLAUDS observations. To derive the SFR of individual galaxies, we used an original method based on the UV-infrared energy budget that does not rely on SED fitting or uncertain assumptions regarding the dust attenuation laws. This energy budget is calibrated by analyzing the behavior of the IRX, $\text{IRX} = L_{\text{IR}}/L_{\text{UV}}$, in the rest-frame color-color diagram, $(\text{NUV}-r)$ versus $(r-K_s)$, for galaxies observed in the UV and FIR in the COSMOS2020 catalog. We show that the use of a single vector, $\text{NrK}$, is effective in predicting the IRX with a small scatter ($\sigma < 0.2-0.3$ dex) independently of the stellar mass. We further extended the IRX calibration to higher-redshift and lower-mass galaxies by stacking the FIR Spitzer and Herschel data.

This approach allowed us to derive the individual SFRs of one million galaxies in the HSC-CLAUDS sample down to $M_\star \sim 10^9 M_\odot$. The main results are as follows:

- Our estimated SFR nicely reproduces the evolution of the MS (SFR versus $M_\star$) up to $z \sim 2$ and the behavior of the attenuation (or (IRX)) with stellar mass reported in the literature.

- We reconstruct the SFR functions over a large range of SFRs ($10^{-2}-10^3 M_\odot \text{yr}^{-1}$) and redshifts (up to $z = 2$), providing a constraint on the faint-end slope for the first time. The SFRFs are well fitted by a Schechter function after correcting for the Eddington bias. The high-SFR tails are in good agreement with previous FIR SFRFs and show a strong evolution of the Schechter parameter, $SFR_\star$, with redshift. On the other hand, the slope of the SFRFs shows almost no evolution up to $z \sim 1.5-2$ with $\alpha = -1.3 \pm 0.1$.

- By integrating the SFRFs, we derive the cosmic SFRD from the current time up to $z = 2$. Despite a relatively good agreement with Madau & Dickinson (2014), we find that our SFRD is a better fit with a plateau between $1 \le z \le 2$ and a steeper drop below $z \sim 1$ than a gradual decline since $z \sim 2$. This has already been suggested in other studies (Katsianis et al. 2021; Gruppioni et al. 2020).

- The contributions to the total SFRD of galaxy populations with different SFR regimes vary rapidly with redshift:

Galaxies with moderate SFRs, $3 \le \text{SFR} \le 30 M_\odot \text{yr}^{-1}$, dominate the SFRD over most of cosmic time, $0.5 \le z \le 2$, while low-SFR galaxies, $0.3 \le \text{SFR} \le 3 M_\odot \text{yr}^{-1}$, dominate at lower $z$. The contribution of high-SFR galaxies, $\text{SFR} \ge 30 M_\odot \text{yr}^{-1}$, sharply increases since $z \sim 1$ to contribute equally at $z \sim 2$ with the intermediate SFR population.

- The contributions to the total SFRD of galaxy populations in different stellar mass regimes vary in a similar way with redshift. Galaxies with moderate stellar masses, $M_\star = 10^{9.5-10.5} M_\odot$, dominate at all redshifts. The decline below $z \sim 1$ affects galaxies of all stellar masses, with a faster pace observed for the highest stellar masses. While feedback and outflows are well-observed phenomena and can be efficient at quenching the SFR in galaxies, there may be a more global origin for the decline at all masses, such as the decline of the cosmological accretion rate as a consequence of an expanding Universe that becomes dominant after cosmic noon.

- We compared the observed SFR functions with four hydrodynamical simulations. Significant differences in the SFR functions are observed between the simulations, and none of them can reproduce the observations at all redshifts. They currently struggle to form high-SFR systems at high redshifts, with only one simulation able to reproduce the evolution of the density of high-SFR galaxies up to $z = 2$. Large differences are also observed at intermediate and low SFRs. This reflects the fact that the SFR functions provide a powerful diagnostic in addition to the more commonly used integrated quantities, such as the stellar mass functions. The SFR functions that give an instantaneous view of the distribution of the in situ star formation at different epochs remain a challenge for the simulations despite the incorporation of diverse, physically motivated prescriptions for the star formation and feedback processes.

All these results emphasize the interest of multiwavelength imaging surveys. While wider and deeper data sets could allow our results to be refined, the use of machine learning algorithms could also improve the calibration of the total infrared luminosity determination. In addition, the use of longer wavelength mid-infrared data would allow the analysis to be extended to higher redshifts.

*Acknowledgements.* V.P. and S.A. wish to thank Véronique Buat for fruitful discussions. We are grateful to the referees for their constructive input. V.P. work was partly funded by NASA APRA grant 80NSSC20K0396. This work was supported by the Spin(e) ANR project (ANR-13-BS05-0005), the DEEP-DIP ANR project (ANR-19-CE31-0023), and by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. These data were obtained and processed as part of the CFHT Large Area U-band Deep Survey (CLAUDS), which is a collaboration between astronomers from Canada, France, and China described in Sawicki et al. (2019). CLAUDS is based on observations obtained with MegaPrime/Mega-Cam, a joint project of CFHT and CEA/DAPNIA, at the CFHT which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Université of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. CLAUDS uses data obtained in part through the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories, Chinese Academy of Sciences, and the Special Fund for Astronomy from the Ministry of Finance of China. CLAUDS uses data products from CALET and the Canadian Astronomy Data Centre (CADC) and was processed using resources from Compute Canada and Canadian Advanced Network For Astrophysical Research (CANFAR) and the CANDIDE cluster at IAP maintained by Stephane Rouberol. The french members of the COSMOS team acknowledge the funding from the Centre National d'Études Spatiales (CNES). The Hyper Suprime-Cam (HSC) Collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics
and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. The Cosmic Dawn Center is funded by the Danish National Research Foundation under Grant No. 140. We acknowledge the Virgo Consortium for making their simulation data available. The EAGLE simulations were performed using the DiRAC-2 facility at Durham, managed by the ICC, and the PRACE facility Curie based in France at TGCC, CEA, Bruyères-le-Châtel. We have benefited from the publicly available programming language Python, including the numpy (Harris et al. 2020), matplotlib (Hunter 2007) packages, and the Topcat analysis tool (Taylor 2005).

## References

Aihara, H., Alsayyad, Y., Ando, M., et al. 2019, PASJ, 71, 1

Anglés-Alcázar, D., Davé, R., Faucher-Giguère, C.-A., Özel, F., & Hopkins, P. F. 2017, MNRAS, 464, 2840

Appleby, S., Davé, R., Kraljic, K., Anglés-Alcázar, D., & Narayanan, D. 2020, MNRAS, 494, 6053

Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355

Arnouts, S., Walcher, C. J., Le Fèvre, O., et al. 2007, A&A, 476, 137

Arnouts, S., Floc'h, E., Chevallard, J., et al. 2013, A&A, 558, A67

Barro, G., Pérez-González, P. G., Cava, A., et al. 2019, ApJS, 243, 22

Bavouzet, N., Dole, H., Le Floc'h, E., et al. 2008, A&A, 479, 83

Bell, E. F., Papovich, C., Wolf, C., et al. 2005, ApJ, 625, 23

Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bethermin, M., Dole, H., Cousin, M., & Bavouzet, N. 2010, A&A, 516, A43

Bethermin, M., Le Floc'h, E., Ilbert, O., et al. 2012, A&A, 542, A58

Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846

Bogdanoska, J., & Burgarella, D. 2020, MNRAS, 496, 5341

Bondi, H. 1952, MNRAS, 112, 195

Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53

Boquien, M., Calzetti, R., Kennicutt, R., et al. 2009, ApJ, 706, 553

Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5

Bothwell, M. S., Kenicutt, R. C., Johnson, B. D., et al. 2011, MNRAS, 415, 1815

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

Bourne, N., Dunlop, J. S., Merlin, E., et al. 2017, MNRAS, 467, 1360

Bouwens, R. J., Illingworth, G. D., Franx, M., et al. 2009, ApJ, 705, 936

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

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

Buat, V., Boselli, A., Gavazzi, G., & Bonfanti, C. 2002, A&A, 383, 801

Cacciato, M., Dekel, A., & Genel, S. 2012, MNRAS, 421, 818

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

Cattaneo, A., Faber, S. M., Binney, J., et al. 2009, Nature, 460, 213

Chabrier, G. 2003, PASP, 115, 763

Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772

Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718

Chary, R., & Elbaz, D. 2001, ApJ, 556, 562

Chevallard, J., Charlot, S., Wandelt, B., & Wild, V. 2013, MNRAS, 432, 2061

Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, 754, 125

Ciesla, L., Boselli, A., Elbaz, D., et al. 2016, A&A, 585, A43

Conseil, S., Vibert, D., Arnouts, S., et al. 2011, ASP Conf. Proc., 442, 107C

Cortese, L., Boselli, A., Franzetti, P., et al. 2008, MNRAS, 386, 1157

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

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

Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

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

Dale, D. A., & Helou, G. 2002, ApJ, 576, 159

Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140

Davé, R., Finlator, K., & Oppenheimer, B. D. 2011a, MNRAS, 416, 1354

Davé, R., Oppenheimer, B. D., & Finlator, K. 2011b, MNRAS, 415, 11

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

Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Davé, R., Crain, R. A., Stevens, A. R. H., et al. 2020, MNRAS, 497, 146

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451

Desprez, G., Picouet, V., Moutard, T., et al. 2023, A&A, 670, A82

Dole, H., Lagache, G., Puget, J. L., et al. 2006, A&A, 451, 417

Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 489, 3036

Draine, B. T. 2009, ASP Conf. Ser., 414, 453

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662

Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453

Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948

Eddington, A. S. 1913, MNRAS, 73, 359

Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Fabian, A. C. 2012, ARA&A, 50, 455

Faucher-Giguère, C.-A., Lidz, A., Zaldarriaga, M., & Hernquist, L. 2009, ApJ, 703, 1416

Felten, J. E. 1976, ApJ, 207, 700

Fensch, J., & Bournaud, F. 2021, MNRAS, 505, 3579

Fontanot, F., Hirschmann, M., & De Lucia, G. 2017, ApJ, 842, L14

Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661

Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486

Garn, T., & Best, P. N. 2010, MNRAS, 409, 421

Gehrz, R. 1989, in Interstellar Dust, eds. L. J. Allamandola, & A. G. G. M. Tielens, 135, 445

Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175

Goto, T., Takagi, T., Matsuhara, H., et al. 2010, A&A, 514, A6

Goto, T., Arnouts, S., Inami, H., et al. 2011, MNRAS, 410, 573

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

Gruppioni, C., Calura, F., Pozzi, F., et al. 2015, MNRAS, 451, 3419

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

Haardt, F., & Madau, P. 1996, ApJ, 461, 20

Haardt, F., & Madau, P. 2001, in Clusters of Galaxies and the High Redshift Universe Observed in X-rays, eds. D. M. Neumann, & J. T. V. Tran, 64

Haardt, F., & Madau, P. 2012, ApJ, 746, 125

Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124

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

Heckman, T. M., Alexandroff, R. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 809, 147

Heinis, S., Buat, V., Béethermin, M., et al. 2013, MNRAS, 429, 1113

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

Hopkins, P. F. 2015, MNRAS, 450, 53

Hopkins, A., & Beacom, F. 2006, ApJ, 651, 142

Hopkins, P. F., & Quataert, E. 2011, MNRAS, 415, 1027

Hopkins, A. M., Connolly, A. J., Haarsma, D. B., & Cram, L. E. 2001, AJ, 122, 288

Hopkins, P. F., Kereš, D., Óñorbe, J., et al. 2014, MNRAS, 445, 581

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Ilbert, O., Tresse, L., Arnouts, S., et al. 2004, MNRAS, 351, 541

Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841

Ilbert, O., Capak, P., Salvato, M., et al. 2009, ApJ, 690, 1236

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

Ilbert, O., Arnouts, S., Le Floc'h, E., et al. 2015, A&A, 579, A2

Inoue, A. K. 2011, Earth Planets Space, 63, 1027

Jarvis, M. J., Bonfield, D. G., Bruce, V. A., et al. 2013, MNRAS, 428, 1281

Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56

Juneau, S., Glazebrook, K., Crampton, D., et al. 2005, ApJ, 619, L135

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

Kassin, S. 2010, MOIRCS Multi-PA Kinematics of Galaxies at z 2+: A Searchfor Disks, Subaru Proposal, 36

Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, ApJ, 758, 106

Katsianis, A., Blanc, G., Lagos, C. P., et al. 2017a, MNRAS, 472, 919

Katsianis, A., Tescari, E., Blanc, G., & Sargent, M. 2017b, MNRAS, 464, 4977

Katsianis, A., Yang, X., & Zheng, X. 2021, ApJ, 919, 88

Katz, N., Hernquist, L., & Weinberg, D. H. 1992, ApJ, 399, L109

Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739

Kennicutt, R. C., Jr 1998, ARA&A, 36, 189

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2

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

Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36

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

Le Floc'h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169

Le Floc'h, E., Aussel, H., Ilbert, O., et al. 2009, ApJ, 703, 222

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

Li, C., & White, S. D. M. 2009, MNRAS, 398, 2177

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

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

Lovell, C. C., Geach, J. E., Davé, R., Narayanan, D., & Li, Q. 2021, MNRAS, 502, 772

Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90

Ly, C., Lee, J. C., Dale, D. A., et al. 2011, ApJ, 726, 109

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., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86

Malefahlo, E. D., Jarvis, M. J., Santos, M. G., et al. 2022, MNRAS, 509, 4291

Mancuso, C., Lapi, A., Cai, Z.-Y., et al. 2015, ApJ, 810, 72

Marchetti, L., Vaccari, M., Franceschini, A., et al. 2016, MNRAS, 456, 1999

Martin, C., Seibert, M., Buat, V., et al. 2005, ApJ, 619, L59

Martin, C. L. 2005, ApJ, 621, 227

Mashian, N., Oesch, P. A., & Loeb, A. 2016, MNRAS, 455, 2101

McAlpine, S., Helly, J. C., Schaller, M., et al. 2016, Astron. Comput., 15, 72

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

McCracken, H. J., Werbaeke, L. V., Morrison, G. E., et al. 2018, ApJS, 10, 1

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

Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H. W. 2011, ApJ, 731, 113

Moutard, T., Arnouts, S., Ilbert, O., et al. 2016a, A&A, 590, A103

Moutard, T., Arnouts, S., Ilbert, O., et al. 2016b, A&A, 590, A102

Moutard, T., Sawicki, M., Arnouts, S., et al. 2020, MNRAS, 494, 1894

Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2

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

Oke, J. B. 1974, ApJS, 27, 21

Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

Pakmor, R., & Springel, V. 2013, MNRAS, 432, 176

Pannella, M., Carilli, C. L., Daddi, E., et al. 2009, ApJ, 698, L116

Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363

Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Planck Collaboration XVI. 2014, A&A, 571, A16

Planck Collaboration XIII. 2016, A&A, 594, A13

Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13

Rahmati, A., Pawlik, A. H., Raičević, M., & Schaye, J. 2013, MNRAS, 430, 2427

Reddy, N. A., Steidel, C. C., Pettini, M., et al. 2008, ApJS, 175, 48

Reddy, N. A., Erb, D. K., Pettini, M., Steidel, C. C., & Shapley, A. E. 2010, ApJ, 712, 1070

Rodighiero, G., Cimatti, A., Gruppioni, C., et al. 2010, A&A, 518, L25

Rodríguez Montero, F., Davé, R., Wild, V., Anglés-Alcázar, D., & Narayanan, D. 2019, MNRAS, 490, 2139

Rosas-Guevara, Y. M., Bower, R. G., Schaye, J., et al. 2015, MNRAS, 454, 1038

Salim, S., Rich, R. M., Charlot, S., et al. 2007, ApJS, 173, 267

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

Sawicki, M., Arnouts, S., Huang, J., et al. 2019, TBC, 17, 1

Schaye, J. 2004, ApJ, 643, 59

Schaye, J., & Dalla Vecchia, C. 2008, MNRAS, 383, 1210

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Schechter, P. 1976, ApJ, 203, 297

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

Schmidt, M. 1959, ApJ, 129, 243

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

Seibert, M., Martin, D. C., Heckman, T. M., et al. 2005, ApJ, 619, L55

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

Shivaei, I., Darvish, B., Sattari, Z., et al. 2020, ApJ, 903, L28

Shuntov, M., McCracken, H. J., Gavazzi, R., et al. 2022, A&A, 664, A61

Silk, J., & Mamon, G. A. 2012, Res. Astron. Astrophys., 12, 917

Smit, R., Bouwens, R. J., Franx, M., et al. 2012, ApJ, 756, 14

Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217

Snedden, A., Coughlin, J., Phillips, L. A., Mathews, G., & Suh, I.-S. 2016, MNRAS, 455, 2804

Sobral, D., Smail, I., Best, P. N., et al. 2013, MNRAS, 428, 1128

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Springel, V. 2005, MNRAS, 364, 1105

Springel, V. 2010, MNRAS, 401, 791

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289

Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776

Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253

Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

Taylor, M. B. 2005, ASP Conf. Ser., 347, 29

Teyssier, R. 2002, A&A, 385, 337

Thomas, N., Davé, R., Anglés-Alcázar, D., & Jarvis, M. 2019, MNRAS, 487, 5764

Tuffs, R. J., Popescu, C. C., Völk, H. J., Kylafis, N. D., & Dopita, M. A. 2004, A&A, 419, 821

Van de Voort, F., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 2782

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Wang, L., Norberg, P., Bethermin, M., et al. 2016, A&A, 592, L5

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

Weaver, J. R., Kauffmann, O., Shuntov, M., et al. 2021, Am. Astron. Soc. Meeting Abstracts, 53, 21506

Weaver, J. R., Davidzon, I., Toft, S., et al. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202245581

Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291

Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104

Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99

Witt, A. N., & Gordon, K. D. 2000, ApJ, 528, 799

Wuyts, S., Labbé, I., Schreiber, N. M. F., et al. 2008, ApJ, 689, 653

Zamojski, M. A., Schiminovich, D., Rich, R. M., et al. 2007, ApJS, 172, 468
# Appendix A: IRX constraints from stacking analysis

Surveys conducted in the thermal infrared regime are gener-
ally not sensitive enough for detecting low-mass sources indi-
vidually. To explore the relationship between NrK and IRX
at low-to-intermediate stellar masses, we constrained galaxy
total infrared properties using standard stacking procedures. The
stacking technique consists in co-adding the signal arising from
a number of sources, and it is commonly used to characterize the
average emission associated with a given galaxy sample. Con-
straints can be inferred down to flux levels much fainter than the
typical $3\sigma$ sensitivity limit of the initial data, depending on the
number of co-added sources (e.g., Dole et al. 2006; Karim et al.
2011).

Our stacking was performed in the subsample lying in the COSMOS field, which benefits from exquisite MIR and FIR coverage from the *Spitzer* and *Herschel* satellites. Galaxies were stacked in various bins of redshift, mass, and *NrK* colors, using the IDL routines of the library released by Bethermin et al. (2010). The size of the redshift bin was fixed to $\Delta z = 0.25$ up to $z = 1.3$ and to $\Delta z = 0.3$ above. We adopted a color bin size of $\Delta NrK = 0.5$ and a bin size of 0.5 dex in logarithmic scale for the stellar mass. Typically, the *NrK* galaxy colors and the stellar masses span a range of to 4 magnitudes at the redshifts explored in this work. The adopted bin sizes thus ensured a statistically large enough number of sources to be stacked in each bin, while mitigating the effects of evolution within the bin. An illustration of the stacked signals that we obtained is shown in Fig. A.1 for two different redshift bins. The procedures we carried out are further described below.

## A.1. Stacking at 24 µm

Given the high space density of 24 µm sources identified in the COSMOS field (6.5 arcmin⁻² down to $3\sigma$), stacks were produced using the 24 µm *residual* map obtained by subtracting with PSF fitting each individual source found in the field (?; see their Fig. 1). For each bin of redshift, mass, and color, a stack was first created by mean-averaging 60''×60''sub-images centered at the sky position of each object in the bin (see Fig. A.1). The averaged stacked signal was measured using aperture photometry, with the same aperture and aperture correction as used by Le Floc'h et al. (2009) to perform the 24 µm PSF fitting and photometry in COSMOS. To account for the emission arising from sources individually detected at 24 µm and thus missed in the stack from the residual map, we then cross-correlated the galaxy subsamples associated with each bin with the list of sources initially subtracted from the 24 µm mosaic.

This correlation was performed with a matching radius of 2'',
as also assumed by Le Floc'h et al. (2009) in their identification
of the COSMOS 24 µm source optical counterparts. The fluxes
of galaxies individually detected at 24 µm were then weighted by
the inverse of the total number of sources considered in the sub-
sample, and they were finally co-added to the averaged stacked
signal to obtain the average 24 µm flux of the population selected
in the bin.

The uncertainties associated with the stacked signal were
estimated using bootstrap techniques as also described by
Bethermin et al. (2010). For each subsample, we stacked as
many sub-images as the total number of sources to be stacked
in the bin, but we randomly allowed some of these sub-images
to appear several times in the stack. This process was repeated
200 times for each bin of mass, redshift, and color, and the uncer-

tainty was estimated from the dispersion of the averaged sig-
nal measured in each of the 200 stacked images. We co-added in
quadrature this uncertainty with the $1\sigma$ flux uncertainties of the
other sources already detected at 24 µm, which gave us the final
uncertainty associated with the mean 24 µm flux characterizing
the galaxy subsample of the bin.

To gain further confidence in the results obtained with this
first approach, we inferred another estimate of the 24 µm emis-
sion associated with each bin of redshift, mass, and color, using
a method based on median stacking and applied to the initial
24 µm image of COSMOS. Here, each pixel in the final stack
represents the median value of the distribution arising from the
same pixel taken from the sub-images of the sample. The sig-
nal in the final stack thus corresponds to the median property of
the stacked population. Systematic offsets may be expected with
the mean estimate discussed above, depending on the underly-
ing distribution of properties associated with the source popula-
tion in the bin. This method has yet the advantage of making the
contribution of neighboring contaminants almost negligible. We
measured the median stacked fluxes and their associated uncer-
tainties with aperture photometry and bootstrapping, similar to
the procedure already used for the mean stacking described ear-
lier. Comparison between the two methods revealed a remark-
able agreement down to faint fluxes (i.e., ~10 µJy), with a small
systematic offset mostly noticeable at low redshift.

This offset between the median and the mean flux densities
can be explained by the generally skewed distribution of galaxy
luminosities in a given bin of redshift and stellar mass. The same
effect can actually be seen by comparing in each bin the average
and the median fluxes of 24 µm sources individually reported
in the COSMOS 24 µm catalog. In the following analysis, we
decided to use the flux estimates obtained with the mean stack-
ing, mostly to remain consistent with the mean values of the IRX
considered for the population of sources individually detected in
the infrared and discussed in Sect. 4.1.

## A.2. Stacking at FIR wavelengths

While the 24 µm emission of galaxies allows us to probe their
component of hot dust heated by young massive stars and/or
radiations produced by BH accretion activity, the FIR wave-
length regime is key for constraining the peak of their SED and
estimating their total infrared luminosity with minimal uncer-
tainty. We stacked the COSMOS imaging performed at 250,
350, and 500 µm with the SPIRE instrument on board the *Her-
schel* satellite (Oliver et al. 2012), to determine the average FIR
emission of sources for each bin of redshift, mass, and color as
achieved at 24 µm. We followed the same approach except that
the mean stacking was performed directly on the initial SPIRE
images instead of the residual mosaics. This choice was moti-
vated by the larger PSF characterizing the FIR *Herschel* data, as
well as the difficulty to associate individual FIR detection with
their counterparts at shorter wavelengths, hence preventing in
each bin of stacked sources a reliable control of the contribution
of galaxies already subtracted from the residual map. Photome-
try on the stacks was estimated in different ways (aperture pho-
tometry, PSF deconvolution, and Gaussian fitting), and we also
compared median and mean stacking with each other. All meth-
ods led to consistent results. Flux uncertainties were finally esti-
mated with bootstrapping techniques, similar to the procedure
employed at 24 µm.

Given the size of the beam characterizing the *Herschel*
mosaics, the stacked emission at FIR wavelength is often
affected by galaxy clustering, which results in a systematic
<figure><img src="image_19.png" /><figcaption>**Fig. A.1.** Results from MIR and FIR stacking in bins of *NrK* and stellar mass in two redshift intervals. Each stamp includes the mean fluxes for the 24 µm data from *Spitzer* and the 250, 350, and 500 µm data from *Herschel*, as indicated in the inset</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. A.2.** Averaged MIR and FIR SED fitting of galaxy subsamples at 0.3 < z < 0.55 (left panel) and 1.3 < z < 1.6 (right panel), based on fluxes derived from the mean stacking technique described in Appendix A. For each redshift range, SED fits are illustrated for two different bins of stellar mass and two bins of *NRK* color. The best SED fits are shown with the solid blue curves, while the red curves represent SEDs extrapolated from the 24 µm flux. Downward filled triangles depict upper limits for the stacked fluxes at 500 µm.</figcaption></figure>

overestimate of the associated signal due to contamination by
neighboring sources. For each SPIRE band, the final flux was
therefore corrected for this bias, which was inferred by decom-
posing the stacked signal into the contribution of the PSF pro-
duced by the stacked sources, and that of a wider underlying
component characterized by a Gaussian profile arising from
the contaminating neighbors (Bethermin et al. 2012). As already
noticed in the literature, we found that the bias can be neglected
at the highest stellar masses, but typically reaches 30 to 40% for
galaxies at M\* ~10\* M⊙.

A.3. Total infrared luminosities and mean IRX estimates

For each bin of redshift, stellar mass, and *NrK* color, a fit-
ting of SEDs at MIR and FIR wavelengths was performed
with LePhare so as to estimate the average Infrared luminos-
ity of galaxies in the bin. This step was carried out following
the same procedure as the one adopted earlier for the individ-
ual detections (see Sect. 3.2). We fixed the 24 µm flux to the
one obtained by combining the stacked signal from the resid-
ual 24 µm map and the contribution of each individual detection
(see Sect. A.1). At FIR wavelengths, we used the peak flux of the mean stacked data, corrected from the clustering bias (see Sect. A.2). In several bins, especially at low stellar mass, we note that the stacked signal is only seen at 24 µm (see Fig. A.1). For these cases, the total infrared luminosity was inferred by extrapolating the average 24 µm flux with the set of templates from Dale & Helou (2002), which adopts the locally observed dust temperature-luminosity relationship. Such extrapolations usually provide luminosity estimates consistent with the more accurate constraints obtained in the FIR where the peak of galaxy infrared SEDs is located, except for galaxies at very high luminosities in the starburst regime. We also verified the validity of this approach using our own galaxy sample, considering the bins with detections in both the MIR and FIR stacked data. For these cases, an overall agreement was found between the luminosity extrapolated from the 24 µm flux and the one obtained with the combination of the *Spitzer* and *Herschel* fluxes, except at high luminosity where the $L_{IR}$ derived from the 24 µm tend to be overestimated (Sect. 3.2).

The average IRX for each bin of redshift, mass, and color was finally derived as the ratio between the average infrared luminosity estimated above and the mean of the UV luminosities measured individually for all the galaxies in the bin. The final uncertainty on the IRX was obtained by combining the uncertainties in quadrature.

# Appendix B: Description of the simulations

## B.1. SIMBA

SIMBA<fnref n="6" /> (Davé et al. 2019) was run with a modified version of the gravity and hydrodynamics solver GIZMO (Hopkins 2015), relying on the GADGET-3 gravity solver (Springel 2005). The SIMBA run used in this work follows the evolution of 1024³ DM and gas particles within a comoving volume of (100 $h^{-1}$ Mpc)³. The simulation assumes a standard ΛCDM cosmology compatible with Planck Collaboration XIII (2016) ($\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$, $\Omega_b = 0.048$, $H_0 = 68$ km $s^{-1}$ Mpc$^{-1}$, $\sigma_8 = 0.82$ and $n_s = 0.97$). The minimum gravitational softening length is 0.5 comoving $h^{-1}$ kpc, and the initial gas and DM particle mass is $1.82 \times 10^7 M_\odot$ and $9.6 \times 10^7 M_\odot$, respectively.

Photoionization heating and radiative cooling models (GRACKLE-3.1 library; Smith et al. 2017) account for metal cooling and the nonequilibrium evolution of primordial elements. A spatially uniform UV ionizing background model (Haardt & Madau 2012) is modified to account for self-shielding based on the Rahmati et al. (2013) prescription.

Star formation is H₂-based and is only allowed to occur in gas with the hydrogen density $n_H \ge 0.1$ cm⁻³. The SFR is computed from the H₂ density and the dynamical time following the Schmidt (1959) relation with 2 percent efficiency. The H₂ fraction follows the prescription of Krumholz & Gnedin (2011), based on the local column density and metallicity (Davé et al. 2016).

The stellar feedback is modeled using metal-enriched, two-phase galactic winds. In addition, the winds are decoupled, that is, hydrodynamics in the winds is turned off until they leave the ISM (Springel & Hernquist 2003). Therefore, they do not deposit energy in the ISM on their way out.

The chemical enrichment model tracks 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) from Type Ia and II SNe and asymptotic giant branch (AGB) stars. In addition, dust growth and destruction are tracked for each individual element on the fly.

The growth of BH particles follows a two-mode accretion model. Hot gas ($T > 10^5$ K) is accreted in a spherically symmetric way following Bondi (1952), while cold gas accretion follows a torque-limited sub-grid prescription capturing the response of gas inflows near the BH to angular momentum loss due to dynamical instabilities (Hopkins & Quataert 2011; Anglés-Alcázar et al. 2017). This combination of BH accretion modes determines the implementation of feedback from AGNs in the form of two-mode kinetic feedback, the so-called radiative mode and jet mode feedback. Ten percent of the material accreted into the central region is assumed to fall onto the BH. These gas particles are immediately ejected in a purely kinetic and bipolar way (i.e., with zero opening angle w.r.t. the angular momentum of the inner disk) in two modes. At high accretion rates (above 0.2 times Eddington rate) and mass above $10^{7.5} M_\odot$, BHs eject material in ~1000 km s⁻¹ winds without changing its temperature (radiative mode). As the BH accretion rate drops below 0.2 of the Eddington rate, jet feedback mode starts to turn on and is fully achieved below 0.02. Gas is ejected with a velocity increment proportional to the logarithm of the inverse of the accretion rate and capped at 7000 km s⁻¹. The temperature of the ejected particles is increased, consistently with observations (Fabian 2012).

In addition, SIMBA implements the X-ray radiation pressure feedback activated in galaxies with low cold gas content and when the jet mode is active. Gas with the hydrogen density $n_H < 0.13$ cm⁻³ is heated by increasing its temperature, while for gas above this density, one-half of the X-ray energy is applied kinetically in the form of radial outward kick, and the second half is added as heat. The main observational quantities used to calibrate the model include the observed $z = 0$ GSMF (Davé et al. 2019).

## B.2. HORIZON-AGN

HORIZON-AGN<fnref n="7" /> (Dubois et al. 2014) was run with the adaptive-mesh refinement code RAMSES (Teyssier 2002) within a comoving volume of (100 $h^{-1}$ Mpc)³ and assuming a standard ΛCDM cosmology compatible with WMAP (Komatsu et al. 2011) ($\Omega_m = 0.272$, $\Omega_\Lambda = 0.728$, $\Omega_b = 0.045$, $H_0 = 70.4$ km $s^{-1}$ Mpc$^{-1}$, $\sigma_8 = 0.81$ and $n_s = 0.967$). It contains 1024³ DM particles (i.e., a mass resolution of $M_{DM,\text{res}} = 8 \times 10^7 M_\odot$), and the initially coarse 1024³ grid (initial gas resolution is $M_{\text{gas,\text{res}}} = 1 \times 10^7 M_\odot$) is refined down to 1 physical kpc. The refinement is triggered when the number of particles becomes greater than 8 (or if the total baryonic mass reaches eight times the initial DM mass resolution in a cell). Heating of the gas from a uniform UV background is activated at $z_{\text{reion}} = 10$ following Haardt & Madau (1996) and gas is allowed to cool to $10^4$ K via H, He, and metals (Sutherland & Dopita 1993).

Star formation only proceeds in regions with hydrogen number density $n_H \ge 0.1$ cm⁻³ (the stellar mass resolution is $\approx 2 \times 10^6 M_\odot$), following a Schmidt relation with 2 percent efficiency, that is to say, 2 percent of gas above the threshold density is converted into stars local free-fall time. HORIZON-AGN implements sub-grid feedback from stellar winds and SN (both type Ia and II) with mass, energy, and metal release.

HORIZON-AGN follows galactic BH formation and growth. BHs can grow by gas accretion at a Bondi-Hoyle-Lyttleton rate capped at the Eddington accretion rate when they form a tight enough binary. The AGN feedback is implemented as a combination of two different modes, the so-called radio mode

<footnote n="6">⁶ http://simba.roe.ac.uk/</footnote>

<footnote n="7">⁷ https://www.horizon-simulation.org/</footnote>
operating when the accretion rate is below 1% of Eddington ratio and the quasar mode active otherwise. The quasar mode consists of an isotropic, spherically symmetric, injection of thermal energy. At low accretion rates, the radio mode deposits AGN feedback energy into a bipolar outflow with a jet velocity of $10^4$ km s$^{-1}$. The efficiency of the radio mode is larger than the quasar mode with efficiencies tuned to match the BH-galaxy scaling relations at $z=0$ (see Dubois et al. 2012, for details).

## B.3. EAGLE

EAGLE<fnref n="8" /> was run using a modified version of the N-body Tree-Particle-Mesh smoothed particle hydrodynamics code GADGET-3 (Springel 2005). The EAGLE run used in this work follows the evolution of $1504^3$ DM particles and an initially equal number of baryonic particles within a comoving volume of (100 Mpc)$^3$, yielding DM and baryonic mass of $9.7 \times 10^6 M_\odot$ and $1.81 \times 10^6 M_\odot$, respectively. The Plummer-equivalent gravitational softening length is 2.66 comoving kpc, limited to a maximum length of proper 0.7 kpc. The simulation assumes a standard $\Lambda$CDM cosmology compatible with Planck Collaboration XVI (2014) ($\Omega_m = 0.307$, $\Omega_\Lambda = 0.693$, $\Omega_b = 0.048$, $H_0 = 67.77$ km s$^{-1}$ Mpc$^{-1}$, $\sigma_8 = 0.8288$ and $n_s = 0.9611$). Radiative cooling and photo-heating are implemented element-by-element for 11 species (Wiersma et al. 2009) exposed to the cosmic microwave background and evolving UV/X-ray background radiation (Haardt & Madau 2001). Star formation is implemented within the gas modeled as a single-phase fluid with a polytropic pressure floor (Schaye & Dalla Vecchia 2008) with a metallicity-dependent density threshold (Schaye 2004), reproducing by construction the observed Kennicutt-Schmidt relation (Kennicutt 1998). The seeding low-mass galaxies with BHs and their growth via gas accretion and merging are based on the method introduced by Springel et al. (2005) and substantially modified by Booth & Schaye (2009) and Rosas-Guevara et al. (2015).

Both stellar and AGN feedback are implemented as stochastic heating of gas particles (Dalla Vecchia & Schaye 2012) with a temperature increase of $10^{7.5}$ K and $10^{8.5}$ K, respectively, chosen to minimize numerical radiative losses and to allow for self-regulation. Therefore, only a single mode of AGN feedback is implemented with a fixed efficiency and injected energy proportional to the gas accretion rate, following a scheme close to the so-called quasar-mode feedback.

Thermal energy is injected into the gas without turning off radiative cooling and without decoupling hydrodynamical forces. The main observational quantities used to calibrate the model include the observed $z=0$ GSMF, galaxy sizes, and stellar to BH mass (see Schaye et al. 2015; Crain et al. 2015, for more details).

## B.4. ILLUSTRISTNG

ILLUSTRISTNG<fnref n="9" /> (Pillepich et al. 2018; Nelson et al. 2019) is a suite of cosmological magnetohydrodynamic simulations run with the moving-mesh code AREPO (Springel 2010), assuming a $\Lambda$CDM cosmology compatible with Planck Collaboration XIII (2016). In this work, we use the TNG100 simulation, with the box length of $75 \ h^{-1}$Mpc (≈ 110 comoving Mpc), with $1820^3$ DM particles and $1820^3$ initial gas cells corresponding to a mass resolution of $7.5 \times 10^6 M_\odot$ for DM and $1.4 \times 10^6 M_\odot$ for baryons. The minimum gravitational softening length is 0.75 kpc for DM and stars, and 190 comoving pc for gas.

ILLUSTRISTNG employs a galaxy formation model built upon the original ILLUSTRIS simulation (Genel et al. 2014; Vogelsberger et al. 2014). It includes radiative gas cooling, both primordial and from metal lines, in the presence of a time-variable, spatially uniform, ionizing UV background instantaneously switched on at $z=6$, with corrections for self-shielding in the dense ISM (Katz et al. 1992; Faucher-Giguère et al. 2009). In addition, cooling is further modulated by the radiation field of nearby AGNs.

Star formation occurs within the gas with the hydrogen density $n_H \ge 0.1 \ cm^{-3}$ following the empirically defined Kennicutt-Schmidt relation. Pressurization of the multiphase ISM from unresolved SNe is modeled for star-forming gas with a two-phase effective equation of state (Springel & Hernquist 2003).

Stellar populations evolve and return mass and metals to their ambient ISM via Type Ia and II SNe and AGB stars following tabulated mass and metal yields. In practice, the model tracks the production and evolution of nine elements: H, He, C, N, O, Ne, Mg, Si, and Fe (Pillepich et al. 2018).

Feedback associated with star formation drives galactic scale outflows implemented with a kinetic wind scheme. Wind particles are hydrodynamically decoupled until they leave the dense ISM. Once hydrodynamically recoupled outside the local ISM, they deposit their mass, momentum, metals, and thermal energy content.

ILLUSTRISTNG follows the formation of massive BHs in sufficiently massive halos, which accrete gas from surrounding gas and inject feedback energy into their environment. The two modes of AGN feedback are implemented: at low accretion rates ILLUSTRISTNG employs a kinetic AGN feedback model producing BH-driven winds, while at high accretion rates thermal energy is injected into the gas surrounding the BH (Weinberger et al. 2017). ILLUSTRISTNG includes the magnetic fields that are followed with ideal magnetohydrodynamic and are dynamically coupled to the gas via the magnetic pressure (Pakmor & Springel 2013). The main quantities used to calibrate the ILLUSTRISTNG model include the global SFRD as a function of cosmic time, the GSMF at $z=0$, and the current stellar-to-halo mass relation (see Pillepich et al. 2018, for more details).

<footnote n="8">⁸ https://icc.dur.ac.uk/Eagle/index.php (Crain et al. 2015; Schaye et al. 2015; McAlpine et al. 2016)</footnote>

<footnote n="9">⁹ https://www.tng-project.org/</footnote>