The evolving far-IR galaxy luminosity function and dust-obscured star-formation rate density out to z~5
M.P. Koprowski, J.S. Dunlop, M.J. Micha?owski, K.E.K. Coppin, J.E. Geach, R.J. McLure, D. Scott, P.P. van der Werf
MMNRAS , 1–17 (2017) Preprint 5 June 2017 Compiled using MNRAS L A TEX style file v3.0
The evolving far-IR galaxy luminosity function anddust-obscured star-formation rate density out to z (cid:39) M. P. Koprowski, , (cid:63) J. S. Dunlop, M. J. Micha(cid:32)lowski, K. E. K. Coppin, J. E. Geach, R. J. McLure, D. Scott, and P. P. van der Werf Centre for Astrophysics Research, Science & Technology Research Institute, University of Hertfordshire, Hatfield AL10 9AB, UK SUPA † , Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh, EH9 3HJ, UK Department of Physics and Astronomy, 6224 Agricultural Road, University of British Columbia, Vancouver V6T 1Z1, Canada Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present a new measurement of the evolving galaxy far-IR luminosity function(LF) extending out to redshifts z (cid:39) , with resulting implications for the level ofdust-obscured star-formation density in the young Universe. To achieve this we haveexploited recent advances in sub-mm/mm imaging with SCUBA-2 on the James ClerkMaxwell Telescope (JCMT) and the Atacama Large Millimeter/Submillimeter Array(ALMA), which together provide unconfused imaging with sufficient dynamic range toprovide meaningful coverage of the luminosity-redshift plane out to z > . Our resultssupport previous indications that the faint-end slope of the far-IR LF is sufficientlyflat that comoving luminosity-density is dominated by bright objects ( (cid:39) L ∗ ). However,we find that the number-density/luminosity of such sources at high redshifts has beenseverely over-estimated by studies that have attempted to push the highly-confused Herschel
SPIRE surveys beyond z (cid:39) . Consequently we confirm recent reports thatcosmic star-formation density is dominated by UV-visible star formation at z > . Us-ing both direct ( / V max ) and maximum likelihood determinations of the LF, we findthat its high-redshift evolution is well characterized by continued positive luminosityevolution coupled with negative density evolution (with increasing redshift). This ex-plains why bright sub-mm sources continue to be found at z > , even though theirintegrated contribution to cosmic star-formation density at such early times is verysmall. The evolution of the far-IR galaxy LF thus appears similar in form to thatalready established for active galactic nuclei, possibly reflecting a similar dependenceon the growth of galaxy mass. Key words: dust, extinction – galaxies: evolution, high-redshift, luminosity function,star formation – cosmology: observations
A key challenge in modern astrophysical cosmology is tocomplete our knowledge of cosmic star-formation history,taking proper account of dust-obscured activity (e.g. Madau& Dickinson 2014). Achieving this requires a reliable mea-surement of the form and evolution of both the rest-frame UV and rest-frame far-IR galaxy luminosity func-tions (LFs) out to the highest redshifts. This is because afair census of both unobscured and dust-obscured comovingstar-formation rate density ( ρ SFR ) requires the luminosity-weighted integration of the relevant LFs over sufficient dy- (cid:63)
E-mail: [email protected] † Scottish Universities Physics Alliance namic range to properly account for the contributions of thebrightest sources, while at the same time quantifying theimpact of the adopted lower luminosity limit for the inte-gration (i.e. reliably establishing the faint-end slope of theLF).In recent years this goal has been largely achieved atrest-frame UV wavelengths over nearly all of cosmic his-tory, 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, meaning-ful disagreements over the comoving UV luminosity density( ρ UV ) produced by the evolving galaxy population are now c (cid:13) a r X i v : . [ a s t r o - ph . GA ] J un M. Koprowski et al. largely confined to extreme redshifts z > (Ellis et al. 2013;Oesch et al. 2014; McLeod et al. 2015), and even here muchof the claimed disagreement can be removed by the adop-tion of consistent limits to the LF integration (McLeod et al.2016; Ishigaki et al. 2017). The reason that the adopted faintintegration limit becomes an issue for calculating ρ UV at ex-treme redshifts (and the resulting contribution of the emerg-ing early galaxy population to reionization; Robertson et al.2013, 2015) is that the faint-end slope ( α ) of the galaxy UVLF steepens with increasing redshift, from relatively modestvalues at intermediate redshifts (e.g. α (cid:39) − . at z (cid:39) ; Parsaet al. 2016; Mehta et al. 2017) to α (cid:39) − at z > (McLureet al. 2013; Bouwens et al. 2015; McLeod et al. 2015, 2016).As a result, despite the discovery of significant numbers ofUV-bright galaxies at z > (Bowler et al. 2014, 2015), ρ UV within the first (cid:39) Gyr of cosmic time is dominated by thecontributions of the numerous faintest galaxies, and so thederived value depends more critically on how far down inluminosity the LF extends than is the case at later times(see Parsa et al. 2016).Despite the long-established importance of the far-IRbackground (Dole et al. 2006), similar progress in our knowl-edge of the far-IR LF has been hampered by the heightenedobservational challenges at mid/far-IR and sub-mm/mmwavelengths (i.e. high background and poor angular reso-lution) and a resulting lack of survey dynamic range (cou-pled with uncertainties in redshift content). Nevertheless,important progress has been made over the past decade, firstwith mid-IR observations using NASA’s
Spitzer Space Tele-scope , and more recently through far-IR imaging with ESA’s
Herschel Space Observatory . Together, these facilities haveenabled the far-IR LF and its evolution to be successfullytraced out to z (cid:39) .First, Spitzer
MIPS 24 µ m imaging was used to studythe mid-IR LF out to z (cid:39) (Le Floc’h et al. 2005; Caputiet al. 2007; Rodighiero et al. 2010), albeit extrapolation from24 µ m to far-IR luminosity becomes increasingly dangerouswith increasing wavelength (although see Elbaz et al. 2010).Attempts were also made to exploit the 70 µ m imaging pro-vided by Spitzer (Magnelli et al. 2009; Patel et al. 2013), butsensitivity/resolution limitations largely restricted the use-fulness of this work to z < (although Magnelli et al. 2011pushed out to z (cid:39) via stacking).Over the last five years, Herschel
PACS and SPIRE sur-veys have enabled this work to be developed through objectselection at more appropriate far-IR wavelengths (i.e. closerto the rest-frame peak of the far-IR emission). Magnelli et al.(2013) used deep
Herschel
PACS imaging from the PEP andGOODS surveys to determine the bright-end of the far-IRLF out to z (cid:39) , while Gruppioni et al. (2013) used the PEPPACS imaging and SPIRE HerMES imaging (at 250, 350and 500 µ m ) to try to extend this work out to z (cid:39) (see alsoBurgarella et al. 2013). The inclusion of the 250, 350 and500 µ m data allowed Gruppioni et al. (2013), in principle,to determine far-IR luminosities without recourse to largeextrapolations, but object selection was still undertaken atthe shorter (PACS) wavelengths. As a result, this study ismostly sensitive to the lower-redshift/warmer sources, andhence, at z > , allows only the detection of the most extremesources. One consequence of the resulting lack of dynamicrange is that the faint-end slope of the far-IR LF could notbe measured at high redshift, and so Gruppioni et al. (2013) simply adopted the z = value at all higher redshifts. Mostrecently, Rowan-Robinson et al. (2016) attempted to expandon the Gruppioni et al. (2013) study, and to extend it outto even higher redshifts ( z (cid:39) ) by including object selectionat 500 µ m (the longest-wavelength Herschel
SPIRE imag-ing band). This study yielded surprisingly high estimates offar-IR luminosity density at high redshifts; however, util-ising the longest-wavelength
Herschel data in this way isfraught with danger due to the large beamsize ( (cid:39) arcsecFWHM) and consequent issues regarding blending, sourcemis-identification, potential AGN contamination and grav-itational lensing (e.g. star-formation rates, SFR, as high as20,000 M (cid:12) yr − are reported in this work).Thus, despite this impressive progress, it is clear thatattempting to reliably measure the far-IR galaxy LF basedon object selection in Spitzer or Herschel surveys becomesincreasingly problematic beyond z (cid:39) . It is now possi-ble to overcome these difficulties at high redshift by mov-ing to ground-based sub-mm/mm object selection, using acombination of wide-area imaging surveys as produced bySCUBA-2 (Holland et al. 2013) on the JCMT (e.g. Geachet al. 2013, 2017; Roseboom et al. 2013; Chen et al. 2016;Micha(cid:32)lowski et al. 2016; Cowie et al. 2017), and higher-resolution smaller-area mapping as can now be achieved withALMA (e.g. Hatsukade et al. 2015, 2016; Umehata et al.2015; Walter et al. 2016; Dunlop et al. 2017). While attemptshave previously been made to explore the basic high-redshiftevolution of the far-IR population based on ground-basedsub-mm data (e.g. Wall, Pope & Scott 2008), until now the‘wedding-cake’ of surveys providing unconfused imaging overa reasonable dynamic range was not of sufficient quality toenable a detailed investigation of the form and evolution ofthe far-IR LF.The merits of moving to ground-based sub-mm (i.e.450 and 850 µ m ) and mm (i.e. 1.1–1.3 mm) selection be-come increasingly obvious as one moves to higher redshifts.First, the smaller beam sizes at these wavelengths offeredby large ground-based single-dish telescopes such as theJCMT, or interferometric arrays such as ALMA, substan-tially reduce/remove the confusion limitations of the longer-wavelength Herschel
SPIRE imaging, minimizing problemsof source blending, false counterpart identification, andhence potentially erroneous redshift information. Second,certainly by z > , object selection at sub-mm/mm wave-lengths is much less susceptible to AGN contamination than Spitzer or Herschel
PACS detections (which sample the spec-tral energy distribution (SED) shortward of λ rest (cid:39) µ m atthese redshifts), and extrapolation to estimate total far-IRluminosities also becomes less problematic. Third, object se-lection in the higher-resolution ground-based sub-mm/mmimaging can be used to help deconfuse and deblend the exist-ing Herschel far-IR imaging, hence enabling its exploitationfor SED determination in a less biased way.A number of recent studies have already provided indi-cations that the high values of ρ FIR inferred from pushingthe
Herschel surveys beyond z (cid:39) . are incorrect. First, es-timates of dust-obscured star-formation activity in knownhigh-redshift galaxies, either derived from analyses of theUV slope (e.g. Dunlop et al. 2013) or from ALMA detec-tions/limits (individual or stacked; e.g. Capak et al. 2015;Bouwens et al. 2016b; Koprowski et al. 2016b) suggest thatthe dust-obscured contribution to ρ SFR at the highest red-
MNRAS , 1–17 (2017) alaxy far-IR LF out to z (cid:39) shifts is relatively small. However, it can reasonably be ar-gued that the most dust-obscured galaxies will not featurein rest-frame UV-selected samples. More significantly, theresults from the deep ALMA imaging of the Hubble UltraDeep Field (HUDF; Dunlop et al. 2017) and the results ofstacking in the deepest SCUBA-2 Cosmology Legacy Sur-vey (S2CLS) images (Bourne et al. 2017) both yield muchlower values of dust-obscured ρ SFR than derived by Grup-pioni et al. (2013) and Rowan-Robinson et al. (2016), es-pecially at z > . Nonetheless, it could still be argued thatthese small-field experiments, while complete to relativelylow dust-obscured SFRs ( (cid:39) (cid:12) yr − ), cover insufficientarea to reveal the contributions from the most extreme ob-jects.In this study we aimed the clarify this situation, andmoreover to properly determine the form and high-redshiftevolution of the far-IR galaxy LF by analysing the resultsfrom the HUDF ALMA 1.3-mm and deep S2CLS 850- µ m imaging in tandem with results from the wider-area − µ m maps recently completed within the S2CLS (Geach et al.2017; Micha(cid:32)lowski et al. 2016). Together these data provide asub-mm/mm survey ‘wedding-cake’ with sufficient dynamicrange and sufficiently-complete redshift content to enable ameaningful measurement of the form and evolution of therest-frame 250- µ m galaxy LF from the available coverage ofthe luminosity-redshift plane. Crucially, the bright tier ofS2CLS covers over 2 deg , and yields a sample of > luminous sources (SFR > (cid:12) yr − ), providing excellentsampling of the bright end of the far-IR LF, while the deeperS2CLS imaging and the ALMA imaging enable the first di-rect measurement of the slope of the faint-end of the far-IRLF at high redshift. Together these data have enabled usto determine the form and evolution of the far-IR LF, andhence the dust-obscured ρ SFR out to z (cid:39) . .This paper is structured as follows. The sub-mm/mmimaging utilised in this work, along with the supportingmulti-wavelength data, are described in the Section 2. Themulti-wavelength identification process, together with themethods used to establish redshifts for the entire sample arethen presented in Section 3. Next, the procedure for deter-mining the IR LFs, using both the / V max and the maximum-likelihood methods, is explained in Section 4. The resultingcalculation of the far-IR and total (UV+far-IR) ρ SFR is pre-sented in Section 5. We discuss our results in the contextof other recent studies in Section 6, and conclude by sum-marising our findings in Section 7.Throughout the paper we use a Chabrier (2003) stellarinitial mass function (IMF) and assume a flat cosmologywith Ω m = . , Ω Λ = . and H = 70 km s − Mpc − . We used the data collected as a part of the SCUBA-2 Cos-mology Legacy Survey (S2CLS). The map-making processand the resulting derived source catalogues are described inGeach et al. (2017). The fields utilised here are the UKIDSS-UDS, where the 850- µ m imaging covers (cid:39) . deg with a1- σ noise of 0.9 mJy (revealing 1085 sources with a signal-to-noise ratio SNR > . ), and the COSMOS field, where Table 1.
The refined SCUBA-2 850- µ m survey fields and sourcesamples used in this work (see text for details). The columns showthe field name, the area, limiting flux density, limiting signal-to-noise ratio and the resulting number of detections. Area S lim SNR N/deg − /mJyCOSMOS deep .
08 4 . .
79 6 . .
71 3 . µ m imaging covers (cid:39) . deg with the σ noiseof 1.6 mJy (revealing 719 sources with SNR > . ). Thesesource catalogues are discussed further in Chen et al. (2016)and Micha(cid:32)lowski et al. (2016). For simplicity we reduced theeffective area of these maps to regions of uniform noise, andfor our analysis retained only sources with a simulated com-pleteness > . . The resulting refined effective survey areas,flux-density limits, SNR, and sample sizes are summarizedin Table 1. To help inform the measurement of the faint-end slope ofthe LF, we used the ALMA 1.3-mm imaging of the HUDFundertaken by Dunlop et al. (2017). A mosaic of 45 ALMApointings was created to cover the full (cid:39) . arcmin areapreviously imaged with WFC3/IR on HST . The ALMA mapreached a noise level of σ . (cid:39) µ Jybeam − , and 16 sourceswere detected with flux densities S . > µ Jy . 13 of the16 sources have spectroscopic redshifts and the remainingthree have accurate photometric redshifts derived from theoptical–near-IR photometry (for the sources of the spectro-scopic redshifts see table 2 in Dunlop et al. 2017). For the S2CLS COSMOS field the optical to mid-IR dataconsist of imaging from the Canada-France-Hawaii Tele-scope Legacy Survey (CFHTLS; Gwyn 2012), as describedby Bowler et al. (2012), Subaru (Taniguchi et al. 2007),the
HST
Cosmic Assembly Near-infrared Deep Extragalac-tic Legacy Survey (CANDELS; Grogin et al. 2011), UltraV-ISTA Data Release 2 (McCracken et al. 2012; Bowler et al.2014) and
Spitzer (S-COSMOS; Sanders et al. 2007). At ra-dio wavelengths the VLA-COSMOS Deep (Schinnerer et al.2010) catalogues were utilised.For the S2CLS UDS field the optical data were ob-tained with Subaru/SuprimeCam (Miyazaki et al. 2002),as described in Furusawa et al. (2008), the near-IR imag-ing was provided by the UKIRT Infrared Deep Sky Survey(UKIDSS; Lawrence et al. 2007; Cirasuolo et al. 2010), themid-IR data are from the
Spitzer
Public Legacy Survey ofthe UKIDSS Ultra Deep Survey (SpUDS; PI: J. Dunlop), asdescribed in Caputi et al. (2011), and the radio (VLA) dataare from Ivison et al. (2005, 2007) and Arumugam et al. (inpreparation).For both the UDS and COSMOS fields, far-IR imagingfrom
Herschel (Pilbratt et al. 2010) was utilised, as provided
MNRAS000
MNRAS000 , 1–17 (2017)
M. Koprowski et al. by the public releases of the HerMES (Oliver et al. 2012)and PEP (Lutz et al. 2011) surveys undertaken with theSPIRE (Griffin et al. 2010) and PACS (Poglitsch et al. 2010)instruments.For the extraction of the far-IR flux densities and lim-its, the
Herschel maps at 100, 160, 250, 350 and 500 µ m were utilised, with beam sizes of 7.4, 11.3, 18.2, 24.9, and36.3 arcsec, and 5 σ sensitivities of 7.7, 14.7, 8.0, 6.6, and9.5 mJy, respectively. The Herschel flux densities of (or up-per limits for) the SCUBA-2 sources were obtained in the fol-lowing way. Square image cut-outs of width 120 arcsec wereextracted from each
Herschel map around each SCUBA-2source, and the PACS (100 and 160 µ m ) maps were used tosimultaneously fit Gaussians with the FWHM of the respec-tive imaging, centred at all radio and 24- µ m sources locatedwithin these cut-outs, and at the positions of the SCUBA-2optical identifications (IDs, or just submm positions if no IDswere selected). Then, to deconfuse the SPIRE (250, 350 and500 µ m ) maps in a similar way, the positions of the 24 µ m sources detected with PACS (at > σ ) were retained, alongwith the positions of all radio sources, and the SCUBA-2IDs (or, again, the SCUBA-2 position in the absence of anyoptical or radio ID). Because of the beam size of the JCMT SCUBA-2 imagingat 850 µ m (FWHM (cid:39) arcsec), we cannot simply adoptthe closest optical/near-IR neighbour as the 850- µ m sourcegalaxy counterpart. Instead, we used the method outlined inDunlop et al. (1989) and Ivison et al. (2007), where we adopta . σ search radius around the SCUBA-2 position based onthe signal-to-noise ratio (SNR): r s = . × . × FWHM / SNR .In order to account for systematic astrometry shifts (causedby pointing inaccuracies and/or source blending; e.g. Dun-lop et al. 2010) we enforced a minimum search radius of4.5 arcsec. Within this radius we calculated the correctedPoisson probability, p , that a given counterpart could havebeen selected by chance.Three imaging wavebands were used when searching forgalaxy counterparts: the VLA 1.4-GHz imaging, the Spitzer
MIPS 24- µ m imaging, and the Spitzer
IRAC 8- µ m imaging.The radio band traces recent star formation via synchrotronradiation from relativistic electrons produced within super-novae (SNe; Condon 1992), whereas the 24- µ m waveband issensitive to the emission from warm dust. Therefore, sincesubmm-selected galaxies are dusty, highly-star-forming ob-jects, they are expected to be very luminous in these bands.Also, at the redshifts of interest, the 8- µ m waveband tracesthe rest-frame near-IR light coming from the older, mass-dominant stellar populations in galaxies, and thus providesa proxy for stellar mass. Given the growing evidence thatsub-mm galaxies are massive, it is expected that 850- µ m sources will have significant 8- µ m fluxes (e.g. Dye et al. 2008;Micha(cid:32)lowski et al. 2010; Biggs et al. 2011; Wardlow et al.2011). Moreover, the surface density of sources in these threewavebands is low enough for chance positional coincidencesto be rare (given a sufficiently small search radius). Oncethe counterparts were found in each of these bands, they . . . . . . . . . . . . . .
60 1 2 3 40 . . . . . . . f r a c t i o np e r b i n COSMOS deepCOSMOS wideUDS
Figure 1.
The redshift distributions of the refined SCUBA-2source samples used in this work. The black histogram depictsthe distribution for all the sources, and yields a mean redshiftof ¯ z = . ± . . From the top, the colour plots show the COS-MOS deep, COSMOS wide and the UDS redshift distributionswith ¯ z = . ± . , ¯ z = . ± . and ¯ z = . ± . respectively. were matched with the optical/near-IR catalogues using asearch radius of r = . arcsec and the closest object takento be the galaxy counterpart (for a more detailed descrip-tion of this process and the relevant identification tables seeMicha(cid:32)lowski et al. 2016). MNRAS , 1–17 (2017) alaxy far-IR LF out to z (cid:39) We used the available multi-wavelength data to derivethe optical/near-IR photometric redshifts for the SCUBA-2source galaxy counterparts using a χ minimisation method(Cirasuolo et al. 2007, 2010) with a code based on theH YPER
Z package (Bolzonella et al. 2000). To create tem-plates of galaxies, the stellar population synthesis models ofBruzual & Charlot (2003) were applied, using the Chabrier(2003) stellar initial mass function (IMF), with a lower andupper mass cut-off of . and (cid:12) respectively. Singleand double-burst star-formation histories with a fixed solarmetallicity were used. Dust reddening was taken into ac-count using the Calzetti et al. (2000) law within the range ≤ A V ≤ . The HI absorption along the line of sight wasapplied according to Madau (1995). The accuracy of thephotometric catalogue of Cirasuolo et al. (2010) is excellent,with a mean | z phot − z spec | / ( + z spec ) = . ± . . Also, forevery source the ‘long-wavelength’ photometric redshift (e.g.Aretxaga et al. 2007; Koprowski et al. 2014, 2016a) was cal-culated using the SCUBA-2 and Herschel data by fitting theaverage sub-mm galaxy SED template of Micha(cid:32)lowski et al.(2010).In addition, as explained in detail in Koprowski et al.(2014, 2016a), we further tested the robustness of our op-tical identifications by comparing the optical photometricredshifts with the ‘long-wavelength’ photometric redshifts.Sources with optical redshifts that transpired to be signifi-cantly lower than their ‘long-wavelength’ ones were consid-ered to be incorrectly identified in the optical (foregroundgalaxies, possibly lenses) and hence their ‘long-wavelength’photometric redshifts were adopted (see Micha(cid:32)lowski et al.2016 for details).The final redshift distributions (100% complete) for therefined fields used in this work (Table 1) are shown in Fig-ure 1. Each coloured histogram shows the redshift distribu-tion for the relevant field (as depicted in the legend), withthe background black histogram showing the results for allthe fields combined. The mean redshifts are ¯ z = . ± . , ¯ z = . ± . and ¯ z = . ± . for the COSMOS deep,COSMOS wide and UDS fields respectively. The mean red-shift for the whole sample used here is ¯ z = . ± . . To derive the evolving far-IR LF, Φ IR , we use two indepen-dent methods. One is the standard / V max method (Schmidt1968), which allows the calculation of the LFs to be per-formed directly from the data, without any assumptions re-garding the functional shape. To derive the functional formwe then fit a set of Schechter functions, where the best-fittingparameter values are found by minimising χ . In order tofind the continuous form of the redshift evolution of the far-IR LF, we additionally use the maximum-likelihood methodpresented in Marshall et al. (1983), where we again take theLFs to be of a Schechter form. . . . . . . . COSMOS deepCOSMOS wide UDS redshift l og ( L / W H z − ) Figure 2.
The coverage of the luminosity-redshift plane providedby the source samples used in this work. The grey solid linesdepict the redshift and luminosity bins used in the LF analysis.The solid colour lines show the corresponding luminosity limits,resulting from the detection limits in each field (see Section 2.1and Table 1). These limits are crucial for determining what isthe lowest luminosity in each redshift bin at which our sampleis complete. Only the luminosity bins for which the minimumluminosity is higher than the luminosity limit are included in theanalysis. / V max / V max / V max method The LF in a given luminosity and redshift bin is calculatedusing: Φ ( L , z ) = ∆ L ∑ i − FDR w i × V max , i , (1)where ∆ L is the width of the luminosity bin, FDR is the falsedetection rate, w i is the completeness for the i -th galaxy and V max , i is the co-moving volume available to the i -th source.The errors on the LFs were calculated using the Poissonianapproach. The available co-moving volume is the volumebetween z min and z max , where z min is just the lower boundaryof a given redshift bin and z max is the maximum redshift atwhich the i -th source would still be visible in a given S2CLSmap, or simply the upper boundary of a given redshift bin,whichever is lower. Therefore: V max , i = ∑ j Ω j π V max , j , (2)where we sum over all the available S2CLS fields and Ω j isthe solid angle subtended by the j -th field on the sky. Thecontribution from the j -th field to the V max , i , in a given red-shift bin, V max , j , is only counted if the maximum redshift atwhich the i -th source would still be visible in that field ishigher than the lower boundary of that redshift bin, other-wise V max , j for that field is simply equal to 0.Since in practice we are working with 850 µ m -detectedsources with a mean redshift of ∼ . at 850 µ m , we decidedto initially calculate the LF at a rest-frame wavelength of250 µ m . In order to do so, the luminosity at µ m / ( + z ) iscalculated and interpolated to λ rest = µ m using the aver-age SMG template from Micha(cid:32)lowski et al. (2010). The range MNRAS , 1–17 (2017)
M. Koprowski et al. − − − − − − − . < z < . . < z < . − − − − − − − . < z < . . < z < . L / W Hz − Φ / M p c − d e x − Figure 3.
The far-IR (rest-frame 250- µ m ) galaxy luminosity functions (LFs) for the four redshift bins studied in this work. The pointswith error bars show the LF values determined using the / V max method. The two faintest points in the . < z < . redshift bin depict theLF values found using the ALMA data. These allowed us to determine the faint-end slope α = − . , which was then adopted for the otherredshift bins. The coloured solid lines show the best-fitting Schechter functions to these data points. The black solid lines (almost perfectlyaligned with the coloured solid lines) depict the results of the maximum-likelihood method, with the derived uncertainty indicated bythe shaded grey region. − − − . < z < . . < z < . − − − . < z < . . < z < . L ? / W Hz − Φ ? / M p c − d e x − Figure 4.
The best-fitting parameter values of the Schechter functions found using the / V max method for the four redshift bins studiedin this work. The contours show ∆ χ = and ∆ χ = above the best-fitting solution, and hence the corresponding σ and σ uncertaintiesin the individual parameter values. MNRAS , 1–17 (2017) alaxy far-IR LF out to z (cid:39) Table 2.
Rest-frame 250- µ m luminosity functions.log( L / WHz − ) log( Φ / Mpc − dex − ) . < z < . . < z < . . < z < . . < z < . . ... − . ± . ... ... . ... − . ± . ... ... . − . ± . − . ± . ... ... . − . ± . − . ± . − . ± . − . ± . . − . ± . − . ± . − . ± . − . ± . . − . ± . − . ± . − . ± . − . ± . . − . ± . − . ± . − . ± . − . ± . . ... − . ± . − . ± . − . ± . . ... − . ± . − . ± . − . ± . . ... ... − . ± . ... of luminosities and redshifts for which the LFs were calcu-lated were determined from the coverage of the luminosity-redshift plane (grey rectangles in Fig. 2). The luminositybins used for each field at a given redshift depend on theluminosity lower limit (corresponding to the flux-density de-tection limit), below which no sources can be detected (solidcoloured lines in Fig. 2). Only the luminosity bins with com-plete luminosity coverage are included in the analysis.The results, with Poissonian errors, are listed in Table2, and plotted in Fig. 3. The coloured solid lines are thebest-fitting Schechter functions: Φ Sch ( L , z ) = Φ (cid:63) (cid:18) LL (cid:63) (cid:19) α exp (cid:18) − LL (cid:63) (cid:19) , (3)with Φ (cid:63) being the normalisation parameter, α the faint-end slope and L (cid:63) the characteristic luminosity that roughlymarks the border between the power-law fit, ( L / L (cid:63) ) α , andthe exponential fit. The black solid lines (almost perfectlyaligned 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 findthe faint-end slope, α , we utilised the ALMA data (Sec-tion 2.2). Since the majority of ALMA sources lie at red-shifts . < z < . , it was decided to determine α only forthat redshift bin and then use that value in the remain-ing bins (two faintest points in . < z < . redshift bin ofFig. 3). Fitting the Schechter function to the . < z < . data yielded a faint-end slope of α = − . . The remainingSchechter-function parameters were determined by minimis-ing χ . In Fig. 4 the 68% and 95% confidence intervals areshown, corresponding to ∆ χ = and 4 respectively; the er-rors on the best-fiting Schechter-function parameters canbe determined by projecting the contours onto the relevantaxis. The best-fitting values for the Schechter-function pa-rameters are given in Table 3, and plotted in Fig. 5 as bluecircles. The likelihood function used here (Marshall et al. 1983) isdefined as a product of the probabilities of observing exactlyone 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 ob-serving zero sources in all the other differential elements in Table 3.
The best-fitting parameter values for the Schechter func-tions as determined by the / V max method. The faint-end slope, α , is fixed here at the z (cid:39) value (see text for details). z α log( Φ (cid:63) / Mpc − dex − ) log( L (cid:63) / WHz − ) . < z < . − . − . + . − . . + . − . . < z < . − . − . + . − . . + . − . . < z < . − . − . + . − . . + . − . . < z < . − . − . + . − . . + . − . the luminosity-redshift plane. Using Poisson probabilities,the likelihood is: L = N ∏ i λ i e − λ i ∏ j e − λ j , (4)where j runs over all differential elements in which no sourceswere observed, and λ is the expected number of galaxies in dzdL at z , L : λ = Φ ( z , L ) Ω ( z , L ) dVdz dzdL , (5)with Ω being the fractional area of the sky in which a galaxywith a given z and L can be detected in our fields. Since forlarge N the test statistic − ( L ) will be χ distributed wedefine: S = − ( L ) = − N ∑ i ln [ Φ ( z i , L i )]+ (cid:90)(cid:90) Φ ( z , L ) Ω ( z , L ) dVdz dzdL , (6)where we dropped terms independent of the model param-eters. This step transforms both exponents in Equation 4into the integral, which represents the model − the expectednumber of sources within the integral limits. This meansthat, as compared to the / V max method, this technique cananalyse the entire luminosity-redshift plane (including itsempty patches) and therefore give a result with higher sta-tistical significance. In addition, it does not bin the data,and also allows a user to define the parameters of a luminos- MNRAS000
The best-fitting parameter values for the Schechter func-tions as determined by the / V max method. The faint-end slope, α , is fixed here at the z (cid:39) value (see text for details). z α log( Φ (cid:63) / Mpc − dex − ) log( L (cid:63) / WHz − ) . < z < . − . − . + . − . . + . − . . < z < . − . − . + . − . . + . − . . < z < . − . − . + . − . . + . − . . < z < . − . − . + . − . . + . − . the luminosity-redshift plane. Using Poisson probabilities,the likelihood is: L = N ∏ i λ i e − λ i ∏ j e − λ j , (4)where j runs over all differential elements in which no sourceswere observed, and λ is the expected number of galaxies in dzdL at z , L : λ = Φ ( z , L ) Ω ( z , L ) dVdz dzdL , (5)with Ω being the fractional area of the sky in which a galaxywith a given z and L can be detected in our fields. Since forlarge N the test statistic − ( L ) will be χ distributed wedefine: S = − ( L ) = − N ∑ i ln [ Φ ( z i , L i )]+ (cid:90)(cid:90) Φ ( z , L ) Ω ( z , L ) dVdz dzdL , (6)where we dropped terms independent of the model param-eters. This step transforms both exponents in Equation 4into the integral, which represents the model − the expectednumber of sources within the integral limits. This meansthat, as compared to the / V max method, this technique cananalyse the entire luminosity-redshift plane (including itsempty patches) and therefore give a result with higher sta-tistical significance. In addition, it does not bin the data,and also allows a user to define the parameters of a luminos- MNRAS000 , 1–17 (2017)
M. Koprowski et al. − − − . . . . . . . . . redshift L ? / W H z − Φ ? / M p c − d e x − Figure 5.
Best-fitting Schechter-function parameters. The bluecircles depict the best-fitting values determined using the / V max method (Table 3), with the σ errors derived by projecting thecontours from Fig. 4 onto the relevant axis. The blue solid lineshows the results from the maximum-likelihood method using thefunctional form given in Equation 7 with the parameter valuesgiven in Table 4. The grey-shaded area depicts the σ errors onthe functional form based on the maximum-likelihood ratios, cor-responding to ∆ χ = . ity function, Φ ( z , L ) , to be continuous functions of redshift,which can be then determined by minimising S .For consistency, we chose Φ ( z , L ) to have the form ofa Schechter function (Equation 3). As before, the faint-endslope, α , was fixed at the . < z < . value of − . . As wesought an acceptable description of the data, we exploredvarious functional forms for the redshift dependence of Φ (cid:63) and log( L (cid:63) ), and in the end found that we could adopt asimple normalised Gaussian function and a linear functionof redshift respectively: Φ (cid:63) ( z | A , σ , µ ) = A σ √ π e − ( z − µ ) σ log ( L (cid:63) ( z | a , b )) = az + b . (7)The redshift limits of the integral in Equation 6 are thesame as in the / V max method, z = . − . . The luminosity The normalised Gaussian was chosen here instead of the linearfunction to ensure that Φ (cid:63) does not become negative at high- z . Table 4.
Best-fit values for the continuous functions of theSchechter function parameters from Equation 7, with σ errorsbased on the maximum-likelihood ratios and corresponding to ∆ χ = . Parameter Value A . + . − . × − σ . + . − . µ . + . − . a . + . − . b . + . − . limits were chosen in order to cover the whole useful luminos-ity range (including the ALMA data), log ( L / WHz − ) = . − . . Minimising S gave the best-fitting parameter val-ues, as summarized in Table 4, which were then used to de-termine the redshift evolution of the best-fitting Schechter-function model parameters, Φ (cid:63) and L (cid:63) , as depicted in Fig. 5by the blue solid lines. Since, for the large number of sourceswe have here, S = − ( L ) is χ distributed, the σ errors(grey area in Fig. 5) are based on the likelihood ratios andcorrespond to ∆ χ = . For comparison the best-fitting val-ues (with σ errors) for four redshift bins, calculated usingthe / V max method, are also plotted here as blue points.With the redshift evolution of the Schechter-functionparameters as described in Equation 7, it is straightforwardto calculate the luminosity function at any redshift between z = . and z = . . We therefore plot the LFs determined us-ing the maximum-likelihood method alongside the results ofthe previous subsection in Fig. 3 as black solid lines. The σ errors (grey area) are again based on the maximum-likelihood ratios and correspond to ∆ χ = . Having constructed the rest-frame 250- µ m LFs, it is nowpossible to establish the redshift evolution of the star for-mation rate density ( ρ SFR ). This is achieved by first in-tegrating the LFs weighted by the total IR luminosity( − µ m ); for consistency this was performed again as-suming the SMG SED of Micha(cid:32)lowski et al. (2010) (which,although peaking at a wavelength consistent with a tempera-ture of 35 K, yields a bolometric luminosity (cid:39) times greaterthan a simple 35 K modified blackbody template), and theLFs were integrated down to a lower luminosity limit of . × L (cid:63) . The resulting inferred total IR luminosity density, ρ IR was then converted into dust-obscured star-formationrate density, ρ SFR , at each redshift using the scaling factorof K IR = . × − M (cid:12) year − erg − sHz from Madau & Dick-inson (2014), with an additional multiplicative factor of 0.63to convert from a Salpeter to a Chabrier IMF.The results are shown in Fig. 6. The red filled squaresdepict the values of ρ SFR derived from the four LFs shownin Fig. 3, established using the / V max method. The red solidline shows the evolution of ρ SFR as determined from thecontinuous form of the redshift evolution of the LF foundusing the maximum-likelihood method, with the grey areashowing the σ errors. Due to the fact we have chosen toparametrize the Schechter function normalisation parame- MNRAS , 1–17 (2017) alaxy far-IR LF out to z (cid:39) . . . . . . . . . − − Madau & Dickinson (2014)Behroozi , Wechsler & Conroy (2013)this work ( L max ) − IR onlyUV + IR ( L max )this work (1 /V max ) − IR onlyUV data (Parsa et al . z l og ( ρ S F R / M (cid:12) y r − M p c − ) Figure 6.
The star-formation rate density, ρ SFR , as a function of redshift. The black dotted and dashed lines represent the recentparametric descriptions of the redshift evolution of the ρ SFR provided by Madau & Dickinson (2014) and Behroozi et al. (2013) respectively(both converted to a Chabrier IMF). The red filled squares show our estimates of the IR SFRDs as determined using the / V max method,converting from IR luminosity density to SFRD using the conversion factor of K IR = . × − M (cid:12) year − erg − sHz given in Madau &Dickinson (2014), multiplied by a factor of 0.63 to convert from a Salpeter to a Chabrier IMF. The red solid line (with the σ uncertaintydepicted by the grey-shaded area) shows our analogous estimate of the redshift evolution of the IR SFRD, as determined from theluminosity-weighted integration of the LF derived through the maximum-likelihood method (with the functional form given in Equation8). The blue squares are the UV SFRD estimates based on the results of Parsa et al. (2016), converting from the rest-frame UV (1500˚A)luminosity to UV-visible SFR using the factor K UV = . × − M (cid:12) year − erg − sHz from Madau & Dickinson (2014) (again with a furthercorrection factor of × ρ SFR given by adding the UV and IR results. ter, Φ (cid:63) as a Gaussian function , and have chosen a simplelinear dependence of redshift for log ( L (cid:63) ) (see Equation 7),the redshift evolution of ρ SFR , depicted as a red solid line, isalso a Gaussian: ρ SFR IR ( z ) = . . √ π e − ( z − . ) × . . (8)The blue squares give the values of UV-visible ρ SFR derived from the Parsa et al. (2016) ρ UV results (inte-grated down to M UV = − ) by converting to UV-visible ρ SFR using the Madau & Dickinson (2014) scaling of K UV = . × − M (cid:12) year − erg − sHz , again also multiplied by afactor of 0.63 to convert to a Chabrier IMF. The blue solidline is a best-fitting Gaussian function to the UV results: ρ SFR UV ( z ) = . . √ π e − ( z − . ) × . . (9)The black solid line shows the total ρ SFR calculated bysimply adding the IR and UV estimates (Equations 8 and 9respectively). The black dotted and dashed lines depict the alternative functional forms of the redshift evolution of total ρ SFR provided by Madau & Dickinson (2014) and Behrooziet al. (2013) respectively (after conversion to a ChabrierIMF).
It is interesting to compare our inferred total IR LFs withthose that have been derived from the
Herschel surveys.First, in Fig. 7, we compare our LFs with those produced byMagnelli et al. (2013) from the deepest PACS surveys in theGOODS fields. Magnelli et al. (2013) utilised the 70, 100 and160- µ 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 blindPACS source extraction, but also created a source catalogueusing positional priors, with the positions of
Spitzer
IRAC3.6- µ m sources used to extract sources from the Spitzer
MIPS 24- µ m imaging, which were then in turn used as MNRAS000
MIPS 24- µ m imaging, which were then in turn used as MNRAS000 , 1–17 (2017) M. Koprowski et al. − − − − − . < z < . . < z < . . < z < . − − − − − . < z < . . < z < . L IR / L (cid:12) Φ / M p c − d e x − Figure 7.
A comparison of the IR LFs found in this work with those of Magnelli et al. (2013), depicted as black solid lines and coloureddashed lines respectively. The solid lines are our IR LFs plotted using Equation 3 (scaled to IR luminosity using the average SMG SEDof Micha(cid:32)lowski et al. 2010), with the parameters found at the median redshift of each bin using Equation 7. It can be seen that ourLFs and those of Magnelli et al. (2013) generally agree reasonably well near the break luminosity, but differ substantially at both faintand bright luminosities in every redshift bin. As discussed further in the text, much of this apparent difference is driven by the differentparameterizations adopted. At the bright end, the Magnelli et al. (2013) sample in fact contains only one object with an estimated L > L (cid:12) , and their adopted bright-end fixed power-law fit results in a severe over-prediction (by a factor (cid:39) − ) of the number ofbright sub-mm sources actually found in the degree-scale SCUBA-2 surveys. At the faint end, Magnelli et al. (2013) fixed all their LFsto have a faint-end power-law slope of − . , whereas we have adopted a slightly shallower slope of − . , derived by incorporating theALMA HUDF results at z (cid:39) . 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 photomet-ric redshift information, and the total infrared (8–1000- µ m)luminosities were inferred by fitting the 70, 100 and 160- µ mphotometry with the SED template library of Dale & Helou(2002). Magnelli et al. (2013) used this information to con-struct LFs in six redshift bins using the / V max method, andfitted the binned LF data with a double power-law functionas used previously by Magnelli et al. (2009, 2011). Becausethese LFs are based on the PACS data, which do not extendto wavelengths longer than 160 µ m, Magnelli et al. (2013)confined their analysis to redshifts z < . .Our results are compared with the Magnelli et al. (2013)LFs in Fig.7, using the five redshift bins adopted by Magnelliet al. (2013). The coloured dashed lines show the Magnelliet al. (2013) LFs, while the solid line in each panel showsour LF at the medium redshift of each bin. With the pos-sible exception of the highest redshift bin, the two sets ofLFs agree fairly well near the break luminosity. As a re-sult, it transpires that, as shown in the next section, the de-rived IR luminosities densities are not very different (albeitthe Magnelli et al. (2013) results are systematically higher).Nevertheless, it is clear that the LFs diverge at both lowerand higher luminosities, with the Magnelli et al. (2013) LFsindicating larger numbers of both faint and bright sources.However, to some extent this difference is exaggerated bythe different parameterisations adopted. Specifically, the difference at the bright end is due, atleast in part, to the fact that Magnelli et al. (2013) adopted adouble power-law function, while we have fitted a Schechterfunction (with an exponential cutoff at bright luminosities).In reality, the Magnelli et al. (2013) sample contains only oneobject with an estimated L > L (cid:12) and only (cid:39) sourceswith L > × L (cid:12) across the entire redshift range (theirfigure 8). As a consequence, the bright end of the Magnelliet al. (2013) LFs is fairly unconstrained, and much of thedifference seen in Fig. 7 is actually driven by their adoptionof a fixed bright-end power-law slope. In fact, as discussedbelow, integrated over the redshift range . < z < . theMagnelli et al. (2013) LFs predict (cid:39) times as many bright850- µ m sources in the wide-area SCUBA-2 sources than areactually observed, and so their LFs are clearly much toohigh at the bright end.At the faint end, Magnelli et al. (2013) again adopteda fixed power-law slope, with φ ∝ L − . . This was based onthe local value derived by Sanders et al. (2003), and sim-ply fixed at higher redshifts because the PACS data didnot really enable meaningful constraints to be placed on thefaint-end slope at z > . . Thus, in Fig. 7, our LFs divergefrom those of Magnelli et al. (2013) towards fainter luminosi-ties because we have adopted a slightly shallower faint-endslope ( φ ∝ L − . ), based on the ALMA measurements at . < z < . . As discussed in more detail next, some other Herschel -based estimates of the faint end of the LF in factuse an even flatter faint-end slope ( φ ∝ L − . ; Gruppioni et al.2013) and, as shown in figure 11 of Magnelli et al. (2013), MNRAS , 1–17 (2017) alaxy far-IR LF out to z (cid:39) − − − − − . < z < . . < z < . . < z < . − − − − − . < z < . . < z < . . < z < . − − − − − . < z < . . < z < . . < z < . L IR / L (cid:12) Φ / M p c − d e x − Figure 8.
A comparison of the IR LFs found in this work with those of Gruppioni et al. (2013), depicted as black solid lines andcoloured areas respectively. The dashed black lines are the best-fitting modified Schechter functions fitted by Gruppioni et al. (2013) totheir binned LF data. The faint-end slope assumed by Gruppioni et al. (2013) was fixed at − z (cid:39) value of − . . In most redshift bins, there is reasonable agreement at the faint end of the LF. However,certainly at z > , the Gruppioni et al. (2013) LFs are much higher (or the sources are much brighter) at the bright end. This differencemanifests itself clearly in the inferred values of IR ρ SFR , where the Gruppioni et al. (2013) LFs yield significantly higher values than thepresent study, especially at z (cid:39) (Fig. 10). In addition, the predicted cumulative number counts at 850 µ m are very different (see Fig. 9and Table 5); the Gruppioni et al. (2013) LFs predict > µ m sources with flux densities S > mJy in the maps used in thiswork, whereas in fact only 17 such objects are actually detected. even using the Spitzer
MIPS data to extend estimates of theLF to fainter luminosities cannot really distinguish betweena faint-end power-law slope of − . or − . . It is thereforeboth interesting and reassuring that our own measured valueof the faint-end slope ( − . ) lies midway between the valuesadopted in the Herschel studies.Next, in Fig. 8, we compare our LFs with those pro-duced by Gruppioni et al. (2013), which were based onPACS far-IR sources extracted from the
Herschel
PEP sur-vey data at 70, 100 and 160 µ m within the COSMOS,ECDFS, GOODS-N and GOODS-S fields (covering a to-tal area (cid:39) . deg ). As well as covering a larger area thanMagnelli et al. (2013) (albeit of course to shallower depths),Gruppioni et al. (2013) extended the SED fitting to includethe 250, 350 and 500- µ m imaging provided by the Herschel
SPIRE imaging in the same fields (by the HerMES sur-vey; Oliver et al. 2012), and attempted to extend their LFanalysis out to redshifts z (cid:39) . Gruppioni et al. (2013) de-tected 373, 7176 and 7376 sources at 70, 100 and 160 µ m ,respectively, used cross-matching (via 24- µ m MIPS imaging)with the shorter-wavelength (optical+near-IR) data in thefields to provide redshift information, and used SED fittingwith a range of templates to the PACS+SPIRE photome-try to derive total IR luminosities. Again, Gruppioni et al.(2013) constructed the LFs using the / V max method, but in their study they fitted the binned LF data with a modifiedSchechter function.Our results are compared with the Gruppioni et al.(2013) LFs in Fig. 8, using the nine redshift bins adoptedby Gruppioni et al. (2013). The colour-shaded areas showthe Gruppioni et al. (2013) results, with the dashed linesrepresenting their best-fitting modified Schechter functions.The solid black lines show our own estimates of the LFat these redshifts. While at the faint-end the LFs are infairly good agreement (as mentioned above, Gruppioni et al.(2013) in fact adopt a slightly shallower faint-end slope), atthe bright-end there is again significant disagreement. More-over, at high redshifts the Gruppioni et al. (2013) LFs aresignificantly higher around the break luminosity, resulting insubstantially larger values of IR luminosity density at z (cid:39) (see next section).One way to quantify how much our new IR LFs dif-fer from those produced from the aforementioned Herschel -based studies is to calculate how many bright sub-mmsources each evolving LF predicts should be present in theSC2LS survey data utilised here. We show the results ofthis in Fig. 9, and quantify the differences at flux densities S > mJy in Table 5. Along with a range of data derivedfrom the SCUBA-2 and ALMA-based literature, the solidblack and blue lines in the left-hand panel of Fig. 9 show the MNRAS000
SPIRE imaging in the same fields (by the HerMES sur-vey; Oliver et al. 2012), and attempted to extend their LFanalysis out to redshifts z (cid:39) . Gruppioni et al. (2013) de-tected 373, 7176 and 7376 sources at 70, 100 and 160 µ m ,respectively, used cross-matching (via 24- µ m MIPS imaging)with the shorter-wavelength (optical+near-IR) data in thefields to provide redshift information, and used SED fittingwith a range of templates to the PACS+SPIRE photome-try to derive total IR luminosities. Again, Gruppioni et al.(2013) constructed the LFs using the / V max method, but in their study they fitted the binned LF data with a modifiedSchechter function.Our results are compared with the Gruppioni et al.(2013) LFs in Fig. 8, using the nine redshift bins adoptedby Gruppioni et al. (2013). The colour-shaded areas showthe Gruppioni et al. (2013) results, with the dashed linesrepresenting their best-fitting modified Schechter functions.The solid black lines show our own estimates of the LFat these redshifts. While at the faint-end the LFs are infairly good agreement (as mentioned above, Gruppioni et al.(2013) in fact adopt a slightly shallower faint-end slope), atthe bright-end there is again significant disagreement. More-over, at high redshifts the Gruppioni et al. (2013) LFs aresignificantly higher around the break luminosity, resulting insubstantially larger values of IR luminosity density at z (cid:39) (see next section).One way to quantify how much our new IR LFs dif-fer from those produced from the aforementioned Herschel -based studies is to calculate how many bright sub-mmsources each evolving LF predicts should be present in theSC2LS survey data utilised here. We show the results ofthis in Fig. 9, and quantify the differences at flux densities S > mJy in Table 5. Along with a range of data derivedfrom the SCUBA-2 and ALMA-based literature, the solidblack and blue lines in the left-hand panel of Fig. 9 show the MNRAS000 , 1–17 (2017) M. Koprowski et al. this work (0 . < z < . . (2013) (0 . < z < . . < z < . . (2013) (0 . < z < . . (2017)Hatsukade et al . (2016)Karim et al . (2013)Simpson et al . (2015)Fujimoto et al . (2016)Oteo et al . (2016)Geach et al . (2017) this work (0 . < z < . . (2013) , T d = 35 KGruppioni et al . (2013) , T d = 45 KGruppioni et al . (2013) , T d = 55 KGruppioni et al . (2013) , T d = 65 KGruppioni et al . (2013) , T d = 75 KGruppioni et al . (2013) , T d = 85 K S µ m / mJy N ( > S ) / d e g − Figure 9.
Cumulative number counts as a function of 850- µ m flux density. In the left-hand panel we show a range of publishedcumulative counts from the SCUBA-2 and ALMA literature (Karim et al. 2013; Simpson et al. 2015; Hatsukade et al. 2016; Fujimotoet al. 2016; Oteo et al. 2016; Dunlop et al. 2017; Geach et al. 2017), along with the predicted number counts calculated from our newIR LFs (presented and discussed in Section 6.1), after integration over the redshift range of . < z < . (black solid line, with the σ uncertainty shown by the grey-shaded area), and after integration over the redshift range . < z < . (red dashed line). For comparison,the 850- µ m number counts predicted by integrating the Gruppioni et al. (2013) evolving LF over . < z < . are plotted as the bluesolid line, while the counts predicted by integrating the Magnelli et al. (2013) LFs over . < z < . are indicated by the green dashedline. These predictions have all been made assuming the SED of Micha(cid:32)lowski et al. (2010) as used throughout this study. Unsurprisingly,since our study is based primarily on 850- µ m data, the number counts predicted by our new IR LFs agree well with the data, whereas the Herschel -based LFs over predict the counts at the bright end by over an order-of-magnitude (as expected, given the dramatic differencesat the bright end of the LFs, as shown in Fig. 7 and Fig. 8). These over-predictions are quantified in Table 5, for 850- µ m sources withflux densities higher than in the combined area of all the S2CLS fields used in this work ( (cid:39) . deg ). The actual number of S > mJysources in this area in the real data is 17 sources in the redshift range . < z < . , or 4 sources for . < z < . . In the right-hand panelwe explore the extent to which such discrepancies can potentially be solved by varying the adopted SED template. To illustrate this weshow the effect of converting the Gruppioni et al. (2013) LFs into predicted 850- µ m number counts for adopted modified black-bodySEDs of increasing temperature. As expected, things agree well at the faint end for T (cid:39) K, but a change to characteristic temperaturesof T (cid:39) K is required to suppress the predicted counts at the bright end. There is no evidence to support such a rapid change in typicalSED effective temperature over such a short range in 850- µ m flux density of IR luminosity. predicted cumulative 850- µ m number counts derived fromour new LFs and from those produced by Gruppioni et al.(2013) (see Fig. 8) after integration over the redshift range . < z < . . The dashed red and green lines show the cor-responding predictions from our own LFs and those pro-duced by Magnelli et al. (2013) (see Fig. 7) after integrationover the redshift range . < z < . . Unsurprisingly, sinceour study is based primarily on 850- µ m data, the numbercounts predicted by our new IR LFs agree well with the ob-served counts, whereas the Herschel -based LFs over predictthe counts at the bright end by over an order-of-magnitude(as expected, given the dramatic differences at the brightend of the LFs, as shown in Fig. 7 and Fig. 8).The extent of this difference is quantified further in Ta-ble 5, where we tabulate how many 850- µ m sources withflux densities S > mJy in the combined area of all theS2CLS fields used here (1.5 deg ) a given LF predicts afterintegration over the appropriate redshift range. The actualdetected number of 850- µ m sources with S > mJy is 17 over the redshift range . < z < . (consistent with ourpredicted 20) of which 4 have . < z < . (consistent withour predicted 8). By contrast the Gruppioni et al. (2013)LFs predict > such galaxies at . < z < . (higher bya factor of (cid:39) − ), while the Magnelli et al. (2013) LFspredict 140 such sources at . < z < . (again high by afactor of − ).As used consistently throughout this study, conversionbetween total IR luminosity, and observed 850- µ m flux den-sity was performed using the average sub-mm source SED ofMicha(cid:32)lowski et al. (2010). In the right-hand panel of Fig. 9we explore what sort of change in the adopted templatewould be required to bring the bright sub-mm number-countpredictions of the Herschel -derived IR LFs into agreementwith reality. Here we show the effect of converting the Grup-pioni et al. (2013) LFs into predicted 850- µ m number countsfor adopted modified black-body SEDs of increasing tem-perature. As expected, things agree well at the faint endfor T (cid:39) K, but a change to characteristic temperatures of
MNRAS , 1–17 (2017) alaxy far-IR LF out to z (cid:39) . . . . . . . . . − − Madau & Dickinson (2014)Behroozi , Wechsler & Conroy (2013)UV + IR ( L max )Gruppioni et al . (2013) + UVMagnelli et al . (2013) + UVRowan − Robinson et al . (2016) + UV z l og ( ρ S F R / M (cid:12) y r − M p c − ) Figure 10.
Total ρ SFR as a function of redshift. The blue solid line shows the results of this work (see Fig. 6), while the data points showthe estimates derived from the
Herschel -based work of Magnelli et al. (2013) (green circles), Gruppioni et al. (2013) (red squares) andRowan-Robinson et al. (2016) (magenta triangles) after addition of the UV estimates of unobscured ρ SFR from Parsa et al. (2016). Theblack dotted and dashed lines represent the recent parametric descriptions of the redshift evolution of ρ SFR found by Madau & Dickinson(2014) and Behroozi et al. (2013) respectively, adopting a Chabrier IMF.
Table 5.
The predicted number of 850- µ m sources with flux den-sities S > mJy in the combined area of all the fields used inthis study ( (cid:39) . deg ) produced by integrating the three IR LFsdiscussed in Section 6.1 and shown in Fig. 7 and Fig. 8. As dis-cussed in the caption to Fig. 9, and indicated here, these predic-tions were based on integration over the redshift range . < z < . or . < z < . as appropriate. The actual number of sources in oursample with S > mJy is 17 in the redshift range . < z < . ,of which four have . < z < . . The predictions are based on theMicha(cid:32)lowski et al. (2010) average sub-mm source SED. Clearly,the Herschel -derived LFs wildly over-predict the actual numberof bright 850- µ m sources.Study z range N ( > ) / . − This Work . < z < . . < z < . . < z < . . < z < . T (cid:39) K is required at the bright end to adequately sup-press the predicted number counts. There is no evidence tosupport such a rapid change in typical SED temperatureover such a short range in 850- µ m flux density or IR lumi-nosity. Indeed, while several studies have reported a corre-lation between IR luminosity and the effective temperature of the typical galaxy SED, the derived slope of this corre-lation is relatively modest, with average dust temperaturereported to rise from T (cid:39) K for L IR (cid:39) L (cid:12) to T (cid:39) Kfor L IR (cid:39) L (cid:12) (Amblard et al. 2010; Hwang et al. 2010;Smith et al. 2012).We conclude that the steep drop-off at the bright endof the sub-mm number counts does simply reflect a near ex-ponential decline at the bright end of the IR LF, and the Herschel results have been contaminated and biased highat the bright end by a mix of blending issues, source iden-tification (and hence redshift) errors, and possibly also theadoption of a double power-law LF (see also Bethermin et al.2017; Liu et al. 2017).
In Fig. 10 we compare the redshift evolution of cosmic star-formation density derived here with the evolution inferredfrom the
Herschel -based studies undertaken by Magnelliet al. (2013), Gruppioni et al. (2013) and Rowan-Robinsonet al. (2016). All results have been converted to a ChabrierIMF, and the values of obscured ρ SFR derived from the sub-mm/far-IR studies have been added to the UV-derived ρ SFR from Parsa et al. (2016), so that all values plotted here rep-resent total ρ SFR .To derive the Magnelli et al. (2013) results (green pointsin Fig. 10) we performed the luminosity-weighted integral
MNRAS000
MNRAS000 , 1–17 (2017) M. Koprowski et al. − − − − − − − z = 0 . z = 1 . z = 1 . z = 2 . z = 2 . z = 3 . z = 3 . z = 4 . z = 4 . L / W Hz − Φ / M p c − d e x − Figure 11.
Our derived rest-frame 250- µ m luminosity functions (LFs) for a range of redshifts, determined from the maximum-likelihoodmethod (plotted using Equations 3 and 7, using the best-fit parameter values listed in Table 4). The combined impact of rising-then-fallingdensity evolution and positive luminosity evolution with redshift is clear. The result is that, as with AGN, many of the most luminoussources are to be found at high redshifts even though overall IR luminosity density is in decline beyond z (cid:39) . The implications for theexpected redshift distribution of 850- µ m sources as a function of flux density are illustrated in Fig. 12. of their IR LFs (shown in Fig. 7) from a lower limit of . × L (cid:63) , where L (cid:63) is the characteristic luminosity foundfor our sample (Table 3). We then derived ρ SFR using theconversion factor given by Madau & Dickinson (2014), be-fore applying the IMF correction factor of 0.63 and addingthe UV-based results. To derive the corresponding valuesfrom Gruppioni et al. (2013), we simply converted their ownIR luminosity density estimates using the same procedure.Rowan-Robinson et al. (2016) attempted to extend the ex-ploitation of the
Herschel
HerMES survey to the highestredshifts via SPIRE-based source selection (over an area of (cid:39) deg ), including sources detected only at 500 µ m. Wehave adopted their derived estimates of obscured ρ SFR , con-verted to a Chabrier (2003) IMF, and again added the UV-based estimates to produce the (magenta) data points shownin Fig. 10.Despite the differences in the LFs discussed extensivelyabove, with the exception of the rather high Gruppioni et al.(2013) results around z (cid:39) , all three Herschel -based studiesyield estimates of ρ SFR consistent with those derived here(as indicated by the blue curve) up to redshifts z (cid:39) . , al-beit the uncertainties in the Herschel -derived estimates aregenerally very large. At higher redshifts the uncertaintiesin the
Herschel -based estimates start to approach 0.5 dex,but even so the smooth decline in ρ SFR beyond z (cid:39) seenin the present study (and consistent with that reported byBehroozi et al. 2013; Madau & Dickinson 2014; Dunlop et al.2017; Bourne et al. 2017; Liu et al. 2017) is inconsistent withthe 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 identificationand redshift estimation stemming from the large-beam long-wavelength SPIRE data, as well as potential blending issues,and extrapolation of ρ SFR from the most extreme sources. In-terestingly, the Rowan-Robinson et al. (2016) results seemsomewhat low near the peak of activity at z (cid:39) . This maysuggest that some subset of numerous far-IR sources whichshould lie at z (cid:39) have been erroneously placed at higherredshifts, where the resulting boost in inferred luminos-ity coupled with the smaller cosmological volumes will in-evitably result in artificially high estimates of ρ SFR . Finally, we discuss how the form of the LF evolution uncov-ered here naturally explains the apparent ‘down-sizing’ ofthe sub-mm source population. There is tentative but per-sistent evidence that the median redshift of sub-mm selectedgalaxies 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(cid:32)lowski et al. 2016), and hence increasing far-IR lumi-nosity. In Fig. 11 we plot our evolving LFs as a functionof redshift from z (cid:39) . to z (cid:39) . . At low redshifts, the in-crease in both Φ (cid:63) and L (cid:63) produces an increase in both thenumber density of sources, and luminosity density, ρ IR , withincreasing redshift (cf Fig. 10). At high redshifts the declinein Φ (cid:63) progressively overcomes the continued positive evolu-tion of L (cid:63) to produce a decline in both these quantities, but MNRAS , 1–17 (2017) alaxy far-IR LF out to z (cid:39) . . . . . . . . . . . . . . . . . . . . . . . . . . . . f r a c t i o np e r b i n S > S > S > S >
10 mJy
Figure 12.
The redshift distributions for 850- µ m source samples selected at different limiting flux densities, as predicted from theanalytic form of the evolving LF shown in Fig. 11 (Equations 3 and 7, with best-fit parameter values listed in Table 4), derived usingthe maximum-likelihood method. It can be seen that the brightest sources tend to lie at higher redshifts, with the peak redshift shiftingfrom z (cid:39) . to z (cid:39) over the flux-density range explored here. This is a direct consequence of the redshift evolution of the IR LF (Figure11), in particular the continued increase of the characteristic luminosity with redshift. The broad range of redshifts found in each panelexplains why it has proved difficult to derive a statistically-significant correlation between S and z , but the predicted relative lack ofbright sub-mm sources at z < . accords well with the latest observations (e.g. Micha(cid:32)lowski et al. 2016). Although bright 850 µ m sourcesat z > . clearly exist, their relatively low number density, and the overall decline of the LF as seen in Fig. 11, means that their discoverydoes not conflict with the result that ρ SFR declines beyond z (cid:39) − . (as seen in Fig. 10). it can be seen that the continued positive evolution of L (cid:63) means that the most luminous sources persist, or indeed arepreferentially found at the highest redshifts explored here.To better connect with observables, we have used ourevolving LF to calculate the predicted redshift distributionof 850- µ m sources as a function of flux density. The resultsare shown for four different flux-density thresholds in Fig. 12.Here it can be seen that the peak in the redshift distributionis expected to naturally increase gradually from z (cid:39) . for S > mJy to z (cid:39) for S > mJy. This is in excellentaccord with what has been reported in the literature (e.g.Koprowski et al. 2014; Micha(cid:32)lowski et al. 2016) and clari-fies why, although dust-enshrouded star-formation is glob-ally less important than UV-visible star-formation activityat z > , bright sub-mm sources will continue to be discov-ered out to high redshifts. We have analysed the coverage of the far-IR luminosity–redshift plane provided by (sub-)mm-selected galaxy sam-ples extracted from the SCUBA-2 Cosmology Legacy Sur-vey and the ALMA imaging of the HUDF to make a newmeasurement of the evolving galaxy far-IR luminosity func-tion (LF) extending out to redshifts z (cid:39) . Using both di-rect ( / V max ) and maximum-likelihood methods we have de-termined the form and evolution of the rest-frame 250- µ m galaxy LF. This LF is well described by a Schechter functionwith a faint-end slope α (cid:39) − . (derived using the ALMAdata at z (cid:39) ) which displays a combination of rising-then-falling density evolution, and positive luminosity evolution.We converted our 250- µ m results to total IR luminosityusing the average sub-mm galaxy SED of Micha(cid:32)lowski et al.(2010), and have then compared our determination of theevolving IR LF with those derived by Magnelli et al. (2013)and Gruppioni et al. (2013) from the Herschel
PEP and Her-MES surveys. Our faint-end slope lies approximately mid-way between the values adopted in these studies, and the LFnormalization near the break luminosity L (cid:63) is also compara-ble at most redshifts. However, both of the Herschel -derivedLFs indicate a much larger number of very bright sources atall redshifts than is found in our JCMT/ALMA-based study.To check which result is correct we derive the predicted 850- µ m number counts from the Herschel -based LFs, and findthat, at bright flux densities ( S > mJy) they predictan order-of-magnitude more sources than are observed. Weexplore whether this discrepancy can be explained/resolvedby adoption of different SED templates, and conclude thatit cannot (without extreme temperatures, and more impor-tantly without a contrived and much more rapid dependenceof T on IR luminosity than has been reported in the litera-ture).We have utilised our measurement of the evolving IRLF to derive comoving IR luminosity density, and hence ob- MNRAS000
PEP and Her-MES surveys. Our faint-end slope lies approximately mid-way between the values adopted in these studies, and the LFnormalization near the break luminosity L (cid:63) is also compara-ble at most redshifts. However, both of the Herschel -derivedLFs indicate a much larger number of very bright sources atall redshifts than is found in our JCMT/ALMA-based study.To check which result is correct we derive the predicted 850- µ m number counts from the Herschel -based LFs, and findthat, at bright flux densities ( S > mJy) they predictan order-of-magnitude more sources than are observed. Weexplore whether this discrepancy can be explained/resolvedby adoption of different SED templates, and conclude thatit cannot (without extreme temperatures, and more impor-tantly without a contrived and much more rapid dependenceof T on IR luminosity than has been reported in the litera-ture).We have utilised our measurement of the evolving IRLF to derive comoving IR luminosity density, and hence ob- MNRAS000 , 1–17 (2017) M. Koprowski et al. scured star-formation rate density, which we then combinewith UV-estimates of unobscured activity (from Parsa et al.2016), to derive the evolution of total ρ SFR . Again we com-pare with values derived from the
Herschel -based studies,and find reasonable agreement out to z (cid:39) , but increas-ing disagreement at higher redshift. Specifically, consistentwith several other recent studies (e.g. Bourne et al. 2017;Dunlop et al. 2017; Liu et al. 2017) we find that ρ SFR de-clines beyond z (cid:39) − . and is dominated by UV-visiblestar-formation activity beyond z (cid:39) . In contrast, Gruppioniet al. (2013) and Rowan-Robinson et al. (2016) report essen-tially no decline in dust-obscured ρ SFR from z (cid:39) to z (cid:39) .Given the severe over-prediction of the 850- µ m counts pro-duced by the Herschel
IR LFs, we conclude that the highvalues reported from these studies most likely reflect prob-lems in source identification and redshift estimation arisingfrom the large-beam long-wavelength SPIRE data, as wellas potential blending issues, and extrapolation of ρ SFR fromthe most extreme sources.Finally, we show how the evolution of the IR LF as de-rived here (with its combination of rising-then-falling charac-teristic density ( Φ (cid:63) ), and positive evolution of characteristicluminosity density ( L (cid:63) )) with redshift, produces a decline ininferred ρ SFR beyond z (cid:39) − . while at the same time pre-dicting that the most luminous sub-mm sources will continueto be found out to very high redshifts ( z (cid:39) − ). Specifically,our evolving LF, with its combined luminosity+density evo-lution, predicts that the median redshift of sub-mm sourcesshould increase with increasing flux density, consistent withseveral reports in the recent literature. Such evolution is con-sistent with many studies of AGN evolution, suggesting thatboth dust-enshrouded star formation and AGN activity arestrongly linked to the growth of stellar mass in galaxies. ACKNOWLEDGEMENTS
KEKC and MPK acknowledge the support of the UK Sci-ence and Technology Facilities Council through grant num-ber ST/M001008/1. JEG is supported by the Royal Soci-ety. JSD, MJM and RJM acknowledge the support of theUK Science and Technology Facilities Council through grantnumber ST/M001229/1.The James Clerk Maxwell Telescope is now operated bythe East Asian Observatory on behalf of The National As-tronomical Observatory of Japan, Academia Sinica Instituteof Astronomy and Astrophysics, the Korea Astronomy andSpace Science Institute, the National Astronomical Observa-tories of China and the Chinese Academy of Sciences (grantno. XDB09000000), with additional funding support fromthe Science and Technology Facilities Council of the UnitedKingdom and par- ticipating universities in the United King-dom and Canada. The data utilised in this paper were takenas part of Program ID MJLSC02.This paper makes use of the following ALMA data:ADS/JAO.ALMA
REFERENCES
Amblard A., et al., 2010, A&A, 518, L9Aretxaga I., et al., 2007, MNRAS, 379, 1571Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57Bethermin M., et al., 2017, preprint, ( arXiv:1703.08795 )Biggs A. D., et al., 2011, MNRAS, 413, 2314Bolzonella M., Miralles J.-M., Pell´o R., 2000, A&A, 363, 476Bourne N., et al., 2017, MNRAS, 467, 1360Bouwens R. J., et al., 2015, ApJ, 803, 34Bouwens R. J., et al., 2016a, ApJ, 830, 67Bouwens R. J., et al., 2016b, ApJ, 833, 72Bowler R. A. A., et al., 2012, MNRAS, 426, 2772Bowler R. A. A., et al., 2014, MNRAS, 440, 2810Bowler R. A. A., et al., 2015, MNRAS, 452, 1817Bruzual G., Charlot S., 2003, MNRAS, 344, 1000Burgarella D., et al., 2013, A&A, 554, A70Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J.,Storchi-Bergmann T., 2000, ApJ, 533, 682Capak P. L., et al., 2015, Nat, 522, 455Caputi K. I., et al., 2007, ApJ, 660, 97Caputi K. I., Cirasuolo M., Dunlop J. S., McLure R. J., FarrahD., Almaini O., 2011, MNRAS, 413, 162Chabrier G., 2003, ApJ, 586, L133Chen C.-C., et al., 2016, preprint, ( arXiv:1601.02630 )Cirasuolo M., et al., 2007, MNRAS, 380, 585Cirasuolo M., McLure R. J., Dunlop J. S., Almaini O., FoucaudS., Simpson C., 2010, MNRAS, 401, 1166Condon J. J., 1992, ARA&A, 30, 575Cowie L. L., Barger A. J., Hsu L.-Y., Chen C.-C., Owen F. N.,Wang W.-H., 2017, ApJ, 837, 139Cucciati O., et al., 2012, A&A, 539, A31Dale D. A., Helou G., 2002, ApJ, 576, 159Dole H., et al., 2006, A&A, 451, 417Dunlop J. S., Peacock J. A., Savage A., Lilly S. J., Heasley J. N.,Simon A. J. B., 1989, MNRAS, 238, 1171Dunlop J. S., et al., 2010, MNRAS, 408, 2022Dunlop J. S., et al., 2013, MNRAS, 432, 3520Dunlop J. S., et al., 2017, MNRAS, 466, 861Dye S., et al., 2008, MNRAS, 386, 1107Elbaz D., et al., 2010, A&A, 518, L29Elbaz D., et al., 2011, A&A, 533, A119Ellis R. S., et al., 2013, ApJ, 763, L7Finkelstein S. L., et al., 2015, ApJ, 810, 71Fujimoto S., Ouchi M., Ono Y., Shibuya T., Ishigaki M., NagaiH., Momose R., 2016, ApJS, 222, 1Furusawa H., et al., 2008, ApJS, 176, 1Geach J. E., et al., 2013, MNRAS, 432, 53Geach J. E., et al., 2017, MNRAS, 465, 1789Griffin M. J., et al., 2010, A&A, 518, L3Grogin N. A., et al., 2011, ApJS, 197, 35Gruppioni C., et al., 2013, MNRAS, 432, 23Gwyn S. D. J., 2012, AJ, 143, 38Hatsukade B., Ohta K., Yabe K., Seko A., Makiya R., AkiyamaM., 2015, ApJ, 810, 91Hatsukade B., et al., 2016, PASJ, 68, 36Holland W. S., et al., 2013, MNRAS, 430, 2513Hwang H. S., et al., 2010, MNRAS, 409, 75Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K.,2017, preprint, ( arXiv:1702.04867 )Ivison R. J., et al., 2002, MNRAS, 337, 1Ivison R. J., et al., 2005, MNRAS, 364, 1025Ivison R. J., et al., 2007, MNRAS, 380, 199Karim A., et al., 2013, MNRAS, 432, 2Koprowski M. P., Dunlop J. S., Micha(cid:32)lowski M. J., Cirasuolo M.,Bowler R. A. A., 2014, MNRAS, 444, 117Koprowski M. P., et al., 2016a, MNRAS, 458, 4321Koprowski M. P., et al., 2016b, ApJ, 828, L21MNRAS , 1–17 (2017) alaxy far-IR LF out to z (cid:39) Lawrence A., et al., 2007, MNRAS, 379, 1599Le Floc’h E., et al., 2005, ApJ, 632, 169Liu D., et al., 2017, preprint, ( arXiv:1703.05281 )Lutz D., et al., 2011, A&A, 532, A90Madau P., 1995, ApJ, 441, 18Madau P., Dickinson M., 2014, ARA&A, 52, 415Magnelli B., Elbaz D., Chary R. R., Dickinson M., Le Borgne D.,Frayer D. T., Willmer C. N. A., 2009, A&A, 496, 57Magnelli B., Elbaz D., Chary R. R., Dickinson M., Le Borgne D.,Frayer D. T., Willmer C. N. A., 2011, A&A, 528, A35Magnelli B., et al., 2013, A&A, 553, A132Marshall H. L., Tananbaum H., Avni Y., Zamorani G., 1983, ApJ,269, 35McCracken H. J., et al., 2012, A&A, 544, A156McLeod D. J., McLure R. J., Dunlop J. S., Robertson B. E., EllisR. S., Targett T. A., 2015, MNRAS, 450, 3032McLeod D. J., McLure R. J., Dunlop J. S., 2016, MNRAS, 459,3812McLure R. J., et al., 2013, MNRAS, 432, 2696Mehta V., et al., 2017, ApJ, 838, 29Micha(cid:32)lowski M., Hjorth J., Watson D., 2010, A&A, 514, A67Micha(cid:32)lowski M. J., et al., 2016, preprint, ( arXiv:1610.02409 )Miyazaki S., et al., 2002, PASJ, 54, 833Oesch P. A., et al., 2014, ApJ, 786, 108Oliver S. J., et al., 2012, MNRAS, 424, 1614Oteo I., Zwaan M. A., Ivison R. J., Smail I., Biggs A. D., 2016,ApJ, 822, 36Parsa S., Dunlop J. S., McLure R. J., Mortlock A., 2016, MNRAS,456, 3194Patel H., Clements D. L., Vaccari M., Mortlock D. J., Rowan-Robinson M., P´erez-Fournon I., Afonso-Luis A., 2013, MN-RAS, 428, 291Pilbratt G. L., et al., 2010, A&A, 518, L1Poglitsch A., et al., 2010, A&A, 518, L2Pope A., Borys C., Scott D., Conselice C., Dickinson M.,Mobasher B., 2005, MNRAS, 358, 149Robertson B. E., et al., 2013, ApJ, 768, 71Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015,ApJ, 802, L19Rodighiero G., et al., 2010, A&A, 515, A8Roseboom I. G., et al., 2013, MNRAS, 436, 430Rowan-Robinson M., et al., 2016, MNRAS, 461, 1100Sanders D. B., Mazzarella J. M., Kim D.-C., Surace J. A., SoiferB. T., 2003, AJ, 126, 1607Sanders D. B., et al., 2007, ApJS, 172, 86Schinnerer E., et al., 2010, ApJS, 188, 384Schmidt M., 1968, ApJ, 151, 393Simpson J. M., et al., 2015, ApJ, 807, 128Smith A. J., et al., 2012, MNRAS, 419, 377Taniguchi Y., et al., 2007, ApJS, 172, 9Umehata H., et al., 2015, ApJ, 815, L8Vieira J. D., et al., 2013, Nat, 495, 344Wall J. V., Pope A., Scott D., 2008, MNRAS, 383, 435Walter F., et al., 2016, ApJ, 833, 67Wardlow J. L., et al., 2011, MNRAS, 415, 1479Weiß A., et al., 2013, ApJ, 767, 88This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000