© The Authors 2026

1. Introduction

One of the main goals of modern extragalactic astronomy is to describe the evolution of star formation (SF) throughout cosmic history. The most direct way of achieving this goal is through the calibration of various relationships that characterize how SF in galaxies evolves with time. The most important examples of such relationships are the ultraviolet (UV) and infrared (IR) luminosity function (LF; e.g., Gruppioni et al. 2013; Koprowski et al. 2017; Wang et al. 2019; Bouwens et al. 2021; Fujimoto et al. 2024) and the star formation rate density (ρ_SFR; e.g., Madau & Dickinson 2014; Dunlop et al. 2017; Traina et al. 2024; Liu et al. 2026; Barrufet et al. 2025).

Since the discovery of the cosmic infrared background (Puget et al. 1996; Hauser et al. 1998), the IR contribution to the total ρ_SFR, determined through the integration of the corresponding luminosity function, has been the primary focus of numerous studies (e.g., Magnelli et al. 2013; Gruppioni et al. 2013; Madau & Dickinson 2014; Rowan-Robinson et al. 2016; Koprowski et al. 2017; Dunlop et al. 2017; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024; Magnelli et al. 2024; Fujimoto et al. 2024; Sun et al. 2025; Liu et al. 2026). While most of the results seem to be consistent out to z ∼ 2 − 3, where the SF activity of the Universe reaches its peak, at ≳4 the derived values disagree by over an order of magnitude.

The accurate assessment of the evolution of SF in the early Universe is essential for evaluating current galaxy evolution models. Thanks to the sensitivity of the rest-frame UV observations, the unobscured portion of the SF in galaxies has been investigated out to redshifts as high as ∼11 (e.g., Bouwens et al. 2014; Laporte et al. 2016; Harikane et al. 2024; McLeod et al. 2024; Donnan et al. 2024). Due to the time required for the formation of the interstellar dust, the unobscured SF is thought to dominate the total budget at such high redshifts. Since it is now well known that by cosmic noon (z ∼ 2 − 3) the vast majority of stellar emission in the UV gets absorbed by interstellar dust and reemitted in the IR, we expect ρ_SFR to transition from being UV dominated to IR dominated. However, because of the observational limitations in the infrared (see Casey et al. 2014 for details), the exact redshift of this transition is yet to be determined.

The contribution of the dust-enshrouded stellar emission to the total SF budget in the early Universe has primarily been assessed through the most extreme IR sources, known as submillimeter galaxies (e.g., Smail et al. 1997; Hughes et al. 1998) or more broadly as dusty star-forming galaxies (DSFGs; Casey et al. 2014), which have very large IR luminosities (L_IR ≳ 10¹² L_⊙; e.g., Chapman et al. 2005; Michałowski et al. 2017) and extremely high SF rates (SFR ≳ 100 M_⊙ yr⁻¹; e.g., Swinbank et al. 2014; Koprowski et al. 2016). Observations conducted using NASA’s Spitzer Space Telescope and ESA’s Herschel Space Observatory enabled the determination of the infrared luminosity functions out to z ≃ 2 − 4 (e.g., Le Floc’h et al. 2005; Caputi et al. 2007; Rodighiero et al. 2010; Patel et al. 2013; Magnelli et al. 2013; Gruppioni et al. 2013; Rowan-Robinson et al. 2016; Wang et al. 2019). However, only the most extreme DSFGs have been detected at z ≳ 2, which necessitated significant extrapolations to account for fainter objects, leading to uncertain and often inconsistent results.

Koprowski et al. (2017) utilized a sample of the James Clerk Maxwell Telescope (JCMT) 850-μm-selected sample to construct the IR luminosity function out to z ∼ 4. While the submillimeter wavelengths benefit from the negative K-correction, allowing for the observation of IR galaxies out to very high redshifts, surprisingly low numbers of sources were detected at z ≳ 2. Similarly, due to the very limited areas observed, the high-resolution ALMA blank-sky surveys conducted in recent years also struggled to detect any sources at very high redshifts (e.g., Walter et al. 2016; Dunlop et al. 2017; Bouwens et al. 2020). Gruppioni et al. (2020) made use of the ALMA sources serendipitously detected as a part of the ALMA Large Program to INvestigate CII at Early Times (ALPINE) in order to derive the IR LF out to z ∼ 6. Fujimoto et al. (2024) reached similar redshifts utilizing the ALMA Lensing Cluster Survey (ALCS) data, with faint objects detected through lensing. Most recently, Magnelli et al. (2024), Traina et al. (2024), and Liu et al. (2026) determined the contribution of IR sources to the total ρ_SFR in the early Universe using the somewhat inhomogeneous sample of individual ALMA pointings collected as a part of the A³COSMOS survey (Liu et al. 2019; Adscheid et al. 2024). Due to the limited sensitivity of the IR data, the corresponding faint-end slope of the IR LF reported in these works ranges between –0.2 and –1.0, with errors approaching 0.5. In addition, Dudzevičiūtė et al. (2020) analyzed a high-resolution ALMA follow-up study (AS2UDS; Stach et al. 2019) of the JCMT SCUBA-2 submillimeter survey of the UKIDSS UDS field (Geach et al. 2017). As a part of their analysis, Dudzevičiūtė et al. (2020) derived photometric redshifts, IR luminosities, IR luminosity functions, and the corresponding star formation rate densities down to the 870 μm flux limit of 3.6 mJy, later extrapolating their results down to a flux limit of 1 mJy using the ALMA number counts of Hatsukade et al. (2018). They found that the contribution of submillimeter galaxies to the total star formation rate density in the Universe peaks at z ∼ 3, with a contribution of ∼15% for sources with S₈₇₀ > 3.6 mJy and ∼60% for S₈₇₀ > 1 mJy, indicating that roughly half of the ρ_SFR at that redshift comes from ultra luminous IR galaxies (ULIRGs)-luminosity sources.

The aim of this study is to add to the existing knowledge by deriving the evolution of the IR LF, including its faint-end slope, up to z ∼ 6. In order to include the contribution from the faintest galaxies, we adopted an indirect approach, where the faint-end portion of the LF is determined through stacking of the optically selected mass-complete sample of star-forming galaxies in the far-IR (FIR) maps of Herschel and JCMT maps, with the stellar mass used as a proxy for the IR luminosity. The paper is structured as follows. We present the data in Section 2. In Section 3.1.1 we explain the details behind our derivation of the IR LF faint-end slope. In Section 3.1.2 we describe how we used data from AS2UDS, a high-resolution ALMA follow-up study of all the SCUBA-2 sources detected within the S2CLS map of the UKIDSS UDS field (Geach et al. 2017; Stach et al. 2019), to construct the bright end of our IR LF. In Section 3.1.3 we construct the functional form of the IR LF and quantify its evolution with redshift out to z ∼ 6, discuss the results, and compare them with recent literature. The corresponding star formation rate density is derived, discussed, and compared with other works in Section 3.2. We provide a summary of this work in Section 4. Throughout the paper we assume a flat cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7.

2. Data

The faint end of the infrared luminosity function was established via stacking the optical/near-IR catalogs of the UKIDSS Ultra Deep Survey (UDS) and COSMOS fields (McLeod et al. 2021) following the methodology outlined in Koprowski et al. (2024), where the stellar masses of galaxies were used as a proxy for the IR luminosity. For the purpose of extracting the stacked FIR flux densities, 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, where Herschel maps at 100, 160, 250, 350, and 500 μm were utilized. In order to constrain the Rayleigh–Jeans tail of the dust emission curve, we also include the data collected as part of the SCUBA-2 Cosmology Legacy Survey (S2CLS; Geach et al. 2017).

The bright end of the IR LF (Section 3.1.2) was derived using high-resolution ALMA follow-up data (AS2UDS; Stach et al. 2019; Dudzevičiūtė et al. 2020) of the S2CLS UKIDSS UDS field (Geach et al. 2017). The original SCUBA-2 map encompasses an area of 0.96 square degrees, with a noise level of less than 1.3 mJy and a median sensitivity of σ₈₅₀ = 0.88 mJy beam⁻¹, where 716 sources were detected. As explained in Geach et al. (2017), the SCUBA-2 source detection completeness for the UDS field was determined by the recovery rate of artificial sources injected into the jackknife maps as a function of input flux.

A detailed description of the ALMA observations and data reduction and the construction of the catalog can be found in Stach et al. (2019). In summary, all 716 SCUBA-2 sources were observed using ALMA Band 7 during Cycles 1, 3, 4, and 5. All the maps were tapered to 0.5 arcsec for source detection purposes, with the resulting AS2UDS catalog consisting of 708 individual ALMA-identified submillimeter galaxies (σ₈₇₀ > 4.3). Based on 60 000 simulated ALMA observations, the sample was found to be complete at S₈₇₀ ≳ 4 mJy.

The associated photometric redshifts and FIR luminosities for all AS2UDS sources were adopted from Dudzevičiūtė et al. (2020). In short, the available photometry spans a wavelength range from the u-band to radio, with low-resolution Herschel SPIRE data deblended using the procedure explained in Swinbank et al. (2014). Photometric redshifts as well as IR luminosities were found using the updated MAGPHYS spectral energy distribution modeling code (da Cunha et al. 2015; Battisti et al. 2019), along with the stellar models from Bruzual & Charlot (2003), the initial mass function of Chabrier (2003), and the dust attenuation model of Charlot & Fall (2000). The resulting photometric redshifts were tested against the available spectroscopic data, yielding a median offset – defined as (z_spec − z_phot)/(1 + z_spec) – of −0.005 ± 0.003, with a dispersion of 0.13.

The corresponding dust temperatures were found by Dudzevičiūtė et al. (2020) using MAGPHYS and the modified blackbody curve:

$$\begin{array}{c}\hfill S\propto (1-e^{-{\tau }_{\mathrm{rest}}})\times B({\nu }_{\mathrm{rest}},T),\end{array}$$(1)

where B(ν_rest, T) is the Planck function and the optical depth defined is as τ_rest ≡ (ν_rest/ν₀)^β, with the emissivity index β set to 1.8. As explained in Dudzevičiūtė et al. (2020), the corresponding temperatures for both MAGPHYS and the modified blackbody fits were found to be in a very good agreement, with a typical fractional difference of (T_d^MBB − T_d^MAGPHYS)/T_d^MBB = −0.28 ± 0.01.

3. Analysis and discussion

3.1. Infrared luminosity function

3.1.1. Faint end

The IR luminosities for the galaxies constituting the faint-end portion of the luminosity function were extracted indirectly by stacking the optical/near-IR samples in the FIR Herschel and JCMT maps and using the stellar mass as proxies. The rationale behind this approach is that the IR luminosities of the star-forming galaxies are tightly correlated with their stellar masses via the so-called main sequence (e.g., Speagle et al. 2014; Tomczak et al. 2016; Daddi et al. 2022; Popesso et al. 2023; Koprowski et al. 2024). The stacking procedure was performed by Koprowski et al. (2024), and we adopted their results in this work (Figure 1 with the individual L_IR values listed in Table B.1). The redshift bins were designed to encompass the approximately 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.2, 1.6, 2.2, 3.2, 5.7]).

Fig. 1. Infrared luminosity-stellar mass relation determined through stacking the optical/near-IR mass-complete samples of McLeod et al. (2021) in the FIR Herschel and JCMT maps, adopted from Koprowski et al. (2024). The stacked LIR values were used to determine the faint-end portion of the IR LF, as explained in Section 3.1.1.

As we consider a sample of star-forming galaxies here, our assessment of the faint end of the IR luminosity function in effect excludes both quiescent and starburst galaxies (see Koprowski et al. 2024 for details). This is warranted, as passive sources are much less numerous than the active galaxies, and they exhibit significantly lower IR luminosities at a given stellar mass, rendering their contribution to the resulting LF negligible. Starbursts, on the other hand, are much more luminous in the IR, and their contribution to the LF is therefore assessed using the ALMA data, as explained in Section 3.1.2.

Since the L_IR–M_* relation, similarly to the star-forming main sequence, exhibits a turnover at high masses, the functional form of Lee et al. (2015) was adopted:

$$\begin{array}{c}\hfill L_{\mathrm{IR}}=L_{\mathrm{IR}}^{\mathit{max}}/(1+{(M_{\ast }/M_0)}^{-\gamma }),\end{array}$$(2)

with the slope, γ, set to one. As explained in Koprowski et al. (2024), the most accurate fits to the L_IR–M_* relation are produced when the logarithm of the free parameters of Eq. (2) are assumed to follow an exponential evolution with redshift, where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \log (L_{\mathrm{IR}}^{\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}\end{array}$$(3)

Eq. (2) was fit to the stacked data (Table B.1) using nonlinear least square fitting. The best-fit curves are shown in Figure 1, with the best-fit parameters of Eq. (3) listed in Table 2.

Because of the lack of ALMA detections (IR LF bright end) at low redshifts (Section 3.1.2), we limited our analysis to the redshift bins of [1.2, 1.6, 2.2, 3.2, 5.7]. In addition, we only considered z-M_* bins with stellar masses below the main-sequence bending mass, M₀, at which point the L_IR-M_* relationship leaves the power-law regime (M₀/M_⊙ ∼ 10^11.0 at 1.2 < z < 2.2 rising to M₀/M_⊙ ∼ 10^11.2 at z > 2.2). The IR luminosity for each z-M_* bin was extracted from the L_IR–M_* best-fit relation (Equations (2) and (3) with the best-fit parameters from Table 2). In order to determine the corresponding IR luminosity function, a standard 1/V_max method (Schmidt 1968) was used:

$$\begin{array}{c}\hfill \mathrm{\Phi }(L,z)=\frac{1}{\mathrm{\Delta }\log L}{\displaystyle \sum _i}\frac{1}{V_i}.\end{array}$$(4)

The width of the L_IR bin in log space, Δlog L, was taken from the best-fit L_IR–M_* functional form of Equation (2). Given that the width of the stellar mass bin is Δlog(M_*/M_⊙) = 0.25 and the slope, γ, in Equation (2) was set to one, Δlog L below the bending mass, M₀, is also ≃0.25. The number of sources in each bin, i, is listed at the bottom panel of Table B.1, and V_i is the comoving volume available to the ith source.

The resulting values for the faint-end portion of the IR LF are listed in Table 1 and depicted as black points in Figure 2. To assess the errors on the number of sources in each bin, the bootstrapping method was used. 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 process was repeated 1000 times, and the uncertainties on the number of sources in each z–M_* bin were taken to be the standard deviation of the resulting simulated values. The errors on the average stacked values of L_IR, listed in Table B.1, are from Koprowski et al. (2024). The errors on the IR LF faint-end data points are driven by the errors on the number of sources in each bin (Bottom panel of Table B.1), and as can be seen in Figure 2 and Table 1, they are significantly smaller than the errors on the L_IR.

Fig. 2. Infrared luminosity function found in this work at different redshift bins, as indicated in the plot. The faint-end data (Table 1), determined through stacking (see Section 3.1.1 for details), are depicted by black points with horizontal error bars. The remaining bright-end data (Table 3) found using the ALMA follow-up data of the S2CLS UKIDSS UDS sources (AS2UDS; Stach et al. 2019) are also shown (Section 3.1.2). The best-fit Schechter functions of Equation (6) are represented by dashed black lines, with the gray area representing 1σ uncertainties. The faint-end slope was determined at two low redshift bins, and the more accurate result, α = −0.26 ± 0.11, was adopted at all the remaining redshifts, as explained in Section 3.1.3.

3.1.2. Bright end

The bright end of the IR LF was determined using high-resolution ALMA follow-up data of all the S2CLS UKIDSS UDS sources from AS2UDS (Stach et al. 2019; Dudzevičiūtė et al. 2020). The main advantage of the ALMA data is the high resolution of the observations, which allow for robust optical counterpart identification. Each source has two values of completeness associated with it – one from the original S2CLS survey (Geach et al. 2017) and one from the ALMA follow-up observations (Stach et al. 2019). To make sure the data used here are flux complete, we limited the ALMA sample to sources with fluxes S₈₇₀ > 4 mJy, for which the completeness in the ALMA maps was established to be close to 100%. Since, as stated in Stach et al. (2019), for SCUBA-2 sources above S₈₅₀ ∼ 3.5 mJy all the single-dish flux in ALMA maps, on average, has been recovered, we can treat ALMA S₈₇₀ > 4 mJy sources as a flux-complete (> 80% completeness in the original S2CLS SCUBA-2 UDS map of Geach et al. 2017) homogeneous sample.

As discussed briefly in Section 2, the IR luminosities were taken from Dudzevičiūtė et al. (2020), where in order to get reliable values of L_IR, we limited our sample to sources with at least one Herschel SPIRE detection. We therefore required our sample to be 870 μm flux complete (S₈₇₀ > 4 mJy) and have at least one SPIRE detection. The median dust temperature of the sources in the AS2UDS sample was estimated by Dudzevičiūtė et al. (2020) to be ∼30 K (with a 68th percentile range of 25.7–37.3 K). Adopting the modified blackbody curve of Equation (1), with the 870 μm flux of 4 mJy and the dust temperature of 37.3 K (upper 68th percentile), we arrived at the conservative IR luminosity completeness limit of log(L_IR/L_⊙) = 12.5 at z < 1.6, respectively rising to 12.75 and 13.00 at z < 3.2 and < 5.7 (Table 3).

The IR luminosity function was found using the 1/V_max method (Schmidt 1968), where

$$\begin{array}{c}\hfill \mathrm{\Phi }(L,z)=\frac{1}{\mathrm{\Delta }\log L}{\displaystyle \sum _i}\frac{1-\mathrm{FDR}}{w_i\times V_i},\end{array}$$(5)

with Δlog L being the width of the luminosity bin in log space and V_i as the comoving volume available to the ith source, while the false detection rate (FDR) and the completeness, w_i, are from the original SCUBA-2 survey (Geach et al. 2017). As explained in Geach et al. (2017), a fraction of low signal-to-noise ratio sources may be missed by the source detection algorithm due to them residing in the noise valleys. For the purpose of the sample completeness, the fraction of missed galaxies was therefore quantified in the form of the source completeness, w_i. Similarly, false sources may be identified due to statistical fluctuations expected from Gaussian noise, as quantified by FDR. In both cases, as described in Geach et al. (2017), the jackknife noise-only map was populated with artificial sources, after which the source extraction algorithm was applied. The ratio between the recovered and the injected numbers were then used to determine w_i and FDR for each source in the parent SCUBA-2 sample of Geach et al. (2017), which were then adopted in our calculations of the IR LF of Equation (5). The resulting values are summarized in Table 3 and depicted as color points with gray error bars in Figure 2. In order to assess the errors, we conducted Monte Carlo simulations in which redshifts and infrared luminosities were randomly sampled from normal distributions, with the means and standard deviations corresponding to the catalog values and their associated errors. The 1σ uncertainties associated with the stacked IR LFs were then obtained by computing the standard deviations of the simulated data.

We note that our derived form of the bright end of the IR LF closely resembles that presented in Dudzevičiūtė et al. (2020). The only notable distinctions are that we have assumed slightly different redshift binning and included the completeness and false detection rate values of the original S2CLS sample of Geach et al. (2017) in our 1/V_max analysis. However, given that the initial SCUBA2 sample has been determined to be nearly flux complete for S₈₅₀ > 4 mJy, incorporating them did not have a substantial impact on our results. The comparison with the IR LF and the ρ_SFR values found in Dudzevičiūtė et al. (2020) is presented in Sections 3.1.3 and 3.2.

3.1.3. IR LF functional form

Following previous works, we fit the IR LF with the Schechter function:

$$\begin{array}{c}\hfill \mathrm{\Phi }(L,z)={\mathrm{\Phi }}_{★}{\left(\frac{L}{L_{★}}\right)}^{\alpha }\exp \left(\frac{-L}{L_{★}}\right),\end{array}$$(6)

where Φ_★ is the normalization parameter, α is the faint-end slope, and L_★ is the characteristic luminosity that marks the border between the power law and the exponential fits (for comparison purposes, the modified Schechter fits are presented in Appendix A). Due to the relatively shallow depth of IR observations, the faint-end slope is usually determined only at redshifts z ≪ 1 (e.g., Gruppioni et al. 2013) and then fixed at higher redshifts, with the exception of Koprowski et al. (2017), where α was found using the ALMA survey of Hubble Ultra Deep Field (Dunlop et al. 2017) at z ≃ 2. Here, we fit the faint-end slope at the two lowest redshift bins (1.2 ≤ z < 1.6 and 1.6 ≤ z < 2.2), finding α = −0.26 ± 0.11 and −0.23 ± 0.18, respectively, and we adopted a more accurate 1.2 ≤ z < 1.6 value of −0.26 at all redshift bins.

The faint- and bright-end data from Figure 2 were fit at each redshift using Equation (6) (solid black curves in Figure 2), with α set to −0.26 and using the nonlinear least squares analysis with the errors estimated from the corresponding covariance matrix. The resulting values of Φ_★ and L_★ are summarized in Table 4 and depicted as black points with error bars in Figure 3. The functional form for the redshift dependence for each parameter (black solid curves in Figure 3) was then found, where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\Phi }}_{★}/{\mathrm{Mpc}}^{-3}{\mathrm{dex}}^{-1}& ={\mathrm{\Phi }}_{★,0}{[{\left(\frac{z}{z_0}\right)}^{a_1}+{\left(\frac{z}{z_0}\right)}^{a_2}]}^{-1}\hfill \\ \hfill L_{★}/{\mathrm{L}}_{\odot }& =b_1\times z^{b_2}\hfill \end{array}\end{array}$$(7)

Fig. 3. Redshift evolution of the Schechter function parameters found in this work. The values determined at each redshift (Table 4) are shown as black points, with the best-fit functions (Equations (7) and (8)) depicted by black solid lines. (For details see Section 3.1.3.)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\Phi }}_{★,0}& =(1.22\pm 0.20)\times {10}^{-3}\hfill \\ \hfill z_0& =2.24\pm 0.23\hfill \\ \hfill a_1& =\mathrm{fixed}\hfill \\ \hfill a_2& =4.95\pm 0.73\hfill \\ \hfill b_1& =(5.26\pm 0.30)\times {10}^{11}\hfill \\ \hfill b_2& =1.38\pm 0.07,\hfill \end{array}\end{array}$$(8)

with a₁ fixed at 0.0 following Casey et al. (2018).

In the top panel of Figure 3 we present the redshift evolution of the characteristic number density, Φ_★, as a solid black line. The steep evolution found here, where log Φ_★ ∝ (1 + z)^−4.95 is close to the dust-poor model of Casey et al. (2018), who assumed log Φ_★ ∝ (1 + z)^−5.9. The redshift evolution of the characteristic luminosity, L_★, is shown in the bottom panel of Figure 3 with a solid black line. The increasing evolution of L_★ with redshift is consistent with the downsizing scenario (Thomas et al. 2010) in which the more luminous (and more massive) galaxies formed earlier than their fainter counterparts (this is also apparent in Figure 4).

Fig. 4. Redshift evolution of the functional form of the IR LF found in this work (Equations (6), (7) and (8)).

In Figures 5 and 6, we present a comparison of our functional form of the IR LF with other data and their corresponding best-fit functions from the literature (Magnelli et al. 2013; Gruppioni et al. 2013; Koprowski et al. 2017; Casey et al. 2018; Dudzevičiūtė et al. 2020; Traina et al. 2024; Fujimoto et al. 2024). Figure 5 shows that our results agree (within the errors) with the majority of the studies shown here, with notable exceptions being the Herschel-selected samples from Gruppioni et al. (2013) and the inhomogeneous ALMA COSMOS data presented by Traina et al. (2024). Previous works include discussion of the discrepancies between submillimeter-selected and Herschel-selected samples (e.g., Koprowski et al. 2017), where the low-resolution Herschel data from Gruppioni et al. (2013) is affected by blending issues, which likely leads to an overestimation of the infrared luminosity values, where for the brightest Herschel sources the 250 μm flux may be overestimated by around 150% (Scudder et al. 2016) and is also prone to active galactic nucleus contamination. The ALMA data used in Traina et al. (2024) consists of the publicly available COSMOS data and is therefore inhomogeneous in terms of observing wavelength and depth. In addition, ALMA observations tend to target the most luminous IR sources, and hence, due to clustering, the IR LFs constructed from such data are likely biased toward larger number densities. Traina et al. (2024) aimed at fixing the clustering bias of their sample by excluding central sources in all of their ALMA pointings. However, the ALMA data, by selection, target relatively dense regions in the field, likely boosting the resulting values of the IR LF. The works of Magnelli et al. (2013) and Fujimoto et al. (2024) largely agree with our results (within errors), with the exception of the brightest L_IR bin at z ∼ 2. In the case of Magnelli et al. (2013), the sample was selected at the Herschel PACS wavelength, with the Herschel fluxes affected by blending and the PACS selection bands likely contaminated by the active galactic nuclei. In the case of Fujimoto et al. (2024), the sample studied was based on the 180 dust continuum sources identified in 33 massive cluster fields and hence can also likely be biased toward larger number densities at the bright end. On the other hand, the 850 μm selection of this work will potentially miss low-dust mass, high-dust temperature sources. Although we tried to account for this effect in our completeness calculations, some sources with a low dust mass and relatively high dust temperature may still be missed in our sample.

Fig. 5. Functional form of the infrared luminosity function found in this work (black solid line). For comparison we show the individual data found in other literature studies (Magnelli et al. 2013; Gruppioni et al. 2013; Wang et al. 2019; Dudzevičiūtė et al. 2020; Traina et al. 2024, and Fujimoto et al. 2024). (For a detailed discussion, see Section 3.1.3.)

Fig. 6. Functional form of the infrared luminosity function found in this work (black solid line). Other IR LFs found in the literature are also shown for comparison (Magnelli et al. 2013; Gruppioni et al. 2013; Koprowski et al. 2017; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024, and Fujimoto et al. 2024). It can be seen that while at the faint end our results seem to be mostly consistent with other works (see Figure 5), the assumed value of the faint-end slope causes the corresponding functions to differ significantly. The shaded region represents the regime bounded at the bottom by the dust-poor models and at the top by the dust-rich models of the early Universe postulated by Casey et al. (2018). (For a detailed discussion, see Section 3.1.3.)

In Figure 5 we also plot the results of Dudzevičiūtė et al. (2020) with green crosses, where the IR LF values were determined at IR luminosities above the limit corresponding to the S₈₇₀ flux of 3.6 mJy and the median temperature of the AS2UDS sample of ∼30 K without the requirement of at least one SPIRE detection. Since we assumed a more conservative dust temperature value of 37.3 K, we effectively limited our analysis to the most luminous IR bins (Section 3.1.2). As the submillimeter flux roughly traces the dust mass of the galaxies (M_d; e.g., Dudzevičiūtė et al. 2020), our lower-limit IR luminosity bins will statistically miss a small fraction of sources with log(M_d/M_⊙)∼9.0 (S₈₇₀ ∼ 4 mJy) and T_d > 37.3 K (increasing to higher dust temperatures at z ≳ 4 due to SPIRE detection requirement). Assuming a lower dust temperature produces less luminous IR lower limits, which become increasingly less complete. The small differences at the highest L_IR bins between both works can mainly be attributed to the different redshift bins adopted but also to the fact that the sample used in Dudzevičiūtė et al. (2020) included sources without any Herschel SPIRE detections and did not adopt the completeness corrections of Equation (5) from the original SCUBA-2 sample of Geach et al. (2017).

The IR LF functional forms corresponding to the data from Figure 5 presented in Figure 6 show the impact of the assumed faint-end slope on the resulting form of the IR LF. Even though, as can be seen in Figure 5, the literature data are mostly consistent with our findings at the faint end, the assumed sα impacts the resulting functional form significantly, which in turn slightly affects the corresponding values of the star formation rate density, ρ_SFR, which we discuss in the next subsection.

3.2. Cosmic star formation history

In order to find the comoving volume density of the IR luminosity, ρ_IR, at redshift z_i, we followed

$$\begin{array}{c}\hfill {\rho }_{\mathrm{IR}}(z_i)={\displaystyle {\int }_{\log (L_{\min })}^{\log (L_{\max })}}\mathrm{\Phi }(L,z_i)\times L\mathrm{dlog}(L),\end{array}$$(9)

where Φ is the Schechter function of Equation (6), with the best-fit parameters at z_i listed in Table 4. Following previous works (e.g., Koprowski et al. 2017; Gruppioni et al. 2013; Traina et al. 2024), we integrate Equation (9) between 8.0 < log(L_IR/L_⊙) < 14.0. To find the star formation rate density, we then followed Kennicutt (1998):

$$\begin{array}{c}\hfill {\rho }_{\mathrm{SFR}}={\rho }_{\mathrm{IR}}\times {\mathcal{K}}_{\mathrm{IR}}\times 0.63,\end{array}$$(10)

with 𝒦_IR = 1.73 × 10⁻¹⁰ M_⊙ year⁻¹ L_⊙⁻¹ and an additional multiplicative factor of 0.63 to convert from a Salpeter to a Chabrier IMF. The resulting values are summarized in Table 5 and depicted as red circles with error bars in Figure 7. The red solid line depicts the redshift evolution of the star formation rate density derived from the best-fit functions of Figure 3 between the redshifts considered in this work (1.2 < z < 5.7). In order to estimate the errors on ρ_SFR at each redshift listed in Table 4, we performed Monte Carlo simulations, where in each of the 1000 realizations, the IR luminosities and their corresponding IR LF values listed in Tables 1 and 3 were randomly sampled from a Gaussian distribution with the mean and the variance equal to the calculated values and their errors, respectively. In each of the realizations, the best-fit Schechter function was found (Equation (6)), and the corresponding ρ_SFR was determined following Equations (9) and 10. The final errors on ρ_SFR at each redshift were then taken to be equal to the standard deviations of all the individual values.

Fig. 7. Star formation rate density for four redshift bins studied in this work (Table 5), depicted with red circles with error bars, found by integrating the corresponding IR LFs between 8.0 < LIR/L⊙ < 14.0. The solid red line represents the redshift evolution of the star formation rate density derived from the best fits to the Schechter function parameters of Figure 3 between redshifts considered in this work (1.2 < z < 5.7). The solid black line represents the IR ρSFR evolution found in Madau & Dickinson (2014), while the dotted blue line shows the evolution of the UV ρSFR of Dunlop (2016). Top: Comparison of our results with those presented in Dudzevičiūtė et al. (2020), where the blue and orange lines depict the contribution to ρSFR from 870 μm sources with fluxes above 3.6 mJy and 1.0 mJy (extrapolated using ALMA number counts of Hatsukade et al. 2018), respectively. The shaded areas represent 1σ errors. The green line shows the evolution down to the IR luminosity detection limits of our IR faint sample presented in Section 3.1.1. The differences in the derived values of the ρSFR at a specific redshift across the samples arise from the different LIR limits adopted (see Section 3.2). At z ≳ 4, these differences decrease, which can be attributed to the growing contribution of bright submillimeter galaxies to the ρSFR. Bottom: Comparison of ρSFR found in this work with the most recent results from the literature plotted with color symbols (Gruppioni et al. 2013; Rowan-Robinson et al. 2016; Koprowski et al. 2017; Dunlop et al. 2017; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024; Magnelli et al. 2024; Fujimoto et al. 2024; Sun et al. 2025; Liu et al. 2026). The shaded region shows the results of Casey et al. (2018) bounded at the bottom by the dust-poor models and at the top by the dust-rich models of the early Universe. It can be seen that our results are in a good agreement with the results of Dunlop et al. (2017), Wang et al. (2019), Magnelli et al. (2024), and Liu et al. (2026). Also, some inconsistencies with other works are apparent, the discussion of which is presented in Section 3.2. Considering the shaded region of Casey et al. (2018), it can be seen that the results of this work point toward the dust-poor early Universe scenario, where the IR sources are expected to dominate the total ρSFR out to z ∼ 4.

In the top panel of Figure 7, we compare our results with those presented in Dudzevičiūtė et al. (2020). The blue line (with 1σ errors represented by the shaded area) depicts the results presented in Dudzevičiūtė et al. (2020) down to the S₈₇₀ flux limit of 3.6 mJy, while the orange line shows the derived values extrapolated down to the S₈₇₀ flux limit of 1.0 mJy, where the ALMA number counts of Hatsukade et al. (2018) were used. The green line gives the evolution of ρ_SFR found in this work down to the IR luminosity detection limits of our IR faint sample presented in Section 3.1.1 (Table 1), and the red circles with 1σ error bars represent the integrated values down to log(L_IR/L_⊙) = 8.0, with the red solid line produced from the evolution of the best-fit Schechter function parameters, presented in Section 3.1.3. The increasing values of ρ_SFR at z ≲ 3 for each of the samples are simply a consequence of different lower L_IR limits adopted, where S₈₇₀ fluxes from Dudzevičiūtė et al. (2020) of 3.6 and 1.0 mJy correspond to a log(L_IR/L_⊙) of ∼12.4 and ∼11.7, respectively. For the sample of this work represented by the green line in the figure, the IR LF lower limit is set to a log(L_IR/L_⊙) of ∼10.9 at z ∼ 1.4, rising to ∼12.0 at z > 4 (Table 1), while the red circles were produced by integrating the best-fit Schechter functions down to log(L_IR/L_⊙) = 8.0. Interestingly, the differences in the derived ρ_SFR values for each sample decrease toward higher redshifts, which indicates the growing contribution of the most IR luminous sources to the total star formation rate density toward earlier epochs (which is also apparent in Figure 4).

In the bottom panel of Figure 7, we compare our findings with a range of recent results from the literature, where whenever possible, the conversion from L_IR to SFR was recalibrated using Equation (10). It can be seen that the numbers found in this work are in a good agreement with those determined in Dunlop et al. (2017), Wang et al. (2019), Magnelli et al. (2024), and Liu et al. (2026). In the case of Dunlop et al. (2017) and Liu et al. (2026), the cosmic star formation rate density was derived from individual ALMA detections. In order to account for sources that went undetected by ALMA, Dunlop et al. (2017) performed stacking down to the stellar mass limit of ∼10^9.3 M_⊙, while Liu et al. (2026) included the contribution from the ALMA-undetected galaxies down to M_* > 10¹⁰ M_⊙. Magnelli et al. (2024) estimated the ρ_SFR at z ∼ 4.5 by uv-plane stacking COSMOS sources in the available ALMA data down to M_* = 10^9.5 M_⊙. These lower limits are consistent with the recent findings in which no dust was detected in high-redshift galaxies with stellar masses below M_* ∼ 10^9.5 M_⊙ (e.g., Pannella et al. 2015; McLure et al. 2018; Koprowski et al. 2018).

It can also be seen in Figure 7 that our results are somewhat lower than those of Traina et al. (2024), Fujimoto et al. (2024), Sun et al. (2025) and significantly lower than the results found in Rowan-Robinson et al. (2016) and Gruppioni et al. (2020) at z ≳ 3. At high redshifts, the selected 350 μm and 500 μm Herschel sample of Rowan-Robinson et al. (2016) consists only of the sources with the most extreme S/N (S/N > 5), and the sources have IR luminosities considerably higher than the knee of the corresponding LF. Gruppioni et al. (2020), as can be seen in Figure 6, predicts bright-end number densities that are significantly larger than what was found in this work. In fact, their IR LF suggests approximately ten times more sources at the bright end than what is predicted by the somewhat extreme case of the dust-rich early Universe scenario of Casey et al. (2018).

The shaded region in Figure 7 depicts the ρ_SFR theoretical regime bounded from the bottom by the dust-poor scenario and from the top by the dust-rich scenario of the early Universe, as calibrated by Casey et al. (2018). In the dust-poor model, the UV-bright sources (dotted blue line in Figure 7; Dunlop 2016) dominate star formation at z ≳ 3.5, where the dust-formation timescale, driven primarily by asymptotic giant branch (AGB) stars and supernovae, is expected to be longer than the time it takes to form the first UV-bright galaxies. The dust-rich model, on the other hand, predicts the cosmic star formation to be dominated by the dusty star forming galaxies between 1.5 < z < 6.5. As noted by Casey et al. (2018), the dust-rich model is somewhat extreme, as very few DSFGs have been found at z ≳ 5, mainly due to observational limitations (see Casey et al. 2014 for details). Indeed, from Figure 7 it can be seen that the IR LF found in this work predicts number densities of DSFGs that are similar to the dust-poor model and consistent with the low-detection rates at high redshifts, placing our results in a dust-poor early Universe regime, where the IR sources seem to dominate the star formation rate density at least out to z ∼ 4.

4. Summary

We have probed the low-luminosity regime of the IR luminosity function out to z = 5.7 by stacking the optical/near-IR catalogs of McLeod et al. (2021) in the FIR Herschel and JCMT maps, using stellar mass as a proxy, where the L_IR–M_* relationship of Koprowski et al. (2024) was used. Together with the ALMA follow-up data of the S2CLS UKIDSS UDS sources (AS2UDS; Stach et al. 2019; Dudzevičiūtė et al. 2020), we have established the evolution of the functional form of IR LF between redshifts 1.2 and 5.7. We fit the data at four redshift bins (Table 4) with the Schechter function of Equation (6), and we established a value of the faint-end slope; the redshift evolution of the characteristic luminosity, L_*; and the characteristic number density, Φ_*. By integrating the best-fit functions, we have determined the values of the comoving volume density of star formation, ρ_SFR. The main results of this work can be summarized as follows:

1. The faint end of the IR luminosity function values of α = −0.26 ± 0.11 and −0.23 ± 0.18 were determined at two low-redshift bins of 1.2 ≤ z < 1.6 and 1.6 ≤ z < 2.2, respectively, and the more accurate 1.2 ≤ z < 1.6 value of –0.26 was adopted at all remaining redshift bins. The characteristic number density at z ≳ 2 is a steeply declining function of redshift, where Φ_* ∝ z^−4.95, and it is very similar to the predicted evolution for the dust-poor model of the early Universe investigated in Casey et al. (2018). The redshift evolution of the characteristic luminosity was found to be a power-law function, where L_★ ∝ z^1.38, supporting the downsizing scenario (Thomas et al. 2010), where most luminous galaxies form before their less luminous counterparts.

2. When comparing the functional form of the IR LF with the results from the recent literature, we found our faint-end slope to be in agreement with most of other calibrations. At the bright end, we predict relatively small number densities of sources, which aligns with recent studies struggling to detect high-redshift galaxies in the FIR (Dunlop et al. 2017; Aravena et al. 2020). A small number of previous studies found significantly more high-L_IR sources, which we attribute to the low resolution of the FIR observations and/or the methodology adopted when dealing with inhomogeneous FIR samples.

3. Finally, the redshift evolution of the star formation rate density, ρ_SFR, was established by integrating the corresponding IR LFs at the four redshift bins investigated in this work. The resulting ρ_SFR peaks at z ≃ 2 and later declines with redshift from z ∼ 2 out to z ∼ 6, with the high-redshift values significantly lower than those presented in Madau & Dickinson (2014). Our data are in very good agreement with the recent works of Dunlop et al. (2017), Wang et al. (2019), Magnelli et al. (2024), and Liu et al. (2026), but they are somewhat below the numbers derived in Traina et al. (2024), Fujimoto et al. (2024) and Sun et al. (2025). The most striking differences are between our results and those established by Rowan-Robinson et al. (2016) and Gruppioni et al. (2020), who significantly larger number densities have been found at the bright end. Considering the two extreme cases of a dust-rich and a dust-poor early Universe presented in Casey et al. (2018), the ρ_SFR redshift evolution found in this work places our findings in the dust-poor scenario, where we predict that the IR sources will dominate the total density of star formation in the Universe out to z ∼ 4.

Acknowledgments

This research was funded in whole or in part by the National Science Center, Poland (grants no. 2020/39/D/ST9/03078, 2023/51/B/ST9/01479, 2023/49/B/ST9/00066 and 2024/53/N/ST9/00350). 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 and DJM acknowledge the support of the Royal Society through a Research Professorship awarded to JSD. KL acknowledges the support of the National Science Centre, Poland, through the PRELUDIUM grant UMO-2023/49/N/ST9/00746. Supported by the Foundation for Polish Science (FNP).

References

Appendix A: Modified-Schechter fits

Fig. A.1. Infrared LF binned data of Figure 2 (black circles with error bars) with the best-fit functional forms overlayed. The modified Schechter fits with α and σ parameters (Equation A.1) allowed to vary freely are depicted with blue solid line. Fixing α and σ at the recent literature values of -0.2 and 0.5, respectively (Gruppioni et al. 2013; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024), produces the modified Schechter function presented with the red dashed line, while the best-fit classical Schechter curve (Equation 6) derived in this work is shown with the black dotted line.

Fig. A.2. Star formation rate density found by integrating the best-fit curves of Figure A.1 between 8.0 < LIR/L⊙ < 14.0. Similarly to Figure 7, the black solid line represents the IR ρSFR evolution found in Madau & Dickinson (2014), while the blue dotted line shown the evolution of the UV ρSFR of Dunlop (2016). The shaded region shows the results of Casey et al. (2018), bounded at the bottom by the dust-poor models and at the top by the dust-rich models of the early Universe. The black circles with the solid curve represent this work’s results found in Section 3.2. Values of ρSFR corresponding to the modified Schechter fits (Eqaution A.1) with the best-fit α and σ parameters of −0.26 and 0.31, respectively, are shown with blue squares. Data found from modified Schechter fits with α and σ fixed at recent literature z ≪ 1 values of −0.2 and 0.5, respectively, are depicted with red crosses. It can be seen that, while using the modified Schechter function with all the parameters allowed to vary freely produces results almost identical to the classical Schechter results of this work, setting α and σ at the z ≪ 1 best-fit values if Gruppioni et al. (2013), produces significantly larger numbers at highest redshift bins.

Since a number of recent works have adopted a modified version of the Schechter function (e.g., Gruppioni et al. 2013; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024), in this appendix we present an alternative IR LF functional fits and the corresponding star formation rate density calculations following the original work of Saunders et al. (1990), where

$$\begin{array}{c}\hfill \mathrm{\Phi }(L)={\mathrm{\Phi }}_{\ast }{\left(\frac{L}{L_{\ast }}\right)}^{\alpha }\exp [-\frac{1}{2{\sigma }^2}{\log }_{10}^2(1+\frac{L}{L_{\ast }})].\end{array}$$(A.1)

The free parameters Φ_* and L_* mark the so-called knee of the luminosity function, with L ≪ L_* behaving as a power-law and L ≫ L_* following a Gaussian curve. This function has proven to be useful, since it allows the control of the shape of the bright-end portion of the curve through the σ parameter. While the classical Schechter function adopted in this work does not require any modifications at the bright end, it is interesting to investigate whether the modified version produces different results in terms of ρ_SFR. In most of the high-redshift studies (e.g., Gruppioni et al. 2013; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024), both α and σ are fixed at z ≪ 1 best-fit values of −0.2 and 0.5, respectively. In Figure A.1 we show the binned data of Figure 2 with black circles with the best-fit classical Schechter function found in this work depicted with black dotted lines. The modified Schechter fits, were α and σ are fixed at −0.2 and 0.5, respectively, are presented with red dashed lines, while the similar fits with α and σ set as free parameters are shown with blue solid lines. It can be seen that fixing α and σ at z ≪ 1 values produces very poor fits at two lowest redshift bins, where the faint-end portion of the curve is best traced by the data. Allowing those parameters to vary (blue solid lines in Figure A.1) given significantly better fits, with the best-fit values of α = −0.26 and σ = 0.31.

The corresponding star formation rate densities are shown in Figure A.2, with the black circles and the solid black line representing the results of this work presented in Section 3.2. The blue squares, showing the ρ_SFR values derived from the modified Schechter fits with α and σ set as free parameters (blue solid lines in Figure A.1), are in an excellent agreement with the classical Schechter fits adopted in this work, while the data corresponding to the z ≪ 1 low-resolution Herschel modified Schechter parameters found in Gruppioni et al. (2013), give significantly higher values of ρ_SFR at high-redshift bins.

Appendix B: Additional table

All Tables

All Figures

	Fig. 1. Infrared luminosity-stellar mass relation determined through stacking the optical/near-IR mass-complete samples of McLeod et al. (2021) in the FIR Herschel and JCMT maps, adopted from Koprowski et al. (2024). The stacked LIR values were used to determine the faint-end portion of the IR LF, as explained in Section 3.1.1.
In the text

Fig. 2. Infrared luminosity function found in this work at different redshift bins, as indicated in the plot. The faint-end data (Table 1), determined through stacking (see Section 3.1.1 for details), are depicted by black points with horizontal error bars. The remaining bright-end data (Table 3) found using the ALMA follow-up data of the S2CLS UKIDSS UDS sources (AS2UDS; Stach et al. 2019) are also shown (Section 3.1.2). The best-fit Schechter functions of Equation (6) are represented by dashed black lines, with the gray area representing 1σ uncertainties. The faint-end slope was determined at two low redshift bins, and the more accurate result, α = −0.26 ± 0.11, was adopted at all the remaining redshifts, as explained in Section 3.1.3.

In the text

	Fig. 3. Redshift evolution of the Schechter function parameters found in this work. The values determined at each redshift (Table 4) are shown as black points, with the best-fit functions (Equations (7) and (8)) depicted by black solid lines. (For details see Section 3.1.3.)
In the text

	Fig. 4. Redshift evolution of the functional form of the IR LF found in this work (Equations (6), (7) and (8)).
In the text

	Fig. 5. Functional form of the infrared luminosity function found in this work (black solid line). For comparison we show the individual data found in other literature studies (Magnelli et al. 2013; Gruppioni et al. 2013; Wang et al. 2019; Dudzevičiūtė et al. 2020; Traina et al. 2024, and Fujimoto et al. 2024). (For a detailed discussion, see Section 3.1.3.)
In the text

Fig. 6. Functional form of the infrared luminosity function found in this work (black solid line). Other IR LFs found in the literature are also shown for comparison (Magnelli et al. 2013; Gruppioni et al. 2013; Koprowski et al. 2017; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024, and Fujimoto et al. 2024). It can be seen that while at the faint end our results seem to be mostly consistent with other works (see Figure 5), the assumed value of the faint-end slope causes the corresponding functions to differ significantly. The shaded region represents the regime bounded at the bottom by the dust-poor models and at the top by the dust-rich models of the early Universe postulated by Casey et al. (2018). (For a detailed discussion, see Section 3.1.3.)

In the text

Fig. 7. Star formation rate density for four redshift bins studied in this work (Table 5), depicted with red circles with error bars, found by integrating the corresponding IR LFs between 8.0 < LIR/L⊙ < 14.0. The solid red line represents the redshift evolution of the star formation rate density derived from the best fits to the Schechter function parameters of Figure 3 between redshifts considered in this work (1.2 < z < 5.7). The solid black line represents the IR ρSFR evolution found in Madau & Dickinson (2014), while the dotted blue line shows the evolution of the UV ρSFR of Dunlop (2016). Top: Comparison of our results with those presented in Dudzevičiūtė et al. (2020), where the blue and orange lines depict the contribution to ρSFR from 870 μm sources with fluxes above 3.6 mJy and 1.0 mJy (extrapolated using ALMA number counts of Hatsukade et al. 2018), respectively. The shaded areas represent 1σ errors. The green line shows the evolution down to the IR luminosity detection limits of our IR faint sample presented in Section 3.1.1. The differences in the derived values of the ρSFR at a specific redshift across the samples arise from the different LIR limits adopted (see Section 3.2). At z ≳ 4, these differences decrease, which can be attributed to the growing contribution of bright submillimeter galaxies to the ρSFR. Bottom: Comparison of ρSFR found in this work with the most recent results from the literature plotted with color symbols (Gruppioni et al. 2013; Rowan-Robinson et al. 2016; Koprowski et al. 2017; Dunlop et al. 2017; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024; Magnelli et al. 2024; Fujimoto et al. 2024; Sun et al. 2025; Liu et al. 2026). The shaded region shows the results of Casey et al. (2018) bounded at the bottom by the dust-poor models and at the top by the dust-rich models of the early Universe. It can be seen that our results are in a good agreement with the results of Dunlop et al. (2017), Wang et al. (2019), Magnelli et al. (2024), and Liu et al. (2026). Also, some inconsistencies with other works are apparent, the discussion of which is presented in Section 3.2. Considering the shaded region of Casey et al. (2018), it can be seen that the results of this work point toward the dust-poor early Universe scenario, where the IR sources are expected to dominate the total ρSFR out to z ∼ 4.

In the text

Fig. A.1. Infrared LF binned data of Figure 2 (black circles with error bars) with the best-fit functional forms overlayed. The modified Schechter fits with α and σ parameters (Equation A.1) allowed to vary freely are depicted with blue solid line. Fixing α and σ at the recent literature values of -0.2 and 0.5, respectively (Gruppioni et al. 2013; Wang et al. 2019; Gruppioni et al. 2020; Traina et al. 2024), produces the modified Schechter function presented with the red dashed line, while the best-fit classical Schechter curve (Equation 6) derived in this work is shown with the black dotted line.

In the text

Fig. A.2. Star formation rate density found by integrating the best-fit curves of Figure A.1 between 8.0 < LIR/L⊙ < 14.0. Similarly to Figure 7, the black solid line represents the IR ρSFR evolution found in Madau & Dickinson (2014), while the blue dotted line shown the evolution of the UV ρSFR of Dunlop (2016). The shaded region shows the results of Casey et al. (2018), bounded at the bottom by the dust-poor models and at the top by the dust-rich models of the early Universe. The black circles with the solid curve represent this work’s results found in Section 3.2. Values of ρSFR corresponding to the modified Schechter fits (Eqaution A.1) with the best-fit α and σ parameters of −0.26 and 0.31, respectively, are shown with blue squares. Data found from modified Schechter fits with α and σ fixed at recent literature z ≪ 1 values of −0.2 and 0.5, respectively, are depicted with red crosses. It can be seen that, while using the modified Schechter function with all the parameters allowed to vary freely produces results almost identical to the classical Schechter results of this work, setting α and σ at the z ≪ 1 best-fit values if Gruppioni et al. (2013), produces significantly larger numbers at highest redshift bins.

In the text