The evolving far-IR galaxy luminosity function and dust-obscured star-formation rate density out to z~5

TL;DR

This study measures the far-IR galaxy luminosity function up to z~5, revealing that dust-obscured star formation is minor at high redshifts and that bright IR sources are rarer than previously thought.

Contribution

It provides new measurements of the evolving far-IR luminosity function using advanced sub-mm/mm imaging, correcting overestimations from earlier Herschel surveys.

Findings

Luminosity density dominated by bright objects (~L*) at all redshifts.

High-redshift star formation is mainly UV-visible, not dust-obscured.

The LF evolution shows increasing luminosity and decreasing density with redshift.

Abstract

We present a new measurement of the evolving galaxy far-IR luminosity function (LF) extending out to redshifts z~5, with resulting implications for the level of dust-obscured star-formation density in the young Universe. To achieve this we have exploited recent advances in sub-mm/mm imaging with SCUBA-2 on the James Clerk Maxwell Telescope (JCMT) and the Atacama Large Millimeter/Submillimeter Array (ALMA), which together provide unconfused imaging with sufficient dynamic range to provide meaningful coverage of the luminosity-redshift plane out to z>4. Our results support previous indications that the faint-end slope of the far-IR LF is sufficiently flat that comoving luminosity-density is dominated by bright objects (~L*). However, we find that the number-density/luminosity of such sources at high redshifts has been severely over-estimated by studies that have attempted to push the highly-confused Herschel SPIRE surveys beyond z~2. Consequently we confirm recent reports that cosmic star-formation density is dominated by UV-visible star formation at z>4. Using both direct (1/Vmax) and maximum likelihood determinations of the LF, we find that its high-redshift evolution is well characterized by continued positive luminosity evolution coupled with negative density evolution (with increasing redshift). This explains why bright sub-mm sources continue to be found at z>5, even though their integrated contribution to cosmic star-formation density at such early times is very small. The evolution of the far-IR galaxy LF thus appears similar in form to that already established for active galactic nuclei, possibly reflecting a similar dependence on the growth of galaxy mass.

Tables

Table 1. Table 1: The refined SCUBA-2 850- μ ​ m 𝜇 m {\rm\mu m} survey fields and source samples used in this work (see text for details). The columns show the field name, the area, limiting flux density, limiting signal-to-noise ratio and the resulting number of detections.

	Area	$S_{\lim }$	SNR	N
	/deg-2	/mJy
COSMOS deep	$0.08$	$4.60$	4.2	34
COSMOS wide	$0.79$	$6.50$	4.3	89
UDS	$0.71$	$3.75$	4.2	454

Table 2. Table 2: Rest-frame 250- μ ​ m 𝜇 m {\rm\mu m} luminosity functions.

log($L_{250}/W{\mathrm{Hz}}^{-1}$)	log($\Phi /{\mathrm{Mpc}}^{-3}{\mathrm{dex}}^{-1}$)
$24.8$	…	$-2.81\pm 0.16$	…	…
$25.0$	…	$-3.20\pm 0.23$	…	…
$25.6$	$-4.23\pm 0.08$	$-3.86\pm 0.05$	…	…
$25.7$	$-4.68\pm 0.11$	$-4.17\pm 0.05$	$-4.40\pm 0.08$	$-4.47\pm 0.12$
$25.8$	$-5.14\pm 0.13$	$-4.42\pm 0.06$	$-4.56\pm 0.08$	$-5.27\pm 0.18$
$25.9$	$-5.42\pm 0.16$	$-4.76\pm 0.07$	$-4.58\pm 0.07$	$-4.93\pm 0.11$
$26.0$	$-6.11\pm 0.30$	$-5.18\pm 0.11$	$-5.05\pm 0.09$	$-5.37\pm 0.14$
$26.1$	…	$-5.83\pm 0.20$	$-5.46\pm 0.14$	$-5.49\pm 0.15$
$26.2$	…	$-6.31\pm 0.30$	$-5.83\pm 0.20$	$-5.97\pm 0.23$
$26.3$	…	…	$-6.31\pm 0.30$	…

Table 3. Table 3: The best-fitting parameter values for the Schechter functions as determined by the 1 / V max 1 subscript 𝑉 max 1/V_{\rm max} method. The faint-end slope, α 𝛼 \alpha , is fixed here at the z ≃ 2 similar-to-or-equals 𝑧 2 z\simeq 2 value (see text for details).

$z$	$\alpha$	log(${\Phi }^{\star }/{\mathrm{Mpc}}^{-3}{\mathrm{dex}}^{-1}$)	log($L^{\star }/W{\mathrm{Hz}}^{-1}$)
	$-0.4$	$-{2.88}_{-0.27}^{+0.33}$	${25.20}_{-0.09}^{+0.08}$
	$-0.4$	$-{3.03}_{-0.10}^{+0.05}$	${25.40}_{-0.03}^{+0.03}$
	$-0.4$	$-{3.73}_{-0.16}^{+0.13}$	${25.63}_{-0.05}^{+0.05}$
	$-0.4$	$-{4.59}_{-0.31}^{+0.28}$	${25.84}_{-0.14}^{+0.16}$

Table 4. Table 4: Best-fit values for the continuous functions of the Schechter function parameters from Equation 8 , with 1 ​ σ 1 𝜎 1\sigma errors based on the maximum-likelihood ratios and corresponding to Δ ​ χ 2 = 1 Δ superscript 𝜒 2 1 \Delta\chi^{2}=1 .

Parameter	Value
$A$	${2.40}_{-0.48}^{+0.60}\times {10}^{-3}$
$\sigma$	${1.04}_{-0.05}^{+0.06}$
$\mu$	${1.28}_{-0.25}^{+0.20}$
$a$	${0.19}_{-0.02}^{+0.02}$
$b$	${25.03}_{-0.05}^{+0.06}$

Table 5. Table 5: The predicted number of 850- μ ​ m 𝜇 m {\rm\mu m} sources with flux densities S 850 > 10 subscript 𝑆 850 10 S_{850}>10  mJy in the combined area of all the fields used in this study ( ≃ 1.5 similar-to-or-equals absent 1.5 \simeq 1.5 deg 2 ) produced by integrating the three IR LFs discussed in Section 6.1 and shown in Fig. 7 and Fig. 8 . As discussed in the caption to Fig. 9 , and indicated here, these predictions were based on integration over the redshift range 0.5 < z < 4.5 0.5 𝑧 4.5 0.5<z<4.5 or 0.5 < z < 2.3 0.5 𝑧 2.3 0.5<z<2.3 as appropriate. The actual number of sources in our sample with S 850 > 10 subscript 𝑆 850 10 S_{850}>10  mJy is 17 in the redshift range 0.5 < z < 4.5 0.5 𝑧 4.5 0.5<z<4.5 , of which four have 0.5 < z < 2.3 0.5 𝑧 2.3 0.5<z<2.3 . The predictions are based on the Michałowski et al. ( 2010 ) average sub-mm source SED. Clearly, the Herschel -derived LFs wildly over-predict the actual number of bright 850- μ 𝜇 \mu m sources.

Study	z range	$N\left(>10\mathrm{mJy}\right)$
		$/\mathrm{\hspace{0.17em}1.5}{\mathrm{deg}}^{-2}$
This Work		20
Gruppioni et al. (2013)		315
This Work		8
Magnelli et al. (2013)		142

Equations

$$\Phi (L,z)=\frac{1}{\Delta L}\sum _i\frac{1-\mathrm{FDR}}{w_i\times V_{\max ,i}},$$

$$\mathrm{lo}{\mathrm{g}}_{10}(\mathrm{FDR})=2.67-0.97\times \mathrm{SNR}.$$

$$V_{\max ,i}=\sum _j\frac{{\Omega }_j}{4\pi }{V}_{\max ,j},$$

$${\Phi }_{\mathrm{Sch}}(L,z)={\Phi }_{\star }{\left(\frac{L}{L_{\star }}\right)}^{\alpha }exp\left(\frac{-L}{L_{\star }}\right),$$

$$\mathcal{L}=\prod _i^N{\lambda }_i{e}^{-{\lambda }_i}\prod _j{e}^{-{\lambda }_j},$$

$$\lambda =\Phi (z,L)\Omega (z,L)\frac{dV}{dz}dzdL,$$

$$\begin{gathered} S=-2\ln (\mathcal{L})=-2\sum _i^N\ln [\Phi (z_i,L_i)]+ \\ 2\iint \Phi (z,L)\Omega (z,L)\frac{dV}{dz}dzdL, \end{gathered}$$

$${\Phi }_{\star }(z|A,\sigma ,\mu )$$

$$\log (L_{\star }(z|a,b))$$

$${\rho }_{\mathrm{SF}{\mathrm{R}}_{\mathrm{IR}}}(z)=\frac{0.18}{1.04\sqrt{2\pi }}{e}^{-\frac{(z-1.77{)}^2}{2\times 1.04^2}}.$$

$${\rho }_{\mathrm{SF}{\mathrm{R}}_{\mathrm{UV}}}(z)=\frac{0.11}{1.48\sqrt{2\pi }}{e}^{-\frac{(z-2.75{)}^2}{2\times 1.48^2}}.$$

Full text

††pubyear: 2017††pagerange: The evolving far-IR galaxy luminosity function and dust-obscured star-formation rate density out to $\mn@boldsymbol{z\simeq 5}$–The evolving far-IR galaxy luminosity function and dust-obscured star-formation rate density out to $\mn@boldsymbol{z\simeq 5}$

1 Introduction

A key challenge in modern astrophysical cosmology is to complete our knowledge of cosmic star-formation history, taking proper account of dust-obscured activity (e.g. Madau & Dickinson 2014). Achieving this requires a reliable measurement of the form and evolution of both the rest-frame UV and rest-frame far-IR galaxy luminosity functions (LFs) out to the highest redshifts. This is because a fair census of both unobscured and dust-obscured comoving star-formation rate density ($\rho_{\rm SFR}$) requires the luminosity-weighted integration of the relevant LFs over sufficient dynamic range to properly account for the contributions of the brightest sources, while at the same time quantifying the impact of the adopted lower luminosity limit for the integration (i.e. reliably establishing the faint-end slope of the LF).

In recent years this goal has been largely achieved at rest-frame UV wavelengths over nearly all of cosmic history, through observations with the Hubble Space Telescope (HST) and wider-area ground-based imaging from Subaru, CFHT, VLT and VISTA (Cucciati et al., 2012; McLure et al., 2013; Bowler et al., 2014, 2015; Bouwens et al., 2015, 2016a; Finkelstein et al., 2015; Parsa et al., 2016). Indeed, meaningful disagreements over the comoving UV luminosity density ($\rho_{\rm UV}$) produced by the evolving galaxy population are now largely confined to extreme redshifts $z>8$ (Ellis et al., 2013; Oesch et al., 2014; McLeod et al., 2015), and even here much of the claimed disagreement can be removed by the adoption of consistent limits to the LF integration (McLeod et al., 2016; Ishigaki et al., 2017). The reason that the adopted faint integration limit becomes an issue for calculating $\rho_{\rm UV}$ at extreme redshifts (and the resulting contribution of the emerging early galaxy population to reionization; Robertson et al. 2013; Robertson et al. 2015) is that the faint-end slope ($\alpha$) of the galaxy UV LF steepens with increasing redshift, from relatively modest values at intermediate redshifts (e.g. $\alpha\simeq-1.3$ at $z\simeq 2$; Parsa et al. 2016; Mehta et al. 2017) to $\alpha\simeq-2$ at $z>7$ (McLure et al., 2013; Bouwens et al., 2015; McLeod et al., 2015; McLeod et al., 2016). As a result, despite the discovery of significant numbers of UV-bright galaxies at $z>6$ (Bowler et al., 2014, 2015), $\rho_{\rm UV}$ within the first $\simeq$ Gyr of cosmic time is dominated by the contributions of the numerous faintest galaxies, and so the derived value depends more critically on how far down in luminosity the LF extends than is the case at later times (see Parsa et al. 2016).

Despite the long-established importance of the far-IR background (Dole et al., 2006), similar progress in our knowledge of the far-IR LF has been hampered by the heightened observational challenges at mid/far-IR and sub-mm/mm wavelengths (i.e. high background and poor angular resolution) and a resulting lack of survey dynamic range (coupled with uncertainties in redshift content). Nevertheless, important progress has been made over the past decade, first with mid-IR observations using NASA’s Spitzer Space Telescope, and more recently through far-IR imaging with ESA’s Herschel Space Observatory. Together, these facilities have enabled the far-IR LF and its evolution to be successfully traced out to $z\simeq 2$.

First, Spitzer MIPS 24 ${\rm\mu m}$ imaging was used to study the mid-IR LF out to $z\simeq 2$ (Le Floc’h et al., 2005; Caputi et al., 2007; Rodighiero et al., 2010), albeit extrapolation from 24 ${\rm\mu m}$ to far-IR luminosity becomes increasingly dangerous with increasing wavelength (although see Elbaz et al. 2010). Attempts were also made to exploit the 70 ${\rm\mu m}$ imaging provided by Spitzer (Magnelli et al., 2009; Patel et al., 2013), but sensitivity/resolution limitations largely restricted the usefulness of this work to $z<1$ (although Magnelli et al. 2011 pushed out to $z\simeq 2$ via stacking).

Over the last five years, Herschel PACS and SPIRE surveys have enabled this work to be developed through object selection at more appropriate far-IR wavelengths (i.e. closer to the rest-frame peak of the far-IR emission). Magnelli et al. (2013) used deep Herschel PACS imaging from the PEP and GOODS surveys to determine the bright-end of the far-IR LF out to $z\simeq 2$, while Gruppioni et al. (2013) used the PEP PACS imaging and SPIRE HerMES imaging (at 250, 350 and 500 ${\rm\mu m}$) to try to extend this work out to $z\simeq 4$ (see also Burgarella et al. 2013). The inclusion of the 250, 350 and 500 ${\rm\mu m}$ data allowed Gruppioni et al. (2013), in principle, to determine far-IR luminosities without recourse to large extrapolations, but object selection was still undertaken at the shorter (PACS) wavelengths. As a result, this study is mostly sensitive to the lower-redshift/warmer sources, and hence, at $z>2$, allows only the detection of the most extreme sources. One consequence of the resulting lack of dynamic range is that the faint-end slope of the far-IR LF could not be measured at high redshift, and so Gruppioni et al. (2013) simply adopted the $z=0$ value at all higher redshifts. Most recently, Rowan-Robinson et al. (2016) attempted to expand on the Gruppioni et al. (2013) study, and to extend it out to even higher redshifts ($z\simeq 6$) by including object selection at 500 ${\rm\mu m}$ (the longest-wavelength Herschel SPIRE imaging band). This study yielded surprisingly high estimates of far-IR luminosity density at high redshifts; however, utilising the longest-wavelength Herschel data in this way is fraught with danger due to the large beamsize ($\simeq 36$ arcsec FWHM) and consequent issues regarding blending, source mis-identification, potential AGN contamination and gravitational lensing (e.g. star-formation rates, SFR, as high as 20,000 ${\rm M_{\odot}yr^{-1}}$ are reported in this work).

Thus, despite this impressive progress, it is clear that attempting to reliably measure the far-IR galaxy LF based on object selection in Spitzer or Herschel surveys becomes increasingly problematic beyond $z\simeq 2$. It is now possible to overcome these difficulties at high redshift by moving to ground-based sub-mm/mm object selection, using a combination of wide-area imaging surveys as produced by SCUBA-2 (Holland et al., 2013) on the JCMT (e.g. Geach et al. 2013, 2017; Roseboom et al. 2013; Chen et al. 2016; Michałowski et al. 2016; Cowie et al. 2017), and higher-resolution smaller-area mapping as can now be achieved with ALMA (e.g. Hatsukade et al. 2015; Hatsukade et al. 2016; Umehata et al. 2015; Walter et al. 2016; Dunlop et al. 2017). While attempts have previously been made to explore the basic high-redshift evolution of the far-IR population based on ground-based sub-mm data (e.g. Wall, Pope & Scott 2008), until now the ‘wedding-cake’ of surveys providing unconfused imaging over a reasonable dynamic range was not of sufficient quality to enable a detailed investigation of the form and evolution of the far-IR LF.

The merits of moving to ground-based sub-mm (i.e. 450 and 850 ${\rm\mu m}$) and mm (i.e. 1.1–1.3 mm) selection become increasingly obvious as one moves to higher redshifts. First, the smaller beam sizes at these wavelengths offered by large ground-based single-dish telescopes such as the JCMT, or interferometric arrays such as ALMA, substantially reduce/remove the confusion limitations of the longer-wavelength Herschel SPIRE imaging, minimizing problems of source blending, false counterpart identification, and hence potentially erroneous redshift information. Second, certainly by $z>3$, object selection at sub-mm/mm wavelengths is much less susceptible to AGN contamination than Spitzer or Herschel PACS detections (which sample the spectral energy distribution (SED) shortward of $\lambda_{\rm rest}\simeq 40\,{\rm\mu m}$ at these redshifts), and extrapolation to estimate total far-IR luminosities also becomes less problematic. Third, object selection in the higher-resolution ground-based sub-mm/mm imaging can be used to help deconfuse and deblend the existing Herschel far-IR imaging, hence enabling its exploitation for SED determination in a less biased way.

A number of recent studies have already provided indications that the high values of $\rho_{\rm FIR}$ inferred from pushing the Herschel surveys beyond $z\simeq 2.5$ are incorrect. First, estimates of dust-obscured star-formation activity in known high-redshift galaxies, either derived from analyses of the UV slope (e.g. Dunlop et al. 2013) or from ALMA detections/limits (individual or stacked; e.g. Capak et al. 2015; Bouwens et al. 2016b; Koprowski et al. 2016b) suggest that the dust-obscured contribution to $\rho_{\rm SFR}$ at the highest redshifts is relatively small. However, it can reasonably be argued that the most dust-obscured galaxies will not feature in rest-frame UV-selected samples. More significantly, the results from the deep ALMA imaging of the Hubble Ultra Deep Field (HUDF; Dunlop et al. 2017) and the results of stacking in the deepest SCUBA-2 Cosmology Legacy Survey (S2CLS) images (Bourne et al., 2017) both yield much lower values of dust-obscured $\rho_{\rm SFR}$ than derived by Gruppioni et al. (2013) and Rowan-Robinson et al. (2016), especially at $z>3$. Nonetheless, it could still be argued that these small-field experiments, while complete to relatively low dust-obscured SFRs ($\simeq 10\,{\rm M_{\odot}yr^{-1}}$), cover insufficient area to reveal the contributions from the most extreme objects.

In this study we aimed to clarify this situation, and moreover to properly determine the form and high-redshift evolution of the far-IR galaxy LF by analysing the results from the HUDF ALMA 1.3-mm and deep S2CLS 850-${\rm\mu m}$ imaging in tandem with results from the wider-area $850-{\rm\mu m}$ maps recently completed within the S2CLS (Geach et al., 2017; Michałowski et al., 2016). Together these data provide a sub-mm/mm survey ‘wedding-cake’ with sufficient dynamic range and sufficiently-complete redshift content to enable a meaningful measurement of the form and evolution of the rest-frame 250-${\rm\mu m}$ galaxy LF from the available coverage of the luminosity-redshift plane. Crucially, the bright tier of S2CLS covers over 2 deg2, and yields a sample of $>1000$ luminous sources (SFR $>300\,{\rm M_{\odot}yr^{-1}}$), providing excellent sampling of the bright end of the far-IR LF, while the deeper S2CLS imaging and the ALMA imaging enable the first direct measurement of the slope of the faint-end of the far-IR LF at high redshift. Together these data have enabled us to determine the form and evolution of the far-IR LF, and hence the dust-obscured $\rho_{\rm SFR}$ out to $z\simeq 4.5$.

This paper is structured as follows. The sub-mm/mm imaging utilised in this work, along with the supporting multi-wavelength data, are described in the Section 2. The multi-wavelength identification process, together with the methods used to establish redshifts for the entire sample are then presented in Section 3. Next, the procedure for determining the IR LFs, using both the $1/V_{\rm max}$ and the maximum-likelihood methods, is explained in Section 4. The resulting calculation of the far-IR and total (UV+far-IR) $\rho_{\rm SFR}$ is presented in Section 5. We discuss our results in the context of other recent studies in Section 6, and conclude by summarising our findings in Section 7.

Throughout the paper we use a Chabrier (2003) stellar initial mass function (IMF) and assume a flat cosmology with $\Omega_{\rm m}=0.3$, $\Omega_{\Lambda}=0.7$ and H0 = 70 km s*-1* Mpc*-1*.

2 Data

2.1 JCMT SCUBA-2 imaging

We used the data collected as a part of the SCUBA-2 Cosmology Legacy Survey (S2CLS). The map-making process and the resulting derived source catalogues are described in Geach et al. (2017). The fields utilised here are the UKIDSS-UDS, where the 850-${\rm\mu m}$ imaging covers $\simeq 0.9$ deg2 with a 1-$\sigma$ noise of 0.9 mJy (revealing 1085 sources with a signal-to-noise ratio SNR $>3.5$), and the COSMOS field, where the 850-${\rm\mu m}$ imaging covers $\simeq 1.3$ deg2 with the $1\sigma$ noise of 1.6 mJy (revealing 719 sources with SNR $>3.5$). These source catalogues are discussed further in Chen et al. (2016) and Michałowski et al. (2016). For simplicity we reduced the effective area of these maps to regions of uniform noise, and for our analysis retained only sources with a simulated completeness $>0.5$. The resulting refined effective survey areas, flux-density limits, SNR, and sample sizes are summarized in Table 1.

2.2 ALMA imaging

To help inform the measurement of the faint-end slope of the LF, we used the ALMA 1.3-mm imaging of the HUDF undertaken by Dunlop et al. (2017). A mosaic of 45 ALMA pointings was created to cover the full $\simeq 4.5$ arcmin2 area previously imaged with WFC3/IR on HST. The ALMA map reached a noise level of $\sigma_{1.3}\simeq 35\,{\rm\mu Jy\,beam^{-1}}$, and 16 sources were detected with flux densities $S_{1.3}>120\,{\rm\mu Jy}$. 13 of the 16 sources have spectroscopic redshifts and the remaining three have accurate photometric redshifts derived from the optical–near-IR photometry (for the sources of the spectroscopic redshifts see table 2 in Dunlop et al. 2017).

2.3 Ancillary data

For the S2CLS COSMOS field the optical to mid-IR data consist of imaging from the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS; Gwyn 2012), as described by Bowler et al. (2012), Subaru (Taniguchi et al., 2007), the HST Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011), UltraVISTA Data Release 2 (McCracken et al., 2012; Bowler et al., 2014) and Spitzer (S-COSMOS; Sanders et al. 2007). At radio wavelengths the VLA-COSMOS Deep (Schinnerer et al., 2010) catalogues were utilised.

For the S2CLS UDS field the optical data were obtained with Subaru/SuprimeCam (Miyazaki et al., 2002), as described in Furusawa et al. (2008), the near-IR imaging was provided by the UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al., 2007; Cirasuolo et al., 2010), the mid-IR data are from the Spitzer Public Legacy Survey of the UKIDSS Ultra Deep Survey (SpUDS; PI: J. Dunlop), as described in Caputi et al. (2011), and the radio (VLA) data are from Ivison et al. (2005, 2007) and Arumugam et al. (in preparation).

For both the UDS and COSMOS fields, far-IR imaging from Herschel (Pilbratt et al., 2010) was utilised, as provided by the public releases of the HerMES (Oliver et al., 2012) and PEP (Lutz et al., 2011) surveys undertaken with the SPIRE (Griffin et al., 2010) and PACS (Poglitsch et al., 2010) instruments.

For the extraction of the far-IR flux densities and limits, the Herschel maps at 100, 160, 250, 350 and 500 ${\rm\mu m}$ were utilised, with beam sizes of 7.4, 11.3, 18.2, 24.9, and 36.3 arcsec, and 5$\sigma$ sensitivities of 7.7, 14.7, 8.0, 6.6, and 9.5 mJy, respectively. The Herschel flux densities of (or upper limits for) the SCUBA-2 sources were obtained in the following way. Square image cut-outs of width 120 arcsec were extracted from each Herschel map around each SCUBA-2 source, and the PACS (100 and 160 ${\rm\mu m}$) maps were used to simultaneously fit Gaussians with the FWHM of the respective imaging, centred at all radio and 24-${\rm\mu m}$ sources located within these cut-outs, and at the positions of the SCUBA-2 optical identifications (IDs, or just submm positions if no IDs were selected). Then, to deconfuse the SPIRE (250, 350 and 500 ${\rm\mu m}$) maps in a similar way, the positions of the 24 ${\rm\mu m}$ sources detected with PACS (at $>3\sigma$) were retained, along with the positions of all radio sources, and the SCUBA-2 IDs (or, again, the SCUBA-2 position in the absence of any optical or radio ID).

3 Redshifts

3.1 Multi-wavelength identifications

Because of the beam size of the JCMT SCUBA-2 imaging at 850 ${\rm\mu m}$ (FWHM $\simeq 15$ arcsec), we cannot simply adopt the closest optical/near-IR neighbour as the 850-${\rm\mu m}$ source galaxy counterpart. Instead, we used the method outlined in Downes et al. (1986) (see their Section 5 for a derivation), where we adopt a $2.5\sigma$ search radius around the SCUBA-2 position based on the signal-to-noise ratio (SNR): $r_{\rm s}=2.5\times 0.6\times\rm{FWHM}/\rm{SNR}$. In order to account for systematic astrometry shifts (caused by pointing inaccuracies and/or source blending; e.g. Dunlop et al. 2010) we enforced a minimum search radius of 4.5 arcsec. Within this radius we calculated the corrected Poisson probability, $p$, that a given counterpart could have been selected by chance.

Three imaging wavebands were used when searching for galaxy counterparts: the VLA 1.4-GHz imaging, the Spitzer MIPS 24-$\mu$m imaging, and the Spitzer IRAC 8-$\mu$m imaging. The radio band traces recent star formation via synchrotron radiation from relativistic electrons produced within supernovae (SNe; Condon 1992), whereas the 24-$\mu$m waveband is sensitive to the emission from warm dust. Therefore, since submm-selected galaxies are dusty, highly-star-forming objects, they are expected to be very luminous in these bands. Also, at the redshifts of interest, the 8-${\rm\mu m}$ waveband traces the rest-frame near-IR light coming from the older, mass-dominant stellar populations in galaxies, and thus provides a proxy for stellar mass. Given the growing evidence that sub-mm galaxies are massive, it is expected that 850-${\rm\mu m}$ sources will have significant 8-${\rm\mu m}$ fluxes (e.g. Dye et al. 2008; Michałowski et al. 2010; Biggs et al. 2011; Wardlow et al. 2011). Moreover, the surface density of sources in these three wavebands is low enough for chance positional coincidences to be rare (given a sufficiently small search radius). Once the counterparts were found in each of these bands, they were matched with the optical/near-IR catalogues using a search radius of $r=1.5$ arcsec and the closest object taken to be the galaxy counterpart (for a more detailed description of this process and the relevant identification tables see Michałowski et al. 2016).

3.2 Redshift distributions

We used the available multi-wavelength data to derive the optical/near-IR photometric redshifts for the SCUBA-2 source galaxy counterparts using a $\chi^{2}$ minimisation method (Cirasuolo et al. 2007; Cirasuolo et al. 2010) with a code based on the HYPERZ package (Bolzonella et al., 2000). To create templates of galaxies, the stellar population synthesis models of Bruzual & Charlot (2003) were applied, using the Chabrier (2003) stellar initial mass function (IMF), with a lower and upper mass cut-off of $0.1$ and $100\,{\rm M_{\odot}}$ respectively. Single and double-burst star-formation histories with a fixed solar metallicity were used. Dust reddening was taken into account using the Calzetti et al. (2000) law within the range $0\leq A_{V}\leq 6$. The HI absorption along the line of sight was applied according to Madau (1995). The accuracy of the photometric catalogue of Cirasuolo et al. (2010) is excellent, with a mean $|z_{\rm phot}-z_{\rm spec}|/(1+z_{\rm spec})=0.008\pm 0.034$. Also, for every source the ‘long-wavelength’ photometric redshift (e.g. Aretxaga et al. 2007; Koprowski et al. 2014; Koprowski et al. 2016a) was calculated using the SCUBA-2 and Herschel data by fitting the average sub-mm galaxy SED template of Michałowski et al. (2010).

In addition, as explained in detail in Koprowski et al. (2014); Koprowski et al. (2016a), we further tested the robustness of our optical identifications by comparing the optical photometric redshifts with the ‘long-wavelength’ photometric redshifts. Sources with optical redshifts that transpired to be significantly lower than their ‘long-wavelength’ ones were considered to be incorrectly identified in the optical (foreground galaxies, possibly lenses) and hence their ‘long-wavelength’ photometric redshifts were adopted (see Michałowski et al. 2016 for details).

The final redshift distributions (100% complete) for the refined fields used in this work (Table 1) are shown in Figure 1. Each coloured histogram shows the redshift distribution for the relevant field (as depicted in the legend), with the background black histogram showing the results for all the fields combined. The mean redshifts are $\bar{z}=2.30\pm 0.23$, $\bar{z}=3.05\pm 0.17$ and $\bar{z}=2.70\pm 0.07$ for the COSMOS deep, COSMOS wide and UDS fields respectively. The mean redshift for the whole sample used here is $\bar{z}=2.73\pm 0.06$.

4 IR Luminosity Function

To derive the evolving far-IR LF, $\Phi_{\rm IR}$, we use two independent methods. One is the standard $1/V_{\rm max}$ method (Schmidt, 1968), which allows the calculation of the LFs to be performed directly from the data, without any assumptions regarding the functional shape. To derive the functional form we then fit a set of Schechter functions, where the best-fitting parameter values are found by minimising $\chi^{2}$. In order to find the continuous form of the redshift evolution of the far-IR LF, we additionally use the maximum-likelihood method presented in Marshall et al. (1983), where we again take the LFs to be of a Schechter form.

4.1 $\boldsymbol{1/V_{\rm max}}$ method

The LF in a given luminosity and redshift bin is calculated using:

[TABLE]

where $\Delta L$ is the width of the luminosity bin, FDR is the false detection rate, $w_{\rm i}$ is the completeness for the $i$-th galaxy and $V_{\rm max,i}$ is the co-moving volume available to the $i$-th source. For S2CLS sources the false detection rate is (from Geach et al. 2017):

[TABLE]

The errors on the LFs were calculated using the Poissonian approach. The available co-moving volume is the volume between $z_{\rm min}$ and $z_{\rm max}$, where $z_{\rm min}$ is just the lower boundary of a given redshift bin and $z_{\rm max}$ is the maximum redshift at which the $i$-th source would still be visible in a given S2CLS map, or simply the upper boundary of a given redshift bin, whichever is lower. Therefore:

[TABLE]

where we sum over all the available S2CLS fields and $\Omega_{j}$ is the solid angle subtended by the $j$-th field on the sky. The contribution from the $j$-th field to the $V_{{\rm max},i}$, in a given redshift bin, $V_{{\rm max},j}$, is only counted if the maximum redshift at which the $i$-th source would still be visible in that field is higher than the lower boundary of that redshift bin, otherwise $V_{{\rm max},j}$ for that field is simply equal to 0.

Since in practice we are working with 850 ${\rm\mu m}$-detected sources with a mean redshift of $\sim 2.5$ at 850 ${\rm\mu m}$, we decided to initially calculate the LF at a rest-frame wavelength of 250 ${\rm\mu m}$. In order to do so, the luminosity at $850\,{\rm\mu m}\,/\,(1+z)$ is calculated and interpolated to $\lambda_{\rm rest}=250$ ${\rm\mu m}$ using the average SMG template from Michałowski et al. (2010). The range of luminosities and redshifts for which the LFs were calculated were determined from the coverage of the luminosity-redshift plane (grey rectangles in Fig. 2). The luminosity bins used for each field at a given redshift depend on the luminosity lower limit (corresponding to the flux-density detection limit), below which no sources can be detected (solid coloured lines in Fig. 2). Only the luminosity bins with complete luminosity coverage are included in the analysis.

The results, with Poissonian errors, are listed in Table 2, and plotted in Fig. 3. The coloured solid lines are the best-fitting Schechter functions:

[TABLE]

with $\Phi_{\star}$ being the normalisation parameter, $\alpha$ the faint-end slope and $L_{\star}$ the characteristic luminosity that roughly marks the border between the power-law fit, $(L/L_{\star})^{\alpha}$, and the exponential fit. The black solid lines (almost perfectly aligned with the coloured solid lines) with errors (grey-shaded area) depict the LFs as determined by the maximum-likelihood method (see next subsection). In order to find the faint-end slope, $\alpha$, we utilised the ALMA data (Section 2.2). Since the majority of ALMA sources lie at redshifts $1.5<z<2.5$, it was decided to determine $\alpha$ only for that redshift bin and then use that value in the remaining bins (two faintest points in $1.5<z<2.5$ redshift bin of Fig. 3). Fitting the Schechter function to the $1.5<z<2.5$ data yielded a faint-end slope of $\alpha=-0.4$. The remaining Schechter-function parameters were determined by minimising $\chi^{2}$. In Fig. 4 the 68% and 95% confidence intervals are shown, corresponding to $\Delta\chi^{2}=1$ and 4 respectively; the errors on the best-fiting Schechter-function parameters can be determined by projecting the contours onto the relevant axis. The best-fitting values for the Schechter-function parameters are given in Table 3, and plotted in Fig. 5 as blue circles.

4.2 Maximum-likelihood method

The likelihood function used here (Marshall et al., 1983) is defined as a product of the probabilities of observing exactly one source in $dzdL$ at the position of the $i$-th galaxy ($z_{i},L_{i}$) for $N$ galaxies in our sample and of the probabilities of observing zero sources in all the other differential elements in the luminosity-redshift plane. Using Poisson probabilities, the likelihood is:

[TABLE]

where $j$ runs over all differential elements in which no sources were observed, and $\lambda$ is the expected number of galaxies in $dzdL$ at $z,L$:

[TABLE]

with $\Omega$ being the fractional area of the sky in which a galaxy with a given $z$ and $L$ can be detected in our fields. Since for large $N$ the test statistic $-2\,{\rm ln}(\mathcal{L})$ will be $\chi^{2}$ distributed we define:

[TABLE]

where we dropped terms independent of the model parameters. This step transforms both exponents in Equation 5 into the integral, which represents the model $-$ the expected number of sources within the integral limits. This means that, as compared to the $1/V_{\rm max}$ method, this technique can analyse the entire luminosity-redshift plane (including its empty patches) and therefore give a result with higher statistical significance. In addition, it does not bin the data, and also allows a user to define the parameters of a luminosity function, $\Phi(z,L)$, to be continuous functions of redshift, which can be then determined by minimising $S$.

For consistency, we chose $\Phi(z,L)$ to have the form of a Schechter function (Equation 4). As before, the faint-end slope, $\alpha$, was fixed at the $1.5<z<2.5$ value of $-0.4$. As we sought an acceptable description of the data, we explored various functional forms for the redshift dependence of $\Phi_{\star}$ and log($L_{\star}$), and in the end found that we could adopt a simple normalised Gaussian function111The normalised Gaussian was chosen here instead of the linear function to ensure that $\Phi_{\star}$ does not become negative at high-$z$. and a linear function of redshift respectively:

[TABLE]

The redshift limits of the integral in Equation 7 are the same as in the $1/V_{\rm max}$ method, ${\rm z}=0.5-4.5$. The luminosity limits were chosen in order to cover the whole useful luminosity range (including the ALMA data), ${\rm log_{10}}(L_{250}\,/\,{\rm W\,Hz^{-1}})=24.6-26.6$. Minimising $S$ gave the best-fitting parameter values, as summarized in Table 4, which were then used to determine the redshift evolution of the best-fitting Schechter-function model parameters, $\Phi_{\star}$ and $L_{\star}$, as depicted in Fig. 5 by the blue solid lines. Since, for the large number of sources we have here, $S=-2\,{\rm ln}(\mathcal{L})$ is $\chi^{2}$ distributed, the $1\sigma$ errors (grey area in Fig. 5) are based on the likelihood ratios and correspond to $\Delta\chi^{2}=1$. For comparison the best-fitting values (with $1\sigma$ errors) for four redshift bins, calculated using the $1/V_{\rm max}$ method, are also plotted here as blue points.

With the redshift evolution of the Schechter-function parameters as described in Equation 8, it is straightforward to calculate the luminosity function at any redshift between $z=0.5$ and $z=4.5$. We therefore plot the LFs determined using the maximum-likelihood method alongside the results of the previous subsection in Fig. 3 as black solid lines. The $1\sigma$ errors (grey area) are again based on the maximum-likelihood ratios and correspond to $\Delta\chi^{2}=1$.

5 Star Formation Rate Density

Having constructed the rest-frame 250-${\rm\mu m}$ LFs, it is now possible to establish the redshift evolution of the star formation rate density ($\rho_{\rm SFR}$). This is achieved by first integrating the LFs weighted by the total IR luminosity ($8-1000\,{\rm\mu m}$); for consistency this was performed again assuming the SMG SED of Michałowski et al. (2010) (which, although peaking at a wavelength consistent with a temperature of 35 K, yields a bolometric luminosity $\simeq 2$ times greater than a simple 35 K modified blackbody template), and the LFs were integrated down to a lower luminosity limit of $0.01\times L^{\star}$. The resulting inferred total IR luminosity density, $\rho_{\rm IR}$ was then converted into dust-obscured star-formation rate density, $\rho_{\rm SFR}$, at each redshift using the scaling factor of $\mathcal{K}_{\rm IR}=4.5\times 10^{-44}\,{\rm M_{\odot}\,year^{-1}\,erg^{-1}\,s}$ from Kennicutt (1998), with an additional multiplicative factor of 0.63 to convert from a Salpeter to a Chabrier IMF.

The results are shown in Fig. 6. The red filled squares depict the values of $\rho_{\rm SFR}$ derived from the four LFs shown in Fig. 3, established using the $1/V_{\rm max}$ method. The red solid line shows the evolution of $\rho_{\rm SFR}$ as determined from the continuous form of the redshift evolution of the LF found using the maximum-likelihood method, with the grey area showing the $1\sigma$ errors. Due to the fact we have chosen to parametrize the Schechter function normalisation parameter, $\Phi^{\star}$ as a Gaussian function , and have chosen a simple linear dependence of redshift for ${\rm log}(L_{\star})$ (see Equation 8), the redshift evolution of $\rho_{\rm SFR}$, depicted as a red solid line, is also a Gaussian:

[TABLE]

The blue squares give the values of UV-visible $\rho_{\rm SFR}$ derived from the Parsa et al. (2016) $\rho_{\rm UV}$ results (integrated down to $M_{\rm UV}=-10$) by converting to UV-visible $\rho_{\rm SFR}$ using the Madau & Dickinson (2014) scaling of $\mathcal{K}_{\rm UV}=1.3\times 10^{-28}\,{\rm M_{\odot}\,year^{-1}\,erg^{-1}\,s\,Hz}$, again also multiplied by a factor of 0.63 to convert to a Chabrier IMF. The blue solid line is a best-fitting Gaussian function to the UV results:

[TABLE]

The black solid line shows the total $\rho_{\rm SFR}$ calculated by simply adding the IR and UV estimates (Equations 9 and 10 respectively). The black dotted and dashed lines depict the alternative functional forms of the redshift evolution of total $\rho_{\rm SFR}$ provided by Madau & Dickinson (2014) and Behroozi et al. (2013) respectively (after conversion to a Chabrier IMF).

6 Discussion

6.1 Comparison of luminosity functions

It is interesting to compare our inferred total IR LFs with those that have been derived from the Herschel surveys. First, in Fig. 7, we compare our LFs with those produced by Magnelli et al. (2013) from the deepest PACS surveys in the GOODS fields. Magnelli et al. (2013) utilised the 70, 100 and 160-$\mu$m imaging produced by the PACS Evolutionary Probe (PEP; Lutz et al. 2011) and GOODS-Herschel (GOODS-H; Elbaz et al. 2011) programmes. They performed blind PACS source extraction, but also created a source catalogue using positional priors, with the positions of Spitzer IRAC 3.6-$\mu$m sources used to extract sources from the Spitzer MIPS 24-$\mu$m imaging, which were then in turn used as positional priors for source extraction in the PACS maps. The latter (MIPS-PACS) source catalogues were then cross-matched with the shorter-wavelength GOODS catalogues (optical+near-infrared) to provide the required photometric redshift information, and the total infrared (8–1000-$\mu$m) luminosities were inferred by fitting the 70, 100 and 160-$\mu$m photometry with the SED template library of Dale & Helou (2002). Magnelli et al. (2013) used this information to construct LFs in six redshift bins using the $1/V_{\rm max}$ method, and fitted the binned LF data with a double power-law function as used previously by Magnelli et al. (2009, 2011). Because these LFs are based on the PACS data, which do not extend to wavelengths longer than 160 $\mu$m, Magnelli et al. (2013) confined their analysis to redshifts $z<2.3$.

Our results are compared with the Magnelli et al. (2013) LFs in Fig.7, using the five redshift bins adopted by Magnelli et al. (2013). The coloured dashed lines show the Magnelli et al. (2013) LFs, while the solid line in each panel shows our LF at the medium redshift of each bin. With the possible exception of the highest redshift bin, the two sets of LFs agree fairly well near the break luminosity. As a result, it transpires that, as shown in the next section, the derived IR luminosities densities are not very different (albeit the Magnelli et al. (2013) results are systematically higher). Nevertheless, it is clear that the LFs diverge at both lower and higher luminosities, with the Magnelli et al. (2013) LFs indicating larger numbers of both faint and bright sources. However, to some extent this difference is exaggerated by the different parameterisations adopted.

Specifically, the difference at the bright end is due, at least in part, to the fact that Magnelli et al. (2013) adopted a double power-law function, while we have fitted a Schechter function (with an exponential cutoff at bright luminosities). In reality, the Magnelli et al. (2013) sample contains only one object with an estimated $L>10^{13}\,{\rm L_{\odot}}$ and only $\simeq 6$ sources with $L>5\times 10^{12}\,{\rm L_{\odot}}$ across the entire redshift range (their figure 8). As a consequence, the bright end of the Magnelli et al. (2013) LFs is fairly unconstrained, and much of the difference seen in Fig. 7 is actually driven by their adoption of a fixed bright-end power-law slope. In fact, as discussed below, integrated over the redshift range $0.5<z<2.3$ the Magnelli et al. (2013) LFs predict $\simeq 20$ times as many bright 850-$\mu$m sources in the wide-area SCUBA-2 sources than are actually observed, and so their LFs are clearly much too high at the bright end.

At the faint end, Magnelli et al. (2013) again adopted a fixed power-law slope, with $\phi\propto L^{-0.6}$. This was based on the local value derived by Sanders et al. (2003), and simply fixed at higher redshifts because the PACS data did not really enable meaningful constraints to be placed on the faint-end slope at $z>0.5$. Thus, in Fig. 7, our LFs diverge from those of Magnelli et al. (2013) towards fainter luminosities because we have adopted a slightly shallower faint-end slope ($\phi~{}\propto~{}L^{-0.4}$), based on the ALMA measurements at $1.5<z<2.5$. As discussed in more detail next, some other Herschel-based estimates of the faint end of the LF in fact use an even flatter faint-end slope ($\phi\propto L^{-0.2}$; Gruppioni et al. 2013) and, as shown in figure 11 of Magnelli et al. (2013), even using the Spitzer MIPS data to extend estimates of the LF to fainter luminosities cannot really distinguish between a faint-end power-law slope of $-0.6$ or $-0.2$. It is therefore both interesting and reassuring that our own measured value of the faint-end slope ($-0.4$) lies midway between the values adopted in the Herschel studies.

Next, in Fig. 8, we compare our LFs with those produced by Gruppioni et al. (2013), which were based on PACS far-IR sources extracted from the Herschel PEP survey data at 70, 100 and 160 ${\rm\mu m}$ within the COSMOS, ECDFS, GOODS-N and GOODS-S fields (covering a total area $\simeq 3.3$ deg2). As well as covering a larger area than Magnelli et al. (2013) (albeit of course to shallower depths), Gruppioni et al. (2013) extended the SED fitting to include the 250, 350 and 500-${\rm\mu m}$ imaging provided by the Herschel SPIRE imaging in the same fields (by the HerMES survey; Oliver et al. 2012), and attempted to extend their LF analysis out to redshifts $z\simeq 4$. Gruppioni et al. (2013) detected 373, 7176 and 7376 sources at 70, 100 and 160 ${\rm\mu m}$, respectively, used cross-matching (via 24-$\mu$m MIPS imaging) with the shorter-wavelength (optical+near-IR) data in the fields to provide redshift information, and used SED fitting with a range of templates to the PACS+SPIRE photometry to derive total IR luminosities. Again, Gruppioni et al. (2013) constructed the LFs using the $1/V_{\rm max}$ method, but in their study they fitted the binned LF data with a modified Schechter function.

Our results are compared with the Gruppioni et al. (2013) LFs in Fig. 8, using the nine redshift bins adopted by Gruppioni et al. (2013). The colour-shaded areas show the Gruppioni et al. (2013) results, with the dashed lines representing their best-fitting modified Schechter functions. The solid black lines show our own estimates of the LF at these redshifts. While at the faint-end the LFs are in fairly good agreement (as mentioned above, Gruppioni et al. (2013) in fact adopt a slightly shallower faint-end slope), at the bright-end there is again significant disagreement. Moreover, at high redshifts the Gruppioni et al. (2013) LFs are significantly higher around the break luminosity, resulting in substantially larger values of IR luminosity density at $z\simeq 3$ (see next section).

One way to quantify how much our new IR LFs differ from those produced from the aforementioned Herschel-based studies is to calculate how many bright sub-mm sources each evolving LF predicts should be present in the SC2LS survey data utilised here. We show the results of this in Fig. 9, and quantify the differences at flux densities $S_{850}>10$ mJy in Table 5. Along with a range of data derived from the SCUBA-2 and ALMA-based literature, the solid black and blue lines in the left-hand panel of Fig. 9 show the predicted cumulative 850-$\mu$m number counts derived from our new LFs and from those produced by Gruppioni et al. (2013) (see Fig. 8) after integration over the redshift range $0.5<z<4.5$. The dashed red and green lines show the corresponding predictions from our own LFs and those produced by Magnelli et al. (2013) (see Fig. 7) after integration over the redshift range $0.5<z<2.3$. Unsurprisingly, since our study is based primarily on 850-$\mu$m data, the number counts predicted by our new IR LFs agree well with the observed counts, whereas the Herschel-based LFs over predict the counts at the bright end by over an order-of-magnitude (as expected, given the dramatic differences at the bright end of the LFs, as shown in Fig. 7 and Fig. 8).

The extent of this difference is quantified further in Table 5, where we tabulate how many 850-${\rm\mu m}$ sources with flux densities $S_{850}>10$ mJy in the combined area of all the S2CLS fields used here (1.5 deg2) a given LF predicts after integration over the appropriate redshift range. The actual detected number of 850-${\rm\mu m}$ sources with $S_{850}>10$ mJy is 17 over the redshift range $0.5<z<4.5$ (consistent with our predicted 20) of which 4 have $0.5<z<2.3$ (consistent with our predicted 8). By contrast the Gruppioni et al. (2013) LFs predict $>300$ such galaxies at $0.5<z<4.5$ (higher by a factor of $\simeq 15-20$), while the Magnelli et al. (2013) LFs predict 140 such sources at $0.5<z<2.3$ (again high by a factor of $15-20$).

As used consistently throughout this study, conversion between total IR luminosity, and observed 850-$\mu$m flux density was performed using the average sub-mm source SED of Michałowski et al. (2010). In the right-hand panel of Fig. 9 we explore what sort of change in the adopted template would be required to bring the bright sub-mm number-count predictions of the Herschel-derived IR LFs into agreement with reality. Here we show the effect of converting the Gruppioni et al. (2013) LFs into predicted 850-$\mu$m number counts for adopted modified black-body SEDs of increasing temperature. As expected, things agree well at the faint end for $T\simeq 45$ K, but a change to characteristic temperatures of $T\simeq 80$ K is required at the bright end to adequately suppress the predicted number counts. There is no evidence to support such a rapid change in typical SED temperature over such a short range in 850-$\mu$m flux density or IR luminosity. Indeed, while several studies have reported a correlation between IR luminosity and the effective temperature of the typical galaxy SED, the derived slope of this correlation is relatively modest, with average dust temperature reported to rise from $T\simeq 30$ K for $L_{\rm IR}\simeq 10^{11}\,{\rm L_{\odot}}$ to $T\simeq 40$ K for $L_{\rm IR}\simeq 10^{13}\,{\rm L_{\odot}}$ (Amblard et al., 2010; Hwang et al., 2010; Smith et al., 2012).

We conclude that the steep drop-off at the bright end of the sub-mm number counts does simply reflect a near exponential decline at the bright end of the IR LF, and the Herschel results have been contaminated and biased high at the bright end by a mix of blending issues, source identification (and hence redshift) errors, and possibly also the adoption of a double power-law LF (see also Bethermin et al. 2017; Liu et al. 2017).

6.2 Comparison of star-formation rate density

In Fig. 10 we compare the redshift evolution of cosmic star-formation density derived here with the evolution inferred from the Herschel-based studies undertaken by Magnelli et al. (2013), Gruppioni et al. (2013) and Rowan-Robinson et al. (2016). All results have been converted to a Chabrier IMF, and the values of obscured $\rho_{\rm SFR}$ derived from the sub-mm/far-IR studies have been added to the UV-derived $\rho_{\rm SFR}$ from Parsa et al. (2016), so that all values plotted here represent total $\rho_{\rm SFR}$.

To derive the Magnelli et al. (2013) results (green points in Fig. 10) we performed the luminosity-weighted integral of their IR LFs (shown in Fig. 7) from a lower limit of $0.01\times L^{\star}$, where $L^{\star}$ is the characteristic luminosity found for our sample (Table 3). We then derived $\rho_{\rm SFR}$ using the conversion factor given by Madau & Dickinson (2014), before applying the IMF correction factor of 0.63 and adding the UV-based results. To derive the corresponding values from Gruppioni et al. (2013), we simply converted their own IR luminosity density estimates using the same procedure. Rowan-Robinson et al. (2016) attempted to extend the exploitation of the Herschel HerMES survey to the highest redshifts via SPIRE-based source selection (over an area of $\simeq 20$ deg2), including sources detected only at 500 $\mu$m. We have adopted their derived estimates of obscured $\rho_{\rm SFR}$, converted to a Chabrier (2003) IMF, and again added the UV-based estimates to produce the (magenta) data points shown in Fig. 10.

Despite the differences in the LFs discussed extensively above, with the exception of the rather high Gruppioni et al. (2013) results around $z\simeq 1$, all three Herschel-based studies yield estimates of $\rho_{\rm SFR}$ consistent with those derived here (as indicated by the blue curve) up to redshifts $z\simeq 2.2$, albeit the uncertainties in the Herschel-derived estimates are generally very large. At higher redshifts the uncertainties in the Herschel-based estimates start to approach 0.5 dex, but even so the smooth decline in $\rho_{\rm SFR}$ beyond $z\simeq 2$ seen in the present study (and consistent with that reported by Behroozi et al. 2013; Madau & Dickinson 2014; Dunlop et al. 2017; Bourne et al. 2017; Liu et al. 2017) is inconsistent with the high values reported by Rowan-Robinson et al. (2016). We suggest that the high values derived by Rowan-Robinson et al. (2016) may reflect problems in source identification and redshift estimation stemming from the large-beam long-wavelength SPIRE data, as well as potential blending issues, and extrapolation of $\rho_{\rm SFR}$ from the most extreme sources. Interestingly, the Rowan-Robinson et al. (2016) results seem somewhat low near the peak of activity at $z\simeq 2$. This may suggest that some subset of numerous far-IR sources which should lie at $z\simeq 2$ have been erroneously placed at higher redshifts, where the resulting boost in inferred luminosity coupled with the smaller cosmological volumes will inevitably result in artificially high estimates of $\rho_{\rm SFR}$.

6.3 Evolution of the LF and the redshift distribution of sub-mm sources

Finally, we discuss how the form of the LF evolution uncovered here naturally explains the apparent ‘down-sizing’ of the sub-mm source population. There is tentative but persistent evidence that the median redshift of sub-mm selected galaxies increases with increasing flux density (Ivison et al., 2002; Pope et al., 2005; Biggs et al., 2011; Vieira et al., 2013; Weiß et al., 2013; Koprowski et al., 2014; Chen et al., 2016; Michałowski et al., 2016), and hence increasing far-IR luminosity. In Fig. 11 we plot our evolving LFs as a function of redshift from $z\simeq 0.5$ to $z\simeq 4.5$. At low redshifts, the increase in both $\Phi^{\star}$ and $L^{\star}$ produces an increase in both the number density of sources, and luminosity density, $\rho_{\rm IR}$, with increasing redshift (cf Fig. 10). At high redshifts the decline in $\Phi^{\star}$ progressively overcomes the continued positive evolution of $L^{\star}$ to produce a decline in both these quantities, but it can be seen that the continued positive evolution of $L^{\star}$ means that the most luminous sources persist, or indeed are preferentially found at the highest redshifts explored here.

To better connect with observables, we have used our evolving LF to calculate the predicted redshift distribution of 850-$\mu$m sources as a function of flux density. The results are shown for four different flux-density thresholds in Fig. 12. Here it can be seen that the peak in the redshift distribution is expected to naturally increase gradually from $z\simeq 1.8$ for $S_{850}>1$ mJy to $z\simeq 3$ for $S_{850}>10$ mJy. This is in excellent accord with what has been reported in the literature (e.g. Koprowski et al. 2014; Michałowski et al. 2016) and clarifies why, although dust-enshrouded star-formation is globally less important than UV-visible star-formation activity at $z>4$, bright sub-mm sources will continue to be discovered out to high redshifts.

7 Conclusion

We have analysed the coverage of the far-IR luminosity–redshift plane provided by (sub-)mm-selected galaxy samples extracted from the SCUBA-2 Cosmology Legacy Survey and the ALMA imaging of the HUDF to make a new measurement of the evolving galaxy far-IR luminosity function (LF) extending out to redshifts $z\simeq 5$. Using both direct ($1/V_{\rm max}$) and maximum-likelihood methods we have determined the form and evolution of the rest-frame 250-$\mu$m galaxy LF. This LF is well described by a Schechter function with a faint-end slope $\alpha\simeq-0.4$ (derived using the ALMA data at $z\simeq 2$) which displays a combination of rising-then-falling density evolution, and positive luminosity evolution.

We converted our 250-$\mu$m results to total IR luminosity using the average sub-mm galaxy SED of Michałowski et al. (2010), and have then compared our determination of the evolving IR LF with those derived by Magnelli et al. (2013) and Gruppioni et al. (2013) from the Herschel PEP and HerMES surveys. Our faint-end slope lies approximately midway between the values adopted in these studies, and the LF normalization near the break luminosity $L^{\star}$ is also comparable at most redshifts. However, both of the Herschel-derived LFs indicate a much larger number of very bright sources at all redshifts than is found in our JCMT/ALMA-based study. To check which result is correct we derive the predicted 850-$\mu$m number counts from the Herschel-based LFs, and find that, at bright flux densities ($S_{850}>10$ mJy) they predict an order-of-magnitude more sources than are observed. We explore whether this discrepancy can be explained/resolved by adoption of different SED templates, and conclude that it cannot (without extreme temperatures, and more importantly without a contrived and much more rapid dependence of $T$ on IR luminosity than has been reported in the literature).

We have utilised our measurement of the evolving IR LF to derive comoving IR luminosity density, and hence obscured star-formation rate density, which we then combine with UV-estimates of unobscured activity (from Parsa et al. 2016), to derive the evolution of total $\rho_{\rm SFR}$. Again we compare with values derived from the Herschel-based studies, and find reasonable agreement out to $z\simeq 2$, but increasing disagreement at higher redshift. Specifically, consistent with several other recent studies (e.g. Bourne et al. 2017; Dunlop et al. 2017; Liu et al. 2017) we find that $\rho_{\rm SFR}$ declines beyond $z\simeq 2-2.5$ and is dominated by UV-visible star-formation activity beyond $z\simeq 4$. In contrast, Gruppioni et al. (2013) and Rowan-Robinson et al. (2016) report essentially no decline in dust-obscured $\rho_{\rm SFR}$ from $z\simeq 3$ to $z\simeq 6$. Given the severe over-prediction of the 850-$\mu$m counts produced by the Herschel IR LFs, we conclude that the high values reported from these studies most likely reflect problems in source identification and redshift estimation arising from the large-beam long-wavelength SPIRE data, as well as potential blending issues, and extrapolation of $\rho_{\rm SFR}$ from the most extreme sources.

Finally, we show how the evolution of the IR LF as derived here (with its combination of rising-then-falling characteristic density ($\Phi^{\star}$), and positive evolution of characteristic luminosity density ($L^{\star}$)) with redshift, produces a decline in inferred $\rho_{\rm SFR}$ beyond $z\simeq 2-2.5$ while at the same time predicting that the most luminous sub-mm sources will continue to be found out to very high redshifts ($z\simeq 5-6$). Specifically, our evolving LF, with its combined luminosity+density evolution, predicts that the median redshift of sub-mm sources should increase with increasing flux density, consistent with several reports in the recent literature. Such evolution is consistent with many studies of AGN evolution, suggesting that both dust-enshrouded star formation and AGN activity are strongly linked to the growth of stellar mass in galaxies.

Acknowledgements

KEKC and MPK acknowledge the support of the UK Science and Technology Facilities Council through grant number ST/M001008/1. MPK acknowledges the support of the Carnegie Trust Research Incentive Grant (PI: M. Michałowski). JEG is supported by the Royal Society. JSD, MJM and RJM acknowledge the support of the UK Science and Technology Facilities Council through grant number ST/M001229/1. MJM acknowledges the support of the National Science Centre, Poland through the POLONEZ grant 2015/19/P/ST9/04010. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 665778.

The James Clerk Maxwell Telescope is now operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observa- tories of China and the Chinese Academy of Sciences (grant no. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and par- ticipating universities in the United Kingdom and Canada. The data utilised in this paper were taken as part of Program ID MJLSC02.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00173.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), to- gether with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.