© The Authors 2024

1. Introduction

The evolution of the star formation rate density (ρ_SFR) over time follows a gradual growth from the Big Bang to cosmic noon (i.e. z ∼ 2), followed by a continuous decline by a factor ∼10 to the present day (e.g. Madau & Dickinson 2014; Bouwens et al. 2015b; Liu et al. 2018; Leslie et al. 2020 and Schreiber et al. 2015, hereafter S15). This evolution alone raises many questions about the growth and death of galaxies. For example, the contribution to the ρ_SFR of galaxies of different stellar masses (M_⋆), or the ρ_SFR at high redshift (i.e. z ≥ 4) are still largely uncertain. Our understanding of the high-redshift part (i.e. z ≥ 4) is mainly built from the ultraviolet (UV), and then deduced by correcting these UV measurements from dust attenuation (e.g. Bouwens et al. 2012b, 2012a; Schenker et al. 2013; Bouwens et al. 2015a, 2015b; Oesch et al. 2018). However, it has recently been claimed, through studies using dust-unbiased measurements (i.e. radio or far-infrared emission), that ρ_SFR is actually higher, by a factor of two to six, at high redshift (i.e. z ≥ 4, e.g. Rowan-Robinson et al. 2016; Novak et al. 2017; Lagache 2018; Gruppioni et al. 2020; Khusanova et al. 2021) compared to Madau & Dickinson (2014).

In order to deduce the cosmic star formation history, it is necessary to correctly infer the star formation rate (SFR) of galaxies. It has been shown that the SFR of star-forming galaxies (SFGs) is positively correlated with their stellar mass (e.g. Elbaz et al. 2007; Daddi et al. 2007; Whitaker et al. 2012, 2014; Speagle et al. 2014; S15; Lee et al. 2015; Delvecchio et al. 2021; Leslie et al. 2020; Popesso et al. 2023) with a small scatter of ∼0.2 − 0.3 dex (e.g. Elbaz et al. 2007; Whitaker et al. 2012; Speagle et al. 2014; S15). This correlation is called the main sequence (MS; Noeske et al. 2007) of SFGs. Although the specifics of this correlation remain a matter of debate, most recent studies tend to favour a linear MS in logarithmic space with a bending appearing at high M_⋆ and for z ≤ 2 − 3 (e.g. Whitaker et al. 2014; Magnelli et al. 2014; S15; Delvecchio et al. 2021; Leslie et al. 2020), rather than a strictly linear correlation between SFR and M_⋆ (e.g. Whitaker et al. 2012; Speagle et al. 2014).

The fraction of star formation obscured by dust is still highly uncertain at high redshift (i.e. z ≥ 3), as it could be higher up to a factor of ten above unobscured formation (Casey et al. 2018). As the correction factors for dust extinction are quite large (e.g. Calzetti et al. 1994; Madau et al. 1998; Steidel et al. 1999), this leads studies to combine SFR from the UV (uncorrected for dust extinction) and SFR deduced directly from the infrared (IR) (e.g. S15; Delvecchio et al. 2021) when seeking the total SFR and ρ_SFR. Although most of them are based on Hubble Space Telescope (HST) detected catalogues to infer the evolution of the SFR, they do not take into account the contribution of so-called HST-dark galaxies in the study of the ρ_SFR. H-dropout (Wang et al. 2019), HST-dark (Zhou et al. 2020) or optically dark and faint (Gómez-Guijarro et al. 2022a; Xiao et al. 2023b) galaxies, represent massive and highly obscured galaxies usually detected with a low significance or not at all in the optical. Although initially thought to be a marginal population with little effect on the ρ_SFR, over time it has been shown that such obscured galaxies could significantly contribute to the ρ_SFR above M_⋆ ∼ 10^10.5 M_⊙ (e.g. Wang et al. 2019; Xiao et al. 2023b).

Since the Atacama Large Millimeter/submillimeter array (ALMA) came online in 2013n it has enabled large and deep surveys, with a better resolution, in the millimetre and submillimetre range (e.g. Franco et al. 2018; Gómez-Guijarro et al. 2022a). Studies can now combine higher-quality millimetre and submillimetre measurements with IR measurements. ALMA has proven to be a powerful tool to probe the gas content of galaxies (e.g. Scoville et al. 2014; Groves et al. 2015; Scoville et al. 2016; Schinnerer et al. 2016; Kaasinen et al. 2019; Liu et al. 2019; Millard et al. 2020; Magnelli et al. 2020; Wang et al. 2022), and developing our global understanding of high-redshift galaxies (see Hodge & da Cunha 2020 for a review). Studies have also been using the discrete ALMA archive pointings to study the ρ_SFR (e.g. Traina et al. 2024).

Spatially resolved studies show that the SFR mostly correlates with the molecular gas (H₂) surface density, and very little with the atomic gas (HI) surface density (e.g. Bigiel et al. 2008; Leroy et al. 2008). The gas content of galaxies is usually probed through CO emission lines because of the difficulty of directly observing the H₂ content (see Bolatto et al. 2013 for a review). However, other techniques have been developed to infer the gas content of galaxies through the study of multi-wavelength dust spectral energy distribution (SED) fits, and by applying the gas-to-dust ratio to infer gas masses (M_gas) from dust masses (M_dust; e.g. Leroy et al. 2011; Magdis et al. 2011, 2012; Magnelli et al. 2012; Rémy-Ruyer et al. 2014; Genzel et al. 2015).

By combining the SFR and the gas content of galaxies, we can infer the star formation efficiency (SFE ≡ SFR/M_gas) of galaxies, which gives a straightforward indicator of how galaxies form stars at a different moment in cosmic time. The SFE has been shown to evolve quickly at a low redshift (i.e. 0 ≤ z ≤ 1), but to be rather constant at higher redhsifts (i.e. z >  1; e.g. Saintonge et al. 2017; Scoville et al. 2017; Tacconi et al. 2018; Liu et al. 2019; Wang et al. 2022). Another way of looking at this is to use the Kennicutt-Schmidt relation (Kennicutt 1998b). This relation has been shown to follow the power-law correlation between the SFR and gas surface density with a slope of ∼1.0 − 1.5 (e.g. Kennicutt 1998b; de los Reyes & Kennicutt 2019; Wang et al. 2022).

In this paper, combining Spitzer, Herschel and ALMA, we study the evolution of galaxy properties from z ∼ 0 to 5 through stacking on four fields: COSMOS, UDS, GOODS-South, and GOODS-North. We try to answer a few open questions about the evolution and global history of galaxies. In particular, we seek to constrain the main sequence, ρ_SFR at a high redshift (i.e. z ≥ 4), the different contributions to ρ_SFR as a function of M_⋆, and the impact of the H-dropout galaxy population on the properties of galaxies as a whole. A similar study using Spitzer and Herschel measurements has been carried out by S15. However, the addition of ALMA measurements to Herschel using stacking has not been carried out yet. The reason is that this requires access to an extensive study programme on ALMA. This has only recently been possible thanks to blind surveys such as the GOODS-ALMA survey (Franco et al. 2018, 2020; Gómez-Guijarro et al. 2022a). A summary of the data used in this study is given in Sect. 2. The stacking method is described in Sect. 3. Section 4 is devoted to the SED fitting procedure. Section 5 reviews the properties of galaxies that can be deduced from this analysis. Sections 6 and 7 present the cosmic star formation history and the cosmic evolution of the gas mass density deduced for this work, respectively. We interpret and discuss our results in a global cosmological context in Sect. 8, and in Sect. 9 we summarise the main results and conclusions.

For this work, we adopted a Salpeter (1955) initial mass function (IMF) and the cosmological parameters (Ω_M, Ω_Λ, h) = [0.30, 0.70, 0.70]. We used a factor of 1.7 to convert M_⋆ and SFR from a Chabrier (2003) to a Salpeter (1955) IMF whenever necessary (e.g. Reddy et al. 2006; Santini et al. 2012; Elbaz et al. 2018; Gómez-Guijarro et al. 2022b). When mentioned, magnitudes are in the AB system, such that $M_{\mathrm{AB}}=23.9-2.5{log}_{10}(S_{\nu }[\mathrm{\mu }\mathrm{Jy}])$.

2. Sample and observations

2.1. Sample

In this study, we worked from catalogues of H-band selected SFGs on four fields: GOODS-South, GOODS North, COSMOS and UDS. We used the ultra-deep H-band catalogue of the CANDELS-HST team (Grogin et al. 2011; Koekemoer et al. 2011) for GOODS-South (Guo et al. 2013), COSMOS (Nayyeri et al. 2017) and UDS (Galametz et al. 2013) fields, while in GOODS-North we use the catalogue from Barro et al. (2019). The 5σ limiting magnitude range from H ∼ 27.4 to 29.7 for GOODS-South, H ∼ 27.4 to 28.8 for COSMOS, H ∼ 27.1 to 27.6 for UDS, and H ∼ 27.8 to 28.7 for GOODS-North. The photometric redshifts and M_⋆ of the galaxies in these catalogues were derived in S15 (GOODS-South, COSMOS and UDS) and Barro et al. (2019) (GOODS-North). Photometric data were fitted up to IRAC 4.5 μm, with EAZY (Brammer et al. 2008), by assuming a delayed exponentially declining star formation history with the Bruzual & Charlot (2003) stellar population synthesis model. The SFGs are UVJ selected, following the definition from Muzzin et al. (2013). The final number of SFGs in the sample is given in Table 1. H-dropout galaxies are by definition not included in these samples, as they are not detected in H-band. However, we used the sample from Wang et al. (2019) to add and discuss their impact in any necessary analyses.

2.2. Observations

Observations with Spitzer-MIPS at 24 μm include maps of the COSMOS field (PI: D. Sanders; LeFloc’h et al. 2009), the GOODS-South and GOODS-North fields (GOODS Legacy programme; PI: M. Dickinson), and the UDS field (SpUDS Spitzer Legacy programme; PI: J. Dunlop). The Herschel PACS and SPIRE maps of the four fields come mainly from the CANDELS-Herschel programme. The PACS GOODS-North and GOODS-South maps are the combined Herschel-PACS data from the PEP (Lutz et al. 2011) and GOODS-Herschel (Elbaz et al. 2011) programmes, as described in Magnelli et al. (2013). For ALMA, we used the $1.13\mathrm{mm}$ low-resolution GOODS-ALMA map (sensitivity reaching an average of $\sigma =95.2\mathrm{\mu }{\mathrm{Jy}\mathrm{beam}}^{-1}$, Gómez-Guijarro et al. 2022a), within the GOODS-South field with the $1\stackrel{″}{.}1$ circularised point spread function (PSF) full width half maximum (FWHM). We note that there is a global and local offset between the position of sources in the ALMA and HST images (Franco et al. 2018, 2020). We corrected them using the offsets provided by Franco et al. (2020).

3. Stacking

Our method consists in stacking several images of galaxies on top of each other in order to increase the overall signal-to-noise ratio (S/N). The result of a stacking procedure is the mean or median flux for all galaxies stacked together. The main advantage of this method is to be able to recover reliable flux measurements for populations of galaxies whose S/N is too low when studied individually. In practice, this means being able to measure the properties of galaxies with lower M_⋆ and higher redshift. An effective way to take advantage of the stacking method is to group galaxies into sub-populations with similar properties. To this end, in this study, we stack our sample of galaxies over different redshift and M_⋆ bins. The total numbers of galaxies per bin are displayed in the Table 2.

In order to access the completeness of the catalogues that we stacked in this study, we followed the method of S15. In summary, we assumed, for each redshift bin, that the observed luminosity at 1.6 μm ($L_{1.6\mathrm{\mu }\mathrm{m}}^{\mathrm{obs}}$) can be related to M_⋆ simply by $M_{\star }=C\times L_{1.6\mathrm{\mu }\mathrm{m}-\mathrm{obs}}^{\alpha }$ with a scatter. The relation and dispersion are fitted from the H-band photometry of the catalogues. We derive the completeness from Monte Carlo simulations by generating, from a uniform redshift distribution in the redshift bin and a given M_⋆, the corresponding $L_{1.6\mathrm{\mu }\mathrm{m}-\mathrm{obs}}$ taking into account the dispersion, and comparing it to the corresponding H-band detection limit for each catalogue. We assume that catalogues are complete when the completeness is above 90% (S15). The corresponding $M_{\star }^{90\%}$, above which all catalogues are assumed to be mass complete, are listed in Table 3. All but two of the bins which yield significant detection (i.e. with S/N ≥ 3σ in at least one band from 24 μm to $1.13\mathrm{mm}$) considered in this study are mass complete: for 9.5 ≤ log10(M_⋆/M_⊙)≤10.0 at 2.3 ≤ z ≤ 3.1, and 10.0 ≤ log10(M_⋆/M_⊙)≤10.5 at 3.1 ≤ z ≤ 3.9. The latter bin is, however, almost mass complete as $M_{\star }^{90\%}={10}^{10.1}{M}_{\odot }$ and is consistent with a $M_{\star }^{90\%}={10}^{10}{M}_{\odot }$ within the uncertainties of the method used to compute the completeness.

The completeness computed in this paper is determined assuming that the H-band is a good proxy for the stellar mass. This hypothesis may however become less robust at high redshifts (i.e. z >  3). To assess whether our derivation of the completeness limit remains relevant for the redshift range 3 <  z <  5, we looked at a deep JWST-NIRCAM survey of the EGS field (Gómez-Guijarro et al. 2023). We found that 89.7% of the galaxies with log₁₀(M_⋆/M_⊙)> 10 (209 galaxies detected, 233 overall) detected in the JWST bands are also detected in the HST H-band (the 5σ limiting magnitude is 27.37 in CEERS – see Gómez-Guijarro et al. 2023, Sect. 2.2 – that is comparable to 27.4, 27.4, 27.1, 27.8 in the GOODS-South, COSMOS, UDS GOODS-North fields respectively). This is consistent with our 90% completeness of ${log}_{10}(M_{\star }^{90\%}/M_{\odot })\sim 10$ for 3.1 ≤ z ≤ 5 and 3.9 ≤ z ≤ 5.0, which gives us confidence that our sample is not strongly biased.

We choose to use mean stacking for our work, rather than median stacking. Indeed, although median stacking better suppresses secondary sources in the stacked image (i.e. bright sources which are close to the stacked target and might appear in the final stacked image), it has been shown that it also yields to systematically biased measurements at low S/N or in the presence of a flux distribution skewed towards low or high values (White et al. 2007; S15). It comes from the fact that the median (⟨.⟩) is not a linear operation: ⟨a + b⟩ ≠ ⟨a⟩+⟨b⟩. As shown in S15, correcting for these systematic biases involves making strong assumptions on the actual flux distribution of the stacked sources, which is non-trivial and very uncertain. On the contrary, the mean stack gives us access to the total flux in the stack, which we currently miss to infer ${\rho }_{\mathrm{SFR}}^{\mathrm{IR}}$.

During the stacking procedure, it is possible to treat detected and undetected galaxies separately. The main way to do this is to stack only the undetected sources on the residual map, and then add them to the fluxes of the detected galaxies via a weighted mean (e.g. Magnelli et al. 2009). Although this method reduces the confusion noise of faint sources and removes most of the contamination from bright neighbours, it can also introduce some biases (S15). Following S15, we have therefore chosen to treat detected and undetected sources in the same way for consistency (i.e. directly from the image and not from the residual). The contamination of bright neighbours is dealt with later in the study (see Sect. 3.2).

We noticed the presence of a global background gradient in most Herschel images. In order to deal with it, we have decided to rotate the stacked postage stamp images successively by 90°.

The resulting flux density in each stacked image was obtained via standard aperture measurements or PSF fitting method (see detail for each wavelength in Sects. 3.1, 3.2 and 3.3). The corresponding S/N was then derived using a simple Monte Carlo approach, that is, S/N = S_ν/σ_MC, where S_ν is the flux density measured at the centre of the stacked image within an aperture of radius, r, while σ_MC is the standard deviation of the signal in 100 circles, of the same radius r, randomly positioned on the edges of the stacked stamp. Depending on the method used to compute the flux, r corresponds to either the radius used for the aperture measurement method, or the radius within which the PSF fitting method was performed.

In addition to this photometric noise, there is an uncertainty in the recovered flux density due to the intrinsic dispersion of the underlying stacked population (i.e. all galaxies in the stack do not have the exact same flux density). This flux dispersion can be quantified using a bootstrap analysis (e.g. S15). The method consists in repeating several time the full stacking and flux density measurement process, picking, for each realisation, the galaxies from the origin sample, with replacement, until the number of galaxies in the original sample is reached. In this study, we compute a 100 runs. The error is then deduced from the standard deviation of the computed flux densities for this bin.

3.1. Spitzer

As the 24 μm MIPS-Spitzer map ($\sim 5\stackrel{″}{.}8$ PSF FWHM) has a better resolution than Herschel-PACS (∼7′′ and $\sim 11\stackrel{″}{.}3$ PSF FWHM at 100 μm and 160 μm, respectively), we do not expect the clustering bias (see Sect. 3.2) to be a dominant effect on the stacked stamp. The fluxes are calculated using a classical aperture photometry method. The aperture is chosen to be 4 pixels (i.e. $4\stackrel{″}{.}8$) in radius, which contains about ∼60% of the total flux for a point source.

3.2. Herschel

When stacking galaxies, a ‘clustering bias’ can occur due to the neighbouring galaxies of the main stacking targets. Such contamination from neighbouring sources can become significant when the size of the PSF becomes comparable to the typical cluster length of SFGs. It has been shown that the clustering bias has a non-negligible impact on the results when stacking galaxies in the Herschel bands (e.g. Bavouzet et al. 2008; Béthermin et al. 2010; Kurczynski & Gawiser 2010; Bourne et al. 2012; Béthermin et al. 2012; Viero et al. 2013; Béthermin et al. 2014; S15; Béthermin et al. 2015; Delvecchio et al. 2021) due to the large PSF of Herschel. Correcting for this bias is crucial to accurately measure the peak of SED of SFGs, as this bias tends to cause the fluxes of Herschel to be increasingly overestimated with increasing wavelength (as the size of the PSF increases with the wavelength).

To mitigate this clustering signal contamination, we choose to follow the method presented in S15. It consists in fitting only a PSF and a local background term:

$$\begin{array}{c}\hfill S(x,y)=\phi \times PSF(x,y)+\epsilon ,\end{array}$$(1)

where φ and ε are the normalisation of the source flux and background, respectively. The fit is performed on a fixed aperture of radius of $0.9\times \mathrm{FWHM}$, as this was found to minimise the clustering contamination to φ (S15). The clustering bias signal is in this case largely included in the background term. Nevertheless, even with this radius the contamination of φ is not null and we still need to apply a correction for what remains of the clustering signal in the mean flux term (φ). The correction factors were calculated by simulating the stacking procedure on mock images (see S15 for more details). The correction factors of S15 are listed in Table 4. For the PACS maps, we used the truncated PSF derived from Vesta, while for Herschel-SPIRE the PSF is assumed to be Gaussian with a FWHM of $18\stackrel{″}{.}15$, $25\stackrel{″}{.}15$ and $36\stackrel{″}{.}3$ respectively at 250, 350 and 500 μm according to Griffin et al. (2010) (see also Shirley et al. 2021).

Finally, for Herschel-PACS, it was shown (Popesso et al. 2012; Magnelli et al. 2013) that the high-pass filter data reduction technique that was used to remove low-frequency noise in the maps could induce an underestimation of the photometric measurements of the unmasked faint sources. It was shown in Popesso et al. (2012) that a correction factor of 17%, and 10% should be taken into account when stacking undetected faint, and detected, sources in the PACS maps (100 μm and 160 μm), respectively.

3.3. ALMA

For the purpose of this study, we decided to work in the image plane when dealing with ALMA. Another option would have been to work in the uv plane (i.e. in Fourier space), as ALMA provides us with measurements in the uv space. Although this may give more robust results, since we do not need to go through the conversion between the uv plane and the image plane first, it is very computationally intensive. Furthermore, on few test bins we only observed a 10% difference between fluxes calculated from stacks in the uv and image planes. We decided not to stack directly in the uv plane but in the image plane in order to save computing time as this should not impact the results of this study.

As the ALMA 1.1$\mathrm{mm}$ map has a much better resolution ($\sim 1\stackrel{″}{.}1$ PSF FWHM; which is comparable to typical sizes of individual SFGs, Suess et al. 2019; Wang et al. 2022) than Herschel-PACS (∼7′′ and $\sim 11\stackrel{″}{.}3$ PSF FWHM at 100 μm and 160 μm, respectively), the clustering bias is negligible on the stacked stamp. Due to the high angular resolution of the ALMA data, and to possible small offsets between the optical centroids (on which our stack positions are based) and millimetre centroids, the stacked images are extended on scales larger than that of the ALMA PSF. However, this has no significant impact of the measure flux density via the aperture photometry method. To measure ALMA flux densities, we applied an aperture photometry within a radius of 1′′, containing 85% of the total flux density (i.e. $S_{\mathrm{AP}}^{\mathrm{PSF}}/S_{\mathrm{tot}}^{\mathrm{PSF}}=0.85$).

4. SED fitting

4.1. SED fitting procedure

One of the objectives of this paper is to retrieve several properties such as SFR, dust temperature (T_dust), dust masses and gas masses from our far-infrared (FIR) to submillimetre stacking analysis, and compare them to the literature. We performed an SED fit from the measured fluxes for each bin of redshift and M_⋆. Because we do not have many points on the FIR at low M_⋆ and high redshifts, we used the library from Schreiber et al. (2018, hereafter S18), which is well suited for our study given its small number of free parameters. Choosing a model with more parameters, such as Draine & Li (2007) and Draine et al. (2014), would imply fixing parameters on some bins. The library of S15 is calibrated on galaxies from z = 0.5 to z = 4, which allows to have a realistic SED while reducing the number of free parameters. We have assumed a form of SED that we would expect for the main sequence galaxies, and no active galactic nuclei (AGN) contribution to the stacked SED.

To fit the SED, we only considered fluxes with S/N ≥ 3. Stacked flux measurements with S/N <  3 are replaced by conservative 5σ upper limits. At z >  4 MIPS-24 stacked fluxes, if available, have been transformed to 5σ upper limits as they are no longer dominated by dust and polycyclic aromatic hydrocarbon (PAH) emissions, but rather by stellar emission. From these SED fits, we measured the corresponding infrared luminosity (L_IR) by integrating the best-fit SED in the range 8–1000 μm rest frame. The error on L_IR was obtained by varying their stacked photometry randomly within their uncertainties.

The T_dust can be defined in different ways (i.e. weighted by mass or luminosity). The luminosity weighted $T_{\mathrm{dust}}^{\mathrm{L}}$, for the library from S18, is calculated from a grey body of effective emissivity β = 1.5 (S18). This means that $T_{\mathrm{dust}}^{\mathrm{L}}$ follows Wien’s law (see Eq. (2)).

$$\begin{array}{c}\hfill T_{\mathrm{dust}}^{\mathrm{L}}[\mathrm{K}]=2.897\times {10}^3/({\lambda }_{\max }[\mathrm{\mu }\mathrm{m}]),\end{array}$$(2)

where λ_max is the wavelength corresponding to the peak of λ^βL_λ. A mass weighted $T_{\mathrm{dust}}^{\mathrm{M}}$, was also calculated for each template by mass weighted averaging each individual template of Galliano et al. (2011) (see S18 for more details). In the cases where we only had one or two points to perform our fit (i.e. the SED peak was not well defined in this case), we chose to restrict the $T_{\mathrm{dust}}^{\mathrm{M}}$ during the fit to the $T_{\mathrm{dust}}^{\mathrm{M}}$ evolution from S18 ±10 K. This is a reasonable way to reduce the error on the deduced L_IR by slightly restricting the $T_{\mathrm{dust}}^{\mathrm{M}}$ to reasonable values. Our best-fits of SEDs are displayed in Fig. 1.

Fig. 1. Best-fit SED for each bin of redshift and stellar mass. Blue dots correspond to the flux measurements, red triangles represent the 5σ upper limits. Blue line is the best-fit SED, the blue shaded area shows the 68% uncertainty of the fit. Red line is the SED maximising LIR in case only upper limits are available.

The templates from the library of S18 are built using the amorphous carbon model from Galliano et al. (2011). This differs from the model from Draine & Li (2007) which takes into account amorphous silicate and graphite grains. This change of model was shown in Galliano et al. (2011) to lower the tension, in the Large Magellanic Cloud and the Milky Way, between the observed dust-to-gas ratio and the stellar abundances. The differences between the two models are mainly reflected in the different emissivity. The choice of a different emissivity does not affect the dust temperature or the L_IR, as these properties correspond to the peak and the area under the SED, respectively. However, the M_dust deduced from Galliano et al. (2011) is about a factor of 2 lower than those deduced from Draine & Li (2007), without affecting the M_gas as the factor cancels out when converting M_dust to M_gas. The debate on the composition of dust grains goes far beyond the scope of this study. We chose to work with M_dust derived from the Draine & Li (2007) model, as it will ease comparison with the literature, which is widely base on the later model. In practice, we have re-fitted the model library of S18 with a Draine & Li (2007) model to associate the M_dust. The impact on M_gas of choosing an amorphous carbon model instead of an amorphous silicate and graphitic grains model is briefly discussed at the end of the Sect. 5.3.

Overall we verified that fitting our stacked flux with the S18 or Draine & Li (2007) libraries has little to no impact on our results, as the differences on some key properties such as L_IR, T_dust and M_dust (once corrected for the emissivity chosen in S18) are quite small: $-7{\%}_{-4\%}^{+7\%}$ for L_IR, $1{\%}_{-6\%}^{+12\%}$ for T_dust and $4{\%}_{-17\%}^{+33\%}$ for M_dust. An essential consequence of this observation is that the models of S18 reproduce well the global SED shape of stacked galaxies. This reinforces the choice we made to use the S18 templates instead of Draine & Li (2007) for this study.

4.2. Adding ALMA to Herschel

The addition of ALMA to this study results in an improvement of the SED at high redshift (i.e, at 3.1 ≤ z ≤ 5.0). More specifically, it allows to get an ALMA measurement at 10.0 ≤ log₁₀(M_⋆)≤10.5 for 3.1 ≤ z ≤ 3.9 and 3.9 ≤ z ≤ 5 that provides information on some of the properties of galaxies (i.e. L_IR and M_gas) at these redshifts and M_⋆ instead of what would have just been an upper limit on these properties. At 10.5 ≤ log₁₀(M_⋆)≤11.0 for 3.1 ≤ z ≤ 3.9 and 3.9 ≤ z ≤ 5, and 11.0 ≤ log₁₀(M_⋆)≤12.0 for 3.9 ≤ z ≤ 5, this also allows the SED peak to be constrained more effectively, as it is only weakly constrained, if at all, by the Herschel data alone. The SED peak is crucial for determining properties such as $T_{\mathrm{dust}}^{\mathrm{M}}$ and M_gas. Elsewhere it also provides an extra point or upper limit that significantly reduces the uncertainty in the properties inferred from the SED fit. These extra ALMA point are thus decisive for this study, which aims to probe the properties of galaxies over a wide dynamic range, extending up to z ∼ 5.

4.3. Simulations to correct for averaging biases

Because a mean stacking procedure is luminosity weighted, it can have non-linear effects on the shape of the resulting SED (e.g. Elbaz et al. 2011; S15; S18). These effects include the widening of the FIR bump and a bias of the peak towards warmer $T_{\mathrm{dust}}^{\mathrm{L}}$. It is mainly a result of mixing galaxies of different redshifts and $T_{\mathrm{dust}}^{\mathrm{L}}$ (e.g. Elbaz et al. 2011; S15; S18). The broadening of the SED increases the difficulty of determining precisely $T_{\mathrm{dust}}^{\mathrm{L}}$ (S18). For the rest of the paper, we chose to work with mass-weighted dust temperatures, and we will refer to $T_{\mathrm{dust}}^{\mathrm{M}}$ simply as T_dust. To ensure that our conclusions are not biased by these effects, we performed simulations to identify any systematic bias due to our stacking procedure.

To this end, for each bin of M_⋆ and redshift, we simulated the biases coming from stacking starting from the distribution of the galaxies in the bin. This has the advantage of taking into account the specificity of the M_⋆ and redshift distributions within our bin. For each set of galaxies, we created a mock counterpart in order to compare the properties resulting from the stacking with those expected.

To each galaxy of mass M_⋆, we assigned the SFR starting from M_⋆ and following SFR = R_SB × SFR_MS. The SFR_MS was calculated from the main sequence trend found in this work (See Eq. (12) and parameters Table 5) as a function of redshift and M_⋆. R_SB represents the starburstiness and is defined by R_SB = SFR/SFR_MS. We want to generate both mock main sequence galaxies and starburst galaxies. It was shown in S15 that both the main sequence width (∼0.3 dex, see S15) and the starburst fraction do not evolve with redshift and M_⋆. We can therefore reasonably assume that the distribution function of R_SB does not vary (S15). This assumption still allows the luminosity functions to be reconstructed properly (e.g. Sargent et al. 2012; S15). Following Sargent et al. (2012), we have modelled the probability density function of R_SB by a double log-normal distribution (see Eq. (3)):

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{R_{\mathrm{SB}}}(x)& =\frac{1-f_{\mathrm{SB}}-f_{\mathrm{miss}}}{\sqrt{2\pi }{\sigma }_{\mathrm{MS}}}exp\left(\frac{-{log}_{10}{(x/x_0)}^2}{2{\sigma }_{\mathrm{MS}}^2}\right)\hfill \\ \hfill & +\frac{f_{\mathrm{SB}}}{\sqrt{2\pi }{\sigma }_{\mathrm{SB}}}exp\left(\frac{-{log}_{10}{(x/B_{\mathrm{SB}})}^2}{2{\sigma }_{\mathrm{SB}}^2}\right),\hfill \end{array}$$(3)

where f_SB is the fraction of starbursts, f_miss is the fraction of galaxies missed by such distribution (neither starburst nor main sequence galaxies), σ_MS and σ_SB are the widths of the main sequence and starburst distributions, B_SB is the median multiplicative boost of star formation that can be expected for a starburst compared to a main sequence galaxy (i.e. the median of starburst galaxies), and x₀ is the median R_SB of main sequence galaxies. We note that with this parametrisation, we expect f_miss and x₀ to be close to 0 and 1 respectively, by construction. We have chosen here to use the parametrisation of S15: σ_MS = σ_SB = 0.31 ± 0.02 dex, f_SB = 3.3%±1.5%, B_SB = 5.3 ± 0.4, f_miss = 0%±2%,and x₀ = 0.87 ± 0.04.

To each galaxy we assign L_IR, deduced from the SFR by subtracting the UV SFR assuming the UV dust attenuation (A_UV) derived from M_⋆ as in Pannella et al. (2015) (see Eqs. (4) and (5)) and Kennicutt (1998a) (see Eq. (6)):

$$\begin{array}{c}\hfill A_{\mathrm{UV}}=1.6\times {log}_{10}(M_{\star })-13.5,\end{array}$$(4)

$$\begin{array}{cc}& SF{R}_{\mathrm{IR}}=SFR-SF{R}_{\mathrm{UV}}\hfill \\ \hfill & \mathrm{where}SF{R}_{\mathrm{UV}}=SFR\times {10}^{-0.4\times A_{\mathrm{UV}}},\hfill \end{array}$$(5)

$$\begin{array}{c}\hfill L_{\mathrm{IR}}[L_{\odot }]=5.8\times {10}^9\times SF{R}_{\mathrm{IR}}[M_{\odot }{\mathrm{yr}}^{-1}].\end{array}$$(6)

Then T_dust was calculated using the best-fit of our work (see Eq. (7)) for the main sequence trend, and we followed Magnelli et al. (2014) (see Eq. (9)) to take into account the impact on T_dust of the distance of the mock galaxy from the main sequence (R_SB). The SED, for each mock galaxy, was then calculated using T_dust, L_IR and redshift with the template library from S18 assuming a contribution of PAH molecules $f_{\mathrm{PAH}}\equiv M_{\mathrm{dust}}^{\mathrm{PAH}}/M_{\mathrm{dust}}$. The value of f_PAH was set to follow a Gaussian distribution with a mean of 0.039 and a scatter 2.5/100 (S18, roughly for main sequence galaxies). Next, the M_dust were obtained from the selected SED template and the M_gas were calculated according to the different methods presented in Sect. 5.3. The individual SEDs were then stacked using a mean stacking method. Fluxes at 24 μm, 100 μm, 160 μm, 250 μm, 350 μm, 500 μm and 1130 μm were deduced from the stacked SED. The fluxes were then fitted with the template library from S18. The properties of the stacked SED were deduced using the method presented in the corresponding section of this paper: see Sect. 5.1 for T_dust, Sect. 5.2 for L_IR and Sect. 5.3 for M_gas. Potential biases were then investigated by comparing the actual average properties with those deduced from our stacking analysis.

The relative differences between the actual average properties with those deduced from our stacking analysis are quite small: $-1{\%}_{-4\%}^{+2\%}$ for L_IR, $-6{\%}_{-2\%}^{+1\%}$ for T_dust and $7{\%}_{-4\%}^{+2\%}$ for M_dust. As a result, we find no clear evidence of significant averaging bias and thus decided not to apply any correction.

4.4. Active galactic nuclei bias

It has been reported that AGN can have a major contribution to the total outgoing light of a galaxy (e.g. Hao et al. 2005; Richards et al. 2006). However, most of this emission is radiated at wavelengths shorter than 24 μm and thus will not affect our FIR measurement. Most extreme AGN may still have an impact on the mid-to-far infrared ratio (in particular 24 μm in our case), but should not affect the FIR colour compared to normal star-forming galaxies (Hatziminaoglou et al. 2010). We checked for any AGN contribution by fitting our SED with a combination of Draine & Li (2007) dust model, and Fritz et al. (2006) AGN model. No conclusive evidence for a major contribution (i.e. $L_{\mathrm{IR}}^{\mathrm{AGN}}/L_{\mathrm{IR}}^{\mathrm{TOT}}\ge 10\%$) from AGNs to the rest-frame FIR was found in any of our M_⋆ and redshift bins.

5. Stacked galaxies properties

5.1. Dust temperature

In this section, we examine the evolution of T_dust. In Fig. 2, we display the corresponding T_dust for our SED as a function of redshift and M_⋆. We do not see a significant dependence of T_dust on M_⋆. It seems that, overall, the T_dust of a main sequence galaxy is mainly determined by its redshift, independently of its M_⋆ (see also Magdis et al. 2012; Magnelli et al. 2014).

Fig. 2. Tdust as a function of M⋆ colour-coded by the redshift bin. Average trend of Tdust as a function redshift from S18, is shown as a reference in faded coloured dashed line.

The temperature of the cosmic microwave background (CMB) ranges from 2.73 K at z = 0, up to 16.4 K at z = 5, and, thus, it could become a significant source of heating at high redshift. We, however, verified, following da Cunha et al. (2013) that this is not the case. Indeed, the observed T_dust of galaxies at high redshift is significantly higher than that of the CMB.

Because our results come from H-band selected galaxies, our stacking analysis does not take into account the H-dropout galaxies (Wang et al. 2019). There are 63 H-dropout galaxies that have been detected in Wang et al. (2019) over ∼600 ${\mathrm{arcmin}}^2$. In this work, we have 1464 galaxies with M_⋆ ≥ 10¹⁰ M_⊙ and z ≥ 3 over ∼1077 ${\mathrm{arcmin}}^2$. The H-dropout would then only account for ∼7.2% of the total sample. We re-fitted the stack of H-dropout from Wang et al. (2019) with S18 templates and deduced a $T_{\mathrm{dust}}^{H\text{-drop}}=37.{34}_{-1.27}^{+1.2}\mathrm{K}$. The contribution of H-dropout from Wang et al. (2019) was then added to our last bin of redshift (i.e. 3.9 ≤ z ≤ 5).

In Fig. 2, our result seems to be globally consistent with the trend of S18. We see only slight evidence that the highest M_⋆ bin (11 ≤ log₁₀(M_⋆)≤12) could be cooler than lower-mass galaxies at a fixed redshift. This effect can be observed for z ≤ 2.5. This could simply show that high-mass galaxies are actually starting to slowly reduce their SFE on their way to become quiescent. This is particularly apparent in our first redshift bin (0.1 ≤ z ≤ 0.4) in which the bending of the main sequence is also the strongest (see Sect. 5.2; S15).

In Fig. 3, we display the mean T_dust over each redshift bin, weighted by the number of galaxies in each bin of M_⋆, as a function of redshift. Our analysis suggests a linear evolution of T_dust as a function of redshift. Our best linear fit, for 0.1 ≤ z ≤ 5, follows:

$$\begin{array}{c}\hfill T_{\mathrm{dust}}[\mathrm{K}]=34.{20}_{-0.15}^{+0.14}+5.{06}_{-0.03}^{+0.03}(z-2).\end{array}$$(7)

Fig. 3. Tdust and ⟨U⟩ as a function of redshift. The blue dots represent the Tdust of this work, the blue line is the best-fit (up to z = 5), the dashed blue line is the best-fit extrapolation (for z ≥ 5), and the blue shaded area represents the 68% uncertainty of the fit. From the literature: S18 (grey solid line), Magnelli et al. (2014) (grey dotted line) and Bouwens et al. (2020) (grey dash-dotted line) converted using Eq. (8). Magdis et al. (2012) (grey, dashed line) and Béthermin et al. (2015) (grey dots) converted using Eq. (9). We also re-fitted the two stacks from Béthermin et al. (2020) (grey dotted error bars) using template library from S18.

We compare our results with those in the literature: S18, Magnelli et al. (2014), Bouwens et al. (2020), Magdis et al. (2012) and Béthermin et al. (2015). We also re-fitted the two stacks of Béthermin et al. (2020) using the template library from S18. For consistency, we converted all $T_{\mathrm{dust}}^{\mathrm{L}}$ to $T_{\mathrm{dust}}^{\mathrm{M}}$ when comparing our results with the literature (i.e. for Magnelli et al. 2014; Bouwens et al. 2020). To do this, we adopt the conversion factor, between mass and light weighted dust temperature, given by S18:

$$\begin{array}{c}\hfill T_{\mathrm{dust}}^{\mathrm{M}}[\mathrm{K}]=0.91\times T_{\mathrm{dust}}^{\mathrm{L}}[\mathrm{K}].\end{array}$$(8)

This conversion factor represents the average conversion factor between $T_{\mathrm{dust}}^{\mathrm{M}}$ and $T_{\mathrm{dust}}^{\mathrm{L}}$ for each individual template. Some studies (i.e. Magdis et al. 2012; Béthermin et al. 2015) consider the mean starlight heating rate (⟨U⟩) (Draine & Li 2007; Draine et al. 2014) instead of T_dust as defined here. To compare our results with those, we have chosen to use the ⟨U⟩-to-T_dust conversion formula of S18 (see Eq. (9)):

$$\begin{array}{c}\hfill T_{\mathrm{dust}}[\mathrm{K}]={\left(⟨U⟩[U_{\odot }]\right)}^{1/5.57}\times 18.2\mathrm{K}.\end{array}$$(9)

As ⟨U⟩ is only a proxy of T_dust, any comparison between the two quantities should be essentially qualitative.

Our T_dust as a function of redshift is consistent within the uncertainties with S18, extended to z = 5. We find no clear evidence for the softening of the T_dust as the redshift increases reported by Magnelli et al. (2014), Magdis et al. (2012). This could stem from the lack of clustering bias correction in their study, as previously reported in S18. Clustering bias correction is a quite important step when stacking in Herschel passbands. Especially in the SPIRE wavelengths, as this effect can account for up to 50%, of the total signal, on average at 500 μm (e.g. S15; Béthermin et al. 2015; Delvecchio et al. 2021). Ignoring this could result in a cooler SED and a lower T_dust. Comparing our trend to other works by Béthermin et al. (2015, 2020), Bouwens et al. (2020), all agree on a linear trend. The differences with our trend may again stem from the way in which the correction for clustering bias is handled in the two works, as it is essential to correctly determine T_dust.

Another way to look at the redshift evolution of the FIR SED is to consider the wavelength of the IR SED bump (λ_peak) as a function of L_IR (see Fig. 4). We observe no dependence of λ_peak as a function of L_IR at fixed redshift, up to z ≃ 4. Only our last two redshift bins (i.e. 3.1 ≤ z ≤ 3.9 and 3.9 ≤ z ≤ 5.0) shows some evolution. However, for the same reason as above, this could result from a bias due to selection effects, as well as the fact that we only have two points in these redshift intervals. We conclude that λ_peak of main sequence galaxies does not depend significantly on L_IR at a given redshift. This property has already been observed by Magnelli et al. (2014), up to z ≃ 2, where no clear dependence of T_dust has been observed along the main sequence in the SFR − M_⋆ plane. However, the normalisation seems to evolve with redshift, reflecting the smooth increase of λ_peak with redshift presented in Drew & Casey (2022). But our results differ from those of Drew & Casey (2022) who studied a different sample built on individual detections and found an evolution of T_dust, as probed by λ_peak, with L_IR. This discrepancy comes from the probable incompleteness of the sample of Drew & Casey (2022). In this paper, we stack galaxies to recover the main sequence, and we look at a mass complete sample, where non-detections are accounted for thanks to the stacking technique. On the other hand, Drew & Casey (2022) worked with individually detected galaxies, and thus might be biased towards the brightest galaxies at all redshifts (i.e. starbursts). In addition, their flux-limited sample suffer from selection effect: at low L_IR, their λ_peak comes from low redshift galaxies, while at high L_IR, their values are coming from distant galaxies. What Drew & Casey (2022) find is an effect of redshift not an effect of L_IR.

Fig. 4. λpeak as a function of LIR. The colour-coded error bars (per redshift bins) represent the results of this work, coloured lines represent the best-fit and shaded area the 68% uncertainty of the fit. The solid black line corresponds to the trend of Drew & Casey (2022), while the dotted black line is an extrapolation of their relation.

5.2. The main sequence of star-forming galaxies

Here we constrain the SFR − M_⋆ correlation, also called the MS of star-forming galaxies (Noeske et al. 2007). The preferred method for estimating the SFR of a galaxy is to study the light from OB stars, because of the close link between their short lifetimes and the instantaneous SFR of galaxies. Although most of their light is emitted in the UV, it can be largely absorbed by dust and then re-emitted as thermal radiation in the IR. To obtain the total SFR associated with a galaxy, or in our case the stack of several galaxies, it is necessary to combine the SFRs deduced from both the UV and the IR as:

$$\begin{array}{c}\hfill SF{R}_{\mathrm{Tot}}=SF{R}_{\mathrm{UV}}+SF{R}_{\mathrm{IR}}.\end{array}$$(10)

The SFR_UV, uncorrected for dust attenuation, is computed from L_UV (rest-frame 1500 Å UV luminosity) following Daddi et al. (2004):

$$\begin{array}{c}\hfill SF{R}_{\mathrm{UV}}[M_{\odot }{\mathrm{yr}}^{-1}]=2.17\times {10}^{-10}{L}_{\mathrm{UV}}[L_{\odot }].\end{array}$$(11)

For our catalogues, the L_UV were calculated for all individual galaxies from EAZY (Pannella et al. 2015). In this study, the SFR_UV was derived for each redshift and M_⋆ bin by averaging the SFR_UV of all individual galaxies in the bin. The L_IR was obtained from the SED fit of our stacks, by integrating the best-fit SED in the range 8–1000 μm rest frame, we deduced the SFR_IR following Kennicutt (1998a) (see Eq. (6)).

The left panel of Fig. 5 displays our SFR for each bin of redshift and M_⋆. The SFR follows a monotonic dependence with increasing M_⋆, at fixed redshift; and with increasing redshift, at fixed M_⋆. The high-mass end presents a bending of the main sequence slope, which is more prominent as the redshift decreases.

Fig. 5. SFRMS as a function of M⋆ over different redshift bins. The dots represent data from this work. The upper limits were used to perform the fit but are not shown here to avoid overloading the figure. Left panel: fitted by the Eq. (12), the shaded area the 68% uncertainty of the fit. Right panel: fitted by the Eq. (13), the shaded area the 68% uncertainty of the fit. The squares represent the best-fit parameter M0. No bending was detected for z >  3.1 and therefore no square error bars are displayed.

Regarding the contribution of H-dropouts to our main sequence estimate, Wang et al. (2019) reported that ALMA-detected H-dropouts mostly fall within the main sequence at z ∼ 4 (the version from S15). Moreover, as they represent only ∼7.1% of the galaxies at equivalent M_⋆, their omission should not alter the overall shape of the main sequence.

The SFR − M_⋆ correlation was then fitted using the formula introduced in S15 as it is able to capture the bending of the main sequence at the high-mass end (see Eq. (12)),

$$\begin{array}{cc}\hfill {log}_{10}(SF{R}_{\mathrm{MS}}& [M_{\odot }{\mathrm{yr}}^{-1}])=m-m_0+a_{0}r\hfill \\ \hfill & -a_1{[max(0,m-m_1-a_{2}r)]}^2,\hfill \end{array}$$(12)

where r ≡ log₁₀(1 + z) and m ≡ log₁₀(M_⋆/10⁹ M_⊙). It was shown in S15 that the main sequence has the shape in logarithmic space of a Gaussian distribution with a homogeneous scatter of σ ≃ 0.3 dex. Because we stacked our galaxies, via mean stacking, we actually recover the ⟨L_IR⟩ of our sample which is different from the mean of a Gaussian distribution in logarithmic space with a dispersion σ_dist, with, ⟨10^X⟩=exp((σ_X × ln(10))²/2) × 10^⟨X⟩, X = log₁₀(SFR) following a Gaussian distribution. To correct for this, we assume that the dispersion of this Gaussian distribution in logarithmic space is σ_dist = 0.3 dex (S15). Our best-fit parameters are given in Table 5, and displayed on the left panel of Fig. 5 (left panel).

The shape chosen for the fit sets a slope of one, which is broadly consistent with our data points. The normalisation of the main sequence increases significantly with redshift. As we have already pointed out in the study of T_dust, at the highest redshift (i.e. 3.9 ≤ z ≤ 5.0) and the two highest M_⋆ bins (i.e. 10.5 ≤ log₁₀(M_⋆)≤11 and 11 ≤ log₁₀(M_⋆)≤12), appear to be highly star-forming and are above our main sequence by a factor of two. This suggests that these bins include a significant fraction of galaxies with high SFR relative to the main sequence (i.e. high R_SB).

We compare our results with several versions of the main sequence from the literature in Fig. 6. Compared to Speagle et al. (2014), who did not fit any bending, our best-fit is consistent with the general trend and evolution of the normalisation. However, our work suggests that the bending is real and becomes stronger at low redshift. Overall, our best-fit remains close to what was found by S15 at z ≥ 0.7. As already confirmed by Delvecchio et al. (2021), the trend of S15 main sequence holds when constructed mainly from FIR data. Nevertheless, we find a stronger curvature of the main sequence at z ≤ 0.7 compared to what we could extrapolate from S15. This is probably due in part to the fact that they did not probe the main sequence for z ≤ 0.3. The version of the main sequence by Leslie et al. (2020) was derived from a study of the radio continuum at 3 GHz, which may explain the differences in slope, normalisation and bending. But both roughly agree on the same evolutionary trend and on a bending of the main sequence with marginal differences. The bending of the main sequence being at the high-mass end, it suggests that it is triggered by some mass driven physical processes.

Fig. 6. SFRMS as a function of M⋆. The red triangles represent upper limits, the blue dots are the observed values from this work, and shaded area the 68% uncertainty of the fit. The blue lines represent the best-fit from this work. The green, cyan and orange lines represents the main sequence from S15, Speagle et al. (2014) and Leslie et al. (2020), respectively. The dashed lines represent extrapolated main sequences to redshifts that were not investigated in their respective studies.

We have also followed the evolution of the stellar mass knee M₀ marking the bending point of the main sequence by following Daddi et al. (2022) (see Eq. (13)):

$$\begin{array}{c}\hfill \frac{\mathit{SFR}}{SF{R}_0}=\frac{1}{1+{(M_0/M_{\star })}^{\gamma }}.\end{array}$$(13)

Following recommendations from Daddi et al. (2022), we set γ = 1.1, which should help to reduce the errors and should not affect the result much. If γ was defined as a free parameter in the fit, the ⟨γ⟩ would be close to 1.1 anyway. Our best-fit parameters are given in Table 6, and displayed on the right panel of Fig. 5. Overall, our results are quite comparable to the evolution of the bending found by Lee et al. (2015) and Daddi et al. (2022) (i.e. a decrease in M₀ as we move to a lower redshift). However, we see little or no evidence of a bending for z >  3. We suggest that this may be because the physical processes that trigger the main sequence bending have not had enough time to impact the main sequence trend at the high-mass end at z >  3. Another possibility is that the bending still occurs at high redshifts, but only for very high-mass galaxies (M_⋆ >  2 − 3 × 10¹¹ M_⊙), which we could not probe with this study because they are extremely rare.

5.3. Gas mass

We chose to calculate M_gas from M_dust deduced from our stacked SEDs. We can link M_gas to M_dust by a gas-to-dust mass ratio (δ_GDR) that depends only on the metallicity Z (logZ = 12 + log₁₀(O/H)) of the galaxy,

$$\begin{array}{c}\hfill M_{\mathrm{gas}}={\delta }_{\mathrm{GDR}}(Z)M_{\mathrm{dust}}.\end{array}$$(14)

Multiples studies (e.g. Leroy et al. 2011; Magdis et al. 2012; Rémy-Ruyer et al. 2014; Genzel et al. 2015) show consistent δ_GDR − Z relations that hold at both low and high redshift. We chose here to use the relation from Magdis et al. (2012):

$$\begin{array}{c}\hfill log{\delta }_{\mathrm{GDR}}=10.54-0.99\times (12+{log}_{10}(\mathrm{O}/\mathrm{H})).\end{array}$$(15)

As for most galaxies we do not have UV or optical spectra, we do not have a direct measurement of their metallicity. To go around this problem, we use the mass-metallicity relation (MZR; Erb et al. 2006) and follow its redshift dependent version from Genzel et al. (2015),

$$\begin{array}{c}\hfill 12+{log}_{10}(\mathrm{O}/\mathrm{H})=a-0.087\times {({log}_{10}(M_{\star }/1.7)-b)}^2,\end{array}$$(16)

where a = 8.74 and b = 10.4 + 4.46log₁₀(1 + z)−1.78(log₁₀(1 + z))². We note that Genzel et al. (2015) use a Chabrier (2003) IMF to define the MZR. In order to use it correctly, our M_⋆ must first be divided by a factor of 1.7. We chose to adopt an uncertainty of 0.2 dex for our metallicity following the recommendation from Magdis et al. (2012). These latter metallicities are calibrated in the PP04 N2 scale ($N_2=[{\mathrm{N}}_{\mathrm{II}}]$$\lambda 6583/{\mathrm{H}}_{\alpha }$; Pettini & Pagel 2004). We note that by using the δ_GDR relation of Leroy et al. (2011), instead of Magdis et al. (2012), would not change our conclusions because for our sample, the relative median difference in linear scale between the two estimates is $(M_{\mathrm{gas}}^{\mathrm{MZR}-\mathrm{Leroy}11}-M_{\mathrm{gas}}^{\mathrm{MZR}-\mathrm{Magdis}12})/M_{\mathrm{gas}}^{\mathrm{MZR}-\mathrm{Magdis}12}=0.{14}_{-0.08}^{+0.05}$.

Another way to calculate the metallicity is the fundamental metallicity relation (FMR; Mannucci et al. 2010, see Eq. (17)). It differs from the MZR by adding some dependence on the SFR,

$$\begin{array}{cc}\hfill 12+{log}_{10}& (\mathrm{O}/\mathrm{H})=8.90+0.37\mathrm{m}-0.14\mathrm{s}\hfill \\ \hfill & -0.{19}^2+0.12\mathrm{ms}-0.054{\mathrm{s}}^2,\hfill \end{array}$$(17)

where m = log₁₀(M_⋆/1.7)−10 and s = log₁₀(SFR/1.7). These metallicities are calibrated for the KD02 photoionisation models (Kewley & Dopita 2002). We used the recipe from Kewley & Ellison (2008) to convert it into a PP04 N2 scale. But once again the differences in metallicities between the MZR and FMR method does not impact our conclusions, as for our sample, the relative median is $(M_{\mathrm{gas}}^{\mathrm{FMR}-\mathrm{Magdis}12}-M_{\mathrm{gas}}^{\mathrm{MZR}-\mathrm{Magdis}12})/M_{\mathrm{gas}}^{\mathrm{MZR}-\mathrm{Magdis}12}=-0.{16}_{-0.06}^{+0.16}$.

The previous methods are only reliable when the SED peak is well defined (i.e. T_dust and M_dust), which is not the case in our bins. The only two bins concerned are 10.0 ≤ log₁₀(M_⋆)≤10.5 for 3.1 ≤ z ≤ 3.9 and 3.9 ≤ z ≤ 5.0. In this situation, we can also estimate M_gas from a single band measurement located in the Rayleigh-Jeans part of the SED (e.g. Scoville et al. 2014; Groves et al. 2015; Schinnerer et al. 2016). The main limitation is that it does not take into account the evolution with redshift of the M_gas metallicity and T_dust (e.g. Genzel et al. 2015; Berta et al. 2016; Schinnerer et al. 2016; Magdis et al. 2017; Harrington et al. 2021). However, this is a reliable way of estimating M_gas at low cost (25% uncertainties; Scoville et al. 2016) and especially when there is only one band measured in the Rayleigh-Jeans tail. For these reasons, while we favoured ${\delta }_{\mathrm{GDR}}-\mathrm{MRZ}$ methods to deduce M_gas as long as the T_dust was well defined. In the opposite case, and when an ALMA measurement is available, we calculated M_gas by the method described in Scoville et al. (2016) (see their Eqs. (16) and (6)), that is

$$\begin{array}{cc}\hfill M_{\mathrm{gas}}& =1.78S_{v_{\mathrm{obs}}}[\mathrm{m}\mathrm{Jy}]{(1+z)}^{-4.8}{\left(\frac{{\nu }_{850\mathrm{\mu }\mathrm{m}}}{{\nu }_{\mathrm{obs}}}\right)}^{3.8}{(d_{\mathrm{L}}[\mathrm{G}pc])}^2\hfill \\ \hfill & \times \frac{6.7\times {10}^{19}}{{\alpha }_{850}}\frac{{\mathrm{\Gamma }}_0}{{\mathrm{\Gamma }}_{\mathrm{RJ}}}{10}^{10}{M}_{\odot }\mathrm{for}{\lambda }_{\mathrm{rest}}>250\mathrm{\mu }\mathrm{m},\hfill \end{array}$$(18)

$$\begin{array}{c}\hfill \mathrm{\Gamma }(z,T_{\mathrm{d}},\nu )=\frac{h{v}_{\mathrm{obs}}(1+z)}{k_{\mathrm{B}}{T}_{\mathrm{d}}}\frac{1}{(exp(\frac{h{v}_{\mathrm{obs}}(1+z)}{k_{\mathrm{B}}{T}_{\mathrm{d}}})-1)},\end{array}$$(19)

where d_L is the luminosity distance, ${\alpha }_{850}=6.2\times {10}^{19}\mathrm{erg}{\sec }^{-1}{\mathrm{Hz}}^{-1}{\mathrm{M}}_{\odot }^{-1}$, ${\mathrm{\Gamma }}_0=\mathrm{\Gamma }(z=0,T_{\mathrm{d}}=25\mathrm{K},\nu ={\nu }_{850\mathrm{\mu }\mathrm{m}})$ and ${\mathrm{\Gamma }}_{\mathrm{RJ}}=\mathrm{\Gamma }(z,T_{\mathrm{d}}=25\mathrm{K},\nu ={\nu }_{1.1\mathrm{mm}})$, h is Planck’s constant, k_B is Boltzmann’s constant, $T_{\mathrm{d}}=25\mathrm{K}$ and ${\nu }_{\lambda }=\frac{c}{\lambda }$. We note that this technique provides, for our sample, consistent results with those of the MZR method, with a relative median $(M_{\mathrm{gas}}^{\mathrm{Scoville}16}-M_{\mathrm{gas}}^{\mathrm{MZR}-\mathrm{Magdis}12})/M_{\mathrm{gas}}^{\mathrm{MZR}-\mathrm{Magdis}12}=0.{19}_{-0.32}^{+0.35}$ when both ${\delta }_{\mathrm{GDR}}-\mathrm{MRZ}$ and Scoville et al. (2016) could be performed.

In summary, we chose to calculate our M_gas using Magdis et al. (2012) for δ_GDR(Z) (Eq. (15)) and Genzel et al. (2015) MZR for the metallicities (Eq. (16)). When the peak of the SED is poorly constrained and an ALMA flux is available, we chose to calculate M_gas from Scoville et al. (2016) method. For the rest of the paper, we will refer to $M_{\mathrm{gas}}^{\mathrm{MZR}-\mathrm{Magdis}12}$, or $M_{\mathrm{gas}}^{\mathrm{Scoville}16}$ simply as M_gas.

Here M_gas represents the total gas budget of a galaxy including both molecular gas (M_H2) and atomic gas (M_HI), that is, M_gas = M_H2 + M_HI. The HI content of galaxies is still poorly known outside the local universe because the 21cm emission line is difficult to detect with current facilities. However, Bauermeister et al. (2010) showed that it is unlikely that the HI content of galaxies varies strongly with redshift. In contrast, the ${\mathrm{H}}_2$ content evolves strongly with redshift (e.g. Daddi et al. 2010; Tacconi et al. 2010, 2013; Lagos et al. 2015; Genzel et al. 2015; Tacconi et al. 2018). In Tacconi et al. (2018), it was concluded that the assumption M_gas ∼ M_H2 should hold for z >  0.4. Therefore, we decided not to consider our M_gas for z <  0.4 (i.e. our first redshift bin), as we could not probe the HI content of the galaxy at these redshifts. In order to add some reliable measurements of M_H2 to our study at z <  0.4, we used data points from Saintonge et al. (2017) that come from very local measurements (0.01 ≤ z ≤ 0.05).

Our M_gas measurements are displayed in Fig. 7. We observe a rapid rise with redshift from Saintonge et al. (2017) data at low redshift, to our first data points at z ≥ 0.4. Then, at fixed M_⋆, M_gas reach a maximum at z ∼ 1 − 2 and remains relatively constant as the redshift increases. At fixed redshift, M_gas gradually increase with M_⋆. All values of our gas masses are provided in Appendix A.

We examined the contribution of the H-dropouts to M_gas in order to have the most unbiased view possible. Here, we have simply calculated the M_gas, associated with each H-dropout in the sample of Wang et al. (2019), using their ALMA measurement at 870 μm and following Scoville et al. (2016) method. The final contribution of the H-dropout to the total M_gas within each bin is $M_{\mathrm{gas}}^{H\text{-dropout}}\times N_{\mathrm{gal}}^{H\text{-dropout}}\times ϵ/(M_{\mathrm{gas}}^{H\text{-dropout}}\times N_{\mathrm{gal}}^{H\text{-dropout}}\times ϵ+M_{\mathrm{gas}}\times N_{\mathrm{gal}})=19{\%}_{-17\%}^{+7\%}$ on average (only the bins with at least one H-dropout are taken into account), where ϵ = Area_This work/Area_{Wang + 19} ∼ 1.8. The contribution of H-dropouts can represent up to ∼32% for some bins, hence we chose to add it to our data. The seven bins that have been corrected for the H-dropout contribution are circled Fig. 7.

Fig. 7. Mgas as a function of redshift for different M⋆. The dots represent Mgas estimates from this work using the ${\delta }_{\mathrm{GDR}}-\mathrm{MRZ}$ method. the crosses represent Mgas estimates from this work following the method from Scoville et al. (2016). The squares are measurements from Saintonge et al. (2017). The bins that have been corrected for the H-dropout contribution are circled. The thick and thin solid colour lines represent the best-fit from this work, using Eq. (21) and Eq. (20), respectively. The dashed lines are the fits extrapolation. The shaded area the 68% uncertainty of the fit.

Tacconi et al. (2018) combined M_gas from stacks in the IR and M_gas from CO emission. We try to fit our data using the formula of Tacconi et al. (2018):

$$\begin{array}{cc}\hfill {log}_{10}& (M_{\mathrm{gas}})=A+B\times {({log}_{10}(1+z)-F)}^{\beta }\hfill \\ \hfill & +D\times ({log}_{10}(M_{\star }/1.7)-10.7)+{log}_{10}(M_{\star }/1.7)\hfill \\ & [+C\times {log}_{10}(R_{\mathrm{SB}}(z,M_{\star }))+E\times {log}_{10}(R_{\mathrm{e}}/R_{\mathrm{e}0}(z,M_{\star }))],\hfill \end{array}$$(20)

where the factor 1.7 represents the conversion of M_⋆ from a Salpeter (1955) to Chabrier (2003) IMF. As our data points represent main sequence galaxies, we cannot probe the C_μ term in Tacconi et al. (2018), that represents the evolution with respect to the distance to the main sequence sSFR/sSFR_MS, where sSFR is the specific SFR (sSFR = SFR/M_⋆). And we do not explore the morphology of galaxies in this study, which makes it impossible to probe the term E_μ, which compares the effective radius of galaxies (R_e) to the mean effective radius of the star forming population R_e0. The results are displayed in Fig. 7, and our best-fit parameters are given in Table 7.

Looking at Fig. 7, we notice that M_gas observed at low redshift (0.4 ≤ z ≤ 1) tends to be higher than the best-fit trend (even though it is in most cases within the error bars), and the trends from literature such as Tacconi et al. (2018). A similar effect has been observed by Tacconi et al. (2018) who measured higher depletion times (τ_dep) deduced from dust observations compared to the ones deduced from CO line fluxes. The effect observed in Tacconi et al. (2018) decreases slowly with redshift (about 0.3 dex at z ∼ 0.4 to 0 at z ∼ 1.5) and roughly matches the effect observed here. We conclude that this effect, which was solved in Tacconi et al. (2018) by matching the zero-point for each method, arises when deducing M_gas from dust observations. It feels that the form use in Tacconi et al. (2018) to fit M_gas fails to properly recover the form of our measurements. In a way to provide for a better fit, we chose to also fit our data using a formula of the form:

$$\begin{array}{c}\hfill {log}_{10}(M_{\mathrm{gas}})=\{\begin{array}{cc}{m}_0+m_{1}m+a_{1}r+a_2{r}^2\hfill & \text{for}z<z_0\hfill \\ m_0+m_{1}m+a_1{r}_0+a_2{r}_0^2\hfill & \text{for}z\ge z_0,\hfill \end{array}\end{array}$$(21)

where m = log₁₀(M_⋆), r = log₁₀(1 + z), r₀ = log₁₀(1 + z₀) and z₀ = −a₁/2a₂. The best-fit parameters, using Eq. (21) are given in Table 8.

We compare our results to the trends observed by Tacconi et al. (2018) and Wang et al. (2022) in Fig. 8. In Fig. 8, the trends are displayed for stellar mass bins in the Salpeter (1955) IMF. Differences between our results and previous works (e.g. Scoville et al. 2016, 2017; Tacconi et al. 2018; Liu et al. 2019; Wang et al. 2022) may result from selection effects and from the method used to derive the gas masses. The present sample was built using an H-band selection that is the best proxy for the stellar mass, although not perfect for the highest redshifts. Other studies generally rely on a CO-line and/or dust continuum selection. This may lead to a potential bias favouring galaxies with higher specific star formation rate, that is, higher sSFR = SFR/M_⋆, hence galaxies above the main sequence. It has been found that such galaxies in general exhibited higher gas fractions. In our analysis, we tried to avoid such bias by checking that our stacked samples followed the median of the SFR − M_⋆ main sequence.

Fig. 8. Mgas as a function of redshift and M⋆. The solid lines represent the best-fit from this work colour-coded by M⋆ (using Eq. (21)), the dashed lines represent the trend from Tacconi et al. (2018), the dotted lines are from Wang et al. (2022). The faded lines represent an extrapolation from their respective laws.

The trends observed by Wang et al. (2022) were obtained as well using a mass-selected sample, but the gas masses were derived from millimetre continuum fluxes using the method from Scoville et al. (2016). We produced a version of our analysis using the same method for consistency check, it is presented in Appendix B. We find that we are able to reproduce the same trends as Wang et al. (2022) if we were to use the same technique based only on the formula of Scoville et al. (2016). Hence, the differences found in the present study can be explained by the use of a different method to derive gas masses. We favoured the present method because it takes into account the variation of metallicity and the associated change in dust-to-gas ratio as a function of stellar mass. This is one of the reasons why we opted for a mass-selection in the first place.

As mentioned in Sect. 4, we chose to work with M_dust from amorphous silicate and graphitic grains (Draine & Li 2007) instead of amorphous carbon (Galliano et al. 2011; S18). We investigated what the impact would have been on M_gas if we had chosen a model based on amorphous carbon such as Galliano et al. (2021). Although the results are quite comparable, using the model from Galliano et al. (2021), would translate into slightly lower M_gas (∼ − 10% at z = 4). These differences would not have changed the conclusions drawn in this paper.

6. Cosmic star formation history

In this section, we present the redshift evolution of ρ_SFR, that is, the cosmic star formation history. To calculate the ρ_SFR, we start from the stellar mass function of SFGs of Davidzon et al. (2017), that is given for different redshift bins. For each redshift bin, we generated galaxies with uniform redshift distribution within the bin, and a M_⋆ distribution following these stellar mass functions. From the redshift and M_⋆, we then assigned a SFR using the same method as in Sect. 4.3. In summary, SFR = R_SB × SFR_MS, where SFR_MS is calculated from the main sequence evolution found in this work (see Eq. (12) and parameters Table 5). And a R_SB is randomly drawn from a double Gaussian distribution representing the position of normal and starbursting galaxies relative to the main sequence (Eq. (3)). Then ρ_SFR is calculated by summing the SFR of galaxies down to $M_{\star }^{\lim }=3\times {10}^9{M}_{\odot }$. Errors were generated by varying the SFR_MS trend of Eq. (12) within the errors of the fit, a 100 times.

The study of Wang et al. (2019) presents the contribution to ρ_SFR of H-dropout galaxies, which can reach up to 10% at z ∼ 4 − 5. We have therefore added the contribution of H-dropout galaxies of Wang et al. (2019) to our ρ_SFR considering the sample of H-dropout galaxies have a median stellar mass of M_⋆ ∼ 10^10.6 M_⊙ (Wang et al. 2019).

We compare the ρ_SFR evolution with that of Madau & Dickinson (2014) who integrate UV and IR luminosity functions down to L_min = 0.03L_*. Deducing ρ_SFR by integrating luminosity functions down to 0.03L_*, or the mass function down to 3 × 10⁹ M_⊙ should yield, to the first order, similar results (Schreiber et al. 2015). We chose for the rest of the paper to derive our ρ_SFR by integrating the mass function down to 3 × 10⁹ M_⊙. Thus, all conclusions on the evolution of the ρ_SFR in this paper should be understood in this framework, as integrating the luminosity function down to lower luminosities, or the mass function down to lower masses, must yield a higher total ρ_SFR.

The total ρ_SFR is displayed in Fig. 9 along with some examples from the literature, and the data points from this work are summarised in Table 9. We observe a rise of ρ_SFR from z ∼ 5 to z ∼ 2, and then it gradually decreases down to z ∼ 0.35. The high masses (i.e. log₁₀(M_⋆/M_⊙)> 10) account for most of the ρ_SFR until z ∼ 4. On the other hand the ρ_SFR associated to low-mass galaxies (i.e. log₁₀(M_⋆/M_⊙)< 10) is roughly constant over 0 ≤ z ≤ 5.

Fig. 9. Cosmic star formation rate density (ρSFR) as a function of redshift. The green line represents the total ρSFR trend from this work (i.e. integrating the stellar mass function down to 3 × 109 M⊙). The red and blue lines represent the contribution of galaxies with log10(M⋆/M⊙)> 10 and log10(M⋆/M⊙)< 10, respectively. The black dotted line represents the ρSFR detected with stacking. The purple squares show the contribution of H-dropout from Wang et al. (2019). The dotted green and red lines show the trend once H-dropout contribution from Wang et al. (2019) has been added to the respective coloured solid lines. The grey line represents the trend from Madau & Dickinson (2014). The cyan and orange squares represent ρSFR estimate from Leslie et al. (2020) and Gruppioni et al. (2020), respectively.

We can see that our estimate of the total ρ_SFR is close to what has been observed by Leslie et al. (2020) and to the evolution of Madau & Dickinson (2014). On the other hand, the measurements from Gruppioni et al. (2020) are mostly in disagreement with our results at high redshift (i.e. z ≥ 2). The ρ_SFR from Gruppioni et al. (2020) is deduced from the integration of the IR-luminosity function (which make it ${\rho }_{\mathrm{SFR}}^{\mathrm{IR}}$) down to 10⁸ L_⊙.

In Fig. 10, we show the contribution to the total ρ_SFR of the full range of stellar masses. We observe that the increase in the total ρ_SFR, from z ∼ 5 to z ∼ 2, comes from the growing number of massive galaxies (i.e. log₁₀(M_⋆/M_⊙)> 10), which can be seen in the evolution of the stellar mass function at these redshifts (Davidzon et al. 2017). Downsizing and the bending of the main sequence explain the fall of the contribution of massive galaxies from z ∼ 2 to z ∼ 0. We again observe that massive galaxies (i.e. log₁₀(M_⋆/M_⊙)> 10) dominate the total ρ_SFR at all redshifts (i.e. they account for more than 70%). In particular, it appears that galaxies with 10.5 ≤ log₁₀(M_⋆/M_⊙)≤11.25 account for more than ∼55% of the total ρ_SFR at z = 2, making them the main driver of the peak in the observed cosmic star formation history at this redshift.

Fig. 10. Contributions over the whole range of stellar masses to the total ρSFR as a function of redshift for SFGs.

7. Cosmic evolution of the gas mass density

In Sect. 6, we have simulated a catalogue of galaxies where for each of them, we calculated M_⋆ and R_SB. From these properties, we now infer M_gas by expanding Eq. (21) as follows:

$$\begin{array}{c}\hfill {log}_{10}(M_{\mathrm{gas}})=\{\begin{array}{c}{m}_0+m_{1}m+a_{1}r+a_2{r}^2\hfill \\ +C\times {log}_{10}(R_{\mathrm{SB}}(z,M_{\star }))\hfill & \text{for}z<z_0\hfill \\ m_0+m_{1}m+a_1{r}_0+a_2{r}_0^2\hfill \\ +C\times {log}_{10}(R_{\mathrm{SB}}(z,M_{\star }))\hfill & \text{for}z\ge z_0,\hfill \end{array}\end{array}$$(22)

where m = log₁₀(M_⋆), r = log₁₀(1 + z), r₀ = log₁₀(1 + z₀) and z₀ = −a₁/2a₂. Here, m₀, m₁, a₁ and a₂ are taken from Table 8, and C = 0.53 is taken from Tacconi et al. (2018). Next, the cosmic evolution of gas density (ρ_gas) is calculated by summing the M_gas of galaxies down to $M_{\star }^{\lim }=3\times {10}^9{M}_{\odot }$. Errors were generated by varying our best-fit of $M_{\mathrm{gas}}^{\mathrm{MS}}$ within its errors, one hundred times. Defined this way, ρ_gas directly represents the gas content of galaxies contributing to the ρ_SFR presented in Fig. 9. We display the evolution of ρ_gas as a function of redshift in Fig. 11. We can see the same kind of features as for the evolution of ρ_SFR: a rise and fall with redshift with a maximum around z ∼ 2; a dominance across all redshifts of the high-mass galaxies (i.e. log₁₀(M_⋆/M_⊙)> 10) contribution; and a relatively flat evolution of the low-mass galaxies (i.e. log₁₀(M_⋆/M_⊙)< 10). This shows that the SFE of both high and low-mass galaxies are not drastically changing with time, and thus that the gas content of galaxies (i.e. the accretion) is the primary driver of their SFRs: high-mass galaxies have higher SFR at fixed redshift because they have more gas. However, comparing the relative evolution of ρ_gas (see Fig. 11) and ρ_SFR (see Fig. 9) shows that ${\rho }_{\mathrm{gas}}^{z=5}/{\rho }_{\mathrm{gas}}^{z=0}$ (and in particular the high-mass contribution) is lower by a factor ∼3 − 4 compared to ${\rho }_{\mathrm{SFR}}^{z=5}/{\rho }_{\mathrm{SFR}}^{z=0}$. This hints that one unit of gas leads to more stars being formed at z ∼ 5 compared to z ∼ 0 (i.e. a higher SFE at z ∼ 5 compared to z ∼ 0).

Fig. 11. ρgas as a function of redshift. The green line represents the total ρgas inferred from this work (i.e. integrating the stellar mass function down to 3 × 109 M⊙). The red and blue lines represent the contribution to the total ρgas of galaxies with log10(M⋆)> 10 and log10(M⋆)< 10, respectively. The purple and cyan lines represent ρgas deduced from the Mgas definition of Wang et al. (2022) and Tacconi et al. (2018), respectively.

We substituted, in our later method to estimate ρ_gas, M_gas from this work, for best-fit of M_gas from the literature and compared it to our estimate of ρ_gas. We can see in Fig. 11, that the redshift evolution of ρ_gas from Tacconi et al. (2018) or Wang et al. (2022) are higher compared to the one from this work. Theses are simply resulting from discrepancies already observed in the respective M_gas trend they were built from, and goes together with a stronger contribution of low-mass galaxies for which the gas mass discrepancy discussed in Fig. 8 is the largest at high redshifts (i.e. z >  3). The latter part is a potential weakness of all studies including ours, since the low stellar mass – high redshift of the cosmic gas mass density range relies on rather uncertain extrapolations. As already pointed out in Liu et al. (2019), the form of the formula chosen to fit M_gas can have a significant impact on the resulting ρ_gas trend, and could be the cause of what is observed here.

We also computed ρ_gas by summing M_gas of all galaxies down to several different $M_{\lim }^{\ast }$ ($M_{\lim }^{\ast }=1.7\times {10}^8$, 1.7 × 10⁹ and 1.7 × 10¹⁰ M_⊙). We compared it to Magnelli et al. (2020) where a ρ_gas was computed for the same $M_{\lim }^{\ast }$ from the stacking of H-band selected galaxies in ALMA and through the method from Scoville et al. (2016). The results are displayed in Fig. 12. Here, our work agrees quite well with that of Magnelli et al. (2020) for the total ρ_gas (i.e. $M_{\star }^{\lim }=1.7\times {10}^8$). However, the mass distribution is not similar: massive galaxies (i.e. $M_{\star }^{\lim }=1.7\times {10}^{10}$) contribute more at low redshift (i.e. z ≤ 0.7) in this work as compared to Magnelli et al. (2020); conversely, their contribution is smaller at higher redshift (lower by a factor ∼3 at z ∼ 2.8 compared to Magnelli et al. 2020).

Fig. 12. ρgas as a function of redshift. The green, orange and red lines represent ρgas derived from this work with $M_{\star }^{\lim }$ of 1.7 × 108, 1.7 × 109 and 1.7 × 1010 M⊙, respectively. The green, orange and red hatches represent estimates of ρgas from Magnelli et al. (2020) with $M_{\star }^{\lim }$ of 1.7 × 108, 1.7 × 109 and 1.7 × 1010 M⊙, respectively. The green dots with black circles depict the estimate of ρgas from Magnelli et al. (2020) for M⋆ >  Mlimit, Mlimit being the stellar mass completeness limit from Magnelli et al. (2020).

8. Discussion

Here, when we discuss the evolution of galaxies, it should be understood that we are discussing the evolution of main sequence galaxies, as all scaling relations of SFR, M_gas and T_dust with M_⋆ and redshift have been deduced for main sequence galaxies.

8.1. Cosmic densities of star formation and gas density as a function of stellar mass

An interesting feature of Fig. 9 is the fact that the contribution to ρ_SFR of galaxies with M_⋆ ≤ 10¹⁰ M_⊙ appears to be constant over 1 ≤ z ≤ 3 and slightly decreasing for z ≥ 3. This means that the balance of the number of these galaxies and their efficiency in producing stars remains constant with cosmic time. This feature can also be seen in Fig. 11, where the contribution to ρ_gas is also quite stable for galaxies with M_⋆ ≤ 10¹⁰ M_⊙.

On the other hand, galaxies with M_⋆ ≥ 10¹⁰ M_⊙ account for most stars formed up to z ∼ 5. This shows that it is galaxies with M_⋆ ≥ 10¹⁰ M_⊙ that are responsible for the observable shape on the ρ_SFR, especially at the cosmic noon at z ∼ 2. In Fig. 13a, we can see that the 10% (of the total number of galaxies contributing to ρ_SFR) most massive galaxies contribute to a large fraction of the ρ_SFR at all redshifts. This contribution goes from ∼24% at z ∼ 0 to ∼48% at z ∼ 5. In particular, this means that by considering only 10% (in number) of the most massive galaxies, we can deduce a relatively good estimation of the total ρ_SFR at intermediate and high redshifts. From these numbers, we can also recognise some hierarchical growth effect, as bins gradually rise in M_⋆ over time.

Fig. 13. Contribution to ρSFR of SFGs as a function of redshift. Catalogues were binned through different methods and a specific property is displayed through colours: see each sub-figure for specifics. (a) Each bin includes 10% of the galaxies in number, picking from the lowest M⋆ onwards. The minimum M⋆ in the bin defines the colour. The numbers give the median M⋆ within the bin along with the lower and maximum extension of the bin. (b) Each bin includes 10% of the total M⋆ of all galaxies, picking from the lowest M⋆ onwards. The median τdep in the bin defines colour. (c) Each bin includes 10% of the total M⋆ of all galaxies, picking from the lowest M⋆ onwards. The minimum M⋆ in the bin defines the colour. Contours of log10(M⋆=M⊙) = 10 and log10(M⋆=M⊙) = 11 are added as black dashed lines.

We find that H-dropout galaxies account for 13% of the total cosmic SFR density (CSFD) at z = 4–6. This is in agreement with the value of 10% found in Wang et al. (2019). We note that H-dropout galaxies account for ∼23% of the stars formed in massive galaxies (i.e. M_⋆ ≥ 10^10.3 M_⊙). This is a large contribution but not a dominant one, as opposed to the claim of Wang et al. (2019) that H-dropouts dominate the SFR density in the most massive galaxies. This discrepancy results from our use of a fully complete sample of galaxies selected in stellar mass, whereas Wang et al. (2019) compared H-dropouts to UV-selected and similarly massive Lyman-break (LBG) galaxies.

We observe a decline in star formation in massive galaxies that mirrors that of the total star formation density. This is an illustration of the bending effect of the MS at high mass (see Fig. 13a). The use of the depletion time in Fig. 13b illustrates this clearly: the depletion time of massive galaxies varies from τ_dep ≲ 200 Myr at z >  4 to τ_dep ≳ 600 Myr at z <  1.

In Fig. 13c, we show that the contribution of the most massive galaxies (e.g. M_⋆ ∼ 10¹¹ M_⊙) to the CSFD is nearly flat (i.e. their contribution to the CSFD is roughly constant), meaning that they tightly follow the global history of cosmic star formation (i.e. the global shape of the CSFD consisting of a rise, a cosmic noon and a fall).

Let us consider a typical galaxy with M_⋆ = M_star, such that 50% of ρ_⋆ (cosmic stellar mass density) is made of galaxies above and below M_star. We can follow the contribution of such typical galaxy to the cosmic SFR density by following the line where the Y-axis on the right of Fig. 13c equals 50%. We find that galaxies with M_⋆ >  M_star contribute to ∼50% of ρ_SFR at z >  3 and drops to only ∼25% at z ∼ 0. This decrease of the relative contribution of massive galaxies to ρ_⋆ illustrates the impact on the cosmic SFR history of the bending of the MS. The fact that this happens continuously with cosmic time, supports a scenario of a slow downfall of star formation rather than a rapid quenching of the most massive galaxies.

In Fig. 14, we track the evolution of the stellar mass above which galaxies contribute to exactly 50% of the CSFD. We can see in Fig. 14 that this particular stellar mass falls around ∼10^10.2 M_⊙ with an evolution with redshift that mirror the one of the CSFD by peaking at the cosmic noon at ∼10^10.4 M_⊙.

Fig. 14. Stellar mass above which galaxies contribute exactly 50% of the CSFD as a function of redshift.

8.2. Comparison with the TNG100 simulation

Here we compare our results with cosmological simulations. In particular, we have investigated whether current simulations are able to reproduce the contribution of the different stellar mass bins to the total ρ_SFR. We chose to examine the TNG100 simulation (Nelson et al. 2019, 2018; Springel et al. 2018; Marinacci et al. 2018; Naiman et al. 2018; Pillepich et al. 2018b, 2018a; Weinberger et al. 2017) of the IllustrisTNG project.

We retrieved the ρ_SFR of TNG100 using the same method as the one with which we calculate the ρ_SFR in this work, that is, by integrating down to 3 × 10⁹ M_⊙. In order to compare the trends with our work, this was also done by integrating only galaxies with log₁₀(M_⋆/M_⊙)> 10 and log₁₀(M_⋆/M_⊙)< 10. We display the evolution of ρ_SFR with redshift retrieved from the TNG100 simulation in Fig. 15. By comparing the trends of the TNG100 simulation with our results, we can see clear discrepancies as also noticed and discussed in Pillepich et al. (2018a). The ρ_SFR of the low masses in TNG100 (i.e. log₁₀(M_⋆/M_⊙)< 10) is not nearly as flat as observed (especially at z ≤ 2). On the other hand, the ρ_SFR of the high masses in TNG100 (i.e. log₁₀(M_⋆/M_⊙)> 10) does not account for as large a part of the total ρ_SFR for z <  3. The high masses account for less than half of the total ρ_SFR at z ∼ 1.7 in the simulation compared to the observations. As a result, the ρ_SFR peak is reached too early (i.e. at z ∼ 3 instead of the observed at z ∼ 1.7) in the simulation. The disparities between the two total ρ_SFR trends can be almost exclusively associated with the contribution of massive galaxies being off in the TNG100 simulation.

Fig. 15. Cosmic star formation rate density (ρSFR) as a function of redshift. The green, red and blue lines represent ρSFR deduced from the TNG100 simulation integrated from 3 × 109 M⊙, including all galaxies, galaxies with log10(M⋆/M⊙)> 10, and galaxies with log10(M⋆/M⊙)< 10, respectively. The green, red and green faded dots represent the total ρSFR trend, ρSFR of galaxies with log10(M⋆/M⊙)> 10, and galaxies with log10(M⋆/M⊙)< 10, respectively (including H-dropout contribution from Wang et al. 2019), from this work. The grey line represents the trend from Madau & Dickinson (2014) for reference.

Another way to show the difference in behaviour between the TNG100 simulation and our observations for the two mass bins, is to look at the evolution over redshift of ρ_SFR^{M_⋆ >  1010}/ρ_SFR^{M_⋆ <  1010} (see Fig. 16). We can see that the ratio decreases continuously with increasing redshift in the simulation, while observations show a clear peak in the ratio around z ∼ 1.7 (i.e. the cosmic noon) where the amount of stars formed in high-mass galaxies (i.e. log₁₀(M_⋆/M_⊙)> 10) exceeds by a factor up to ∼5.6 the one from low-mass galaxies (i.e. log₁₀(M_⋆/M_⊙)< 10).

Fig. 16. Ratio of the ρSFR of high-mass (M⋆ >  1010 M⊙) over low-mass (M⋆ <  1010 M⊙) galaxies, as a function of redshift. The red line represents the ratio deduced from this work (i.e. from observations). The blue line represents the ratio deduced from the TNG100 simulation.

In addition to underestimating the contribution of massive galaxies to ρ_SFR, these simulations tend to recover, for SFGs, lower gas fractions (see Fig. 17) and SFR (see Fig. 18) compared to observations. In TNG100, f_gas <  15% for M_⋆ >  10¹⁰ M_⊙ at all redshifts (Fig. 17), as it was already addressed in Pillepich et al. (2019). Moreover, the SFR of massive individual galaxies is lower in the TNG100 simulation than in observations at, for example, z >  4 (see Fig. 18). Lewis et al. (2024) shows a lack of metal content in IllustrisTNG simulations at z = 1, that might be due to AGN feedback that removes large quantities of metal-rich gas from the centres of massive galaxies. The lower metallicities observed by Lewis et al. (2024) in IllustrisTNG simulations at z = 1 could be a residual memory of how AGN affected star formation at higher z, which matches the lack of star formation we observe in Fig. 18. This could indicate flaws in the way gas accretion, the efficiency of galaxies to form stars from their gas content, or the ejection of gas through galactic winds is treated in the simulations. This shows that the problems encountered in the simulations may be due to the feedback processes used to regulate star formation and the gas reservoir of galaxies. In particular, the impact of AGN feedback may be overestimated to regulate star formation and kill massive galaxies. It has been shown that AGN feedback does not directly expel the gas and quench star formation at large scales. High-resolution simulations (i.e. with maximum spatial resolution $\sim 5-100\mathrm{pc}$) show that the AGN feedback tends to have little or no effect on the dense gas inside the galactic disc, because most of the out-flowing winds escape perpendicular to the galactic disc (Gabor & Bournaud 2014). AGN-driven outflows (i.e. ejective feedback), if sustained, could only quench a galaxy after a long time scale, that is, more than one Gyr (Gabor & Bournaud 2014; Biernacki & Teyssier 2018). These high-resolution simulations favour preemptive feedback (i.e. cutting out inflows of gas into the disc), though AGN-driven winds, to quench a galaxy (Gabor & Bournaud 2014; DeGraf et al. 2017; Biernacki & Teyssier 2018). However, even strong AGN-driven winds (i.e. $\sim 2500\mathrm{km}/\mathrm{s}$) would only reduce star formation in the galaxy by a factor of 2 (DeGraf et al. 2017). The TNG100 simulation has a much lower spatial resolution (i.e. a softening length of $\sim 0.74\mathrm{kpc}$ at z = 1), and is therefore not able to resolve the AGN-feedback interactions correctly. The sub-grid model of AGN feedback used in the IllustrisTNG has already been advocated to be responsible for a lower submillimetre galaxies (SMGs) number counts compared to observations (Hayward et al. 2021). The first James Webb Space Telescope (JWST) results interestingly go in the same direction: the number of bright, possibly massive, galaxies is found to exceed predictions at high z (Finkelstein et al. 2022; Donnan et al. 2023; Mason et al. 2023). This discrepancy on ρ_SFR could be a starting point to the underestimation of the global ρ_SFR in TNG100 for (1 ≤ z ≤ 2) reported early by Donnari et al. (2019) or non-trivial issues in the simulated galaxy and halo populations found in the Illustris simulation (Nelson et al. 2015).

Fig. 17. fgas as a function of M⋆ for SFGs in TNG100. Coloured lines represent the evolution deduced from the best-fit of Mgas from this work (see Table 7).

Fig. 18. SFR as a function of M⋆ for SFGs in TNG100. Coloured lines represent the best-fit from this work (see Table 5).

8.3. The SFE and Kennicutt-Schmidt relation

Through this work, we have obtained a measurement of both SFR and M_gas of main sequence galaxies over various redshift and M_⋆ bins. We can thus deduce the corresponding SFE of main sequence galaxies (i.e. examine the SFR-M_gas plane) and its evolution as a function of redshift and M_⋆. We display the SFE deduced from this work in Fig. 19, where, bins of redshift were regrouped into 3 bins to empathise results discussed later on.

Fig. 19. SFE as a function of stellar mass. Squares linked with a dashed line represent the measured from this work, colour coded by bin of redshift. Dotted lines are the value SFE deduced from our best fit law of the SFR and Mgas main sequences.

The SFE we measured shows an increase with redshift, which is accentuated at the highest masses. To investigate whether the evolution of SFE can be explained by the gas content of galaxies at high redshift we look at the Kennicutt-Schmidt relation (Kennicutt 1998b). The Kennicutt-Schmidt relation illustrates similar properties by examining SFR and M_gas surface density (i.e. Σ_SFR and Σ_gas respectively):

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{SFR}}=SFR/(2\pi R_{\mathrm{e}}^2),\end{array}$$(23)

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{gas}}=M_{\mathrm{gas}}/(2\pi R_{\mathrm{e}}^2).\end{array}$$(24)

It also links them as:

$$\begin{array}{c}\hfill {log}_{10}({\mathrm{\Sigma }}_{\mathrm{SFR}})=A+N\times {log}_{10}({\mathrm{\Sigma }}_{\mathrm{gas}}).\end{array}$$(25)

We do not probe R_e in this study, thus, we chose to use R_e(z, M_⋆), for galaxies with M_⋆ ≥ 10¹⁰, from Wang et al. (2022). We divided galaxies into bins of redshift (0 ≤ z ≤ 6) and M_⋆ (10 ≤ log₁₀(M_⋆)≤12). When fitting all bins, we added the sample with low redshift and low M_⋆ from Kennicutt (1998b). The error on the fit was generated by randomly varying the bin values within the uncertainties of the corresponding bin. We display the Schmidt-Kennicutt relation deduced from this work in Fig. 20.

Fig. 20. Kennicutt-Schmidt relation: ΣSFR as a function of Σgas. The dots represent the data from this work, and are colour coded by redshift. The cyan squared error bar represents the median of the sample distribution from Kennicutt (1998b), with individual galaxies displayed as grey crosses. The dashed black line represents the best-fit of all points. Coloured lines represent the best-fit, by fixing the slope to the all points best-fit value (i.e. $N=1.{18}_{-0.01}^{+0.01}$), at fix redshift.

All studies agree on a linear dependence of log₁₀(Σ_SFR) as a function of log₁₀(Σ_gas) with an overall slope ranging from 1 to 1.5 (e.g. Kennicutt 1998b; Wang et al. 2022), with no evolution of the slope or normalisation with redshift or M_⋆. Our best-fit, over all redshifts and M_⋆ bins, also suggests a power-law scaling with a slope of $N=1.{18}_{-0.01}^{+0.01}$.

We found that Σ_SFR = A × (Σ_SFR)^N with N ∼ 1.18 which means that an increase in M_gas induces some increase in SFE (i.e. $\mathrm{SFE}\sim {({\mathrm{M}}_{\mathrm{gas}})}^{0.18}$). It would not be the case if the slope of the Kennicutt-Schmidt relation was N = 1 as in Tacconi et al. (2013). We can correct Fig. 19 from this effect by plotting SFE/$M_{\mathrm{gas}}^{(N-1)}$ instead of SFE (see Fig. 21).

Fig. 21. SFE/$M_{\mathrm{gas}}^{(N-1)}$ (N = 1.18) as a function of stellar mass. Squares linked with a dashed line represent the measured from this work, colour coded by bin of redshift.

At z ≤ 2, in Fig. 21, the normalisation does not evolve much with redshift. On the other hand, in Fig. 7, M_gas increases with redshift and M_⋆ at z ≤ 2. This means that an increase in SFR as a function of redshift or stellar mass, at z ≤ 2, is mainly due to the variation of M_gas.

At z >  2, in Fig. 21, the normalisation increases with redshift: by a factor ∼3 between z ∼ 2 and z ∼ 4.3. On the other hand, in Fig. 7, M_gas remains roughly constant for z >  2. This means that the SFE plays a non-negligible role in increasing the normalisation of the main sequence SFR between 2 ≤ z ≤ 4.3.

In Fig. 21, we can see that the increase of SFE seen in Fig. 19 cannot be explained by the amount of gas as it is still appearing in Fig. 21.

Overall, it means that, we need to differentiate two effects: an increase of SFE due to the increase of M_gas, and an intrinsic increase in SFE. For example, at fixed stellar mass for 10.5 ≤ log₁₀(M_⋆/M_⊙)≤11 between z = 1.15 and z = 2, SFR increases by a factor ∼4 (see Fig. 5), while M_gas increase by a factor ∼1.6 (see Fig. 7): it means that SFE increases by a factor ∼2.5, which reflects the increase of the normalisation in Fig. 19 by ∼2.5 for the same redshift range. The first effect in an increase of SFE due to the increase of M_gas (i.e. resulting from the Kennicutt-Schmidt relation of slope N ∼ 1.18): ${\alpha }_{\mathrm{SFE}}^{\text{due to}}{M}_{\mathrm{gas}}\sim {\alpha }_{M_{\mathrm{gas}}}^{N_1}={(1.6)}^{0.18}\sim 1.09$. The second effect is an intrinsic increase in SFE between the two redshifts (i.e. independent of the increase of M_gas): ${\alpha }_{\mathrm{SFE}}^{\mathrm{intrinsic}}\sim {\alpha }_{\mathrm{SFE}}/{\alpha }_{\mathrm{SFE}}^{\text{due to}{M}_{\mathrm{gas}}}\sim 2.3$. Finally, we could conclude by saying that a factor ∼2.3 is solely due to an increase of SFE, and a factor $\sim {\alpha }_{M_{\mathrm{gas}}}\times {\alpha }_{\mathrm{SFE}}^{\text{due to}{M}_{\mathrm{gas}}}=1.75$ is due to the increase of M_gas.

We note that this effect is not inherent to our specific data points, or even best-fit evolution laws of SFR and M_gas, as a similar evolution could be observed if we were to generate catalogues from the laws of S15 for SFR, and Wang et al. (2022) or Tacconi et al. (2018) for M_gas (following the method used to generate ρ_SFR and ρ_gas see Sects. 6 and 7). If we were to choose laws from Speagle et al. (2014) for SFR and Tacconi et al. (2018) for M_gas, this effect would disappear completely, and we would end up observing an evolution of the SFE with redshift. We conclude, in this case, that the correct normalisation and the presence of the bending of the main sequence (which are lacking in the study of Speagle et al. 2014) are essential to observe this effect. One could advocate that the observation of an increase in efficiency with redshift when generating catalogues, for example, from S15 for SFR, and Tacconi et al. (2018) for M_gas, comes from the fact that the M_gas fitting form from Tacconi et al. (2018), or the one use in this work, do not allow a bending or an evolution of it. However, we found no evidence of a bending or its evolution at high mass when we tried to fit our data using various forms allowing them. Overall, it suggests that some physical processes limit the star formation efficiency in low-redshift galaxies, as compared to high-redshift galaxies, beyond the simple fact that they have less available gas.

Overall, even though the evolution of the SFE with redshift is not well constrained, this study suggests that galaxies at high-redshift form stars more efficiently at a fixed gas surface density. This excess cannot be explained by a universal Kennicutt-Schmidt relation nor SFE over all redshifts. This may partly explain the early results of JWST: a higher ρ_UV at z ≥ 8 (Finkelstein et al. 2023; Donnan et al. 2023; Mason et al. 2023) compared to pre-JWST extrapolation Oesch et al. (2018); more massive galaxies at z ≥ 5 (Labbe et al. 2023; Xiao et al. 2023a). This may be explained by the model proposed by Dekel et al. (2023) who showed that the environment of the most massive dark-matter haloes at z ≥ 10 favour high star formation efficiencies.

In conclusion, the gradual appearance of the bending of the main sequence with redshift appears to result from a reduction of the SFE, that is, a gradual disappearance of this extra factor affecting the SFE (observed at z ∼ 4.3 and absent at z ∼ 2). Then, at a fixed lower redshift (i.e. z <  2), the bending of the main sequence (as M_⋆ increases) is mainly due to an increase in M_gas with stellar mass, and to second order, to a decrease in SFE with stellar mass. This characteristics is also observed in the local universe (i.e. for 0.01 <  z <  0.05) in Saintonge et al. (2016): when the stellar mass increases by a factor 10, the SFR increases by a factor ∼1.4, which is due to an increase of M_gas by a factor ∼1.8 and a decrease in SFE by a factor ∼1.25. Saintonge et al. (2016) reach the same conclusion in the local universe: at a fixed redshift, M_gas is the main reason for the bending of the main sequence. The slow downfall of the star formation efficiency in the most massive galaxies (i.e. M_⋆ ≥ 10¹¹ M_⊙) observed between 0.1 ≤ z ≤ 4.3, was also observed in Schreiber et al. (2016). It seems that a discontinuation of gas accretion cannot be the only reason for the bending of the SFR main sequence (Daddi et al. 2022). It must be coupled with other effects that prevent the gas in high M_⋆ galaxies from cooling and collapsing into stars (i.e. an evolution of the SFE). The main driver of the bending could be the quenching of galaxies due to environmental effects, as massive galaxies populate denser environments where they suffer from ram pressure stripping, galaxy harassment (e.g. Kalita et al. 2022).

9. Summary and conclusions

For this work we gathered catalogues of H-band selected SFGs over four fields (GOODS-South, GOODS-North, COSMOS, UDS). This sample is stacked over seven wavelengths (24 μm, 100 μm, 160 μm, 250 μm, 350 μm, 500 μm and 1.13$\mathrm{mm}$). The main addition to this work, compared to previous studies, is the use of the GOODS-ALMA survey (Gómez-Guijarro et al. 2022b), which helps to reach lower M_⋆ and higher redshifts, when combined with Herschel. We combined the stacking of H-band selected galaxies with H-dropout galaxy properties to obtain a complete view of the cosmic history of galaxies. We derived the evolution of MS galaxy properties such as SFR, T_dust, and M_gas through a consistent analysis. We generated a detailed view of ρ_SFR and ρ_gas depending on key galaxy properties such as M_⋆, L_IR and τ_dep.

This study allowed us to confirm several features of galaxies, the main ones being summarised as follows:

• The T_dust evolves linearly with redshift. There is no clear evidence for a dependence of T_dust depending on M_⋆ at a fixed redshift, when considering MS galaxies.

• The SFR-M_⋆ MS has a slope close to unity at low masses (i.e. M_⋆ ≲ 10¹⁰ − 10¹¹ M_⊙), with a bending appearing below z ∼ 2 on the high-mass end. The bending of the MS corresponds to the slow downfall of SFE in massive galaxies.

• The total ρ_SFR follows a form close to the one presented in Madau & Dickinson (2014), at least up to z = 5. As a result, the contribution of H-dropout galaxies to the total ρ_SFR is ∼13% at z = 5.

The main new results from this paper can be summarised as follows:

• We measured the contribution of galaxies of different M_⋆ to the total ρ_SFR over cosmic time. We find that massive galaxies (i.e. M_⋆ ≥ 10¹⁰ M_⊙) account for most of the ρ_SFR up to z ∼ 5. Low-mass galaxies (i.e. M_⋆ ≤ 10¹⁰ M_⊙) contribute a roughly constant fraction to the total ρ_SFR at 0.4 ≤ z ≤ 5.

• The global discrepancy between the cosmic SFR density observed and the TNG100 one is discussed in Pillepich et al. (2018b). However, we show that the TNG100 simulation fails to reproduce the contributions of the different bins of M_⋆ to the total ρ_SFR. In particular, in the TNG100 simulation, massive galaxies do not form enough stars around cosmic noon, when the cosmic SFR density peaks (i.e. ρ_{SFR − TNG100}^{M_⋆ ≥ 1010 M_⊙} ≃ 50% ρ_{SFR − observed}^{M_⋆ ≥ 1010 M_⊙} at z = 2). We conclude that the mechanisms used to regulate star formation (in particular feedback from AGNs) are possibly too strong in this simulation, namely too efficient at producing a drop in the gas content in massive galaxies.

• We estimate that H-dropout galaxies (as defined in Wang et al. 2019) contribute to ∼23% of the ρ_SFR of galaxies with M_⋆ ≥ 10^10.3 M_⊙ at z = 4–6.

• We find hints that high-redshift galaxies convert their gas reservoirs more efficiently than local galaxies, that is, suggesting a possible evolution of SFE with redshift.

Acknowledgments

This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. C.G.G. acknowledges support from CNES.

References

Appendix A: Table of gas masses

Appendix B: Using only the Scoville16 method

In this appendix, we look at the results that would be deduced if we only considered gas masses that could be deduced from our ALMA measurements and using the Scoville et al. (2016) method to derive M_gas. This appendix goes through the accordingly updated Figures 7, 11, 12, 19, 20 and 21.

We follow the same method we went through in Sect. 5.3, we add the contribution of the H-dropouts to M_gas. In that case, Eq. 20 allows to get an acceptable fit. Therefore, we chose to fit the final M_gas with Eq. 20. Our best fit parameters are given in Table B.1.

Our new gas masses measurements are deduced from the method from Scoville et al. (2016) are displayed in Fig. B.1. In Fig. B.1, we also compare our results to the trends observed by Tacconi et al. (2018) and Wang et al. (2022). In Fig. B.1, the trends are displayed for stellar mass bins in the Salpeter (1955) IMF. Using this method, our results coincide with the one from Wang et al. (2022), which is to be expected as they are both coming from a mass selected sample, deducing gas masses from millimetre measurements through the method of Scoville et al. (2016).

Next, the resulting gas mass densities are deduced through the method we describe in Sect. 7. We display the results in Fig. B.2. When deducing gas masses using the method from Scoville et al. (2016) only we end up with results much closer to what is found in the literature (e.g. Tacconi et al. 2018; Wang et al. 2022).

We compare our results to Magnelli et al. (2020) in Fig. B.3. Again, results are much more in line with what is found in Magnelli et al. (2020).

Fig. B.1. Mgas as a function of redshift for different M⋆. The crosses represent Mgas estimates from this work following the method from Scoville et al. (2016). The squares are measurements from Saintonge et al. (2017). The bins that have been corrected for the H-dropout contribution are circled. The thick solid colour lines represent the best-fit from this work using Eq. 20.The shaded area the 68% uncertainty of the fit. The faded solid lines represent the best-fit from this work colour-coded by M⋆ (using Eq. 21), the faded dashed lines represent the trend from Tacconi et al. (2018), the faded dotted lines are from Wang et al. (2022).

Fig. B.2. ρgas as a function of redshift. The green dotted lines represent the total ρgas inferred from this work (i.e. integrating the stellar mass function down to 3 × 109 M⊙) using Mgas deduced with the Scoville et al. (2016) method only. The red and blue dotted lines represent the contribution to the total ρgas of galaxies with log10(M⋆)> 10 and log10(M⋆)< 10, respectively. The corresponding result deduced throughout this paper are displayed in faded solid red, blue and green coloured lines. The purple and cyan lines represent ρgas deduced from the Mgas definition of Wang et al. (2022) and Tacconi et al. (2018), respectively.

Fig. B.3. ρgas as a function of redshift. The green, orange and red dotted lines represent ρgas derived from this work with $M_{\star }^{\lim }$ of 1.7 × 108, 1.7 × 109 and 1.7 × 1010 M⊙, respectively, using Mgas deduced with the Scoville et al. (2016) method only. The corresponding result deduced throughout this paper are displayed in faded solid red, orange and green coloured lines. The green, orange and red hatches represent estimates of ρgas from Magnelli et al. (2020) with $M_{\star }^{\lim }$ of 1.7 × 108, 1.7 × 109 and 1.7 × 1010 M⊙, respectively. The green dots with black circles depict the estimate of ρgas from Magnelli et al. (2020) for M⋆ >  Mlimit, Mlimit being the stellar mass completeness limit from Magnelli et al. (2020).

Following the method we describe in Sect. 8.3, we computed the SFE corresponding to the gas masses deduced from the method of Scoville et al. (2016). We display the results in Fig. B.4. In that case, we can still observe an increase of the SFE with redshift.

Fig. B.4. SFE as a function of stellar mass. Squares linked with a dashed line represent the measured from this work when computed using Mgas deduced with the Scoville et al. (2016) method only, colour coded by bin of redshift. Dotted lines are the value SFE deduced from our best fit law of the SFR and Mgas main sequences.

We also investigate whether the evolution of SFE can be explained by the gas content of galaxies at high redshift by looking at the Kennicutt-Schmidt relation (Kennicutt 1998b). We follow the same method presented throughout Sect. 8.3. In that case, our best-fit of the linear dependence of log₁₀(Σ_SFR) as a function of log₁₀(Σ_gas) suggests a power-law scaling with a slope of $N=1.{29}_{-0.01}^{+0.01}$. The results of the Kennicutt-Schmidt relation are shown in Fig. B.5. In Fig. B.5, an increase of the normalisation can still be observed at z >  3. This effect is, however, much more limited compared to the one observed though the main studied (see Fig. 20).

Fig. B.5. Kennicutt-Schmidt relation: ΣSFR as a function of Σgas. The dots represent the data from this work when computed using Mgas deduced with the Scoville et al. (2016) method only, and are colour coded by redshift. The cyan squared error bar represents the median of the sample distribution from Kennicutt (1998b), with individual galaxies displayed as grey crosses. The dashed black line represents the best-fit of all points. Coloured lines represent the best-fit, by fixing the slope to the all points best-fit value (i.e. $N=1.{29}_{-0.01}^{+0.01}$), at fix redshift.

We corrected the SFE from the gas contribution by looking at SFE/$M_{\mathrm{gas}}^{(N-1)}$ instead of SFE (see Fig. B.6). In Fig. B.6, we can still observe an evolution with redshift (i.e. SFE cannot be explain only by the amount of M_gas).

Fig. B.6. SFE/$M_{\mathrm{gas}}^{(N-1)}$ (N = 1.29) as a function of stellar mass. Squares linked with a dashed line represent the measured from this work when computed using Mgas deduced with the Scoville et al. (2016) method only, colour coded by bin of redshift.

Our conclusion on these results (i.e. using the method from Scoville et al. (2016) to deduce gas masses) and the discrepancies observed from the ones of the main work of this paper (i.e. using an MZR method to deduce gas masses from the IR to millimetre) is that the main reason these higher and evolving SFE were not observed is mostly due to the method used to infer the gas masses: IR to millimetre compared to millimetre only methods.

All Tables

All Figures

	Fig. 1. Best-fit SED for each bin of redshift and stellar mass. Blue dots correspond to the flux measurements, red triangles represent the 5σ upper limits. Blue line is the best-fit SED, the blue shaded area shows the 68% uncertainty of the fit. Red line is the SED maximising LIR in case only upper limits are available.
In the text

	Fig. 2. Tdust as a function of M⋆ colour-coded by the redshift bin. Average trend of Tdust as a function redshift from S18, is shown as a reference in faded coloured dashed line.
In the text

Fig. 3. Tdust and ⟨U⟩ as a function of redshift. The blue dots represent the Tdust of this work, the blue line is the best-fit (up to z = 5), the dashed blue line is the best-fit extrapolation (for z ≥ 5), and the blue shaded area represents the 68% uncertainty of the fit. From the literature: S18 (grey solid line), Magnelli et al. (2014) (grey dotted line) and Bouwens et al. (2020) (grey dash-dotted line) converted using Eq. (8). Magdis et al. (2012) (grey, dashed line) and Béthermin et al. (2015) (grey dots) converted using Eq. (9). We also re-fitted the two stacks from Béthermin et al. (2020) (grey dotted error bars) using template library from S18.

In the text

	Fig. 4. λpeak as a function of LIR. The colour-coded error bars (per redshift bins) represent the results of this work, coloured lines represent the best-fit and shaded area the 68% uncertainty of the fit. The solid black line corresponds to the trend of Drew & Casey (2022), while the dotted black line is an extrapolation of their relation.
In the text

Fig. 5. SFRMS as a function of M⋆ over different redshift bins. The dots represent data from this work. The upper limits were used to perform the fit but are not shown here to avoid overloading the figure. Left panel: fitted by the Eq. (12), the shaded area the 68% uncertainty of the fit. Right panel: fitted by the Eq. (13), the shaded area the 68% uncertainty of the fit. The squares represent the best-fit parameter M0. No bending was detected for z >  3.1 and therefore no square error bars are displayed.

In the text

Fig. 6. SFRMS as a function of M⋆. The red triangles represent upper limits, the blue dots are the observed values from this work, and shaded area the 68% uncertainty of the fit. The blue lines represent the best-fit from this work. The green, cyan and orange lines represents the main sequence from S15, Speagle et al. (2014) and Leslie et al. (2020), respectively. The dashed lines represent extrapolated main sequences to redshifts that were not investigated in their respective studies.

In the text

Fig. 7. Mgas as a function of redshift for different M⋆. The dots represent Mgas estimates from this work using the ${\delta }_{\mathrm{GDR}}-\mathrm{MRZ}$ method. the crosses represent Mgas estimates from this work following the method from Scoville et al. (2016). The squares are measurements from Saintonge et al. (2017). The bins that have been corrected for the H-dropout contribution are circled. The thick and thin solid colour lines represent the best-fit from this work, using Eq. (21) and Eq. (20), respectively. The dashed lines are the fits extrapolation. The shaded area the 68% uncertainty of the fit.

In the text

	Fig. 8. Mgas as a function of redshift and M⋆. The solid lines represent the best-fit from this work colour-coded by M⋆ (using Eq. (21)), the dashed lines represent the trend from Tacconi et al. (2018), the dotted lines are from Wang et al. (2022). The faded lines represent an extrapolation from their respective laws.
In the text

Fig. 9. Cosmic star formation rate density (ρSFR) as a function of redshift. The green line represents the total ρSFR trend from this work (i.e. integrating the stellar mass function down to 3 × 109 M⊙). The red and blue lines represent the contribution of galaxies with log10(M⋆/M⊙)> 10 and log10(M⋆/M⊙)< 10, respectively. The black dotted line represents the ρSFR detected with stacking. The purple squares show the contribution of H-dropout from Wang et al. (2019). The dotted green and red lines show the trend once H-dropout contribution from Wang et al. (2019) has been added to the respective coloured solid lines. The grey line represents the trend from Madau & Dickinson (2014). The cyan and orange squares represent ρSFR estimate from Leslie et al. (2020) and Gruppioni et al. (2020), respectively.

In the text

	Fig. 10. Contributions over the whole range of stellar masses to the total ρSFR as a function of redshift for SFGs.
In the text

Fig. 11. ρgas as a function of redshift. The green line represents the total ρgas inferred from this work (i.e. integrating the stellar mass function down to 3 × 109 M⊙). The red and blue lines represent the contribution to the total ρgas of galaxies with log10(M⋆)> 10 and log10(M⋆)< 10, respectively. The purple and cyan lines represent ρgas deduced from the Mgas definition of Wang et al. (2022) and Tacconi et al. (2018), respectively.

In the text

Fig. 12. ρgas as a function of redshift. The green, orange and red lines represent ρgas derived from this work with $M_{\star }^{\lim }$ of 1.7 × 108, 1.7 × 109 and 1.7 × 1010 M⊙, respectively. The green, orange and red hatches represent estimates of ρgas from Magnelli et al. (2020) with $M_{\star }^{\lim }$ of 1.7 × 108, 1.7 × 109 and 1.7 × 1010 M⊙, respectively. The green dots with black circles depict the estimate of ρgas from Magnelli et al. (2020) for M⋆ >  Mlimit, Mlimit being the stellar mass completeness limit from Magnelli et al. (2020).

In the text

Fig. 13. Contribution to ρSFR of SFGs as a function of redshift. Catalogues were binned through different methods and a specific property is displayed through colours: see each sub-figure for specifics. (a) Each bin includes 10% of the galaxies in number, picking from the lowest M⋆ onwards. The minimum M⋆ in the bin defines the colour. The numbers give the median M⋆ within the bin along with the lower and maximum extension of the bin. (b) Each bin includes 10% of the total M⋆ of all galaxies, picking from the lowest M⋆ onwards. The median τdep in the bin defines colour. (c) Each bin includes 10% of the total M⋆ of all galaxies, picking from the lowest M⋆ onwards. The minimum M⋆ in the bin defines the colour. Contours of log10(M⋆=M⊙) = 10 and log10(M⋆=M⊙) = 11 are added as black dashed lines.

In the text

	Fig. 14. Stellar mass above which galaxies contribute exactly 50% of the CSFD as a function of redshift.
In the text

Fig. 15. Cosmic star formation rate density (ρSFR) as a function of redshift. The green, red and blue lines represent ρSFR deduced from the TNG100 simulation integrated from 3 × 109 M⊙, including all galaxies, galaxies with log10(M⋆/M⊙)> 10, and galaxies with log10(M⋆/M⊙)< 10, respectively. The green, red and green faded dots represent the total ρSFR trend, ρSFR of galaxies with log10(M⋆/M⊙)> 10, and galaxies with log10(M⋆/M⊙)< 10, respectively (including H-dropout contribution from Wang et al. 2019), from this work. The grey line represents the trend from Madau & Dickinson (2014) for reference.

In the text

	Fig. 16. Ratio of the ρSFR of high-mass (M⋆ >  1010 M⊙) over low-mass (M⋆ <  1010 M⊙) galaxies, as a function of redshift. The red line represents the ratio deduced from this work (i.e. from observations). The blue line represents the ratio deduced from the TNG100 simulation.
In the text

	Fig. 17. fgas as a function of M⋆ for SFGs in TNG100. Coloured lines represent the evolution deduced from the best-fit of Mgas from this work (see Table 7).
In the text

	Fig. 18. SFR as a function of M⋆ for SFGs in TNG100. Coloured lines represent the best-fit from this work (see Table 5).
In the text

	Fig. 19. SFE as a function of stellar mass. Squares linked with a dashed line represent the measured from this work, colour coded by bin of redshift. Dotted lines are the value SFE deduced from our best fit law of the SFR and Mgas main sequences.
In the text

Fig. 20. Kennicutt-Schmidt relation: ΣSFR as a function of Σgas. The dots represent the data from this work, and are colour coded by redshift. The cyan squared error bar represents the median of the sample distribution from Kennicutt (1998b), with individual galaxies displayed as grey crosses. The dashed black line represents the best-fit of all points. Coloured lines represent the best-fit, by fixing the slope to the all points best-fit value (i.e. $N=1.{18}_{-0.01}^{+0.01}$), at fix redshift.

In the text

	Fig. 21. SFE/$M_{\mathrm{gas}}^{(N-1)}$ (N = 1.18) as a function of stellar mass. Squares linked with a dashed line represent the measured from this work, colour coded by bin of redshift.
In the text