Issue |
A&A
Volume 675, July 2023
|
|
---|---|---|
Article Number | A164 | |
Number of page(s) | 28 | |
Section | Extragalactic astronomy | |
DOI | https://doi.org/10.1051/0004-6361/202245756 | |
Published online | 14 July 2023 |
HSC-CLAUDS survey: The star formation rate functions since z ∼ 2 and comparison with hydrodynamical simulations
1
Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France
e-mail: vincent@picouet.fr
2
Department of Astronomy, Columbia University, 550 W. 120th Street, New York, NY 10027, USA
3
Laboratoire AIM, CEA/DSM-CNRS-Université Paris Diderot, IRFU/Service d’Astrophysique, Bât. 709, CEA-Saclay, 91191 Gif-sur-Yvette Cedex, France
4
Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de Strasbourg, 67000 Strasbourg, France
5
Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary’s University, 923 Robie Street, Halifax, Nova Scotia B3H 3C3, Canada
6
Institut d’Astrophysique de Paris, UMR 7095, CNRS, UPMC Univ. Paris VI, 98 bis boulevard Arago, 75014 Paris, France
7
Herzberg Astronomy and Astrophysics, National Research Council of Canada, 5071 West Saanich Rd., Victoria, BC V9E 2E7, Canada
8
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
9
University of the Western Cape, Bellville, Cape Town 7535, South Africa
10
South African Astronomical Observatories, Observatory, Cape Town 7925, South Africa
11
Cosmic Dawn Center (DAWN), Copenhagen, Denmark
12
Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark
13
Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA
Received:
21
December
2022
Accepted:
5
April
2023
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 = LIR/LUV) 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 σ ∼ 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 ⟨IRX⟩) with stellar mass. In addition to the known lack of evolution of this relation up to z = 2 for galaxies with M⋆ ≤ 1010.3 M⊙, we observe a plateau in ⟨IRX⟩ 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⋆ ≥ 2.109 M⊙ at z = 2). The SFR functions cover a wide range of SFRs (0.01 ≤ SFR ≤ 1000 M⊙ 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 log10(SFR⋆) = 5.8z + 0.76. The slope of the SFR functions, α, shows almost no evolution up to z = 1.5 − 2 with α = −1.3 ± 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 ≤ z ≤ 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
Note to the reader: Reference Sawicki et al. 2019 was incorrectly written. The reference was corrected on 22 August 2023.
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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 an increase from early time up to z ∼ 3 − 4 (Smit et al. 2012; Mashian et al. 2016), the SFRD reaches a maximum at cosmic noon, z ∼ 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 ∼ 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, tdep = Mmol/SFR ∼ 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 ≤ 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α, or far-infrared (FIR) estimators. The simulations are in better agreement with the high-end SFR functions from UV/Hα 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 ∼ 108 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 = LIR/LUV), a measurement of the UV attenuation (AUV), and the slope of the UV continuum (β 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-β 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 (LIR) 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α luminosity (LHα) is produced from the gas ionized by short-lived massive stars (t ∼ 107 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α 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 ∼ 2, all the SFRFs derived from UV luminosities show a shortage of high-SFR sources (SFR > 100 M⊙ yr−1), while FIR selection reveals sources up to SFR = 1000 M⊙ yr−1. In the low-SFR regime, the slope of the SFRFs from Hα and UV luminosities shows a large range of values, −1.4 ≤ α ≤ −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 ∼ 2, based on the contribution of galaxies with M⋆ ≥ 109 M⊙.
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, σ(IRX)∼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 LIR 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 ∼ 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 ∼ 27) and optical imaging (i ∼ 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 (∼5.5 deg2 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 ≤ z ≤ 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 ΛCDM cosmology with Ωm = 0.3, ΩΛ = 0.7 and the Hubble constant H0 = 70 km s−1 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 ∼ 2.
2.1. COSMOS and far-infrared catalogs
2.1.1. COSMOS2020 Catalog
COSMOS is a major extragalactic field with a large multiwavelength photometric coverage over a 2 deg2 field. It collects ground-based optical observations with intermediate and broadband filters, NIR photometry for a total of 31 filter passbands 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 SExtractor 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 SExtractor-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 Ks-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 mAB(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 deg2, 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.
![]() |
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. |
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 SExtractor software (Bertin & Arnouts 1996) is produced using a multiband χ2 image in dual mode. In the SExtractor catalog we also include the far-UV (FUV; λ ∼ 1500 Å) and near-UV (NUV; λ ∼ 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 available1. We restricted our analysis to the E-COSMOS (α ∼ 150°, δ ∼ 2.5°) and XMM-LSS (α ∼ 35.5°, δ ∼ −5°) fields, where the deep NIR photometry is available. Figure 2 shows the layouts of the CLAUDS (with u and u⋆ filters), HSC-SSP, and VIRCAM observations in these two regions. The 5σ 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 iAB ≤ 26.5 and KAB ≤ 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.
![]() |
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). |
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.
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 iAB ≤ 25, with a scatter σ ≤ 0.03 (σ being defined as σ = 1.48 × Median(|zp − zs|/(1 + zs)) and a small fraction of outliers η ≤ 7% (η being defined as the fraction of galaxies with Δz > 0.15 × (1 + z)).
At fainter magnitude, iAB ≥ 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 τ−1e−t/τ with τ varying between 0.1 and 30 Gyr, and two delayed star formation histories; two stellar 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/τ > 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, Ks 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 (NUVabs − rabs) versus (rabs − Kabs) (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.
![]() |
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. |
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 |zCOSMOS − zCLAUDS|/ . This cut will always be used when comparing the physical parameters. Our stellar masses (top panels) appear 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%.
![]() |
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. |
The stellar mass completeness (smallest mass at which most of the objects would still be observable) of our SFG sample (iAB < 26.5 and KAB < 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 Mmin 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 (Mlim) it would have at its redshift if its KAB apparent magnitude was equal to the limiting magnitude of the survey (Klim = 24.5): log10(Mlim) = log10(Mmed)+0.4(K − Klim). At 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 (LIR) is defined as the luminosity integrated from 8 to 1000 μm:
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 LIR 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 () and the infrared luminosity (LIR) measured with the PACS and SPIRE photometry. However, they noticed a bias at high luminosity (
), where the
systematically overestimates the LIR. 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 LIR estimate. We confirm the excellent agreement for the SFGs with low infrared luminosity, with a small scatter (σ ∼ 0.07 dex) in between the two estimates, while at high luminosity we confirm the overestimation of
compared to LIR. We therefore included the Herschel photometry in the analysis to improve the estimate at high infrared luminosity; it does not impact the LIR estimate for the bulk of the 24 μm sample.
![]() |
Fig. 5. Difference between the infrared luminosity, LIR, estimated from the 24 μm photometry alone and by adding the Herschel photometry up to z ∼ 2. |
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:
![]() |
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 ≤ 24.5; see text). The horizontal gray lines in the bottom panel indicate the SFR threshold that corresponds to the luminous (LIRG, LIR ≥ 1011 L⊙) and ultra-luminous (ULIRG, LIR ≥ 1012 L⊙) infrared galaxies. |
where ℒNUV is the monochromatic NUV luminosity: ℒNUV/L⊙ = νLν(2300 Å). As already shown by Le Floc’h et al. (2005), below z = 0.5, the population is composed of moderately SFGs with SFR ≤ 10 M⊙ yr−1.
The fraction of luminous infrared galaxies (LIRGs) gradually increases from z = 0.5 to z = 1 and dominates at z ≥ 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⋆ ≥ 5 × 109 M⊙. 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 ; where
and
are the Vmax weighted comoving volume densities of FIR (f24 μm ≥ 50 μJy) and Ks-selected (K ≤ 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⋆ ∼ 1 − 2 × 1010 M⊙ at z ≥ 0.5. The lower completeness near z ∼ 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 = LIR/ℒNUV, allows for assessing the obscured star-formation contribution and offers an interesting alternative to SED fitting derived SFR, which depends strongly on the adopted attenuation laws2.
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 ∼ 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:
Since the dispersion σ[IRX(ϕ)] reaches a minimum when the vector NRK(ϕ) is perpendicular to the stripes, we minimize σ[IRX(ϕ)] to find the best angle. We obtain an angle of ϕ = 21°, consistent with A13 (ϕ = 18°), 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 (σ ≤ 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 different regions as their NrK distributions shift to higher values at higher stellar mass.
![]() |
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. |
We first adopted the same parametrization for the ⟨IRX⟩ versus NrK relation as A13, with two separated quantities ⟨IRX⟩=f(z)+α ⋅ NrK, where f(z) is a third-order polynomial function describing the redshift evolution, and α 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) = a0 + a1 ⋅ z + a2 ⋅ z2 + a3 ⋅ z3 we derived a0 = −0.73 ± 0.02; a1 = 2.02 ± 0.06; a2 = −1.33 ± 0.06; a3 = 0.33 ± 0.02; for the NrK term, α = 0.63 ± 0.003.
The result of this calibration is added as dashed black lines in Fig. 7, and the predicted LIR () that comes from it are compared to the reference LIR in Fig. 8. The predicted and reference LIR are in good agreement with almost no bias with redshift and a global dispersion σ ∼ 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.
![]() |
Fig. 8. Comparison of the infrared luminosity, LIR, estimated from the COSMOS2020 FIR sample with that estimated with the IRX–NrK calibration ( |
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.
![]() |
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. |
We also observe that at high redshift and stellar mass bins (bottom-right panels), the LIR 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 LIR 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
In each (M⋆, 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 bottom right panel combines the ⟨IRX⟩ versus stellar mass relationship for the four redshift bins.
![]() |
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. |
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⋆ ≤ 1010.3 M⊙), 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⋆ ≥ 1010.3 M⊙), a ⟨IRX⟩ plateau is observed despite the large scatter in the measurements. It appears at different stellar mass with redshift, near M⋆ ∼ 1010.7 M⊙ at z ∼ 2 and M⋆ ∼ 1010.3 M⊙ 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).
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 (
) and
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⋆ ≥ 1010 M⊙), 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 LIR and the β-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-βUV relations (85% of the ∼10 K matched sources). In all the mass bins the different SFR estimates are consistent with each other within their relative uncertainties. At low mass, M⋆ ≤ 1010 M⊙, our results are in overall good agreement with CANDELS SFRs, with a small scatter (σ < 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 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 masses3, 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⋆ > 1011 M⊙, 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.
![]() |
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. |
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 ∼ 0, estimate is in excellent agreement with GALEX (Salim et al. 2007) and HRS Ciesla et al. (2016) surveys.
![]() |
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. |
At higher redshift, the normalization of the MS agrees with previous studies up to z ∼ 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⋆ ∼ 1010 M⊙. 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 ≤ z ≤ 1. One noticeable feature is that the unobscured SFR distributions for the three high-mass bins saturate in SFR around a similar value (LNUV ∼ 109.7 L⊙). 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 LUV ∼ 1010 L⊙, while the LFIR continues to increase by a factor of ∼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.
![]() |
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. |
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 LUV is observed, so that fainter LFUV 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 LUV 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 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 Vmax estimator (Felten 1976; Ilbert et al. 2004) combining the i- and Ks-band selections. Each galaxy is then weighted by 1/Vmax(z = min(zi, zk)), 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 (σpoisson) and the cosmic variance due to large-scale density fluctuations (σcv):
The cosmic variance term, σcv, is estimated by using the Moster et al. (2011) cookbook. The DM cosmic variance, , is given by their Eq. (10) and depends on the redshift and the considered surface. To scale this term to the surface of our survey, we linearly extrapolated their estimate predicted for 2 deg2 to 5.5 deg2. The galaxy bias,
, 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 ≤ z ≤ 2.05 and derived with the IRX–NrK calibration without (black dots) and with (red dots) the stacking correction.
![]() |
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 (α = −1.3). The gray area corresponds at the bright end to the Eddington correction and at the faint end to the slope uncertainty (Δα = 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. |
The orange shaded area shows the Vmax weighted SFR functions derived for the COSMOS2020 FIR-selected catalog (f24 μm ≥ 50 μJy), with the total SFR derived according to Eq. (2).
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):
where α is the faint-end slope of the power-law regime and and 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 σ = 0.15 at z ≤ 1 and σ = 0.2 at z ≥ 1, assuming that the two SFR estimators used in the comparison have similar errors (i.e., by dividing the observed error by ). We then convolve the Schechter function (Eq. (8)) by the SFR uncertainty:
A least square fit procedure is performed between ϕ(SFR)conv and the Vmax 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⋆ 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 α = −1.3 ± 0.1. Since there is no evidence for the evolution of α 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 σ(α) = 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 (δα = ±0.1) at the faint end and the Eddington correction at the bright end.
Parameter values of the SFR function Schechter fit for both individual FIR calibration and the verified calibration.
Finally, we also fit the SFR functions by a double power law with the same fixed faint-end slope:
The double-power-law fits give very similar results to the Schechter functions down to Φ(SFR)∼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 (SFR⋆) exhibits a monotonic decline from ∼60 to ∼2 M⊙ yr−1 from z = 2 to z = 0 and a smooth, though a noisier, increase of the normalization Φ⋆ 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α emission lines, and we apply a correction factor for studies using a different initial mass function than Chabrier4.
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 (≤11 Mpc) to constrain the faint end of the local SFR function down to SFR < 0.01 M⊙ 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 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 ≤ z ≤ 2.7 after correction for incompleteness effects. Marchetti et al. (2016) used Herschel observations to infer the infrared luminosity down to LIR = 109 L⊙ 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α Survey corresponding to Hα 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α emitters at z = 0.40, 0.84, 1.47 and 2.23. The Hα luminosity functions were then corrected for incompleteness, [NII] contamination, and dust extinction assuming an average attenuation of A(Hα) = 1 mag.
Parsa et al. (2016) used the deep fields (HUDF, CANDELS, and UltraVista-COSMOS) to measure the UV luminosity functions at 1.5 ≤ z ≤ 2.5. To convert it into an SFR function, Katsianis et al. (2017b) adopt a luminosity-dependent evolution of the β-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 ≤ z ≤ 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 (α = −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 end5 but their normalizations appear consistent with us around SFR ∼ 1–3 M⊙ 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 α = −1.2. The SFRF from Hα by Sobral et al. (2013, dark green triangles) at z ∼ 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α) = 1 mag) applied (see below).
At intermediate redshifts, 0.5 ≤ z ≤ 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 (log10(M⋆) = 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 (α ∼ −1.3) with no sign of evolution up to z ∼ 1.5. The Hα SFRF from Ly et al. (2011, light green triangles) at z ∼ 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α luminosity function.
At high redshifts, 1.5 ≤ z ≤ 2, our SFRFs are still in good agreement with Gruppioni et al. (2013, 2015) at the bright end, while Hα SFRF slightly underestimate the high star-forming population. At z ∼ 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 β-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 Log10(M⋆/M⊙) = 9.5 at this redshift. We thus can miss the potential contribution of lower-mass galaxies in the SFRF around SFR∼10 M⊙ 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 (Δz ∼ 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 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 ≤ 2, Chapman et al. 2005). It could help us better match our SFR functions above SFR ≥ 100 M⊙ 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 (α ∼ −1.3), with no evolution at least up to z ∼ 1.6 and potentially z ∼ 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α- 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.
Summary of the different simulations used in this work: SIMBA, HORIZON-AGN, EAGLE, and TNG100.
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 log10(M⋆/M⊙) = 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⋆, init/100 Myr) of 0.18 M⊙ yr−1 for SIMBA and 0.02 M⊙ 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⋆ < 1011 M⊙; however, for more massive galaxies the aperture reduces the stellar masses (see, e.g., Schaye et al. 2015). Finally, 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: log10(sSFR/yr−1) = −10.8 + 0.3z. This is equivalent to removing galaxies with Δ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.
![]() |
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). |
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.
![]() |
Fig. 16. Redshift evolution of the Schechter parameters – α (top panel), SFR⋆ (middle panel), and Φ⋆ (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. |
Second, while all simulations reproduce the high-SFR side at low z (z ≤ 0.5), SIMBA is the only one able to reproduce its evolution with redshift up to z ∼ 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 ∼ 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 ≥ 1.
Finally, the comoving density of galaxies with low to intermediate SFRs varies significantly between simulations. At z ∼ 2, they all agree with a slope −1.5 ≤ α ≤ −1.3, consistent with the observations. At later times, while the slope in HORIZON-AGN shows almost no evolution, in EAGLE and TNG100 it gradually flattens, but yet within the range of the observations (except at z ∼ 0, where (α + 1) > 0 for TNG100). In contrast, the evolution is much more pronounced for SIMBA, with (α + 1) = 0 already at z = 1 and ends up with a positive slope ((α + 1)≫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 ≥ 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 ≥ 0.5, but reproduce reasonably well the low to intermediate 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 GSMFs are represented as a gray-shaded region. At z ∼ 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 particularly at z ∼ 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 1011 M⊙. At z ≳ 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 ≥ 1011 M⊙) 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 ≥ 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.
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⋆ = 1011 M⊙. 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 over-predict the low-mass end for the star-forming and the global population at z ≳ 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⋆ relation at high redshifts (z ∼ 2) by a factor of a few, while the agreement improves at low redshift (z ≲ 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⋆ ∼ 109.5 M⊙. 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⋆ ∼ 109 M⊙ and a factor of ∼3 for galaxies with M⋆ ≥ 1010 M⊙. 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⋆ ≤ 1010 M⊙.
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 ≳ 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 ≲ 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⋆ 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 ≳ 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⋆ relation, and the discrepancies at the bright and low end of SFR functions while keeping GSMFs in relatively good agreement with observations.
![]() |
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, α = −1.3 ± 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). |
Another interesting aspect is the scatter of the star-forming MS. This scatter is significantly broader for SIMBA (σ ∼ 0.4) compared to the other simulations (σ ∼ 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, ψ, is derived by integrating the SFR Schechter functions as
The redshift evolution of ψ(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 α = −1.3 (Δα = ±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 ∼ 2, the cosmic noon, up to the present day by a factor of 10. While Madau & Dickinson (2014) predicts a decline scaling as ψ(z)∝(1 + z)2.7 since z = 1.5 − 2, our observations suggest a slightly steeper rate with ψ(z)∝(1 + z)2.8 when using the same parametric form as their Eq. (15).
The redshift evolution of ψ(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 ≤ z ≤ 1 followed by a plateau between 1 ≤ z ≤ 2 rather than a peak at z ∼ 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 ≤ 0.5, the main contributors are galaxies with low star formation activity (0.3 ≤ SFR ≤ 3 M⊙ yr−1) while moderate SFGs (3 ≤ SFR ≤ 30 M⊙ yr−1) take over at higher redshift. The contribution of the high star-forming population (SFR ≥ 30 M⊙ yr−1) shows a sharp increase from 0.5 ≤ z ≤ 2 contributing equally at z ∼ 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 ≤ log(M⋆/M⊙)≤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, ψ(z)∝(1 + z)α, with α = 2.2 for the lowest-mass population and α = 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 ≤ Log(M⋆/M⊙)≤12 never exceeds 15% of the total SFRD. While we find that the two low-mass bins (Log(M⋆/M⊙)≤10 and 10 ≤ Log(M⋆/M⊙)≤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 ψ(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 ∼ 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, IRX = LIR/LUV, in the rest-frame color-color diagram, (NUV − r) versus (r − Ks), for galaxies observed in the UV and FIR in the COSMOS2020 catalog. We show that the use of a single vector, NrK, is effective in predicting the IRX with a small scatter (σ < 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⋆ ∼109 M⊙. The main results are as follows:
-
Our estimated SFR nicely reproduces the evolution of the MS (SFR versus M⋆) up to z∼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 − 103 M⊙ 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⋆, with redshift. On the other hand, the slope of the SFRFs shows almost no evolution up to z ∼ 1.5 − 2 with α = −1.3 ± 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 ≤ z ≤ 2 and a steeper drop below z∼1 than a gradual decline since z ∼ 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 ≤ SFR ≤ 30 M⊙ yr−1, dominate the SFRD over most of cosmic time, 0.5 ≤ z ≤ 2, while low-SFR galaxies, 0.3 ≤ SFR ≤ 3 M⊙ yr−1, dominate at lower z. The contribution of high-SFR galaxies, SFR ≥ 30 M⊙ yr−1, sharply increases since z ∼ 1 to contribute equally at z ∼ 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⋆ = 109.5 − 10.5 M⊙, dominate at all redshifts. The decline below z ∼ 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.
Most of them were already compiled in Katsianis et al. (2017b).
Acknowledgments
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 DEEPDIP 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/ MegaCam, 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 Science de l’Univers 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 [NASA ADS] [CrossRef] [Google Scholar]
- Anglés-Alcázar, D., Davé, R., Faucher-Giguère, C.-A., Özel, F., & Hopkins, P. F. 2017, MNRAS, 464, 2840 [Google Scholar]
- Appleby, S., Davé, R., Kraljic, K., Anglés-Alcázar, D., & Narayanan, D. 2020, MNRAS, 494, 6053 [NASA ADS] [CrossRef] [Google Scholar]
- Arnouts, S., Moscardini, L., Vanzella, E., et al. 2002, MNRAS, 329, 355 [Google Scholar]
- Arnouts, S., Walcher, C. J., Le Fèvre, O., et al. 2007, A&A, 476, 137 [CrossRef] [EDP Sciences] [Google Scholar]
- Arnouts, S., Floc’h, E., Chevallard, J., et al. 2013, A&A, 558, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Barro, G., Pérez-González, P. G., Cava, A., et al. 2019, ApJS, 243, 22 [NASA ADS] [CrossRef] [Google Scholar]
- Bavouzet, N., Dole, H., Le Floc’h, E., et al. 2008, A&A, 479, 83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Bell, E. F., Papovich, C., Wolf, C., et al. 2005, ApJ, 625, 23 [Google Scholar]
- Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697 [Google Scholar]
- Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Bethermin, M., Dole, H., Cousin, M., & Bavouzet, N. 2010, A&A, 516, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Bethermin, M., Le Floc’h, E., Ilbert, O., et al. 2012, A&A, 542, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846 [NASA ADS] [CrossRef] [Google Scholar]
- Bogdanoska, J., & Burgarella, D. 2020, MNRAS, 496, 5341 [NASA ADS] [CrossRef] [Google Scholar]
- Bondi, H. 1952, MNRAS, 112, 195 [Google Scholar]
- Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53 [Google Scholar]
- Boquien, M., Calzetti, D., Kennicutt, R., et al. 2009, ApJ, 706, 553 [NASA ADS] [CrossRef] [Google Scholar]
- Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5 [Google Scholar]
- Bothwell, M. S., Kenicutt, R. C., Johnson, B. D., et al. 2011, MNRAS, 415, 1815 [Google Scholar]
- Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001 [Google Scholar]
- Bourne, N., Dunlop, J. S., Merlin, E., et al. 2017, MNRAS, 467, 1360 [NASA ADS] [Google Scholar]
- Bouwens, R. J., Illingworth, G. D., Franx, M., et al. 2009, ApJ, 705, 936 [Google Scholar]
- Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503 [Google Scholar]
- Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000 [NASA ADS] [CrossRef] [Google Scholar]
- Buat, V., Boselli, A., Gavazzi, G., & Bonfanti, C. 2002, A&A, 383, 801 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Cacciato, M., Dekel, A., & Genel, S. 2012, MNRAS, 421, 818 [NASA ADS] [Google Scholar]
- Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682 [NASA ADS] [CrossRef] [Google Scholar]
- Cattaneo, A., Faber, S. M., Binney, J., et al. 2009, Nature, 460, 213 [NASA ADS] [CrossRef] [Google Scholar]
- Chabrier, G. 2003, PASP, 115, 763 [Google Scholar]
- Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772 [NASA ADS] [CrossRef] [Google Scholar]
- Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718 [Google Scholar]
- Chary, R., & Elbaz, D. 2001, ApJ, 556, 562 [Google Scholar]
- Chevallard, J., Charlot, S., Wandelt, B., & Wild, V. 2013, MNRAS, 432, 2061 [CrossRef] [Google Scholar]
- Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, 754, 125 [NASA ADS] [CrossRef] [Google Scholar]
- Ciesla, L., Boselli, A., Elbaz, D., et al. 2016, A&A, 585, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Conseil, S., Vibert, D., Arnouts, S., et al. 2011, ASP Conf. Proc., 442, 107C [NASA ADS] [Google Scholar]
- Cortese, L., Boselli, A., Franzetti, P., et al. 2008, MNRAS, 386, 1157 [NASA ADS] [CrossRef] [Google Scholar]
- Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839 [Google Scholar]
- Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937 [NASA ADS] [CrossRef] [Google Scholar]
- Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, A&A, 539, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686 [NASA ADS] [CrossRef] [Google Scholar]
- Dahlen, T., Mobasher, B., Dickinson, M., et al. 2007, ApJ, 654, 172 [Google Scholar]
- Dale, D. A., & Helou, G. 2002, ApJ, 576, 159 [Google Scholar]
- Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140 [NASA ADS] [CrossRef] [Google Scholar]
- Davé, R., Finlator, K., & Oppenheimer, B. D. 2011a, MNRAS, 416, 1354 [CrossRef] [Google Scholar]
- Davé, R., Oppenheimer, B. D., & Finlator, K. 2011b, MNRAS, 415, 11 [Google Scholar]
- Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98 [Google Scholar]
- Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265 [Google Scholar]
- Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 [Google Scholar]
- Davé, R., Crain, R. A., Stevens, A. R. H., et al. 2020, MNRAS, 497, 146 [Google Scholar]
- Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451 [Google Scholar]
- Desprez, G., Picouet, V., Moutard, T., et al. 2023, A&A, 670, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Dole, H., Lagache, G., Puget, J. L., et al. 2006, A&A, 451, 417 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 489, 3036 [NASA ADS] [CrossRef] [Google Scholar]
- Draine, B. T. 2009, ASP Conf. Ser., 414, 453 [Google Scholar]
- Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662 [Google Scholar]
- Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453 [Google Scholar]
- Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948 [Google Scholar]
- Eddington, A. S. 1913, MNRAS, 73, 359 [Google Scholar]
- Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Fabian, A. C. 2012, ARA&A, 50, 455 [Google Scholar]
- Faucher-Giguère, C.-A., Lidz, A., Zaldarriaga, M., & Hernquist, L. 2009, ApJ, 703, 1416 [Google Scholar]
- Felten, J. E. 1976, ApJ, 207, 700 [NASA ADS] [CrossRef] [Google Scholar]
- Fensch, J., & Bournaud, F. 2021, MNRAS, 505, 3579 [CrossRef] [Google Scholar]
- Fontanot, F., Hirschmann, M., & De Lucia, G. 2017, ApJ, 842, L14 [Google Scholar]
- Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661 [Google Scholar]
- Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486 [Google Scholar]
- Garn, T., & Best, P. N. 2010, MNRAS, 409, 421 [Google Scholar]
- Gehrz, R. 1989, in Interstellar Dust, eds. L. J. Allamandola, & A. G. G. M. Tielens, 135, 445 [NASA ADS] [CrossRef] [Google Scholar]
- Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175 [Google Scholar]
- Goto, T., Takagi, T., Matsuhara, H., et al. 2010, A&A, 514, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Goto, T., Arnouts, S., Inami, H., et al. 2011, MNRAS, 410, 573 [NASA ADS] [CrossRef] [Google Scholar]
- Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23 [Google Scholar]
- Gruppioni, C., Calura, F., Pozzi, F., et al. 2015, MNRAS, 451, 3419 [Google Scholar]
- Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Haardt, F., & Madau, P. 1996, ApJ, 461, 20 [Google Scholar]
- Haardt, F., & Madau, P. 2001, in Clusters of Galaxies and the High RedshiftUniverse Observed in X-rays, eds. D. M. Neumann, & J. T. V. Tran, 64 [Google Scholar]
- Haardt, F., & Madau, P. 2012, ApJ, 746, 125 [Google Scholar]
- Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124 [Google Scholar]
- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
- Heckman, T. M., Alexandroff, R. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 809, 147 [Google Scholar]
- Heinis, S., Buat, V., Béethermin, M., et al. 2013, MNRAS, 429, 1113 [CrossRef] [Google Scholar]
- Heinis, S., Buat, V., Béthermin, M., et al. 2014, MNRAS, 437, 1268 [Google Scholar]
- Hopkins, P. F. 2015, MNRAS, 450, 53 [Google Scholar]
- Hopkins, A., & Beacom, F. 2006, ApJ, 651, 142 [NASA ADS] [CrossRef] [Google Scholar]
- Hopkins, P. F., & Quataert, E. 2011, MNRAS, 415, 1027 [NASA ADS] [CrossRef] [Google Scholar]
- Hopkins, A. M., Connolly, A. J., Haarsma, D. B., & Cram, L. E. 2001, AJ, 122, 288 [NASA ADS] [CrossRef] [Google Scholar]
- Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581 [NASA ADS] [CrossRef] [Google Scholar]
- Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
- Ilbert, O., Tresse, L., Arnouts, S., et al. 2004, MNRAS, 351, 541 [NASA ADS] [CrossRef] [Google Scholar]
- Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ilbert, O., Capak, P., Salvato, M., et al. 2009, ApJ, 690, 1236 [NASA ADS] [CrossRef] [Google Scholar]
- Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ilbert, O., Arnouts, S., Le Floc’h, E., et al. 2015, A&A, 579, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Inoue, A. K. 2011, Earth Planets Space, 63, 1027 [NASA ADS] [CrossRef] [Google Scholar]
- Jarvis, M. J., Bonfield, D. G., Bruce, V. A., et al. 2013, MNRAS, 428, 1281 [Google Scholar]
- Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56 [Google Scholar]
- Juneau, S., Glazebrook, K., Crampton, D., et al. 2005, ApJ, 619, L135 [NASA ADS] [CrossRef] [Google Scholar]
- Karim, A., Schinnerer, E., Martínez-Sansigre, A., et al. 2011, ApJ, 730, 61 [Google Scholar]
- Kassin, S. 2010, MOIRCS Multi-PA Kinematics of Galaxies at z 2+: A Searchfor Disks, Subaru Proposal, 36 [Google Scholar]
- Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, ApJ, 758, 106 [NASA ADS] [CrossRef] [Google Scholar]
- Katsianis, A., Blanc, G., Lagos, C. P., et al. 2017a, MNRAS, 472, 919 [NASA ADS] [CrossRef] [Google Scholar]
- Katsianis, A., Tescari, E., Blanc, G., & Sargent, M. 2017b, MNRAS, 464, 4977 [NASA ADS] [CrossRef] [Google Scholar]
- Katsianis, A., Yang, X., & Zheng, X. 2021, ApJ, 919, 88 [NASA ADS] [CrossRef] [Google Scholar]
- Katz, N., Hernquist, L., & Weinberg, D. H. 1992, ApJ, 399, L109 [NASA ADS] [CrossRef] [Google Scholar]
- Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739 [NASA ADS] [Google Scholar]
- Kennicutt, R. C., Jr 1998, ARA&A, 36, 189 [NASA ADS] [CrossRef] [Google Scholar]
- Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2 [Google Scholar]
- Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [Google Scholar]
- Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36 [Google Scholar]
- Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24 [Google Scholar]
- Le Floc’h, E., Papovich, C., Dole, H., et al. 2005, ApJ, 632, 169 [Google Scholar]
- Le Floc’h, E., Aussel, H., Ilbert, O., et al. 2009, ApJ, 703, 222 [Google Scholar]
- Leslie, S., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58 [Google Scholar]
- Li, C., & White, S. D. M. 2009, MNRAS, 398, 2177 [NASA ADS] [CrossRef] [Google Scholar]
- Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119 [NASA ADS] [CrossRef] [Google Scholar]
- Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172 [Google Scholar]
- Lovell, C. C., Geach, J. E., Davé, R., Narayanan, D., & Li, Q. 2021, MNRAS, 502, 772 [NASA ADS] [CrossRef] [Google Scholar]
- Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ly, C., Lee, J. C., Dale, D. A., et al. 2011, ApJ, 726, 109 [Google Scholar]
- Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
- Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Malefahlo, E. D., Jarvis, M. J., Santos, M. G., et al. 2022, MNRAS, 509, 4291 [Google Scholar]
- Mancuso, C., Lapi, A., Cai, Z.-Y., et al. 2015, ApJ, 810, 72 [Google Scholar]
- Marchetti, L., Vaccari, M., Franceschini, A., et al. 2016, MNRAS, 456, 1999 [Google Scholar]
- Martin, C., Seibert, M., Buat, V., et al. 2005, ApJ, 619, L59 [NASA ADS] [CrossRef] [Google Scholar]
- Martin, C. L. 2005, ApJ, 621, 227 [NASA ADS] [CrossRef] [Google Scholar]
- Mashian, N., Oesch, P. A., & Loeb, A. 2016, MNRAS, 455, 2101 [CrossRef] [Google Scholar]
- McAlpine, S., Helly, J. C., Schaller, M., et al. 2016, Astron. Comput., 15, 72 [NASA ADS] [CrossRef] [Google Scholar]
- McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- McCracken, H. J., Werbaeke, L. V., Morrison, G. E., et al. 2018, ApJS, 10, 1 [Google Scholar]
- Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64 [Google Scholar]
- Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H. W. 2011, ApJ, 731, 113 [NASA ADS] [CrossRef] [Google Scholar]
- Moutard, T., Arnouts, S., Ilbert, O., et al. 2016a, A&A, 590, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Moutard, T., Arnouts, S., Ilbert, O., et al. 2016b, A&A, 590, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Moutard, T., Sawicki, M., Arnouts, S., et al. 2020, MNRAS, 494, 1894 [NASA ADS] [CrossRef] [Google Scholar]
- Nelson, D., Springel, V., Pillepich, A., et al. 2019, Comput. Astrophys. Cosmol., 6, 2 [Google Scholar]
- Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43 [CrossRef] [Google Scholar]
- Oke, J. B. 1974, ApJS, 27, 21 [Google Scholar]
- Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614 [NASA ADS] [CrossRef] [Google Scholar]
- Pakmor, R., & Springel, V. 2013, MNRAS, 432, 176 [NASA ADS] [CrossRef] [Google Scholar]
- Pannella, M., Carilli, C. L., Daddi, E., et al. 2009, ApJ, 698, L116 [Google Scholar]
- Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194 [Google Scholar]
- Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363 [CrossRef] [Google Scholar]
- Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493 [CrossRef] [Google Scholar]
- Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077 [Google Scholar]
- Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196 [Google Scholar]
- Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Rahmati, A., Pawlik, A. H., Raičević, M., & Schaye, J. 2013, MNRAS, 430, 2427 [NASA ADS] [CrossRef] [Google Scholar]
- Reddy, N. A., Steidel, C. C., Pettini, M., et al. 2008, ApJS, 175, 48 [CrossRef] [Google Scholar]
- Reddy, N. A., Erb, D. K., Pettini, M., Steidel, C. C., & Shapley, A. E. 2010, ApJ, 712, 1070 [NASA ADS] [CrossRef] [Google Scholar]
- Rodighiero, G., Cimatti, A., Gruppioni, C., et al. 2010, A&A, 518, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Rodríguez Montero, F., Davé, R., Wild, V., Anglés-Alcázar, D., & Narayanan, D. 2019, MNRAS, 490, 2139 [CrossRef] [Google Scholar]
- Rosas-Guevara, Y. M., Bower, R. G., Schaye, J., et al. 2015, MNRAS, 454, 1038 [NASA ADS] [CrossRef] [Google Scholar]
- Salim, S., Rich, R. M., Charlot, S., et al. 2007, ApJS, 173, 267 [NASA ADS] [CrossRef] [Google Scholar]
- Sanders, D. B., Mazzarella, J. M., Kim, D. C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607 [Google Scholar]
- Sawicki, M., Arnouts, S., Huang, J., et al. 2019, MNRAS, 489, 5202 [NASA ADS] [Google Scholar]
- Schaye, J. 2004, ApJ, 643, 59 [Google Scholar]
- Schaye, J., & Dalla Vecchia, C. 2008, MNRAS, 383, 1210 [Google Scholar]
- Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521 [Google Scholar]
- Schechter, P. 1976, ApJ, 203, 297 [Google Scholar]
- Schiminovich, D., Ilbert, O., Arnouts, S., et al. 2005, ApJ, 619, L47 [Google Scholar]
- Schmidt, M. 1959, ApJ, 129, 243 [NASA ADS] [CrossRef] [Google Scholar]
- Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Seibert, M., Martin, D. C., Heckman, T. M., et al. 2005, ApJ, 619, L55 [NASA ADS] [CrossRef] [Google Scholar]
- Shapley, A. E., Sanders, R. L., Salim, S., et al. 2022, ApJ, 926, 145 [NASA ADS] [CrossRef] [Google Scholar]
- Shivaei, I., Darvish, B., Sattari, Z., et al. 2020, ApJ, 903, L28 [NASA ADS] [CrossRef] [Google Scholar]
- Shuntov, M., McCracken, H. J., Gavazzi, R., et al. 2022, A&A, 664, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Silk, J., & Mamon, G. A. 2012, Res. Astron. Astrophys., 12, 917 [Google Scholar]
- Smit, R., Bouwens, R. J., Franx, M., et al. 2012, ApJ, 756, 14 [NASA ADS] [CrossRef] [Google Scholar]
- Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217 [NASA ADS] [CrossRef] [Google Scholar]
- Snedden, A., Coughlin, J., Phillips, L. A., Mathews, G., & Suh, I.-S. 2016, MNRAS, 455, 2804 [CrossRef] [Google Scholar]
- Sobral, D., Smail, I., Best, P. N., et al. 2013, MNRAS, 428, 1128 [NASA ADS] [CrossRef] [Google Scholar]
- Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51 [Google Scholar]
- Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15 [Google Scholar]
- Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
- Springel, V. 2010, MNRAS, 401, 791 [Google Scholar]
- Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289 [Google Scholar]
- Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776 [Google Scholar]
- Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [Google Scholar]
- Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781 [Google Scholar]
- Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157 [NASA ADS] [CrossRef] [Google Scholar]
- Taylor, M. B. 2005, ASP Conf. Ser., 347, 29 [Google Scholar]
- Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
- Thomas, N., Davé, R., Anglés-Alcázar, D., & Jarvis, M. 2019, MNRAS, 487, 5764 [NASA ADS] [CrossRef] [Google Scholar]
- Tuffs, R. J., Popescu, C. C., Völk, H. J., Kylafis, N. D., & Dopita, M. A. 2004, A&A, 419, 821 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Van de Voort, F., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 2782 [NASA ADS] [CrossRef] [Google Scholar]
- Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518 [Google Scholar]
- Wang, L., Norberg, P., Bethermin, M., et al. 2016, A&A, 592, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Nature, 572, 211 [Google Scholar]
- Weaver, J. R., Kauffmann, O., Shuntov, M., et al. 2021, Am. Astron. Soc. Meeting Abstracts, 53, 21506 [Google Scholar]
- Weaver, J. R., Davidzon, I., Toft, S., et al. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202245581 [Google Scholar]
- Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291 [Google Scholar]
- Whitaker, K. E., van Dokkum, P. G., Brammer, G., & Franx, M. 2012, ApJ, 754, L29 [Google Scholar]
- Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104 [NASA ADS] [CrossRef] [Google Scholar]
- Wiersma, R. P. C., Schaye, J., & Smith, B. D. 2009, MNRAS, 393, 99 [NASA ADS] [CrossRef] [Google Scholar]
- Witt, A. N., & Gordon, K. D. 2000, ApJ, 528, 799 [Google Scholar]
- Wuyts, S., Labbé, I., Schreiber, N. M. F., et al. 2008, ApJ, 689, 653 [NASA ADS] [CrossRef] [Google Scholar]
- Zamojski, M. A., Schiminovich, D., Rich, R. M., et al. 2007, ApJS, 172, 468 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: IRX constraints from stacking analysis
Surveys conducted in the thermal infrared regime are generally not sensitive enough for detecting low-mass sources individually. 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. Constraints can be inferred down to flux levels much fainter than the typical 3σ sensitivity limit of the initial data, depending on the number of co-added sources (e.g., Dole et al. 2006; Karim et al. 2011).
![]() |
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 |
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 Δz = 0.25 up to z = 1.3 and to Δz = 0.3 above. We adopted a color bin size of Δ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−2 down to 3σ), 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 subsample, 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 uncertainty was estimated from the dispersion of the averaged signal measured in each of the 200 stacked images. We co-added in quadrature this uncertainty with the 1σ 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 emission 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 signal 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 underlying distribution of properties associated with the source population 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 uncertainties with aperture photometry and bootstrapping, similar to the procedure already used for the mean stacking described earlier. Comparison between the two methods revealed a remarkable 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 stacking, 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.
![]() |
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. |
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 wavelength regime is key for constraining the peak of their SED and estimating their total infrared luminosity with minimal uncertainty. We stacked the COSMOS imaging performed at 250, 350, and 500 μm with the SPIRE instrument on board the Herschel 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 motivated 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. Photometry on the stacks was estimated in different ways (aperture photometry, PSF deconvolution, and Gaussian fitting), and we also compared median and mean stacking with each other. All methods led to consistent results. Flux uncertainties were finally estimated 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 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 decomposing the stacked signal into the contribution of the PSF produced 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⋆ ∼109 M⊙.
A.3. Total infrared luminosities and mean IRX estimates
For each bin of redshift, stellar mass, and NrK color, a fitting of SEDs at MIR and FIR wavelengths was performed with LePhare so as to estimate the average Infrared luminosity of galaxies in the bin. This step was carried out following the same procedure as the one adopted earlier for the individual detections (see Sect. 3.2). We fixed the 24 μm flux to the one obtained by combining the stacked signal from the residual 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 LIR 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
SIMBA6 (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 10243 DM and gas particles within a comoving volume of (100 h−1 Mpc)3. The simulation assumes a standard ΛCDM cosmology compatible with Planck Collaboration XIII (2016) (Ωm = 0.3, ΩΛ = 0.7, Ωb = 0.048, H0 = 68 km s−1 Mpc−1, σ8 = 0.82 and ns = 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 ×107 M⊙ and 9.6 ×107 M⊙, 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 H2-based and is only allowed to occur in gas with the hydrogen density nH≥ 0.13 cm−3. The SFR is computed from the H2 density and the dynamical time following the Schmidt (1959) relation with 2 percent efficiency. The H2 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 > 105 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 107.5 M⊙, BHs eject material in ∼1000 km s−1 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−1. 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 nH< 0.13 cm−3 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-AGN7 (Dubois et al. 2014) was run with the adaptive-mesh refinement code RAMSES (Teyssier 2002) within a comoving volume of (100 h−1 Mpc)3 and assuming a standard ΛCDM cosmology compatible with WMAP (Komatsu et al. 2011) (Ωm = 0.272, ΩΛ = 0.728, Ωb = 0.045, H0 = 70.4 km s−1 Mpc−1, σ8 = 0.81 and ns = 0.967). It contains 10243 DM particles (i.e., a mass resolution of MDM, res = 8 × 107 M⊙), and the initially coarse 10243 grid (initial gas resolution is Mgas, res = 1 × 107 M⊙) 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 zreion = 10 following Haardt & Madau (1996) and gas is allowed to cool to 104 K via H, He, and metals (Sutherland & Dopita 1993).
Star formation only proceeds in regions with hydrogen number density nH≥ 0.1 cm−3 (the stellar mass resolution is ≃2 × 106 M⊙), 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 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 104 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
EAGLE8 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 15043 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 ×106 M⊙ and 1.81 ×106 M⊙, 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 ΛCDM cosmology compatible with Planck Collaboration XVI (2014) (Ωm = 0.307, ΩΛ = 0.693, Ωb = 0.048, H0 = 67.77 km s−1 Mpc−1, σ8 = 0.8288 and ns = 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 107.5 K and 108.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
ILLUSTRISTNG9 (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 ΛCDM cosmology compatible with Planck Collaboration XIII (2016). In this work, we use the TNG100 simulation, with the box length of 75 h−1Mpc (≈ 110 comoving Mpc), with 18203 DM particles and 18203 initial gas cells corresponding to a mass resolution of 7.5 × 106 M⊙ for DM and 1.4 × 106 M⊙ 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 nH≥ 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).
All Tables
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.
Parameter values of the SFR function Schechter fit for both individual FIR calibration and the verified calibration.
Summary of the different simulations used in this work: SIMBA, HORIZON-AGN, EAGLE, and TNG100.
All Figures
![]() |
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. |
In the text |
![]() |
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). |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
Fig. 5. Difference between the infrared luminosity, LIR, estimated from the 24 μm photometry alone and by adding the Herschel photometry up to z ∼ 2. |
In the text |
![]() |
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 ≤ 24.5; see text). The horizontal gray lines in the bottom panel indicate the SFR threshold that corresponds to the luminous (LIRG, LIR ≥ 1011 L⊙) and ultra-luminous (ULIRG, LIR ≥ 1012 L⊙) infrared galaxies. |
In the text |
![]() |
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. |
In the text |
![]() |
Fig. 8. Comparison of the infrared luminosity, LIR, estimated from the COSMOS2020 FIR sample with that estimated with the IRX–NrK calibration ( |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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. |
In the text |
![]() |
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 (α = −1.3). The gray area corresponds at the bright end to the Eddington correction and at the faint end to the slope uncertainty (Δα = 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. |
In the text |
![]() |
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). |
In the text |
![]() |
Fig. 16. Redshift evolution of the Schechter parameters – α (top panel), SFR⋆ (middle panel), and Φ⋆ (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. |
In the text |
![]() |
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, α = −1.3 ± 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). |
In the text |
![]() |
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 |
In the text |
![]() |
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. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.