© The Authors 2024

1. Introduction

The increasing amount of multiwavelength data, obtained as part of wide-field and deep surveys, enables a progressively detailed examination of large samples of galaxies. Substantial advancements in determining important physical properties of galaxies (e.g., photometric redshifts, stellar masses, dust attenuation curves) have enabled the identification (in the rest-frame color-color space) of two separate types of galaxy populations: star-forming ones and quiescent ones (SFGs and QGs; e.g., Labbé et al. 2005; Wuyts et al. 2007; Ilbert et al. 2010). The main difference between these two distinct populations of galaxies is their star formation activity; it has been shown that the star formation rate (SFR) and stellar mass (M_*) of the SFGs follow a tight correlation out to the highest redshifts probed – the so-called star-forming main sequence (MS; e.g., Speagle et al. 2014; Tomczak et al. 2016; Daddi et al. 2022; Popesso et al. 2023). Most commonly, the SFR-M_* relation can be parameterized with a simple power law:

$$\begin{array}{c}\hfill \log (\mathrm{SFR}/{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1})=\alpha \log ({\mathrm{M}}_{\ast }/{\mathrm{M}}_{\odot })+\beta ,\end{array}$$(1)

where α and β represent the slope and normalization of the MS, respectively. While the normalization has been shown to increase in redshift, with galaxies exhibiting higher SFRs at earlier epochs, the slope of the MS has been found to be close to unity, mostly independent of redshift, with values ranging between ∼0.8 and 1.3 (e.g., Lee et al. 2015; Tomczak et al. 2016; Popesso et al. 2023). More recently, a number of studies have found that the star-forming MS flattens at higher masses (e.g., Lee et al. 2015; Schreiber et al. 2015; Leslie et al. 2020; Leja et al. 2022; Popesso et al. 2023), reaching an asymptotic value of the SFR.

The exact shape and evolution of the MS is still widely discussed in the literature, with many studies presenting somewhat different perspectives. The sources of the discrepancies are many, including the assumed initial mass function (IMF), stellar population synthesis modeling, star formation histories (SFHs), and dust extinction curves. Most importantly, different techniques used for the estimation of the total SFRs can affect the results significantly. Since part of the UV emission coming from young stars gets absorbed by the interstellar dust, it is crucial to be able to properly quantify this effect. Measuring the dust attenuation from the rest-frame UV/optical data alone produces very uncertain results, because of the degeneracy between the age (different SFHs) and reddening (different metallicities, extinction curves). The advent of the Herschel Space Observatory (Pilbratt et al. 2010) allowed for the direct detection of the energy being re-emitted by dust in the far-infrared (FIR), and therefore a more reliable determination of the total (UV+IR) SFRs. However, since Herschel FIR data suffers from a relatively low resolution, only the most IR-luminous sources can be detected at high redshifts, and hence samples of SFGs studied with Herschel tend to be incomplete, often biasing the results toward large values of the SFRs. Speagle et al. (2014) and, more recently, Popesso et al. (2023) compiled various studies in the literature in order to convert them to the same absolute calibrations in an attempt to solve some of the issues listed above. While they carefully address all the inconsistencies caused by different modeling assumptions, the uncertainties associated with the various total SFR indicators, as well as the completeness of samples compiled, remain a concern.

In this work, we make use of the sample of ∼100k K-band-selected galaxies from the ground-based UV+optical+near-IR ∼2 deg² imaging of the UKIDSS Ultra Deep Survey (UDS) and COSMOS fields (McLeod et al. 2021). Galaxies are stacked in bins of redshift and stellar mass in the FIR Herschel and the James Clerk Maxwell Telescope (JCMT) SCUBA-2 (Holland et al. 2013) maps, in order to determine the shape and redshift evolution of the star-forming MS out to z ≲ 5.7. The stacking approach, used by several authors before (e.g., Tomczak et al. 2016; Leslie et al. 2020; Mérida et al. 2023), has several advantages. The sample of galaxies stacked is mass-complete, assuring that the resulting values of the IR luminosity (and SFR) are not biased due to selection effects. In addition, thanks to the inclusion of the SCUBA-2 850 μm data, our stacked photometry maps the Rayleigh-Jeans portion of the dust emission spectrum, allowing for the precise determination of the dust temperatures, and thus the SFRs.

The paper is structured as follows. In Sect. 2, the data used in this work is described. We outline the stacking procedure and the star-forming sample selection process in Sect. 3. Section 4 explains the methodology used in the determination of the star-forming MS for our sample. We discuss the results of this work in Sect. 5. In particular, the shape and evolution of the star-forming MS and its relation to the so-called Schmidt-Kennicutt relation (SK relation; Schmidt 1959; Kennicutt 1998) is explained in Sect. 5.1, followed by a discussion of the variations in the peak dust temperature with IR luminosity and redshift in Sect. 5.2 and a comparison to other works in Sect. 5.3. In Sect. 5.4, we attempt to interpret some of the discrepancies between the different shapes of the MS found in the literature, by “growing” the observed galaxy stellar mass function (GSMF) in time. Finally, we summarize our conclusions in Sect. 6. Throughout the paper, we use the Chabrier (2003) stellar IMF with an assumed flat cosmology of Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. Data

2.1. Optical/near-IR catalogs

We started with the optical/near-IR catalogs of the UKIDSS Ultra Deep Survey (UDS) and COSMOS fields. The details of the data used, galaxies selection process and determination of photometric redshifts and stellar masses are presented in McLeod et al. (2021), with the shortened description given below.

2.1.1. Imaging

For the UKIDSS UDS field (Lawrence et al. 2007), we included UDS DR11 imaging (Almaini et al., in prep.) in the near-IR JHK bands, with Y-band imaging from the VISTA VIDEO DR4 release over XMM-LSS (Jarvis et al. 2013). We complemented this with UV imaging in u^*-band from CFHT MegaCam, and optical imaging in BVRiz′ from Subaru Suprime-Cam (Furusawa et al. 2008). We also employed imaging in the narrowband NB921 (Koyama et al. 2011) and the refurbished Suprime-Cam z-band (z′_new; Furusawa et al. 2016). To extend our wavelength coverage into the mid-infrared, we utilized Spitzer IRAC imaging in 3.6 μm and 4.5 μm, combining the SEDS (Ashby et al. 2013), S-CANDELS (Ashby et al. 2015), and SPLASH (PI Capak; see Mehta et al. 2018) programs. The effective area of the overlapping UDS imaging was 0.69 deg², after masking for bright stars.

For the COSMOS field, we combined near-IR YJHK_s imaging from UltraVISTA DR4 (McCracken et al. 2012) with UV/optical u^*griz imaging from the CFHTLS-D2 T0007 data release (Hudelot et al. 2012). We confine our study to the 1 deg² overlapping area covered by these data sets. After masking, this reduced to an effective area of 0.86 deg². As with the UDS, we also included Subaru Suprime-Cam $z{\prime }_{\mathrm{new}}$ imaging. There is a wealth of Spitzer IRAC imaging available in 3.6 μm and 4.5 μm over the COSMOS field. The programs we combined include S-COSMOS (Sanders et al. 2007), SPLASH (PI Capak), SEDS (Ashby et al. 2013), S-CANDELS (Ashby et al. 2015) and SMUVS (Ashby et al. 2018).

2.1.2. Photometric redshifts and stellar masses

As is explained in detail in McLeod et al. (2021), the parent multiwavelength catalogs were constructed by running SOURCE-EXTRACTOR (Bertin & Arnouts 1996) in dual-image mode, with the K-band as the detection image. Photometry was measured in 2-arcsec diameter apertures on the PSF-homogenized UV to near-IR images. For the lower resolution Spitzer IRAC 3.6 and 4.5 μm data, the photometry was measured using the software TPHOT (Merlin et al. 2015).

For each object in the catalog, the photometric redshift we assigned was the median of six SED fitting runs, employing three different software codes. With LEPHARE, we included runs with the Bruzual & Charlot (2003), PEGASE (Fioc & Rocca-Volmerange 1999) and COSMOS (Ilbert et al. 2009) template libraries. We also included two EAZY (Brammer et al. 2008) runs, using the PCA and PEGASE templates, and one further run with the BPZ code (Benítez 2000), using the Coleman et al. (1980) templates. The resulting values of σ_z, defined as 1.483 × MAD(dz), where MAD is the median absolute deviation and dz = (z_phot − z_spec)/(1 + z_spec), for COSMOS and UDS fields were found to be equal to 0.019 and 0.022, respectively (McLeod et al. 2021).

Stellar masses were measured for each object by fixing the photometric redshift to z_med, and re-fitting the SED using LEPHARE with the Bruzual & Charlot (2003) library. This template set included a Chabrier (2003) IMF, a Calzetti et al. (2000) dust attenuation law and IGM absorption as in Madau (1995). The values were found to be in excellent agreement with those derived by the CANDELS team (see Table 2 of McLeod et al. 2021), with a tight 1:1 relation and a typical scatter of ±0.05 dex. The resulting catalogs were then cleaned of stars, AGN, artifacts and objects with contaminated photometry. Finally, a further χ² < 40 goodness-of-fit cut was implemented on the catalogs.

2.2. Far-infrared

For the extraction of the stacked FIR flux densities in the COSMOS and UDS fields we used the Herschel (Pilbratt et al. 2010) Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012) and the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) Evolutionary Probe (PEP; Lutz et al. 2011) data obtained with the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) and PACS instruments. We utilized Herschel maps at 100, 160, 250, 350 and 500 μm with beam sizes of 7.39, 11.29, 18.2, 24.9, and 36.3 arcsec, and 5σ sensitivities of 7.7, 14.7, 24.0, 27.5, and 30.5 mJy, respectively.

In order to constrain the Rayleigh-Jeans tail of the dust emission curve, we have also included the data collected as a part of the SCUBA-2 Cosmology Legacy Survey (S2CLS; Geach et al. 2017), with the UKIDSS-UDS and COSMOS 850 μm imaging covering ≃0.9 deg² with a 1σ noise of 0.9 mJy and ≃1.3 deg² with the 1σ noise of 1.6 mJy, respectively.

3. Stacking

3.1. Sample selection

For the purpose of successful stacking (i.e., to maximize the number of sources in each stack), we designed each redshift bin to encompass ∼1 billion years of Universe evolution, ranging from 1 to 9 billion years after the Big Bang (with the corresponding redshift bins of [0.45, 0.6, 0.75, 1.0, 1.25, 1.6, 2.2, 3.2, 5.7]; see Table A.1). We also made sure that the sample we stacked is mass-complete. The procedure for estimating the mass completeness in the UDS and COSMOS fields is explained in McLeod et al. (2021). In brief, the completeness was estimated using simulations, in which K-band imaging was produced by injecting a set of artificial galaxies with a wide range of physical properties. Adopting the original setup, the corresponding photometry catalogs were then constructed using SEXTRACTOR. In order to determine the 90% mass-completeness limit, the relation from Pozzetti et al. (2010) was used to derive the limiting stellar mass for each object in our sample. The 90% mass completeness limit, below which 90% of the limiting stellar masses lie, could then be derived as a function of redshift. The resulting redshift and stellar mass bins adopted in this work are listed in Table A.1.

In order to derive the SFR-M_* relation for the star-forming sample, we first needed to separate it from the quiescent population. Here, we adopted the selection based on the NUVrJ colors (Ilbert et al. 2013):

$$\begin{array}{cc}\hfill (NUV-r)>& 3\times (r-J)+1;\hfill \\ \hfill (NUV-r)>& 3.1,\hfill \end{array}$$(2)

where NUV, r and J rest-frame magnitudes were generated from the best-fit SEDs. The advantage of this method is that it is able to separate sources for which reddening is caused by aging from those for which it is produced by dust extinction, thereby decreasing the chance of contaminating the dust-enshrouded star-forming sample with quiescent galaxies. Because at high redshifts the rest-frame J-band can no longer be traced by the photometry available in this work, the NUVrJ color-color selection was performed out to z ∼ 4, as is depicted in Fig. 1.

Fig. 1. Results of the NUVrJ color-color selection results, defined in Eq. (2) shown for redshift bins adopted in this work. The NUVrJ colors (Ilbert et al. 2013) are used in order to distinguish between sources reddened due to the lack of star formation (quiescent galaxies – upper left section of the plots) from those where red colors are caused by dust extinction (SFGs – main locus in the plots).

To avoid overestimating the SFRs at the high-mass end of the MS, we further detected and removed the starburst galaxy population from our stacking sample. To identify starbursts, we used the ALMA-detected sources in the COSMOS (Liu et al. 2019) and UDS (Stach et al. 2019) fields, with associated stellar masses and SFRs, as has been found by Liu et al. (2019) and Dudzevičiūtė et al. (2020), respectively. Following Elbaz et al. (2018), starbursts are defined as galaxies with SFR/SFR_MS > 3, where SFR_MS is the corresponding MS value, given the redshift and stellar mass of each ALMA-detected source, with the functional form of the star-forming MS of Popesso et al. (2023).

3.2. Far-infrared fluxes

Since current FIR/submillimeter surveys struggle to reach the (sub-)millijansky depths necessary for the detection of the dust emission in typical SFGs, we have resorted here to the well-known method of stacking (e.g., Whitaker et al. 2014; Schreiber et al. 2015; Tomczak et al. 2016; Koprowski et al. 2018). In order to determine the FIR fluxes at a given band, we cut a square stamp around the position of each source, with the size of approximately 3× the full width at half maximum (FWHM) of the corresponding beam. The average FIR fluxes at each pixel were then evaluated following the inverse-variance weighting:

$$\begin{array}{c}\hfill S_{\nu }=\frac{{\sum }_i{S}_{\nu ,i}/{\sigma }_i^2}{{\sum }_{i}1/{\sigma }_i^2},\end{array}$$(3)

where S_ν, i is the flux density of the ith source and σ_i is the 1σ instrumental noise, extracted at each pixel in a given source’s stamp. The error on the mean stacked fluxes for each pixel was calculated according to

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{stack}}^2(S_{\nu })=\frac{1}{{\sum }_{i}1/{\sigma }_i^2},\end{array}$$(4)

with the example 500 μm stamps for each redshift and stellar mass bin given in Data availability.

The stacked average FIR flux densities at each bin were then extracted from the nonlinear least squares best fits to the light profiles derived from corresponding stamps. The reason for this is that, at the Herschel 250, 350, and 500 μm SPIRE bands, due to the relatively large beam sizes, the clustering of objects associated with the stacked sources is seen to significantly contribute to the mean flux. Figure 2 shows the light profile and the best fits for the sample of Herschel SPIRE 500 μm sources at 1.25 < z < 1.75 with stellar masses between 10.25 < log(M_*/M_⊙) < 10.5. As can be seen, the stacked flux (red line) consists of the average flux from the stacked sample (dashed blue line) and the background flux from the associated sources (dotted green line), both of which can be modeled with a Gaussian function. In cases when a double Gaussian (beam + background) gives better fits than a single Gaussian (beam), the background is subtracted and the resulting average flux density of the stacked sample is extracted from a beam component of the best fit.

Fig. 2. 1.25 < z < 1.75, 10.25 < log(M*/M⊙) < 10.5 light profile derived at 500 μm. The red points represent mean values of the flux densities derived from stacked images as a function of distance from the image center. The dashed blue, dotted green, and solid red lines represent the corresponding beam, the contribution from background sources, and the sum of the two, respectively. For Herschel SPIRE bands, the contribution from background sources is significant and must therefore be subtracted from the overall light profile, when deriving stacked flux densities (see Sect. 3.2).

The errors on the stacked FIR fluxes, which are related to the redshifts, stellar masses, and the NUVrJ quiescent galaxy selection procedure uncertainties, were estimated using the bootstrapping method. In each step, a mock catalog was constructed by drawing, at random and with replacement, a sample of sources from the original mass-complete dataset. The stacking procedure, described above, was then applied and the average fluxes were extracted. The process was repeated 1000 times and the uncertainties on the FIR fluxes, σ_bootstrap, taken to be the standard deviation of the resulting simulated values. In addition, the contribution to the final stacked FIR fluxes errors from the uncertainties related to the functional fits to the light profiles, σ_fit, described above, as well as the errors related to the individual observed FIR fluxes from Eq. (4) were included, with the final errors on the stacked values equal to

$$\begin{array}{c}\hfill \sigma (S_{\nu })=\sqrt{{\sigma }_{\mathrm{bootstrap}}^2+{\sigma }_{\mathrm{fit}}^2+{\sigma }_{\mathrm{stack}}^2}.\end{array}$$(5)

4. The star formation rate–stellar mass relations

4.1. Star formation rate measurements

A total SFR consists of the UV and IR components, where the rest-frame UV light is sensitive to the young massive stars and is therefore a good tracer of recent, unobscured star formation, while the IR emission accounts for the stellar emission that has been absorbed by the interstellar dust. Here, we followed the recent calibration of Kennicutt & Evans (2012), where

$$\begin{array}{c}\hfill {\mathrm{SFR}}_{\mathrm{tot}}=1.71\times {10}^{-10}(L_{\mathrm{FUV}}+0.46\times L_{\mathrm{IR}}).\end{array}$$(6)

The far-UV (FUV) part of the SFR was derived from the FUV luminosity measured at 150 nm, where L_FUV = νL_ν, with FUV SFRs calculated for each galaxy and then the mean value evaluated at each redshift or stellar mass bin.

In order to find the IR luminosity, L_IR, the dust emission curve of Casey (2012) was fit to the available FIR photometry using the nonlinear least squares analysis. The adopted dust emission curve, described in Table 1 of Casey (2012), can be approximated as a classic single-temperature gray-body function at FIR wavelengths, with the addition of the power law fit at mid-IR wavelengths. The data was initially fit with both the single gray-body and the power law gray-body curves and it was found that the latter gives better fits to the data at mid-IR bands. The power law gray-body curve has four free parameters – the mid-IR power law slope, α, the emissivity index, β, the effective dust temperature, $T_{\mathrm{d}}^{\mathrm{eff}}$, and the overall normalization, N. Since, at z ≳ 2, the available photometry redward of the dust emission peak consists only of the 850 and 500 μm, with the 500 μm being affected by blending (see Fig. 2), we decided to fix β at 1.96. Similarly, due to the lack of photometry at a rest-frame wavelength of λ ≲ 80 μm, we also fixed α at 2.3. Values of α and β were fixed according to the recommendations of Drew & Casey (2022). The curve was then integrated between 8 and 1000 μm and the errors were estimated from the corresponding covariance matrix. The stacked photometry and the best fits for the whole (star-forming and quiescent) and star-forming subsamples are given in Data availability. The resulting best-fit values of the total SFR at each redshift or stellar mass bin for the whole and star-forming subsamples are listed in the top and bottom panels of Table A.1, as well as being depicted as color points with error bars in the left and right panels of Fig. 3, respectively.

Fig. 3. SFR-M* relations derived in Sect. 4 for the whole (star-forming and quiescent – left panel) and the star-forming (MS – right panel) galaxy samples. The best-fit lines, color-coded with redshift, were derived according to Eq. (8).

4.2. Fitting the star formation rate–stellar mass function

Since the SFR-M_* relation was shown to exhibit a turnover at high masses (e.g., Lee et al. 2015; Popesso et al. 2023), we used the functional form of the SFR-M_* relation from Lee et al. (2015):

$$\begin{array}{c}\hfill \mathrm{SFR}={\mathrm{SFR}}_{\max }/(1+{(M_{\ast }/M_0)}^{-\gamma }),\end{array}$$(7)

where the SFR transitions from a power law with slope γ at masses below M₀ to an asymptotic value, SFR_max, at M_* ≫ M₀. The slope, γ, has been shown to vary with redshift, with values between 0.8 ≲ γ ≲ 1.3 (e.g., Lee et al. 2015; Tomczak et al. 2016). Because of the low dynamic range in M_*, particularly at high redshift, we decided to fix γ. With the best-fit value of γ = 1.0, adopted from Popesso et al. (2023), Eq. (7) becomes

$$\begin{array}{c}\hfill \mathrm{SFR}={\mathrm{SFR}}_{\max }/(1+(M_0/M_{\ast })).\end{array}$$(8)

We find that the most accurate fits to the SFR-M_* relation are produced when the logarithm of the free parameters of Eq. (8) are assumed to follow an exponential evolution with redshift, where

$$\begin{array}{cc}\hfill \log ({\mathrm{SFR}}_{\max }/{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1})& =a_1+a_2\times e^{-a_{3}z}\hfill \\ \hfill \log (M_0/{\mathrm{M}}_{\odot })& =b_1+b_2\times e^{-b_{3}z}.\hfill \end{array}$$(9)

The functional form of Eq. (8) was fit to the stacked data (Table A.1) using the nonlinear least square fitting, where, initially, the fits were performed for the whole sample, including both star-forming and quiescent galaxies. The best-fit curves can be seen in the left panel of Fig. 3 and the best-fit parameters from Eq. (9) for the whole sample are listed in Table 1.

As can be seen in Fig. 3, the SFR-M_* relation shows a clear turnover at high stellar masses, represented in Eq. (8) by the parameter M₀. However, since the whole sample consists of both the star-forming and quiescent galaxies, the apparent turnover may simply be a consequence of the increasing population of quiescent galaxies at low redshifts. We therefore repeated the procedure for the SFGs only. The corresponding best SFR-M_* fits are shown in the right panel of Fig. 3, with the best-fit parameters listed in Table 1. As can be seen, the turnover is still visible at low redshifts, albeit slightly less pronounced.

5. Discussion

5.1. Evolution of the star-forming main sequence

As has been discussed in previous studies (e.g., Santini et al. 2014; Wang et al. 2022), the shape of the star-forming MS is controlled by the relationship between the molecular gas mass and SFR surface densities – the SK relation. After stacking a mass-complete sample of 0.4 ≲ z ≲ 3.6, M_* > 10¹⁰ M_⊙ MS galaxies, Wang et al. (2022) found that all the MS galaxies lie on the stellar mass and redshift independent MS-only SK relation, with the slope of ∼1.13. For comparison, at high ΣM_gas, starburst galaxies exhibit two to three times higher star-forming efficiencies, producing the SK relation with slightly higher normalization (e.g., Daddi et al. 2010; Kennicutt & De Los Reyes 2021). In addition, Wang et al. (2022) showed that the FIR sizes of all the MS galaxies do not vary with stellar mass or redshift, with the mean half-light radius of ∼2.2 kpc, implying the so-called integrated SK law, where $\mathrm{SFR}\propto M_{\mathrm{gas}}^{1.13}$. The star-forming MS stellar mass and redshift evolution is therefore driven by the shape of the integrated SK relation and galaxy molecular gas content (f_gas). The molecular gas fraction, f_gas ≡ M_gas/M_*, on the other hand, is driven by many processes, including the star formation efficiency (SFE ≡ SFR/M_gas), gas cooling, feedback, and the merger rate. As has been shown by Santini et al. (2014), the redshift and stellar mass evolution of the molecular gas content in SFGs, together with the redshift-independent SK relation, give rise to the fundamental f_gas–M_*-SFR three-dimensional surface. Santini et al. (2014) stacked GOODS-S, GOODS-N and the COSMOS SFGs out to z = 2.5 in Herschel bands to show that the molecular gas fraction at given stellar mass and SFR does not depend on redshift. The projection of the f_gas–M_*-SFR surface onto the M_*-SFR plane, given f_gas variations with stellar mass and redshift, produces the evolution of the star-forming MS.

From the best-fit function describing the MS evolution below the bending mass M₀ (Eqs. (7) and (9)) at any given redshift, the SFR dependence on the stellar mass is close to linear (γ ∼ 1), consistent with the constancy of the low-mass power law slope of the stellar mass function of α ∼ 1.5 (e.g., Tomczak et al. 2014, although see Sect. 5.4). Since, from the MS-only integrated SK relation, $\mathrm{SFR}\propto M_{\mathrm{gas}}^{1.13}$, a linear dependence of the SFR on stellar mass implies that at a given redshift and below M₀, the gas fraction should decrease with stellar mass, as was indeed found in previous works (e.g., Schreiber et al. 2016; Wang et al. 2022). Similarly, since SFR ∝ M_* below M₀ and the gas fraction decreases with stellar mass, the star-forming efficiency (≡SFR/M_gas) will increase with M_*. On the other hand, for a given stellar mass, increasing the SFR with redshift (evolution of the MS normalization) should also increase the molecular gas mass (MS-only SK relation). For a given stellar mass, the gas fraction must therefore increase with redshift, which is consistent with previous findings (e.g., Geach et al. 2011; Wang et al. 2022).

At masses > M₀, the MS flattens. As has been suggested (e.g., Daddi et al. 2022; Popesso et al. 2023), the SFR-M_* turnover is most likely caused by the phasing out of cold accretion. Based on theoretical scenarios, the cold gas can accrete freely for dark matter halos with M_DM < M_shock, where M_shock (≃10^11.8 M_⊙ independent of redshift) marks the mass above which shocks can efficiently heat incoming baryons. Above M_shock, the cold accretion can still occur at higher redshifts through cold, collimated streams for halos with M_DM < M_stream, with M_stream boundary evolving with redshift from ∼10^12.5 M_⊙ at z = 2 to ∼10^13.5 M_⊙ at z = 3. Converting from halo to stellar masses using the prescription of Behroozi et al. (2013), we plot the redshift evolution of the M_shock–M_stream boundary in Fig. 4, with the M_shock, M_stream and M_shock + M_stream boundaries shown in dotted, dash-dotted, and solid blue lines, respectively. We also plot the results of this work in black, where it can be seen that the boundary is in good agreement with the redshift evolution of the SFR-M_* turnover mass M₀ (increasing exponentially toward higher redshifts; Eq. (9)). For comparison, the results of Behroozi et al. (2013) and Popesso et al. (2023) are also shown. This suggests that the flattening of the SFR-M_* relation, at a given redshift, is most likely caused by the suppression of the cold-gas feeding of the galaxy, which effectively decreases its molecular gas fraction. At the same time, Daddi et al. (2022) points out that this process is not equivalent to quenching, describing galaxies with very low SFRs, for which additional processes are required, the most likely candidates of which include mergers and AGNs.

Fig. 4. Mshock and Mstream boundaries, above which the star formation is theoretically expected to be quenched, as is described in Sect. 5.1, shown in blue, where halo masses were converted to stellar masses using the prescription of Behroozi et al. (2013). The empirical curve of Behroozi et al. (2013) is also shown with a solid green line. For comparison, the results of this work are plotted in black. Red data represents the M0 redshift evolution found in Popesso et al. (2023).

Flattening of the MS at M_* > M₀ gives a relatively small increase in SFR for a large increase in stellar mass. Provided that MS galaxies lie on the integrated SK relation, a small increase in SFR gives a small increase in molecular gas mass, causing the molecular gas fractions above M₀ to decrease significantly, consistent with the above scenario. However, Schreiber et al. (2016) found that flattening of the star-forming MS at z ∼ 1 is caused not by the decrease in the molecular gas fraction but rather by the significant drop in the star-forming efficiency, with only a small change in f_gas. For galaxy stellar masses above the bending mass, the large increase in M_* gives only a small increase in SFR. A small change in f_gas would imply that the corresponding molecular gas mass should also significantly increase. This is inconsistent with the finding of Wang et al. (2022) that all the SFGs lie on the SK relation, where $\mathrm{SFR}\propto M_{\mathrm{gas}}^{1.13}$, since a small increase in the SFR above M₀ for MS galaxies should give similarly small rise in the molecular gas mass (and hence large decrease in f_gas). A near-constant gas fraction above M₀ would therefore suggest that the MS galaxies at these large masses should lie below the SK relation. More work, with larger samples, is clearly required to resolve this inconsistency.

5.2. Temperature-redshift relation

The evolution of the dust temperature with redshift has been debated numerous times in the literature. There are three main sets of claims that can be distinguished – those that argue that, at fixed L_IR, T_d is lower at higher redshifts (e.g., Chapman et al. 2002; Pope et al. 2006; Magnelli et al. 2014), those that claim T_d to be increasing with z (e.g., Schreiber et al. 2018; Faisst et al. 2020), and those that show no evolution at all (e.g., Dudzevičiūtė et al. 2020; Reuter et al. 2020; Drew & Casey 2022). As has been shown by Drew & Casey (2022), the apparent increase in the dust temperature in lower-redshift galaxies is often caused by the bias, whereby colder systems had been fit with warmer-temperature SEDs due to the lack of long wavelength data (and hence no photometric constraints on the Rayleigh-Jeans tail of the dust emission curve). A number of works that have found no evolution of T_d with redshift are based on brightest sources selected at (sub-)mm wavelengths, not necessarily representative of the bulk of the star-forming galaxy population on the MS (e.g., Dudzevičiūtė et al. 2020; Reuter et al. 2020). Drew & Casey (2022) have found that for the sample of z < 2 galaxies, the IR luminosity anti-correlates with the dust emission rest-frame peak wavelength (and hence correlates with the luminosity-weighted dust temperature). They have also determined that the z < 2 T_d–L_IR relation does not evolve with redshift. Given that the SFRs (and L_IR) for star-forming MS galaxies, at fixed stellar mass, increase with redshift, the redshift-independent correlation of L_IR and T_d produces hotter SEDs at higher z’s.

In the left panel of Fig. 5 we plot the peak dust temperatures (derived from best-fit SEDs given in Data availability using Wien’s law) for the star-forming sample of galaxies as a function of the IR luminosity, color-coded with redshift. We can see that, as in the case of Drew & Casey (2022), T_d clearly correlates with L_IR (even when confining our sample to sources at z < 2). However, at any given redshift, we do not find any evolution of the dust temperature with L_IR. If anything, similarly to Schreiber et al. (2018), we can see a minor decline of T_d at highest IR luminosities at z ≲ 2. As is discussed in Liang et al. (2019), the differences in dust temperatures are caused mainly by the variations in the specific star formation rate (sSFR; see the right panel of Fig. 6). In this scenario, higher sSFR systems have more young (t ≲ 10 Myrs) star clusters, which heat the dust to higher temperatures, producing more emission in mid-IR wavelengths and hence increasing the values of the effective peak temperatures. This is consistent with the slight decrease in T_d at the highest stellar masses and the increase with redshift seen in the left panel of Fig. 5. Since, at any given redshift, the dust temperature variations are small, we plot the mean values of T_d as a function of z (black dots) in the right panel of Fig. 5. We find that the data can be well fit with the quadratic equation, where

Fig. 5. Evolution of peak dust temperature found in this work. Left panel: The peak dust temperatures, Td, as a function of the IR luminosity, LIR, color-coded with redshift. It can be seen that while Td increases with LIR between z = 0.45 − 5.7, it does not seem to significantly vary at any given redshift. Right panel: Since Td variations with LIR at any given redshift remain small, we plot the average peak dust temperature as a function of redshift as black points, with the best-fit solid black line. For comparison, we show the relations of Schreiber et al. (2018), Viero et al. (2022) and Ismail et al. (2023) as dotted red, dash-dotted green, and dashed blue lines, respectively. For discussion see Sect. 5.2.

$$\begin{array}{cc}\hfill T_{\mathrm{d}}^{\mathrm{MS}}/\mathrm{K}& =(1.13\pm 0.17)\times z^2+(0.88\pm 0.78)\hfill \\ \hfill & \times z+(21.74\pm 0.70).\hfill \end{array}$$(10)

Fig. 6. Logarithm of the sSFR, derived from the best-fit functional forms of the star-forming MS found in this work, plotted as a function of redshift and color-coded with stellar mass. The shape of the relations is consistent with the redshift evolution of the molecular gas fraction in SFGs, as is discussed in Sect. 5.3. For comparison, we show the findings of Popesso et al. (2023) at two stellar mass regimes. The inset plot shown the relation between the sSFR and peak dust temperature. The similar shapes of the sSFR evolution with z and Td is consistent with the scenario, presented in Liang et al. (2019), in which different dust temperature are mainly driven by the variations in the sSFR (see Sect. 5.2).

For comparison, the results of Schreiber et al. (2018), Viero et al. (2022) and Ismail et al. (2023) are shown as dotted red, dash-dotted green, and dashed blue lines, respectively. It can be seen that our results, as well as the relation of Schreiber et al. (2018), give slightly higher temperatures at low redshifts than Ismail et al. (2023). This is most likely caused by the fact that both in this work and in the work of Schreiber et al. (2018), the emissivity index, β, has been fixed at constant value (1.96 and 1.5, respectively). As is shown in Ismail et al. (2023), β correlates with dust temperature, where warmer SEDs are characterized by lower emissivity values. At low redshifts, where T_d is ∼20 K, β is expected to vary between 2.5 and 3.5. Fixing β at lower values will naturally increase the best-fit dust temperatures. The effect is stronger for lower assumed values of emissivity, which can be seen in the figure. Comparing to Viero et al. (2022), we can see that our results also point toward slightly hotter dust in the early Universe. As is explained in Viero et al. (2022), however, assuming hotter dust at higher redshifts is necessary in order to explain a number of observations, like the lack of dusty objects in the early Universe (Sommovigo et al. 2020), inconsistencies in the high-z IRX-β relation (Bouwens et al. 2016), or unusually high dust masses (Bakx et al. 2020).

5.3. Comparison to the literature

In Fig. 7, we show the comparison of the SFR-M_* relation for SFGs (star-forming MS) found in this work (solid black lines) with those found by Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023), plotted with dashed blue, dash-dotted green, and dotted red lines, respectively. The lines represent the best-fit functions at the central value of each redshift bin studied in this work. The first feature that can be seen is the lack of the turnover for the SFR-M_* relation of Speagle et al. (2014). As has been pointed out by Popesso et al. (2023), this was likely caused by the fact that the sample studied in Speagle et al. (2014) was mostly confined to sources with stellar masses below 10¹¹ M_⊙.

Fig. 7. Best-fit star-forming MS shown at each redshift bin studied in this work as the solid black line. For comparison, we show the relations of Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023), plotted with dashed blue, dash-dotted green, and dotted red lines, respectively. While all the functions tend to roughly agree with each other, some inconsistencies can clearly be seen, discussion of which is presented in Sect. 5.3.

At low redshifts, for stellar masses log(M_*/M_⊙)≲10.5, our values of the SFR are lower than those of Tomczak et al. (2016) and Popesso et al. (2023). This is most likely caused by fitting our IR data with colder SEDs, where dust temperatures of ∼20 K were assumed (see Fig. 5). At z > 2 and log(M_*/M_⊙)≳11, our values of SFR tend to be slightly larger than those of Tomczak et al. (2016) and Popesso et al. (2023). As is explained in Sect. 3.1, all the starburst galaxies that were detected with ALMA in COSMOS and UDS fields were excluded from the sample of this work. However, because the relatively long wavelength of ALMA observations tends to miss hot dust sources (e.g., Chen et al. 2022), we likely failed to identify a fraction of a high-mass starburst population that, due to its high SFRs, would boost the average stacked values. In addition, the determination of the high-mass end of the MS at high redshifts is based on a handful of sources. Looking at Fig. 3 of Popesso et al. (2023), it can be seen that at the highest redshift both the bending and the normalization of the MS are underpredicted by the adopted functional forms.

In Fig. 6, we show the corresponding evolution of the sSFR with redshift, color-coded with the stellar mass. A clear turnover can be seen at z ∼ 2, which is a direct consequence of the bending of the normalization of the MS with redshift (Eq. (9)). This, on the other hand, as is explained in Sect. 5.1, is likely caused by the increase in the molecular gas fraction, which was found by Wang et al. (2022) to be exponentially increasing with redshift. The sSFR-z relation exhibits lower normalization at higher stellar masses, since at masses ≳M₀ the MS flattens, decreasing the resulting sSFR. In the inset plot, the evolution of the sSFR with the peak dust temperature can be seen. Similar shapes of the sSFR relation with both redshift and dust temperature are consistent with the suggestion (Liang et al. 2019) that the increase in the dust temperature with redshift, found in Sect. 5.2, is indeed caused by the corresponding increase in the sSFRs.

For comparison, the relation of Popesso et al. (2023) for log(M_*/M_⊙)≃10.0 and 11.5 are shown as dotted and dashed black lines, respectively. The inconsistencies at low and high redshifts come from different normalizations of the corresponding MSs, as was explained above. In order to investigate those further, in the following section we explore how the evolution of the MS found in this work and the works of Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023) affects the “growth” of the GSMF.

5.4. Growth of the galaxy stellar mass function

The star-forming MS effectively quantifies the stellar mass growth of a galaxy with time. We can therefore compare the evolution of the GSMF as implied by any given MS with the observed one. The procedure is described in detail in Leja et al. (2015) and can be summarized as follows. At each redshift, the stellar mass for a given GSMF is being grown assuming the SFR from the MS, causing the GSMF to shift toward higher masses. Additionally, mainly due to winds and outflows, a fraction of the mass is ejected from a stellar population, a careful treatment of which has been presented in Moster et al. (2013). The stellar mass growth with time can therefore be expressed as

$$\begin{array}{c}\hfill \dot{M}(M,z)=\mathrm{SFR}(M,z)\times (1-R),\end{array}$$(11)

where R is the fraction of the mass lost during passive stellar evolution, starting at 0 and reaching 0.36 after ∼10 Gyr, assuming the Chabrier IMF. However, since most of the mass loss happens during the first 100 Myr, following Leja et al. (2015), we assume it to be instantaneous and fix R at 0.36.

We start the procedure with the observed star-forming stellar mass function of McLeod et al. (2021). The stellar mass at any given redshift is evolved according to Eq. (11) assuming different star-forming MSs. In Fig. 8 we show the results of this exercise, with the observed star-forming GSMF at the final redshift shown in gray. The stellar mass function, produced by evolving the observed GSMF between redshifts listed in the upper right corner of each panel, for the MS of Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023) are plotted with dashed blue, dash-dotted green, and dotted red lines, respectively. The results from this work are shown with a solid black line. It can be seen that the evolved GSMFs differ at both low and high-mass ends, which is due to the differences in the shapes of the corresponding MSs. At low redshifts and stellar masses the MS found in this work predicts lower values of the SFR, which causes the GSMF to grow more slowly in mass and hence the shift toward higher masses is in this case smaller than those produced by Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023). Similarly, at high redshift and stellar masses, the stellar mass shifts are slightly larger (except for the case of Speagle et al. 2014, where no bending was assumed). The most striking feature of Fig. 8 is the inconsistency of all the evolved functions with the observed GSMF. There are two main reasons for that. Firstly, the mass evolution of the GSMF is going to be affected by galaxy mergers, where mergers are expected to decrease the number density and contribute to the mass growth of galaxies. As is shown by Leja et al. (2015), however, including mergers will change the evolved GSMF at low masses by only 0.1 − 0.2 dex between redshifts 2.25 and 0.5. The second reason is the fact that, as GSMF evolves in time, a fraction of SFGs will transform into passive sources, effectively shifting the evolved stellar mass function down toward lower densities.

Fig. 8. Comparison between the observed and simulated star-forming GSMFs. The simulated GSMFs are constructed by “growing” the observed function of McLeod et al. (2021), between the redshifts listed in the upper right corner of each panel, according to Eq. (11). The mass evolution resulting from assumed star-forming MSs of Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023), are plotted with dashed blue, dash-dotted green, and dotted red lines, respectively. The simulated function constructed using the results of this work are presented in solid black, while the thick solid gray line represents the observed GSMF at the final redshift of each simulation. The inconsistencies at low- and high-mass ends of the simulated GSMFs come from differences in the corresponding shapes of the MSs involved (Fig. 7). The significant shift in densities between the observed and simulated GSMFs is caused by the fraction of SFGs, simulated between the redshifts in question, entering a quiescent mode, as is discussed in Sect. 5.4.

In order to account for this effect, we additionally evolve the GSMF, where both star-forming and passive galaxies have been included. In order to extend this exercise to higher redshifts, we adopt the stellar mass function of Davidzon et al. (2017), where z ≃ 6 has been reached. The SFR-M_* relation used in this case includes all galaxies (star-forming and quiescent) and is shown in the left panel of Fig. 3. The resulting evolved GSMF at different redshifts is compared to the observed one in Fig. 9 (thin black and thick gray lines, respectively). Since Speagle et al. (2014) and Popesso et al. (2023) did not publish the SFR-M_* relations for the galaxy samples that include the quiescent sources, we compare our results with Tomczak et al. (2016) only (dashed green line). It is clear from Fig. 9 that the obvious shift in galaxies number density between the observed and simulated functions, apparent in Fig. 8, is no longer present. This confirms that the inconsistency in the star-forming GSMFs is caused by the decrease in density due to the fraction of SFGs transforming into passive sources. For low masses at z ≳ 2 the low-mass end of the simulated stellar mass functions tends to be steeper than the observed ones. This is caused by the low-mass slope of the observed GSMF undergoing the evolution in redshift toward higher values. In order for our SFR-M_* relation to produce results consistent with this, the low-mass slope, γ, from Eq. (7) would also need to evolve. However, as is explained in Sect. 4.2, since our data does not probe low stellar masses at high redshifts, γ was assumed to be constant. The fact that the constant γ gives steeper low-mass slopes of the GSMF at high z’s indicates that the power law slope of the SFR-M_* relation for all galaxies (star-forming and quiescent) should increase in redshift. Similarly, since the low-mass slope of the star-forming GSMF, found in both Davidzon et al. (2017) and McLeod et al. (2021), increases at z ≳ 2, the slope of the star-forming MS should also evolve. At z ≳ 3 the high-mass end of the simulated GSMF overpredicts the number of galaxies. While the errors on the observed high-redshift stellar mass function at these masses are significantly large (∼1 dex), it is also possible that the high-mass end of the SFR-M_* relations in this work, and hence the growth of the stellar masses, is contaminated by the presence of starburst galaxies.

Fig. 9. Comparison between the observed and simulated GSMFs, where both star-forming and quiescent galaxies are considered. To extend the redshift range, the simulated GSMFs are constructed by “growing” the observed function of Davidzon et al. (2017), between the redshifts listed in the upper right corner of each panel, according to Eq. (11). The mass evolution results from assuming the SFR-M* relation characteristic of galaxy samples including both the star-forming and quiescent sources, where we compare our results, shown in solid black, with those of Tomczak et al. (2016), depicted with a dash-dotted green line. The thick solid gray line represents the observed GSMF at the final redshift of each simulation. It can be seen that the significant gap in densities, apparent in Fig. 8, is no longer present, indicating that, indeed, the inconsistencies were caused by the fraction of sources evolving into quiescent mode, as is discussed in Sect. 5.4. The low-mass slope of simulated GSMFs at z ≳ 2 is inconsistent with the observed function, indicating that the corresponding slope of the SFR-M* relations most likely evolves toward steeper values toward earlier epochs.

6. Summary

We have stacked ∼100k K-band selected COSMOS and UDS galaxies in the IR Herschel and JCMT bands in bins of stellar mass and redshift in order to investigate the SFR evolution out to z ∼ 7. Inverse variance stacking was used, where, for the Herschel SPIRE bands, background sources were found to contribute significantly to the stacked fluxes. After removing the background sources’ contribution (fitted with the Gaussian function), the IR luminosities, and hence the SFRs, were determined and the functional forms of the SFR evolution with stellar mass and redshift found. This was done for both the star-forming and whole (star-forming and quiescent) galaxy samples, where the quiescent sources were selected using the UVJ color-color space. The main findings of this work can be summarized as follows:

• (i)

  The SFR of SFGs populating the MS exhibits a nearly linear growth with stellar mass at any given redshift at low masses. Above the so-called bending mass, M₀, the MS flattens toward an asymptotic value of the SFR. At constant stellar mass, the SFR also increases with redshift – the redshift evolution of the normalization of the MS. This behavior is consistent with the scenario in which the evolution of the star-forming MS is driven by the position of the galaxy on the integrated SK relation, given its molecular mass fraction (e.g., Santini et al. 2014; Wang et al. 2022). We find that the redshift evolution of the bending mass, M₀, is exponential, increasing sharply from z ∼ 0.5 to z ∼ 2 and flattening at z > 2. The similarity to the shape of the relation between the stellar mass, above which shocks can efficiently heat incoming baryons, and redshift, found in Daddi et al. (2022), suggests that the flattening of the MS is likely caused by the phasing out of the cold accretion.

• (ii)

  We investigated the evolution of the peak dust temperature with IR luminosity and redshift, by fitting the stacked photometry with the dust emission curve of Casey (2012). Similarly to Drew & Casey (2022), an increase in T_d with IR luminosity was found. However, when confined to any given redshift, our data does not show any significant evolution with L_IR. We conclude that the increase in T_d with IR luminosity is driven by the evolution of the L_IR with redshift for the MS galaxies. We find that the dust temperature increases quadratically with redshift, consistent with the scenario in which higher values of T_d for SFGs are driven by the increase in their sSFRs (e.g., Liang et al. 2019).

• (iii)

  The evolution of the MS found in this work was compared with the relations of Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023). We find that the normalization of our MS at z ≲ 2 for masses below the bending mass, M₀, is slightly lower, while the SFRs above M₀ at high redshifts exhibit larger values (except for the MS of Speagle et al. 2014, where no MS bending was found). We attribute these inconsistencies to the dust temperatures assumed, when fitting the IR data. As was found in this work, T_d increases from ∼20 K to ∼60 K between redshifts 0.5 and 5.5, and hence assuming dust temperatures at the level of ∼30 K, as is often done in the absence of the FIR data, will correspondingly affect the determined IR luminosities (and hence the SFRs). We also note that the higher values of the SFRs at high redshifts for M_* > 10¹¹ M_⊙ could be, at least in part, caused by the contamination of the MS sample by starburst galaxies.

• (iv)

  Adopting different shapes of the SFR-M_* relation, we simulated the time evolution for the observed GSMFs of Davidzon et al. (2017) and McLeod et al. (2021). For the star-forming sample, we find that the GSMF “grows” too fast, producing simulated relations significantly above the observed ones. By performing the simulations using the SFR-M_* relation, found using both star-forming and quiescent galaxies, we confirm that the overprediction of SFG number densities is caused by a significant fraction of sources evolving into the quiescent phase and thereby leaving the star-forming MS. In addition, we find the simulated low-mass end of the GSMF gives higher number densities at high z, indicating that the low-mass slope of the MS likely evolves toward steeper values with redshift. However, due to the significant uncertainties, more detailed analysis based on larger high-z samples is required.

Data availability

Example stacked image and best-fit SEDs are available here.

Acknowledgments

This research was funded in whole or in part by the National Science Center, Poland (grant no. 2020/39/D/ST9/03078). For the purpose of Open Access, the author has applied a CC-BY public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission. JSD acknowledges the support of the Royal Society through a Research Professorship. MJM acknowledges the support of the National Science Centre, Poland through the SONATA BIS grant 2018/30/E/ST9/00208 and the Polish National Agency for Academic Exchange Bekker grant BPN/BEK/2022/1/00110. DJM acknowledges the support of the Science and Technology Facilities Council. KL acknowledges the support of the Polish Ministry of Education and Science through the grant PN/01/0034/2022 under ‘Perły Nauki’ program.

References

Appendix A: Additional table

All Tables

All Figures

Fig. 1. Results of the NUVrJ color-color selection results, defined in Eq. (2) shown for redshift bins adopted in this work. The NUVrJ colors (Ilbert et al. 2013) are used in order to distinguish between sources reddened due to the lack of star formation (quiescent galaxies – upper left section of the plots) from those where red colors are caused by dust extinction (SFGs – main locus in the plots).

In the text

Fig. 2. 1.25 < z < 1.75, 10.25 < log(M*/M⊙) < 10.5 light profile derived at 500 μm. The red points represent mean values of the flux densities derived from stacked images as a function of distance from the image center. The dashed blue, dotted green, and solid red lines represent the corresponding beam, the contribution from background sources, and the sum of the two, respectively. For Herschel SPIRE bands, the contribution from background sources is significant and must therefore be subtracted from the overall light profile, when deriving stacked flux densities (see Sect. 3.2).

In the text

	Fig. 3. SFR-M* relations derived in Sect. 4 for the whole (star-forming and quiescent – left panel) and the star-forming (MS – right panel) galaxy samples. The best-fit lines, color-coded with redshift, were derived according to Eq. (8).
In the text

Fig. 4. Mshock and Mstream boundaries, above which the star formation is theoretically expected to be quenched, as is described in Sect. 5.1, shown in blue, where halo masses were converted to stellar masses using the prescription of Behroozi et al. (2013). The empirical curve of Behroozi et al. (2013) is also shown with a solid green line. For comparison, the results of this work are plotted in black. Red data represents the M0 redshift evolution found in Popesso et al. (2023).

In the text

Fig. 5. Evolution of peak dust temperature found in this work. Left panel: The peak dust temperatures, Td, as a function of the IR luminosity, LIR, color-coded with redshift. It can be seen that while Td increases with LIR between z = 0.45 − 5.7, it does not seem to significantly vary at any given redshift. Right panel: Since Td variations with LIR at any given redshift remain small, we plot the average peak dust temperature as a function of redshift as black points, with the best-fit solid black line. For comparison, we show the relations of Schreiber et al. (2018), Viero et al. (2022) and Ismail et al. (2023) as dotted red, dash-dotted green, and dashed blue lines, respectively. For discussion see Sect. 5.2.

In the text

Fig. 6. Logarithm of the sSFR, derived from the best-fit functional forms of the star-forming MS found in this work, plotted as a function of redshift and color-coded with stellar mass. The shape of the relations is consistent with the redshift evolution of the molecular gas fraction in SFGs, as is discussed in Sect. 5.3. For comparison, we show the findings of Popesso et al. (2023) at two stellar mass regimes. The inset plot shown the relation between the sSFR and peak dust temperature. The similar shapes of the sSFR evolution with z and Td is consistent with the scenario, presented in Liang et al. (2019), in which different dust temperature are mainly driven by the variations in the sSFR (see Sect. 5.2).

In the text

Fig. 7. Best-fit star-forming MS shown at each redshift bin studied in this work as the solid black line. For comparison, we show the relations of Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023), plotted with dashed blue, dash-dotted green, and dotted red lines, respectively. While all the functions tend to roughly agree with each other, some inconsistencies can clearly be seen, discussion of which is presented in Sect. 5.3.

In the text

Fig. 8. Comparison between the observed and simulated star-forming GSMFs. The simulated GSMFs are constructed by “growing” the observed function of McLeod et al. (2021), between the redshifts listed in the upper right corner of each panel, according to Eq. (11). The mass evolution resulting from assumed star-forming MSs of Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023), are plotted with dashed blue, dash-dotted green, and dotted red lines, respectively. The simulated function constructed using the results of this work are presented in solid black, while the thick solid gray line represents the observed GSMF at the final redshift of each simulation. The inconsistencies at low- and high-mass ends of the simulated GSMFs come from differences in the corresponding shapes of the MSs involved (Fig. 7). The significant shift in densities between the observed and simulated GSMFs is caused by the fraction of SFGs, simulated between the redshifts in question, entering a quiescent mode, as is discussed in Sect. 5.4.

In the text

Fig. 9. Comparison between the observed and simulated GSMFs, where both star-forming and quiescent galaxies are considered. To extend the redshift range, the simulated GSMFs are constructed by “growing” the observed function of Davidzon et al. (2017), between the redshifts listed in the upper right corner of each panel, according to Eq. (11). The mass evolution results from assuming the SFR-M* relation characteristic of galaxy samples including both the star-forming and quiescent sources, where we compare our results, shown in solid black, with those of Tomczak et al. (2016), depicted with a dash-dotted green line. The thick solid gray line represents the observed GSMF at the final redshift of each simulation. It can be seen that the significant gap in densities, apparent in Fig. 8, is no longer present, indicating that, indeed, the inconsistencies were caused by the fraction of sources evolving into quiescent mode, as is discussed in Sect. 5.4. The low-mass slope of simulated GSMFs at z ≳ 2 is inconsistent with the observed function, indicating that the corresponding slope of the SFR-M* relations most likely evolves toward steeper values toward earlier epochs.

In the text