11institutetext:  Institute of Astronomy, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University, Grudzia̧dzka 5, 87-100 Toruń,
Poland<EMAIL_ADDRESS>22institutetext: Institute for Astronomy,
University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK
33institutetext: Astronomical Observatory Institute, Faculty of Physics, Adam
Mickiewicz University, ul. Słoneczna 36, 60-286 Poznań, Poland
44institutetext: National Centre for Nuclear Research, Pasteura 7, 093,
Warsaw, Poland

# Charting the main sequence of star-forming galaxies out to redshifts $z\sim
7$

M. P. Koprowski Charting the main sequence of star-forming galaxies out to
redshifts $z\sim 7$ Charting the main sequence of star-forming galaxies out to
redshifts $z\sim 7$ J. V. Wijesekera Charting the main sequence of star-
forming galaxies out to redshifts $z\sim 7$ Charting the main sequence of
star-forming galaxies out to redshifts $z\sim 7$ J. S. Dunlop Charting the
main sequence of star-forming galaxies out to redshifts $z\sim 7$ Charting the
main sequence of star-forming galaxies out to redshifts $z\sim 7$ D. J. McLeod
Charting the main sequence of star-forming galaxies out to redshifts $z\sim 7$
Charting the main sequence of star-forming galaxies out to redshifts $z\sim 7$
M. J. Michałowski Charting the main sequence of star-forming galaxies out to
redshifts $z\sim 7$ Charting the main sequence of star-forming galaxies out to
redshifts $z\sim 7$ Charting the main sequence of star-forming galaxies out to
redshifts $z\sim 7$ Charting the main sequence of star-forming galaxies out to
redshifts $z\sim 7$ K. Lisiecki Charting the main sequence of star-forming
galaxies out to redshifts $z\sim 7$ Charting the main sequence of star-forming
galaxies out to redshifts $z\sim 7$ Charting the main sequence of star-forming
galaxies out to redshifts $z\sim 7$ Charting the main sequence of star-forming
galaxies out to redshifts $z\sim 7$ R. J. McLure Charting the main sequence of
star-forming galaxies out to redshifts $z\sim 7$ Charting the main sequence of
star-forming galaxies out to redshifts $z\sim 7$

We present a new determination of the star-forming main sequence (MS),
obtained through stacking 100k $K$-band-selected galaxies in the far-IR
Herschel and James Clerk Maxwell Telescope (JCMT) imaging. By fitting the dust
emission curve to the stacked far-IR photometry, we derive the IR luminosities
($L_{\rm IR}$) and, hence, star formation rates (SFR) out to $z\sim 7$. The
functional form of the MS is found, with the linear SFR-$M_{\ast}$ relation
that flattens at high stellar masses and the normalization that increases
exponentially with redshift. We derive the corresponding redshift evolution of
the specific star formation rate (sSFR) and compare our findings with the
recent literature. We find our MS to be exhibiting slightly lower
normalization at $z\lesssim 2$ and to flatten at larger stellar masses at high
redshifts. By deriving the relationship between the peak dust temperature
($T_{\rm d}$) and redshift, where $T_{\rm d}$ increases linearly from $\sim
20$ K at $z=0.5$ to $\sim 50$ K at $z=6$, we conclude that the apparent
inconsistencies in the shapes of the MS are most likely caused by the
different dust temperatures assumed when deriving SFRs in the absence of far-
IR data. Finally, we investigate the derived shape of the star-forming MS by
simulating the time evolution of the observed galaxy stellar mass function
(GSMF). While the simulated GSMF is in good agreement with the observed one,
some inconsistencies persist. In particular, we find the simulated GSMF to be
somewhat overpredicting the number density of low-mass galaxies at $z\gtrsim
2$.

###### Key Words.:

dust, extinction – galaxies: ISM – galaxies: evolution – galaxies: star
formation – galaxies: high-redshift

##  1 Introduction

The increasing amount of multi-wavelength 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
and quiescent (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, where it has been shown that the
SFGs follow a tight correlation between their star formation rate, SFR, and
stellar mass, $M_{\ast}$, 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_{\ast}$
relation can be parametrized with a simple power-law:

${\rm log(SFR/\hbox{$\rm\thinspace M_{\odot}$}\,{\rm yr}^{-1})}=\alpha\,{\rm
log}(M_{\ast}/\hbox{$\rm\thinspace M_{\odot}$})+\beta,$ (1)

where $\alpha$ and $\beta$ represent the slope and normalization of the main
sequence, respectively. While the normalization has been shown to increase in
redshift, where galaxies exhibit higher SFRs at earlier epochs, the slope of
the main sequence was found to be close to unity, mostly independent of
redshift, with values ranging between $\sim 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 found that the star-forming main sequence 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 main sequence 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 (SPS) modeling, star
formation histories (SFHs), dust extinction curves, etc. 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-IR and, therefore, more reliable
determination of the total (UV+IR) SFRs. However, since Herschel far-IR data
suffers from relatively low resolution, only the most IR-luminous sources can
be detected at high redshifts and, hence, samples of star-forming galaxies
studied with Herschel tend to be incomplete, often biasing the results towards
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 modelling 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 $\sim 100$k $K$-band-selected
galaxies from the ground-based UV+optical+near-IR $\sim 2$ deg2 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 far-IR
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 main sequence out to $z\sim 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 ${\rm\mu m}$ data, our stacked photometry maps the
Rayleigh–Jeans portion of the dust emission spectrum, allowing for a precise
determination of the dust temperatures and, thus, star formation rates.

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

##  2 Data

###  2.1 Optical/near-IR catalogs

We start 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 preparation) 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^{\ast}$-band from CFHT MegaCam, and
optical imaging in BVRiz′ from Subaru Suprime-Cam (Furusawa et al. 2008). We
also employed imaging in the narrow-band NB921 (Koyama et al. 2011) and the
refurbished Suprime-Cam $z$-band ($z^{\prime}_{\rm new}$; Furusawa et al.
2016). To extend our wavelength coverage into the mid-infrared, we utilised
Spitzer IRAC imaging in 3.6$\mu$m and 4.5$\mu$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
deg2, after masking for bright stars.

For the COSMOS field, we combined near-IR $YJHK_{\rm s}$ imaging from
UltraVISTA DR4 (McCracken et al. 2012) with UV/optical $u^{\ast}griz$ imaging
from the CFHTLS-D2 T0007 data release (Hudelot et al. 2012). We confine our
study to the 1 deg2 overlapping area covered by these data sets. As with the
UDS, we also included Subaru Suprime-Cam z${}^{\prime}_{new}$ imaging. There
is a wealth of Spitzer IRAC imaging available in 3.6$\mu$m and 4.5$\mu$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

We constructed our parent multi-wavelength catalogs 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 ${\rm\mu 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.

Stellar masses were measured for each object by fixing the photometric
redshift to $z_{\rm med}$, and re-fitting the SED using LePhare with the
Bruzual & Charlot (2003) library. This template set included a Chabrier (2003)
initial mass function, a Calzetti et al. (2000) dust attenuation law and IGM
absorption as in Madau (1995). The resulting catalogs were cleaned of stars,
AGN, artefacts and objects with contaminated photometry. Finally, a further
$\chi^{2}<40$ goodness-of-fit cut was implemented on the catalogs.

###  2.2 Far-IR

For the extraction of the stacked far-IR 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 ${\rm\mu m}$ with beam sizes of 7.39, 11.29,
18.2, 24.9, and 36.3 arcsec, and 5$\sigma$ 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
${\rm\mu m}$ imaging covering $\simeq 0.9$ deg2 with a 1$\sigma$ noise of 0.9
mJy and $\simeq 1.3$ deg2 with the 1$\sigma$ noise of 1.6 mJy, respectively.

##  3 Stacking

###  3.1 Sample selection

Since we use the derived SFR-$M_{\ast}$ relations to ‘grow’ the McLeod et al.
(2021) galaxy stellar mass function in Section 5.4, we decided to follow them
in choosing the mass and redshift bins, while at the same time extending our
analysis to $z=6.75$. At $z>3.75$, in order to be able to reach the detection
threshold, the mass and redshift bins have been accordingly stretched. The
procedure for estimating the mass completeness in the UDS and COSMOS fields is
explained in McLeod et al. (2021). In brief, the mass completeness was
estimated using simulations, where $K$-band imaging was produced by injecting
a set of artificial galaxies with a wide range of physical properties.
Adopting the original set-up, the corresponding photometry catalogs were then
constructed using sextractor and the 90% mass-completeness limit, defined as
the stellar mass below which 90% of the limiting stellar masses lie, at each
redshift was derived. The resulting redshift and stellar mass bins adopted in
this work are listed in Table 2.

In order to derive the SFR-$M_{\ast}$ relation for star-forming sample, we
first need to separate it from the quiescent population. Similarly to McLeod
et al. (2021) we identify passive galaxies using the UVJ color-color selection
(Williams et al. 2009) from Carnall et al. (2020):

$\displaystyle\begin{split}(U-V)>&\,0.88\times(V-J)+0.69;\\\ (U-V)>&\,1.3;\\\
(V-J)<&\,1.6,\end{split}$ (2)

where $U$,$V$ and $J$ rest-frame magnitudes were generated from the best-fit
SEDs as explained in McLeod et al. (2021). The advantage of this method is
that it is able to separate sources where reddening is caused by aging from
those where 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 UVJ color-color selection
was performed out to $z\sim 4$, as depicted in Figure 1.

Figure 1: The results of the $UVJ$ color-color selection results, defined in
Equation 2 shown for redshift bins adopted in this work. The $UVJ$ colors
(Williams et al. 2009) 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 (star-forming galaxies – main locus in the plots).

###  3.2 Far-IR fluxes

Since current far-IR/submm surveys struggle to reach the (sub-)mJy depths,
necessary for the detection of the dust emission in typical star-forming
galaxies, we resort 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 far-IR fluxes at a given band, we cut a
square stamp around the position of each source, with the size of
approximately 3 $\times$ the full width at half maximum (FWHM) of the
corresponding beam. The average far-IR fluxes at each pixel are then evaluated
following the inverse-variance weighting:

$S_{\nu}=\frac{\sum_{i}S_{\nu,i}/\sigma_{i}^{2}}{\sum_{i}1/\sigma_{i}^{2}},$
(3)

where $S_{\nu,i}$ is the flux density of the $i$th source and $\sigma_{i}$ is
the 1$\sigma$ instrumental noise, extracted at each pixel in a given source’s
stamp. The error on the mean stacked fluxes for each pixel is calculated
according to:

$\sigma^{2}(S_{\nu})=\frac{1}{\sum_{i}1/\sigma_{i}^{2}},$ (4)

with the example 500 ${\rm\mu m}$ stamps for each redshift and stellar mass
bin shown in Figure 9.

The stacked average far-IR flux densities at each bin were then extracted from
the non-linear 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 ${\rm\mu 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 ${\rm\mu m}$ sources at
$1.25<z<1.75$ with stellar masses between $10.25<{\rm log}(M_{\ast}/{\rm
M_{\odot}})<10.5$. As can be seen, the stacked flux (red line) consists of the
average flux from the stacked sample (blue dashed line) and the background
flux from the associated sources (green dotted line), both of which can be
modeled with a Gaussian function. In cases when a double Gaussian (${\rm
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. The flux density
errors are then derived from the corresponding covariance matrices.

Figure 2: The $1.25<z<1.75$, $10.25<{\rm log}(M_{\ast}/{\rm M_{\odot}})<10.5$
light profile derived at 500 ${\rm\mu 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 blue dashed, green dotted and red solid
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 Section 3.2).

##  4 The SFR-stellar mass relations

###  4.1 SFR measurements

A total star formation rate 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 the recent, unobscured star formation, while the IR emission
accounts for the stellar emission that has been absorbed by the interstellar
dust, where:

${\rm SFR_{tot}}=\mathcal{K}_{\rm FUV}L_{\rm FUV}+\mathcal{K}_{\rm IR}L_{\rm
IR}.$ (5)

The far-UV (FUV) part of the SFR is derived from the FUV luminosity measured
at 150 nm, where $L_{\rm FUV}=\nu L_{\nu}$ and $\mathcal{K}_{\rm
FUV}=2.5\times 10^{-10}\,{\rm M_{\odot}\,yr^{-1}\,L_{\odot}^{-1}}$ (Kennicutt
1998), with an additional multiplicative factor of 0.63 to convert from a
Salpeter to a Chabrier IMF. The FUV SFRs were calculated for each galaxy and
then the mean value was evaluated at each redshift/stellar mass bin. In order
to find the IR portion of the SFR we follow Madau & Dickinson (2014), where
the IR luminosity is converted to the SFR using $\mathcal{K}_{\rm
IR}=1.73\times 10^{-10}\,{\rm M_{\odot}\,yr^{-1}\,L_{\odot}^{-1}}$, again with
an additional multiplicative factor of 0.63 to convert from a Salpeter to a
Chabrier IMF.

Figure 3: SFR-$M_{\ast}$ relations derived in Section 4 for the whole (star-
forming and quiescent – left panel) and the star-forming (main sequence –
right panel) galaxy samples. The best-fit lines, color-coded with redshift,
are derived according to Equation 7.

In order to find the IR luminosity, $L_{\rm IR}$, the dust emission curve of
Casey (2012) was fit to the available far-IR photometry using the non-linear
least squares analysis. The adopted dust emission curve, described in Table 1
of Casey (2012), can be approximated as a classic single-temperature greybody
function at far-IR wavelengths, with the addition of the power-law fit at mid-
IR wavelengths. The data was initially fit with both the single greybody and
the power-law greybody curves and it was found that the latter gives better
fits to the data at mid-IR bands. The power-law greybody curve has 4 free
parameters – mid-IR power-law slope, $\alpha$, the emissivity index, $\beta$,
the effective dust temperature, $T^{\rm eff}_{\rm d}$, and the overall
normalization, $N$. Since, at $z\gtrsim 2$, the available photometry redwards
of the dust emission peak consists only of the 850 and 500 ${\rm\mu m}$, with
the 500 ${\rm\mu m}$ being affected by blending (see Figure 2), we decided to
fix $\beta$ at 1.96. Similarly, due to the lack of photometry at rest-frame
wavelength $\lambda\lesssim 80\,{\rm\mu m}$, we also fix $\alpha$ at 2.3.
Values of $\alpha$ and $\beta$ were fixed according to the recommendations of
Drew & Casey (2022). The curve is then integrated between 8 and 1000 ${\rm\mu
m}$ and the errors are 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 shown in Figures 10 and 11,
respectively. The resulting best-fit values of the total SFR at each
redshift/stellar mass bin for the whole and star-forming subsamples are listed
in the top and bottom panels of Table 2, as well as depicted as color points
with error bars in the left and right panels of Figure 3, respectively.

###  4.2 Fitting the SFR-stellar mass function

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

${\rm SFR}={\rm SFR_{max}}/(1+(M_{\ast}/M_{0})^{-\gamma}),$ (6)

where the SFR transitions from a power-law with slope $\gamma$ at masses below
$M_{0}$, to an asymptotic value, ${\rm SFR_{max}}$, at $M_{\ast}\gg M_{0}$.
The slope, $\gamma$, has been shown to vary with redshift, with values between
$0.8\lesssim\gamma\lesssim 1.3$ (e.g. Lee et al. 2015; Tomczak et al. 2016).
Because of the low dynamic range in $M_{\ast}$, particularly at high redshift,
we decided to fix $\gamma$. With the best-fit value of $\gamma=1.0$, adopted
from Popesso et al. (2023), Equation 6 becomes:

${\rm SFR}={\rm SFR_{max}}/(1+(M_{0}/M_{\ast})).$ (7)

We find that the most accurate fits to the SFR-$M_{\ast}$ relation are
produced when the logarithm of the free parameters of Equation 7 are assumed
to follow an exponential evolution with redshift, where:

$\begin{split}{\rm log(SFR_{max}/M_{\odot}yr^{-1})}&=a_{1}+a_{2}\times
e^{-a_{3}z}\\\ {\rm log}(M_{0}/{\rm M_{\odot}})&=b_{1}+b_{2}\times
e^{-b_{3}z}.\end{split}$ (8)

Table 1: Best-fit values to parameters from Equation 8 for the star-forming and all galaxies (star-forming and quiescent). Parameter | All galaxies | SF galaxies
---|---|---
a1 | $\phantom{-1}3.69\pm 0.32$ | $\phantom{-1}3.47\pm 0.22$
a2 |  1$-3.81\pm 0.27$ |  1$-3.13\pm 0.16$
a3 | $\phantom{-1}0.47\pm 0.07$ | $\phantom{-1}0.56\pm 0.08$
b1 | $\phantom{-}11.91\pm 0.45$ | $\phantom{-}11.69\pm 0.31$
b2 |  1$-2.48\pm 0.38$ |  1$-1.66\pm 0.24$
b3 | $\phantom{-1}0.44\pm 0.15$ | $\phantom{-1}0.53\pm 0.21$

The functional form of Equation 7 is fitted to the stacked data (Table 2)
using the non-linear least square fitting, where, initially, the fits are
performed for the whole sample, including both star-forming and quiescent
galaxies. The best-fit curves can be seen in the left panel of Figure 3 and
the best-fit parameters from Equation 8 for the whole sample are listed in
Table 1.

As can be seen in Figure 3, the SFR-$M_{\ast}$ relation shows a clear turnover
at high stellar masses, represented in Equation 7 by the parameter $M_{0}$.
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,
repeat the procedure for the star-forming galaxies only. The corresponding
best SFR-$M_{\ast}$ fits are shown in the right panel of Figure 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

Figure 4: Left panel The peak dust temperatures, $T_{\rm d}$, as a function of
the IR luminosity, $L_{\rm IR}$, color-coded with redshift. It can be seen
that while $T_{\rm d}$ increases with $L_{\rm IR}$ between $z\sim 0.5-6$, it
does not seem to significantly vary at any given redshift. Right panel Since
$T_{\rm d}$ variations with $L_{\rm IR}$ 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 black solid line. For comparison, we show the
relations of Schreiber et al. (2018), Viero et al. (2022) and Ismail et al.
(2023) as red dotted, green dash-dotted and blue dashed lines, respectively.
For discussion see Section 5.2.

###  5.1 Evolution of the star-forming main sequence

As discussed in previous studies (e.g. Santini et al. 2014; Wang et al. 2022),
the shape of the star-forming main sequence is controlled by the relationship
between the molecular gas mass and star formation rate surface densities – the
Schmidt-Kennicutt relation. After stacking a mass-complete sample of
$0.4\lesssim z\lesssim 3.6$, $M_{\ast}>10^{10}\,{\rm M_{\odot}}$ main-sequence
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 $\sim
1.13$. For comparison, at high $\Sigma M_{\rm 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 $\sim 2.2$ kpc, implying the so-called
integrated Schmidt-Kennicutt law, where ${\rm SFR}\propto M_{\rm gas}^{1.13}$.
The star-forming main-sequence stellar mass and redshift evolution is,
therefore, driven by the shape of the integrated SK relation and galaxy
molecular gas content ($f_{\rm gas}$). The molecular gas fraction, $f_{\rm
gas}\equiv M_{\rm gas}/M_{\ast}$, on the other hand, is driven by many
processes, including the star formation efficiency (${\rm SFE\equiv
SFR}/M_{\rm gas}$), gas cooling, feedback, merger rate, etc. As shown by
Santini et al. (2014), the redshift and stellar mass evolution of the
molecular gas content in star-forming galaxies, together with the redshift-
independent SK relation, give rise to the fundamental $f_{\rm
gas}$-$M_{\ast}$-SFR three-dimensional surface. Santini et al. (2014) stacked
GOODS-S, GOODS-N and the COSMOS star-forming galaxies 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_{\rm
gas}$-$M_{\ast}$-SFR surface onto the $M_{\ast}$-SFR plane, given $f_{\rm
gas}$ variations with stellar mass and redshift, produces the evolution of the
star-forming main sequence.

From the best-fit function describing the MS evolution below the bending mass
$M_{0}$ (Equations 6 and 8) at any given redshift, the SFR dependence on the
stellar mass is close to linear ($\gamma\sim 1$), consistent with the
constancy of the low-mass power-law slope of the stellar mass function of
$\alpha\sim 1.5$ (e.g. Tomczak et al. 2014, although see Section 5.4). Since,
from the MS-only integrated SK relation, ${\rm SFR}\propto M_{\rm
gas}^{1.13}$, linear dependence of SFR on stellar mass implies that at given
redshift and below $M_{0}$, 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 ${\rm SFR}\propto M_{\ast}$ below $M_{0}$ and
the gas fraction decreases with stellar mass, the star-forming efficiency
($\equiv{\rm SFR}/M_{\rm gas}$) will increase with $M_{\ast}$. On the other
hand, for a given stellar mass, increasing the SFR with redshift (evolution of
the main-sequence 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, consistently with previous findings (e.g.
Geach et al. 2011; Wang et al. 2022).

At masses $>M_{0}$ the main sequence flattens. As suggested (e.g. Daddi et al.
2022; Popesso et al. 2023), the SFR-$M_{\ast}$ 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 (DM) halos with $M_{\rm DM}<M_{\rm
shock}$, where $M_{\rm shock}$ ($\simeq 10^{11.8}\,{\rm M_{\odot}}$
independent of redshift) marks the mass above which shocks can efficiently
heat incoming baryons. Above $M_{\rm shock}$ the cold accretion can still
occur at higher redshifts through cold, collimated streams for halos with
$M_{\rm DM}<M_{\rm stream}$, with $M_{\rm stream}$ boundary evolving with
redshift from $\sim 10^{12.5}\,{\rm M_{\odot}}$ at $z=2$ to $\sim
10^{13.5}\,{\rm M_{\odot}}$ at $z=3$. Converting from halo to stellar masses
using the prescription of Behroozi et al. (2013), Daddi et al. (2022) shows in
their Figure 2 that the $M_{\rm shock}$-$M_{\rm stream}$ boundary is in
remarkable agreement with the redshift evolution of the SFR-$M_{\ast}$
turnover mass $M_{0}$ (increasing exponentially towards higher redshifts;
Equation 8). This suggests that the flattening of the SFR-$M_{\ast}$ 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 amounts of SFR, for
which additional processes are required, most likely candidates of which
include mergers and AGNs. Flattening of the main sequence at $M_{\ast}>M_{0}$
gives a relatively small increase in SFR for large increase in stellar mass.
Provided that main-sequence galaxies lie on the integrated SK relation, small
increase in SFR gives small increase in molecular gas mass, causing the
molecular gas fractions above $M_{0}$ to decrease significantly, consistent
with the above scenario. However, Schreiber et al. (2016) found that
flattening of the star-forming main sequence at $z\sim 1$ is not caused by the
decrease of the molecular gas fraction but rather the significant drop in the
star-forming efficiency, with only a small change in the $f_{\rm gas}$. For
galaxy stellar masses above the bending mass, the large increase in $M_{\ast}$
gives only a small increase in SFR. A small change in $f_{\rm 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 star-forming galaxies lie on the SK relation, where ${\rm SFR}\propto
M_{\rm gas}^{1.13}$, since small increase in SFR above $M_{0}$ for main-
sequence galaxies should give similarly small rise in the molecular gas mass
(and hence large decrease in $f_{\rm gas}$). A near-constant gas fraction
above $M_{0}$ would, therefore, suggest that the main-sequence 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_{\rm IR}$, $T_{\rm d}$ is
lower at higher redshifts (e.g. Chapman et al. 2002; Pope et al. 2006;
Magnelli et al. 2014), those that claim $T_{\rm 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 shown by Drew & Casey (2022), the apparent increase of the
dust temperature at 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 find
no evolution of $T_{\rm 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 main sequence (e.g. Dudzevičiūtė
et al. 2020; Reuter et al. 2020). Drew & Casey (2022) find 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 also determine that the $z<2$ $T_{\rm d}$-$L_{\rm IR}$
relation does not evolve with redshift. Given that the star formation rates
(and $L_{\rm IR}$) for star-forming main-sequence galaxies, at fixed stellar
mass, increase with redshift, the redshift-independent correlation of $L_{\rm
IR}$ and $T_{\rm d}$ produces hotter SEDs at higher $z$’s.

Figure 5: The best-fit star-forming main sequence shown at each redshift bin
studied in this work as the black solid line. For comparison, we show the
relations of Speagle et al. (2014), Tomczak et al. (2016) and Popesso et al.
(2023), plotted with blue dashed, green dash-dotted and red dotted 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
Section 5.3.

In the left panel of Figure 4 we plot the peak dust temperatures (derived from
best-fit SEDs shown in Figure 11 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_{\rm d}$ clearly
correlates with $L_{\rm 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_{\rm IR}$. If anything, similarly to Schreiber et al.
(2018), we can see a minor decline of $T_{\rm d}$ at highest IR luminosities
at $z\lesssim 2$. As discussed in Liang et al. (2019), the differences in dust
temperatures are caused mainly by the variations of the specific star
formation rate (sSFR; see the right panel of Figure 6). In this scenario,
higher sSFR systems have more young ($t\lesssim 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 of $T_{\rm d}$ at
highest stellar masses and the increase with redshift seen in the left panel
of Figure 4. Since, at any given redshift, the dust temperature variations are
small, we plot the mean values of $T_{\rm d}$ as a function of $z$ (black
dots) in the right panel of Figure 4. A clear linear relation can be seen,
with the best-fit empirical formula:

$T^{\rm MS}_{\rm d}/{\rm K}=(6.0\pm 0.3)\times z+(17.8\pm 1.0)$ (9)

shown as a black solid line. For comparison, the results of Schreiber et al.
(2018), Viero et al. (2022) and Ismail et al. (2023) are shown as red dotted,
green dash-dotted and blue dashed 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, $\beta$, has been fixed at constant value
(1.96 and 1.5, respectively). As shown in Ismail et al. (2023), $\beta$
correlates with dust temperature, where warmer SEDs are characterized by lower
emissivity values. At low redshifts, where $T_{\rm d}$ is $\sim 20$ K, $\beta$
is expected to vary between 2.5 and 3.5. Fixing $\beta$ 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. It is
also most obvious for low redshifts, as at higher $z$’s, where the dust
temperature increases, the emissivity approaches the values assumed in this
work. The $T_{\rm d}$-$z$ relation determined by Viero et al. (2022) increases
significantly faster than the remaining curves with redshift. This is simply
caused by the different definition of the dust temperature in Viero et al.
(2022), where Casey (2012) dust emission SEDs were used with effective
temperature $T^{\rm eff}_{\rm d}$. Different temperature definition is
necessary because, due to the interstellar dust in galaxies exhibiting a range
of temperatures, the Casey (2012) dust emission curve has a power law
component at mid-IR wavelengths. As explained in Section 4.1, the power-law
component at mid-IR gives better fits to Herschel PACS photometry. Translating
from $T^{\rm eff}_{\rm d}$ to $T_{\rm d}$ for the relationship found by Viero
et al. (2022) produces a curve which is consistent with the findings of Ismail
et al. (2023).

###  5.3 Comparison to Literature

In Figure 5 we show the comparison of the SFR-$M_{\ast}$ relation for star-
forming galaxies (star-forming main sequence) found in this work (black solid
lines) with those found by Speagle et al. (2014), Tomczak et al. (2016) and
Popesso et al. (2023), plotted with blue dashed, green dash-dotted and red
dotted 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_{\ast}$ relation
of Speagle et al. (2014). As 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^{11}\,{\rm
M_{\odot}}$.

At $z\gtrsim 2$ and ${\rm log}(M_{\ast}/{\rm M_{\odot}})\gtrsim 11$, our
values of SFR tend to be larger than those of Popesso et al. (2023). The main
reason for this is that the determination of the high-mass end of the main
sequence at high redshifts is, in both cases, based on a handful of sources.
Looking at Figure 3 of Popesso et al. (2023), it can be seen that at the
highest redshift both the bending and the normalization of the main sequence
are underpredicted by the adopted functional forms. The errors on the main
sequence parameters at high redshifts are significant enough for our best-fit
function to be consistent with the individual data points of Popesso et al.
(2023). In addition, in most of the works compiled by Popesso et al. (2023),
the IR luminosities, and hence SFRs, were determined from near- and mid-IR
data fit with relatively cold SEDs ($\sim 30$ K as opposed to $\sim 50$ K peak
temperatures found at high-$z$’s in this work – Figure 4), which can bias the
results towards lower values of SFR. On the other hand, the presence of
starburst galaxies can overpredict the high-mass end of the main sequence at
high redshifts. At low redshifts, for ${\rm log}(M_{\ast}/{\rm
M_{\odot}})\lesssim 10.5$, our values of the SFR are lower than those of
Tomczak et al. (2016) and Popesso et al. (2023). This is, again, most likely
caused by fitting the IR data with warm SEDs, where dust temperatures higher
than $\sim 20$ K were assumed.

Figure 6: The logarithm of the specific star formation rate, derived from the
best-fit functional forms of the star-forming main sequence 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 star-forming galaxies, as discussed in Section 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 $T_{\rm d}$
is consistent with the scenario, presented in Liang et al. (2019), in which
different dust temperature are mainly driven by the variations in the specific
star formation rate (see Section 5.2).

In Figure 6 we show the corresponding evolution of the specific star formation
rate with redshift, color-coded with the stellar mass. A clear turnover can be
seen at $z\sim 2$, which is a direct consequence of the bending of the
normalization of the main sequence with redshift (Equation 8). This, on the
other hand, as explained in Section 5.1, is likely caused by the increase in
molecular gas fraction, found by the Wang et al. (2022) to be exponentially
increasing with redshift. The sSFR-$z$ relation exhibits lower normalization
at higher stellar masses, since at masses $\gtrsim M_{0}$ the main sequence
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 of the dust temperature with
redshift, found in Section 5.2, is indeed caused by the corresponding increase
in the sSFRs.

For comparison, relation of Popesso et al. (2023) for ${\rm log}(M_{\ast}/{\rm
M_{\odot}})\simeq 9.5$ and 11.0 are shown as black dotted and dashed lines,
respectively. The inconsistencies at low and high redshifts come from
different normalizations of the corresponding main sequences, as explained
above. In order to investigate those further, in the following section we
explore how the evolution of the main sequence 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 galaxy stellar mass function (GSMF).

###  5.4 Growth of the Galaxy Stellar Mass Function

Figure 7: The comparison between the observed and simulated star-forming
galaxy stellar mass functions (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 Equation 10. The
mass evolution resulting from assumed star-forming main sequences of Speagle
et al. (2014), Tomczak et al. (2016) and Popesso et al. (2023), are plotted
with blue dashed, green dash-dotted and red dotted lines, respectively. The
simulated function constructed using the results of this work are presented in
solid black, while the thick solid grey 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 (Figure 5). The significant shift in densities
between the observed and simulated GSMFs is caused by the fraction of star-
forming galaxies, simulated between the redshifts in question, entering a
quiescent mode, as discussed in Section 5.4. Figure 8: The comparison between
the observed and simulated galaxy stellar mass functions (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 Equation 10. The mass evolution results
from assuming the SFR-$M_{\ast}$ 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 green dash-dotted line. The thick solid grey line represents the observed
GSMF at the final redshift of each simulation. It can be seen that the
significant gap in densities, apparent in Figure 7, is no longer present,
indicating that, indeed, the inconsistencies were caused by the fraction of
sources evolving into quiescent mode, as discussed in Section 5.4. The low-
mass slope of simulated GSMFs at $z\gtrsim 2$ is inconsistent with the
observed function, indicating that the corresponding slope of the
SFR-$M_{\ast}$ relations most likely evolves towards steeper values towards
earlier epochs.

The star-forming main sequence effectively quantifies the stellar mass growth
of a galaxy with time. We can, therefore, compare the evolution of the galaxy
stellar mass function as implied by any given main sequence 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 main sequence, causing the GSMF to
shift towards 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:

$\dot{M}(M,z)={\rm SFR}(M,z)\times(1-R),$ (10)

where $R$ is the fraction of the mass lost during passive stellar evolution,
starting at 0 and reaching 0.36 after $\sim 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 Equation 10 assuming different star-forming main sequences. In
Figure 7 we show the results of this exercise, with the observed star-forming
GSMF at the final redshift shown in grey. 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 blue dashed, green dash-
dotted and red dotted lines, respectively. The results from this work are
shown with black solid 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 main sequences. 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 towards 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 Figure 7 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 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 star-forming galaxies will transform into passive sources,
effectively shifting the evolved stellar mass function down towards lower
densities.

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\simeq 6$ has been reached. The SFR-$M_{\ast}$
relation used in this case includes all galaxies (star-forming and quiescent)
and is shown in the left panel of Figure 3. The resulting evolved GSMF at
different redshifts is compared to the observed one in Figure 8 (black thin
and grey thick lines, respectively). Since Speagle et al. (2014) and Popesso
et al. (2023) did not publish the SFR-$M_{\ast}$ relations for the galaxy
samples that include the quiescent sources, we compare our results with
Tomczak et al. (2016) only (green dashed line). It is clear from Figure 8 that
the obvious shift in galaxies number density between the observed and
simulated functions, apparent in Figure 7, 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 star-forming galaxies transforming into passive
sources. For low masses at $z\gtrsim 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 towards higher values. In order for our SFR-$M_{\ast}$ relation to
produce results consistent with this, the low-mass slope, $\gamma$, from
Equation 6 would also need to evolve. However, as explained in Section 4.2,
since our data does not probe low stellar masses at high redshifts, $\gamma$
was assumed to be constant. The fact that the constant $\gamma$ gives steeper
low-mass slopes of the GSMF at high $z$’s indicates that the power-law slope
of the SFR-$M_{\ast}$ 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\gtrsim 2$, the slope of the star-forming main sequence should
also evolve. At $z\gtrsim 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 ($\sim
1$ dex), it is also possible that the high-mass end of the SFR-$M_{\ast}$
relations in this work, and hence the growth of the stellar masses, is
contaminated by the presence of starburst galaxies. Finally, it can be seen in
Figure 8 that at $z\lesssim 2$ the SFR-$M_{\ast}$ all-galaxies relation of
Tomczak et al. (2016) slightly overpredicts the resulting GSMFs at low masses,
most likely indicating that the normalization of the main sequence at low $z$
is in this case too high.

##  6 Summary

We have stacked $\sim 100$k $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\sim 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 star formation rates 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:

  1. (i)

The SFR of star-forming galaxies populating the main sequence exhibits a
nearly linear growth with stellar mass at any given redshift at low masses.
Above the so-called bending mass, $M_{0}$, the main sequence flattens towards
an asymptotic value of the SFR. At constant stellar mass, the SFR also
increases with redshift – the redshift evolution of the normalization of the
main sequence. This behavior is consistent with the scenario in which the
evolution of the star-forming main sequence is driven by the position of the
galaxy on the integrated Schmidt-Kennicutt 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_{0}$, is exponential, increasing
sharply from $z\sim 0.5$ to $z\sim 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 main sequence is likely caused by the
phasing out of the cold accretion.

  2. (ii)

We investigate the evolution of the peak dust temperature with IR luminosity
and redshift, by fitting the stacked photometry with dust emission curve of
Casey (2012). Similarly to Drew & Casey (2022), an increase of $T_{\rm d}$
with IR luminosity was found. However, when confined to any given redshift,
our data does not show any significant evolution with $L_{\rm IR}$. We
conclude that the increase of $T_{\rm d}$ with IR luminosity is driven by the
evolution of the $L_{\rm IR}$ with redshift for the main-sequence galaxies. We
find that the dust temperature increases linearly with redshift, consistent
with the scenario in which higher values of $T_{\rm d}$ for star-forming
galaxies are driven by the increase in their specific star formation rates
(e.g. Liang et al. 2019).

  3. (iii)

The evolution of the main sequence 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 main sequence at $z\lesssim 2$
for masses below the bending mass, $M_{0}$, is slightly lower, while the SFRs
above $M_{0}$ at high redshifts exhibit larger values (except for the MS of
Speagle et al. 2014, where no main-sequence bending was found). We attribute
these inconsistencies to the dust temperatures assumed, when fitting the IR
data. As found in this work, $T_{\rm d}$ increases from $\sim 20$ K to $\sim
50$ K between redshifts 0 and 7 and hence, assuming dust temperatures at the
level of $\sim 30$ K, as is often done in the absence of the far-IR data, will
correspondingly affect the determined IR luminosities (and hence SFRs). We
also note that the relatively higher values of the SFRs at high redshifts for
$M_{\ast}>10^{11}\,{\rm M_{\odot}}$ can be, at least in part, caused by the
contamination of the main-sequence sample with starburst galaxies.

  4. (iv)

Adopting different shapes of the SFR-$M_{\ast}$ relation, we simulate the time
evolution for the observed galaxy stellar mass functions, 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_{\ast}$ relation, found using both star-forming and quiescent galaxies,
we confirm that the overprediction of star-forming galaxies number densities
is caused by a significant fraction of sources evolving into the quiescent
phase and, thereby, leaving the star-forming main sequence. In addition, we
find the simulated low-mass end of the GSMF to give higher number densities at
high $z$, indicating that the low-mass slope of the main sequence likely
evolves towards steeper values with redshift. However, due to the significant
uncertainties, more detailed analysis based on larger high-$z$ samples is
required.

###### Acknowledgements.

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

  * Ashby et al. (2013) Ashby M.L.N., Willner S.P., Fazio G.G., et al., 2013, ApJ, 769, 80
  * Ashby et al. (2015) Ashby M.L.N., Willner S.P., Fazio G.G., et al., 2015, ApJS, 218, 33
  * Behroozi et al. (2013) Behroozi P.S., Wechsler R.H., Conroy C., 2013, ApJ, 770, 57
  * Benítez (2000) Benítez N., 2000, ApJ, 536, 571
  * Bertin & Arnouts (1996) Bertin E., Arnouts S., 1996, A&AS, 117, 393
  * Brammer et al. (2008) Brammer G.B., van Dokkum P.G., Coppi P., 2008, ApJ, 686, 1503
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R.C., et al., 2000, ApJ, 533, 682
  * Carnall et al. (2020) Carnall A.C., Walker S., McLure R.J., et al., 2020, MNRAS, 496, 695
  * Casey (2012) Casey C.M., 2012, MNRAS, 425, 3094
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chapman et al. (2002) Chapman S.C., Smail I., Ivison R.J., et al., 2002, ApJ, 573, 66
  * Coleman et al. (1980) Coleman G.D., Wu C.C., Weedman D.W., 1980, ApJS, 43, 393
  * Daddi et al. (2010) Daddi E., Elbaz D., Walter F., et al., 2010, ApJ, 714, L118
  * Daddi et al. (2022) Daddi E., Delvecchio I., Dimauro P., et al., 2022, A&A, 661, L7
  * Davidzon et al. (2017) Davidzon I., Ilbert O., Laigle C., et al., 2017, A&A, 605, A70
  * Drew & Casey (2022) Drew P.M., Casey C.M., 2022, ApJ, 930, 142
  * Dudzevičiūtė et al. (2020) Dudzevičiūtė U., Smail I., Swinbank A.M., et al., 2020, MNRAS, 494, 3828
  * Faisst et al. (2020) Faisst A.L., Fudamoto Y., Oesch P.A., et al., 2020, MNRAS, 498, 4192
  * Fioc & Rocca-Volmerange (1999) Fioc M., Rocca-Volmerange B., 1999, arXiv e-prints, astro–ph/9912179
  * Furusawa et al. (2008) Furusawa H., Kosugi G., Akiyama M., et al., 2008, ApJS, 176, 1
  * Furusawa et al. (2016) Furusawa H., Kashikawa N., Kobayashi M.A.R., et al., 2016, ApJ, 822, 46
  * Geach et al. (2011) Geach J.E., Smail I., Moran S.M., et al., 2011, ApJ, 730, L19
  * Geach et al. (2017) Geach J.E., Dunlop J.S., Halpern M., et al., 2017, MNRAS, 465, 1789
  * Griffin et al. (2010) Griffin M.J., Abergel A., Abreu A., et al., 2010, A&A, 518, L3
  * Holland et al. (2013) Holland W.S., Bintley D., Chapin E.L., et al., 2013, MNRAS, 430, 2513
  * Hudelot et al. (2012) Hudelot P., Cuillandre J.C., Withington K., et al., 2012, VizieR Online Data Catalog: The CFHTLS Survey (T0007 release) (Hudelot+ 2012), VizieR On-line Data Catalog: II/317. Originally published in: SPIE Conf. 2012
  * Ilbert et al. (2009) Ilbert O., Capak P., Salvato M., et al., 2009, ApJ, 690, 1236
  * Ilbert et al. (2010) Ilbert O., Salvato M., Le Floc’h E., et al., 2010, ApJ, 709, 644
  * Ismail et al. (2023) Ismail D., Beelen A., Buat V., et al., 2023, arXiv e-prints, arXiv:2307.15747
  * Jarvis et al. (2013) Jarvis M.J., Bonfield D.G., Bruce V.A., et al., 2013, MNRAS, 428, 1281
  * Kennicutt (1998) Kennicutt J. Robert C., 1998, ApJ, 498, 541
  * Kennicutt & De Los Reyes (2021) Kennicutt J. Robert C., De Los Reyes M.A.C., 2021, ApJ, 908, 61
  * Koprowski et al. (2018) Koprowski M.P., Coppin K.E.K., Geach J.E., et al., 2018, MNRAS, 479, 4355
  * Koyama et al. (2011) Koyama Y., Kodama T., Nakata F., Shimasaku K., Okamura S., 2011, ApJ, 734, 66
  * Labbé et al. (2005) Labbé I., Huang J., Franx M., et al., 2005, ApJ, 624, L81
  * Lawrence et al. (2007) Lawrence A., Warren S.J., Almaini O., et al., 2007, MNRAS, 379, 1599
  * Lee et al. (2015) Lee N., Sanders D.B., Casey C.M., et al., 2015, ApJ, 801, 80
  * Leja et al. (2015) Leja J., van Dokkum P.G., Franx M., Whitaker K.E., 2015, ApJ, 798, 115
  * Leja et al. (2022) Leja J., Speagle J.S., Ting Y.S., et al., 2022, ApJ, 936, 165
  * Leslie et al. (2020) Leslie S.K., Schinnerer E., Liu D., et al., 2020, ApJ, 899, 58
  * Liang et al. (2019) Liang L., Feldmann R., Kereš D., et al., 2019, MNRAS, 489, 1397
  * Lutz et al. (2011) Lutz D., Poglitsch A., Altieri B., et al., 2011, A&A, 532, A90
  * Madau (1995) Madau P., 1995, ApJ, 441, 18
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Magnelli et al. (2014) Magnelli B., Lutz D., Saintonge A., et al., 2014, A&A, 561, A86
  * McCracken et al. (2012) McCracken H.J., Milvang-Jensen B., Dunlop J., et al., 2012, A&A, 544, A156
  * McLeod et al. (2021) McLeod D.J., McLure R.J., Dunlop J.S., et al., 2021, MNRAS, 503, 4413
  * Mehta et al. (2018) Mehta V., Scarlata C., Capak P., et al., 2018, ApJS, 235, 36
  * Mérida et al. (2023) Mérida R.M., Pérez-González P.G., Sánchez-Blázquez P., et al., 2023, ApJ, 950, 125
  * Merlin et al. (2015) Merlin E., Fontana A., Ferguson H.C., et al., 2015, A&A, 582, A15
  * Moster et al. (2013) Moster B.P., Naab T., White S.D.M., 2013, MNRAS, 428, 3121
  * Oliver et al. (2012) Oliver S.J., Bock J., Altieri B., et al., 2012, MNRAS, 424, 1614
  * Pilbratt et al. (2010) Pilbratt G.L., Riedinger J.R., Passvogel T., et al., 2010, A&A, 518, L1
  * Poglitsch et al. (2010) Poglitsch A., Waelkens C., Geis N., et al., 2010, A&A, 518, L2
  * Pope et al. (2006) Pope A., Scott D., Dickinson M., et al., 2006, MNRAS, 370, 1185
  * Popesso et al. (2023) Popesso P., Concas A., Cresci G., et al., 2023, MNRAS, 519, 1526
  * Reuter et al. (2020) Reuter C., Vieira J.D., Spilker J.S., et al., 2020, ApJ, 902, 78
  * Santini et al. (2014) Santini P., Maiolino R., Magnelli B., et al., 2014, A&A, 562, A30
  * Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
  * Schreiber et al. (2015) Schreiber C., Pannella M., Elbaz D., et al., 2015, A&A, 575, A74
  * Schreiber et al. (2016) Schreiber C., Elbaz D., Pannella M., et al., 2016, A&A, 589, A35
  * Schreiber et al. (2018) Schreiber C., Elbaz D., Pannella M., et al., 2018, A&A, 609, A30
  * Speagle et al. (2014) Speagle J.S., Steinhardt C.L., Capak P.L., Silverman J.D., 2014, ApJS, 214, 15
  * Tomczak et al. (2014) Tomczak A.R., Quadri R.F., Tran K.V.H., et al., 2014, ApJ, 783, 85
  * Tomczak et al. (2016) Tomczak A.R., Quadri R.F., Tran K.V.H., et al., 2016, ApJ, 817, 118
  * Viero et al. (2022) Viero M.P., Sun G., Chung D.T., Moncelsi L., Condon S.S., 2022, MNRAS, 516, L30
  * Wang et al. (2022) Wang T.M., Magnelli B., Schinnerer E., et al., 2022, A&A, 660, A142
  * Whitaker et al. (2014) Whitaker K.E., Franx M., Leja J., et al., 2014, ApJ, 795, 104
  * Williams et al. (2009) Williams R.J., Quadri R.F., Franx M., van Dokkum P., Labbé I., 2009, ApJ, 691, 1879
  * Wuyts et al. (2007) Wuyts S., Labbé I., Franx M., et al., 2007, ApJ, 655, 51

##  Appendix A Plots & Tables

Figure 9: Example stacked images, shown in bins of stellar mass and redshift,
for the Herschel SPIRE 500 ${\rm\mu m}$ band. Figure 10: Stacked far-IR
photometry for the whole sample of galaxies (star-forming and quiescent),
listed in the top panel of Table 2. The best-fit dust emission curves, adopted
from Casey (2012), are shown in solid green with $1\sigma$ standard deviation
in grey. Figure 11: As in Figure 10 but for star-forming sample only – bottom
panel in Table 2.

Table 2: The stacked values of the SFR (in $\hbox{$\rm\thinspace M_{\odot}$}\,{\rm yr}^{-1}$) for all the galaxies (star-forming and quiescent; top panel) and star-forming galaxies (bottom panel) between redshifts 0.25 and 6.75. The first and eighth columns list the mass bins, where $\mathcal{M}\leavevmode\nobreak\ \equiv\leavevmode\nobreak\ \log_{10}(M_{\star}/\hbox{$\rm\thinspace M_{\odot}$})$. At $z>3.75$, in order to increase the statistics in each mass bin, the width of the bins was increased, as shown in column 8. | $0.25\leq z<0.75$ | $0.75\leq z<1.25$ | $1.25\leq z<1.75$ | $1.75\leq z<2.25$ | $2.25\leq z<2.75$ | $2.75\leq z<3.75$ |  | $3.75\leq z<4.75$ | $4.75\leq z<6.75$
---|---|---|---|---|---|---|---|---|---
| log(SFR) | log(SFR) | log(SFR) | log(SFR) | log(SFR) | log(SFR) |  | log(SFR) | log(SFR)
$\phantom{0}9.25\leq\mathcal{M}<\phantom{0}9.50$ | $-0.13\pm 0.06$ |  |  |  |  |  |  |  |
$\phantom{0}9.50\leq\mathcal{M}<\phantom{0}9.75$ | $\phantom{-}0.18\pm 0.04$ | $0.53\pm 0.05$ |  |  |  |  |  |  |
$\phantom{0}9.75\leq\mathcal{M}<10.00$ | $\phantom{-}0.37\pm 0.05$ | $0.74\pm 0.05$ | $1.04\pm 0.06$ |  |  |  |  |  |
$10.00\leq\mathcal{M}<10.25$ | $\phantom{-}0.54\pm 0.05$ | $0.91\pm 0.05$ | $1.21\pm 0.05$ | $1.39\pm 0.06$ | $1.52\pm 0.07$ |  |  |  |
$10.25\leq\mathcal{M}<10.50$ | $\phantom{-}0.60\pm 0.05$ | $1.05\pm 0.05$ | $1.38\pm 0.05$ | $1.50\pm 0.06$ | $1.68\pm 0.07$ | $1.88\pm 0.07$ |  |  |
$10.50\leq\mathcal{M}<10.75$ | $\phantom{-}0.66\pm 0.05$ | $1.14\pm 0.05$ | $1.46\pm 0.05$ | $1.74\pm 0.05$ | $1.82\pm 0.06$ | $1.98\pm 0.07$ | $10.50\leq\mathcal{M}<11.00$ | $2.34\pm 0.08$ | $2.47\pm 0.14$
$10.75\leq\mathcal{M}<11.00$ | $\phantom{-}0.69\pm 0.05$ | $1.20\pm 0.05$ | $1.59\pm 0.05$ | $1.78\pm 0.05$ | $1.96\pm 0.06$ | $2.19\pm 0.07$
$11.00\leq\mathcal{M}<11.25$ | $\phantom{-}0.74\pm 0.05$ | $1.22\pm 0.05$ | $1.66\pm 0.05$ | $1.98\pm 0.05$ | $2.22\pm 0.06$ | $2.42\pm 0.06$ | $11.00\leq\mathcal{M}<11.50$ | $2.61\pm 0.07$ | $2.80\pm 0.07$
$11.25\leq\mathcal{M}<11.50$ | $\phantom{-}0.60\pm 0.07$ | $1.27\pm 0.05$ | $1.71\pm 0.06$ | $2.13\pm 0.06$ | $2.40\pm 0.06$ | $2.55\pm 0.07$
$\phantom{0}9.25\leq\mathcal{M}<\phantom{0}9.50$ | $-0.05\pm 0.06$ |  |  |  |  |  |  |  |
$\phantom{0}9.50\leq\mathcal{M}<\phantom{0}9.75$ | $\phantom{-}0.27\pm 0.04$ | $0.57\pm 0.05$ |  |  |  |  |  |  |
$\phantom{0}9.75\leq\mathcal{M}<10.00$ | $\phantom{-}0.47\pm 0.05$ | $0.80\pm 0.05$ | $1.06\pm 0.06$ |  |  |  |  |  |
$10.00\leq\mathcal{M}<10.25$ | $\phantom{-}0.68\pm 0.05$ | $1.01\pm 0.05$ | $1.27\pm 0.06$ | $1.42\pm 0.06$ | $1.53\pm 0.08$ |  |  |  |
$10.25\leq\mathcal{M}<10.50$ | $\phantom{-}0.81\pm 0.05$ | $1.22\pm 0.05$ | $1.49\pm 0.05$ | $1.59\pm 0.06$ | $1.76\pm 0.06$ | $1.89\pm 0.07$ |  |  |
$10.50\leq\mathcal{M}<10.75$ | $\phantom{-}0.93\pm 0.05$ | $1.35\pm 0.05$ | $1.63\pm 0.05$ | $1.85\pm 0.05$ | $1.93\pm 0.06$ | $2.03\pm 0.07$ | $10.50\leq\mathcal{M}<11.00$ | $2.38\pm 0.08$ | $2.47\pm 0.14$
$10.75\leq\mathcal{M}<11.00$ | $\phantom{-}1.01\pm 0.05$ | $1.47\pm 0.05$ | $1.78\pm 0.05$ | $1.93\pm 0.05$ | $2.07\pm 0.06$ | $2.27\pm 0.07$
$11.00\leq\mathcal{M}<11.25$ | $\phantom{-}1.12\pm 0.05$ | $1.51\pm 0.05$ | $1.88\pm 0.05$ | $2.11\pm 0.05$ | $2.31\pm 0.06$ | $2.50\pm 0.06$ | $11.00\leq\mathcal{M}<11.50$ | $2.64\pm 0.08$ | $2.80\pm 0.07$
$11.25\leq\mathcal{M}<11.50$ | $\phantom{-}1.12\pm 0.06$ | $1.65\pm 0.05$ | $1.96\pm 0.05$ | $2.27\pm 0.05$ | $2.53\pm 0.06$ | $2.61\pm 0.07$