SCUBA-2 imaging of colour-selected Herschel sources

Red, redder, reddest: SCUBA-2 imaging of colour-selected Herschel sources

S. Duivenvoorden, S. Oliver, J. M. Scudder, J. Greenslade, D. A. Riechers, S. M. Wilkins, V. Buat, S. C. Chapman, D. L. Clements, A. Cooray, K. E. K. Coppin, H. Dannerbauer, G. De Zotti, J. S. Dunlop, S. A. Eales,A. Efstathiou, D. Farrah, J. E. Geach, W. S. Holland, P. D. Hurley, R. J. Ivison, L. Marchetti, G. Petitpas, M. T. Sargent, D. Scott, M. Symeonidis, M. Vaccari, J. D. Vieira, L. Wang, J. Wardlow and M. Zemcov
Affiliations are listed at the end of the paper
E-mail: S.Duivenvoorden@Sussex.ac.uk
Accepted xxx. Received xxx; in original form xxx
Abstract

High-redshift, luminous, dusty star forming galaxies (DSFGs) constrain the extremity of galaxy formation theories. The most extreme are discovered through follow-up on candidates in large area surveys. Here we present 850 m SCUBA-2 follow-up observations of 188 red DSFG candidates from the Herschel Multi-tiered Extragalactic Survey (HerMES) Large Mode Survey, covering 274 deg. We detected 87 per cent with a signal-to-noise ratio 3 at 850 m. We introduce a new method for incorporating the confusion noise in our spectral energy distribution fitting by sampling correlated flux density fluctuations from a confusion limited map. The new 850 m data provide a better constraint on the photometric redshifts of the candidates, with photometric redshift errors decreasing from to . Comparison spectroscopic redshifts also found little bias (). The mean photometric redshift is found to be 3.6 with a dispersion of and we identify 21 DSFGs with a high probability of lying at . After simulating our selection effects we find number counts are consistent with phenomenological galaxy evolution models. There is a statistically significant excess of WISE-1 and SDSS sources near our red galaxies, giving a strong indication that lensing may explain some of the apparently extreme objects. Nevertheless, our sample should include examples of galaxies with the highest star formation rates in the Universe ( Myr).

keywords:
Infrared: galaxies – galaxies: high-redshift – galaxies: starburst.
pubyear: 2018pagerange: Red, redder, reddest: SCUBA-2 imaging of colour-selected Herschel sourcesA

1 Introduction

Over the last few decades, great progress has been made in understanding the star formation history of the Universe (see e.g. review by Madau & Dickinson 2014). It has become apparent that observing at UV and optical wavelengths is insufficient as a large fraction of the star formation is obscured, resulting in dusty star forming galaxies (DSFGs see e.g reviews by Lonsdale et al. 1984; Cesarsky et al. 1996; Smail et al. 1997; Burgarella et al. 2013 and Casey et al. 2014). The most extreme forms of obscured star formation at high-redshift still pose serious challenges to galaxy evolution models (e.g. Baugh et al., 2005; Lacey et al., 2010; Narayanan et al., 2010; Hayward et al., 2013; Béthermin et al., 2017). The discovery and characterisation of the rarest and most extreme galaxies (star formation rates, SFR, Myr, number densities , Gruppioni et al. 2013) is thus an important goal, but requires large volume surveys at long wavelengths.

This is now possible with deep large-area surveys ( deg) at far infrared (FIR) and sub-mm wavelengths with e.g. the South Pole Telescope (SPT, Carlstrom et al., 2011) and the Herschel Space Observatory (Pilbratt et al., 2010).

Follow-up of SPT sources has been very successful in finding high redshift DSFGs (Vieira et al., 2013; Weiß et al., 2013; Strandet et al., 2016, 2017). The SPT source selection at a wavelength of 1.4 mm has however a broader redshift distribution than Herschel detected sources (Greve et al., 2012)

Herschel surveys cover a huge area (the largest being HerMES Oliver et al. 2012 and H-ATLAS Eales et al. 2010) and while most detections are associated with starburst galaxies (e.g, Casey et al. 2012a; Casey et al. 2012b) it has been clearly demonstrated that selecting those with red colours is extremely efficient for identifying a tail extending towards higher redshift () (Cox et al., 2011; Riechers et al., 2013; Dowell et al., 2014; Asboth et al., 2016; Ivison et al., 2016; Riechers et al., 2017). The challenge now is using these very large Herschel surveys to find and systematically study, large, homogeneous samples of rare, extremely luminous, sources.

Asboth et al. (2016) probed this high-redshift population in the largest HerMES field, the HerMES Large Mode Survey (HeLMS, covering approximately 300 deg) by selecting all bright “500 m riser” () DSFGs candidates. This sample was selected over an area a factor of 13 times larger than previous 500 m riser HerMES surveys (Dowell et al., 2014). The number of sources that fulfilled these criteria (477) is an order of magnitude higher than predicted by galaxy evolution models (Béthermin et al., 2011, 2012; Dowell et al., 2014)

Another large 600 deg red DSFGs search in the H-ATLAS survey (Ivison et al., 2016) used a (30 mJy) detection threshold at in combination with and colour selection criteria to obtain a sample of 7961 candidate high-redshift DSFGs. After a visual inspection (Ivison et al., 2016) a sub-sample of 109 DSFGs, candidates were selected for follow up at longer wavelengths with SCUBA-2 or LABOCA.

All these red sources are candidates for high-luminosity sources. Some, particularly those with a flux density at 100 mJy, are likely to be strongly gravitationally lensed (Negrello et al., 2010; Conley et al., 2011; Nayyeri et al., 2016; Negrello et al., 2017) others may be blends (e.g. Scudder et al., 2016). Nevertheless, they are extremely interesting because, those that are not lensed, blended or otherwise boosted may represent the most active galaxies in cosmic history.

In this work, we present a follow-up study of 188 of the brightest 200 ( mJy), of the 477 Asboth et al. (2016) objects using SCUBA-2 (Holland et al., 2013) on the James Clerk Maxwell Telescope (JCMT). With the addition of the data provided by SCUBA-2 we have a better constraint on both the FIR luminosities and the redshifts of these DSFGs and prepare the way for high resolution follow-up.

With our sample of 188 galaxies observed by SCUBA-2 we roughly double the number of 500 m riser galaxies possessing longer submm wavelength data.

The format of this paper is as follows. We describe the data in Section 2. We describe our methods for determining the photometric redshifts, FIR luminosities and star formation rates (SFRs) in Section 3. The results are described in Section 4, and the discussion and conclusions in Sections 5 and 6, respectively. We use a standard flat cosmology with and .

2 Data

2.1 Selecting high-redshift dusty galaxies in HeLMS

We use the red HeLMS sample identified in Asboth et al. (2016) and below follows a short summary of their selection. The area mapped by HeLMS is a 300 deg equatorial field which is part of the HerMES project. The observations were performed using the SPIRE instrument (Griffin et al., 2010) on board the Herschel Space Observatory. Some parts of the HeLMS field were masked. Edge effects, along with a “seagull-shaped” region of strong Galactic cirrus were removed, leaving a useful area of 274 deg.

Sources were detected using a map-based search method described in Asboth et al. (2016), similar to what was used in Dowell et al. (2014), instead of sources from the HerMES catalogue derived directly from the 250 m map (Clarke et al. in prep.). For a description of how the sources were selected and the exact spatial filters adopted we refer the reader to Asboth et al. (2016), but we give a brief description here for completeness.

The SPIRE 250, 350, 500 m maps are created with the same pixel size (6 arcsec) and (for source detection only) smoothed to the same resolution using an optimal filter for easy comparison between wavebands (Chapin et al., 2011). The local background is removed by smoothing the maps with a 2D median boxcar filter on 3 arcmin scales to remove any cirrus contamination. The filters are also applied to the error map to find the typical instrumental noise in the smoothed map. The instrumental noise values are 7.56, 6.33 and 7.77 mJy, in the 250, 350, and 500 m SPIRE bands.

The confusion noise () in the SPIRE map is caused by sources which emit at all three SPIRE wavelengths. This causes the confusion noise to be correlated between wavelengths. This information is used to construct a difference map (D) from the SPIRE 500 m () and SPIRE 250 m () maps with a reduced confusion limit (Dowell et al., 2014);

 D=√1−k2M500−kM250 (1)

with a k value of 0.392 to maximize the D/. This D-map has a confusion noise of 3.50 mJy, which is much lower than in the three smoothed SPIRE bands (13.66, 11.21, 6.98 mJy at 250, 350 and 500 m, respectively).

The bright peaks in the D-map are selected with a 4 cut-off at 34 mJy. At these positions the SPIRE flux densities are determined from the (higher resolution) nominal resolution map while taking into account the positional uncertainty of 6 arcsec (as measured with simulations in Asboth et al. 2016). From these flux densities a catalogue of sources is created. There is no requirement for a detection in both 250 and 350 m, in order to avoid biassing the selection against the reddest objects.

The smoothed and raw images are compared with each other within a 30 30 arcsec region around each source to find cosmic rays. All candidate sources with are removed. The final catalogue is selected to have mJy in order to minimize the effect of faint cosmic rays which are not found by the described technique. All 17 sources with radio fluxes in excess of 1 mJy are removed using the the 21-cm radio catalogues from the NRAO VLA Sky Survey (NVSS, Becker et al., 1995) and the Radio Sky at Twenty-cm (FIRST) survey (Condon et al., 1998) to avoid contamination by flat spectrum quasars at z 1. The rejection of NVSS/FIRST sources means that we potentially miss some genuine red sources that are lensed by radio-loud galaxies (Haas et al., 2014; Leung & Riechers, 2016). The final Asboth et al. (2016) catalogue contains a total of 477 sources.

2.2 Scuba-2

We selected the 200 brightest galaxies i.e. 63 mJy, of the 477 Asboth et al. (2016) sources, and we observed a random sub-set of 188 of them for 15 minutes each using the DAISY pattern with the SCUBA-2 camera at the JCMT (Holland et al., 2013). The observations were taken in semester 15B between 31-7-2015 and 15-11-2015 with an opacity at 225GHz between 0.05 and 0.12.

Our integration times were based on the previous observations of 28 red objects from Dowell et al. (2014) with almost identical selection criteria as our sample. Those observations were 12.5 minute DAISY observations and 27 out of the 28 where detected. Using the colour distribution from these data to simulate the 850 m fluxes of the HeLMS sample we estimated that a 1 RMS = 4.5mJy would detect 70 per cent of our targets at 3.

We explored several data reduction methods including the data reduction used for the SCUBA-2 Cosmology Legacy Survey (Geach et al., 2017, S2CLS), and the quick pipeline reduction using REDUCE SCAN FAINT POINT SOURCES. We found that the “zero-mask” (Holland et al., 2017) data reduction used in Ivison et al. (2016) provided us with the highest signal-to-noise values and a RMS ranging between 3.2 mJy and 6.4 mJy with a mean of 4.3 mJy where the S2CLS method reaches an average RMS of 4.9 mJy. The flux densities obtained with the zero-mask method are on average 2.64.0 mJy higher than the S2CLS method. We decided to use the zero-mask data reduction technique for all our observations because of its effectiveness in suppressing large scale noise (Ivison et al., 2016; Holland et al., 2017).

The zero-mask data reduction uses the Dynamic Iterative Map Maker within the SMURF package (Chapin et al., 2013). This algorithm assumes that the image is free of significant emission except for a 60 arcsec diameter region centred on our target. Since the positions of our targets are in the centres of our DAISY observations this algorithm is very effective in suppressing large scale noise. This has an advantage over the S2CLS pipeline (Geach et al., 2017), which can make no prior assumptions about the positions of the targets. The maps are generated with 1 arcsec 1 arcsec pixels.

We use the same data reduction technique for the SCUBA-2 flux calibrators to get accurate flux conversion factors (FCF). These FCFs, ranging between 658 and 777 Jy pW beam, are used to convert our reduced image to units of Jy/beam. The FCFs are expected to be accurate to within 5 per cent (Dempsey et al., 2013).

Our prior positions are derived from the Herschel data and have a typical positional uncertainty () of 6 arcsec (Asboth et al., 2016). Another positional uncertainty arises from the JCMT 2-3 arcsec rms pointing accuracy (). We combine both uncertainties to obtain the final positional uncertainty ():

 σp=√σ2H+σ2J. (2)

We apply our source extraction by taking the flux density of the brightest pixel within a 20 arcsec radius of our prior position in the beam convolved image. This 20 arcsec radius corresponds roughly to the 3 positional uncertainty of our prior source in the SCUBA-2 map. We obtained an average noise level of mJy for our point source extraction.

For the purpose of analysis we divide our sample into three sub-groups with fairly arbitrary signal-to-noise ratio boundaries. Group 1 contains objects that have a clear detection, . Group 2 consists of detections between . Finally, Group 3 are galaxies for which we do not have a clear detection, . (Due to the large uncertainty in position we are unable to obtain a significant detection in the stacked signal for the Group 3 galaxies.) The three groups contain 64, 99 and 25 objects respectively.

As we are considering SCUBA-2 measurements of Herschel detected galaxies we are concerned about the accuracy of the flux measurement, rather than the reality of a catalogued source (as we would be with a blank field survey). Nevertheless we would expect random noise fluctuations and confusion noise from galaxies not associated with our original target. Furthermore we are using the brightest pixel, so our flux measurements are biased high (Coppin et al., 2008). We quantify this bias using the simulation shown in Figure 1. This simulation takes all deep S2CLS fields as the “truth”.

We add noise to the S2CLS maps by adding extra Gaussian noise to reach a total noise of  mJy, similar to those of our observations, we call this new maps the noise-added map. We then add positional errors to the S2CLS catalogue with a mean of zero and a standard deviation of 7 arcsec to the S2CLS positions to simulate the positional uncertainty of our DSFGs. We then apply our photometric measurement at the original S2CLS position and compare with the original S2CLS flux. We repeat this process 5 times to get the results from different random noise simulations.

The comparison shows that for sources with below 13 mJy (3) we are (on average) overestimating the flux density, but this overestimation is on average lower than 4.3 mJy (1). We also tested the sources extraction method of picking the nearest pixel to our prior positions and find that this method underestimates the flux density significantly for sources with 13 mJy. We decided to use our brightest pixel sources extraction because we expect that a significant percentage of our sources will lie above 13 mJy given that 63 mJy.

2.3 Ancillary data

It is unlikely that our high-redshift galaxy sample will be directly detected in any shallow large-field surveys at optical/NIR wavelengths which are not likely to contain galaxies without an AGN (Section 4.2). However, low-redshift galaxies can significantly magnify a higher redshift source behind them via gravitational lensing.

Therefore it is possible to identify a lens using the available low-redshift galaxies from the Wide-field Infrared Survey Explorer (Wright et al., 2010, WISE) and the Sloan Digital Sky Survey (York et al., 2000, SDSS). We examined the SDSS images for possible contamination from large extended nearby galaxies and we found none. However, we do find several SDSS galaxies nearby and within the FWHM area of the SPIRE beam. Due to the large SPIRE/SCUBA-2 beam it will not be possible to unambiguously identify which of the several galaxies within the beam is potentially lensing the DSFG or is the optical/NIR counterpart of the DSFG.

For all our sources (excluding HELMS_RED_80 and HELMS_RED_421, see AGN Section 4.2) we find a total of 400 WISE detected sources (Cutri et al., 2013) within a 20 arcsec radius. Of those sources only one is detected () in WISE-4 and this source is located 19.8 arcsec away from the SPIRE detection, additionally we find four WISE-3 detections () near other sources which are all located 11.2 arcsec away from the SPIRE detection. For the numerous detections in the WISE-1 band it is not clear if the WISE source is a random aligned nearby galaxy, associated with our source, is an AGN or is lensing the background DSFG. We therefore did not use WISE data in our SED fit. We can, however, study the statistical excess of galaxies nearby to our sources (Wang et al., 2011), where we only use WISE-1 sources as all but two WISE-2 galaxies are detected in WISE-1. We use SDSS DR9 (Ahn et al., 2012) and the Cutri et al. (2013) WISE catalogue to select all detected galaxies near the line-of-sight of our targets (see Section 5.1).

Strong gravitational lensing, with a lensing magnification factor () larger than 2, could provide an explanation for our high flux densities. Wide field Herschel surveys show that galaxies with a flux density at 100 mJy are likely to be strongly gravitationally lensed (Negrello et al., 2010; Conley et al., 2011; Nayyeri et al., 2016; Negrello et al., 2017). This 100 mJy limit comes from the steep slope in the FIR luminosity function, which causes intrinsically luminous ( 100 mJy) sources to be extremely rare. Our sample of 500 m riser galaxies contains 9 galaxies with 100 mJy, of which we expect 80 per cent to be strongly lensed (Negrello et al., 2010; Wardlow et al., 2013). The probability that a DSFG is strongly lensed declines for 100 mJy, but for galaxies around 70 mJy at there is still a significant (20 per cent) chance that they are lensed (Bussmann et al., 2015; Nayyeri et al., 2016).

Other follow-up programs have observed part of our sample:

• Four of the sources (HELMS_RED_3, HELMS_RED_4, HELMS_RED_6 and HELMS_RED_7) were observed at the CSO using MUSIC (Sayers et al., 2014) at four wavelengths, 2.09, 1.4, 1.1 and 0.92 mm. The resulting flux densities can be found in section 6.2 and Table 4 of Asboth et al. (2016).

• Two sources (HELMS_RED_4, HELMS_RED_31) have spectroscopic follow up with the Atacama Large Millimeter Array (ALMA). The resulting spectra can be found in Asboth et al. (2016). The redshift of HELMS_RED_4 is 5.162 and the redshift of HELMS_RED_31 is 3.798 or 4.997 depending on the line detection being the CO(5-4) or the CO(4-3) line.

• Two sources (HELMS_RED_1, HELMS_RED_2) have spectroscopic follow up by the Combined Array for Research in Millimeter-wave Astronomy (CARMA). The detected redshifts are 4.163 and 4.373, respectively (Riechers et al. in prep., Leung et al. in prep).

• Five sources (HELMS_RED_1, 2, 4, 10, 13) have been observed with the Submillimeter Array (SMA), and will be discussed in detail in Greenslade et al., in prep.

• Two sources (HELMS_RED_1, 3) are detected in the Atacama Cosmology Telescope (ACT) equatorial survey (Su et al., 2017). The measured flux densities at 148, 218 and 278 GHz are 12.491.74, 35.112.62 and 72.326.26 mJy for HELMS_RED_1 and 6.141.76, 19.502.56 and 35.326.24 mJy for HELMS_RED_3.

MUSIC and ACT provide even more data points in the Rayleigh-Jeans part of the spectrum. These additional long wavelength data will improve our SED-fitting process. The spectroscopic redshifts from CARMA and ALMA will be used to help validate our SED-fitting process and to confirm that our selection process does indeed pre-select high-redshift galaxies. We use the preliminary SMA results to get accurate information about the source positions and to determine if any sources are blended.

3 Modeling the DSFGs

3.1 SED fitting for photometric redshifts

Fits to the FIR/submm spectral energy distributions (SED) to obtain photometric redshifts and integrated properties are performed using the EAZY code (Brammer et al., 2008) using a sample of representative FIR/submm templates (e.g. Aretxaga et al., 2003).

The FIR peak of luminous infrared galaxies () can be, crudely, characterized by cool dust with average temperatures in the 25-45 K range (e.g. Soifer et al., 1984; Klaas et al., 1997). The lack of strong features means it is difficult to distinguish between either very cold dust or high-redshift galaxies using only SPIRE photometry. The addition of the data enables us to estimate the peak of the FIR emission, and therefore able to place far tighter constraints on the redshift (Section 3.3). However, since temperature and redshift are degenerate the choice of templates is a critical factor in photometric redshift estimation and so our templates have been carefully chosen to cover a broad range of temperatures.

Our six templates consist of the broad star-forming galaxy (BSFG) derived by Berta et al. (2013), cosmic Eyelash and three warm starburst galaxies M82, Arp220 (Polletta et al., 2007) and HFLS3 (Riechers et al., 2013). However, these templates have a gap at an effective temperature 37 K so we create an extra SED template from a modified black body (MBB) with a temperature of 37 K, a dust emissivity index () of 1.5 and a MIR power law component () of 2.0 (Casey, 2012). These templates are illustrated in Figure 2. With EAZY we fit all possible linear combinations of our templates set.

In Figure 3 we show the colour-colour plot of our observations. We overlay the redshift tracks from our sample of SED-templates. Our template set thus contains a wide range of representative DSFGs over a large redshift range. We can exclude very cold ( K) galaxies at as they would not be a 500 m riser. Such galaxies at higher redshift could potentially contaminate our sample. But this type of galaxies are very rare between (Symeonidis et al., 2013). Such a cold galaxy would furthermore have a higher colour than any of our measured colours at .

We only use broad-band FIR data, and we neglect the contribution of emission lines. At redshifts of FIR lines have a per cent effect at 250 , however, they have a negligible effect at 350, and 500 ; at 850 they have a per cent contribution at though this rises to per cent at (Smail et al., 2011).

We adjust EAZY to allow for 10 per cent systematic error for the data. This 10 per cent incorporates both the 5 per cent error in the FCF for SCUBA-2 and our use of a different algorithm to reduce the data for SCUBA-2 and SPIRE. The advantage of using this extra 10 per cent systematic error is that it dominates unrealistically small statistical errors for very bright () sources.

In Section 3.4 we directly compare our method with other methods, other template choices and with spectroscopic redshifts.

3.2 Noise estimates

The SPIRE and SCUBA-2 maps contain both confusion and instrumental noise. Both have to be included in the SED fitting to ensure that the errors on fitted parameters, e.g. photometric redshift are assigned the appropriate errors. The confused background in the SPIRE band is caused by coincident sources; these contribute in all three wavelength bands. The instrumental noise can be assumed to be uncorrelated and included straightforwardly in the calculations within EAZY. However, to incorporate the confusion noise we need to consider that this is correlated noise.

The confusion noise at from SCUBA-2 is significantly lower than the confusion in the SPIRE bands (1 mJy vs. 6-7 mJy; Geach et al., 2017; Nguyen et al., 2010) due to the smaller beams size of SCUBA-2 and lower number counts. The SCUBA-2 confusion noise is subdominant to the instrumental noise we obtained in the images. We can therefore safely neglect the effects of confusion in our SCUBA-2 flux density estimates. We can simulate possible values for the SPIRE contribution in the following way.

In a confusion limited map, where the instrumental noise is negligible compared with the confusion noise, the fluctuations in that map can be considered to be caused by confusion noise alone. We can randomly sample such a confusion limited map at the same position in all three bands drawing a 3-tuple of flux density values that represent the confusion noise. These samples automatically include the correlation between the bands111An alternative, would be to estimate the covariance matrix between the maps, and synthesise correlated flux density values from this assuming Gaussian fluctuations. However, by sampling directly from the map we skip this step and get a more direct model of the correlated confusion noise.

The HELMS field is not confusion limited so we sample the confusion limited COSMOS (Scoville et al., 2007) field. COSMOS has a 1 instrumental noise 2.5 mJy, though small, this residual instrumental noise means we will slightly overestimate the confusion noise values. We perturb the 3-tuple flux of each object in our catalogue by one of the sample 3-tuples drawn from COSMOS. We then run EAZY on the perturbed catalogue. We do this simulation exercise 1000 times (however, due to the finite size of the field these are not independent).

We average the redshift probability distribution function (PDF) over all simulation runs to obtain the final PDF for each galaxy. The results of 1000 runs for a single representative galaxy are shown in Figure 4. The resulting PDF is slightly broader than the PDF from the traditional method of not using the confusion noise. This effect would be larger if the noise in HeLMS had been dominated by confusion noise.

In Figure 5 we show the improvement in photometric redshift by adding the longer wavelength SCUBA-2 data. The average uncertainty (calculated from the variance of the estimated PDF from EAZY) when we only use the SPIRE flux densities is larger, , than the uncertainty with the additional SCUBA-2 data . This Figure also shows that we overestimate the photometric redshift when we only use the SPIRE data.

3.3 Physical parameters

Using EAZY we obtain the full PDF and the best fitting SED template for every galaxy. With this template we compute the total infrared luminosity, . The FIR luminosity is defined as the integral over the rest frame spectrum between 8m and 1000m, i.e. . In practice we lack a good measurement of the flux in the rest frame mid-infrared (MIR) from 8 - 30m. We therefore integrate between 30m and 1000m and use a correction factor for the potentially large amounts of missed flux in the MIR. We calculate the correction factor from the average fraction of the FIR luminosity contained in the MIR regime for 5 of our 6 templates. We exclude the HFLS3 template for this measurement due to a lack of constraints in the MIR. We obtain a correction factor of 1.17 and we multiply our measured integral by this factor to obtain the resulting . We also obtain an error on using both the errors on our flux density estimates and the scatter from our 1000 EAZY runs.

The negative K-correction (for galaxies measured at longer wavelengths than the peak of their SED) counteracts (to some extent) the dimming with distance, and so these galaxies are relatively constant in brightness (e.g. Casey et al., 2014). Therefore our estimates of can be tightly constrained even with a large uncertainty in the redshift.

Our can be translated into SFR estimates using Kennicutt (1998) for a Salpeter IMF

 SFRM⊙yr−1=1.96×10−10LIRL⊙. (3)

Here the fraction of ultraviolet energy absorbed by dust has been assumed to be , for which we have no constraint. Our estimates for the SFR would be the same if we had used the Rowan-Robinson et al. (1997) calibration factor with a . We assume no gravitational lensing (Section 5.1) and no contamination by AGN (Section 4.2) in our calculation of the SFR. The resulting SFRs should be multiplied by a factor 0.63 or 0.67 if assuming a Chabrier or Krupa IMF (Madau & Dickinson, 2014).

Our final catalogue is presented in Appendix A, where we list the positions, flux densities, redshifts and of all our galaxies observed with SCUBA-2.

3.4 Testing the photometric redshifts

Ivison et al. (2016) made a similar assumption with the selection of their templates, and tested their photometric redshift code against 25 red high-redshift DSFGs with spectroscopic redshift. Their photometric redshifts where found by finding the lowest value for their set of three templates. The main difference between our method is that EAZY not only fits the provided templates but also any linear combination of those templates. The results from Ivison et al. (2016) show only a small offset in with a scatter of 0.14.

We compare our photometric redshift method () directly with Ivison et al. (2016), by running our code on their sample. We obtain a mean () offset in of 0.11 and a median () offset of 0.12. We note that this offset is smaller than the mean estimated error in our redshift ().

The main difference between our method and that of Ivison et al. (2016) is that they tested a set of 6 templates individually with a sample of available spectroscopic redshifts, and discarded the ones with the poorest fit in . Two of the poorest fitting templates in their analysis were the Arp 220 and HFLS3, which are on the “blue” end of the range of FIR SEDs. If we discard our “blue” templates (M82, HFLS3 and Arp 220) we find that our photometric redshift estimates are very close to the Ivison et al. (2016) estimates ( = 0.024 and = 0.035). However, we keep these “blue” templates in our analysis, to ensure conservative errors, noting that EAZY produces a full redshift PDF using all our templates (and all linear combinations of them) simultaneously.

We can see how our results would change if we made a different choice of templates. Strandet et al. (2016) used a Monte Carlo method to sample a range of MBB from Greve et al. (2012) with dust temperature parameter sampled from a Gaussian with mean and standard deviation K. We use a similar full MCMC approach to fit using the FITIR module of the INTERROGATOR code (Wilkins et al. in prep). With this method we can specify prior information about all free parameters. We consider both the MBB parameterisation of Greve et al. (2012) (which has two free parameters, the temperature , and the emissivity ) and the parameterisation of Casey (2012) (which has three free parameters: the temperature, emissivity , and the slope of the near-IR power law ). For the Greve et al. (2012) parameterisation we fix the emissivity and consider 3 different priors on the temperature : fixed to K, a normal distribution centred at K with K, and a uniform prior . For the Casey (2012) we assume uniform prior on the temperature of and consider cases where both and are fixed (to and respectively) and where they have a uniform prior: and .

The results are shown in Table 1 where we compare the output of each different template set to our chosen templates when applied to our sample. We compute a number of comparison statistics, the mean offset (, the rms scatter in () and the in comparison with our three spectroscopic redshifts. For the normal distributed (T = 39 K) method we find a = 0.056 and a = 1.35, for the uniform prior () is 0.011 and the = 0.67 and for the single temperature model we find a = 54. From these results we can see that the Gaussian prior produces very similar results as our method, and that the flat 20-60 K prior models are consistent with the spectroscopic redshifts, but overestimate the size of the error bars (). The single temperature model is insufficient in fitting photometric redshifts.

The ultimate test is the comparison against spectroscopic redshifts. We obtain a good total of 3.07 for our three spectroscopic redshifts. But due to the limited number of spectroscopic redshifts in our sample we also use the SPT detected DSFGs which fulfil our colour selection criteria (Weiß et al., 2013; Strandet et al., 2016, 2017), HFLS3 and the H-atlas 500 m risers Fudamoto et al. (2017). The results are shown in Figure 6. We estimate a bias of with a rms of 0.19 and a reduced of 1.4. The rms scatter in the bias (0.19), our average uncertainty per galaxy () and all have comparable values.

There is a visible trend in Figure 6 that is decreasing with redshift, the reduced for linear decreasing model is 0.9 compared to 1.4 for the non-evolving model. This result indicates that we underestimate the redshift of high-redshift galaxies due to a rising dust temperature of our spectroscopic sample towards higher redshift (Ivison et al., 2016). However, this same result could also arise from selection effects, where a warm HFLS3 type galaxy would not have made our selection criteria at as it would not be a 500 m riser (Figure 7). Another possible explanation for this trend is that higher redshift galaxies need to be brighter to fulfil our flux density selection criteria, and these brighter galaxies tend to be warmer (e.g. Symeonidis et al., 2013; Kirkpatrick et al., 2017).

Any fitting methods with a range of temperatures and no explicit prior on the temperatures is effectively assigning a uniform prior to the temperatures. This is what our method does as do most photometric redshift fitting methods. In the low signal-to-noise regime the prior has a stronger influence on the posterior and so there will be a trend to fit mid-range temperatures rather than high or low temperatures. This naturally tempers the extremes of redshifts distributions based on the best redshift. However, the redshift PDFs are a reasonable representation of the information available.

4 Results

4.1 Statistical properties

In Figure 7 we show the SFR vs redshift distribution of our sources. Our sources have a median redshift of 3.6 and a median SFR (uncorrected for flux boosting or the possible presence of gravitational lensing) of Myr. All our galaxies could be classified as distant hyper-luminous infrared galaxies (HyLIRGS), i.e. with exceeding L and a mean of 2.7 L.

We find that 31.4 4.7 per cent lie between redshifts of 4 and 6. This finding is consistent with Ivison et al. (2016) who found 33 6 per cent of their sample to lie within this redshift range. The inferred space density in this redshift range is Mpc. Due to the predicted short lifetime for the starburst () phase we need to apply a duty-cycle correction to the observed space density to infer the actual underlying space density () for these type of galaxies

 ρ=tobs/tburst×ρobs, (4)

were is the time between . For we assume 100 Myr, which is in agreement with their expected gas depletion times (Ivison et al., 2011; Bothwell et al., 2013). The final inferred space density estimate is thus Mpc. The assumption of 100 Myr is the same as used by Ivison et al. (2016) and while longer timescales (0.5-1.0 Gyr) have been postulated (e.g. Lapi et al., 2014; Aversa et al., 2015) these would result in an even lower space-density.

The primary difference between the Ivison et al. (2016) sample and our sample is that Ivison et al. (2016) used a mJy selection where we use a mJy sample. Therefore our sample has a space density of about a factor of 10 lower than the Ivison et al. (2016) estimate of Mpc.

We use our sample to calculate the SFRD for bright 500 m risers in the SPIRE bands as shown in Figure 7. The contribution to the overall SFRD is below 1 per cent at any redshift. For comparison we also show the SFRD results from the S2CLS mJy selected sources, which is complete for galaxies with a SFR 300 M yr (Michałowski et al., 2017). The Michałowski et al. (2017) result comes from 2 deg blank fields, which observe the more common population of DSFGs and contribute more to the overall SFRD at any epoch.

4.1.1 Luminosity function

The SPIRE sources luminosity function and its evolution to z has been reported in Gruppioni et al. (2013). We can use this luminosity function as a basis to predict the number of galaxies we expect in our sample. To get an accurate estimate for our incompleteness we need to know the relative distribution of different galaxy types at these high luminosities and redshifts. The intrinsic colours of different galaxy types can be used to determine whether or not they fulfil our selection criteria as a function of redshift.

Due to the lack of information on the distribution of galaxy types at high-redshift we have to extrapolate what we know about the distribution of SED shapes at lower redshift and luminosity to the redshifts and luminosities of our sample. We do this using the results from Symeonidis et al. (2013), who measured the correlation between average dust temperatures and infrared luminosities. They characterised the rising dust temperature with luminosity for a sample of galaxies, and , to provide a simple phenomenological characterisation of this, we apply a linear fit in temperature vs. to predict the average temperature for galaxies. We also use the average value for the variance in the temperature for galaxies.

Using this temperature-luminosity-redshift distribution we draw 200 galaxies at every redshift between 1.5 and 8 () and luminosities between () and then each galaxy is assigned a temperature drawn from a Gaussian with mean from the temperature-luminosity and a sigma of 6 K. This produces a mock catalogue of 325,000 galaxies, for which we have mock , and values. We use the Casey (2012) MBB to calculate the expected flux densities at SPIRE and SCUBA-2 wavelengths for each galaxy. The upper limit of is used for practical reasons to simplify the drawing of a random luminosity. It was not intended to indicate a realistic physical limit. However, the number density is dropping off very steeply at high luminosity so exactly where this cut is made makes little difference to the outcome.

We add Gaussian noise with a mean of zero and a sigma of the mean instrumental error of our observations to simulate the variations caused by instrumental noise. On top of the Gaussian noise we also draw a correlated confusion noise estimate for every source using the COSMOS map (see Section 3.2), and we add this correlated confusion to our mock observed flux density estimates. Our novel way of adding the correlated confusion noise is crucial as it partly conserves the colour of the source. The standard deviation of the confusion noise we added is 6.7, 7.1 and 6.8 mJy at 250, 350 and 500 m, respectively and together with the instrumental noise of order 7 mJy this leads to 1 fluctuations of 10 mJy. It will therefore not be uncommon that sources of order 30 mJy at 500 m will be boosted to the selection criteria of 63 mJy due to the noise and the steepness of luminosity function.

We multiply the fraction of mock galaxies in every luminosity and redshift bin which fulfil our selection criteria by the expected space density for such galaxies (Gruppioni et al., 2013) to obtain the number of galaxies we would expect in the HeLMS field. This results in a total sample of galaxies in our mock catalogue over an area of 274 square degrees. This is mildly larger than, but consistent with, the 200 galaxies we observed in the HeLMS field. The error bars are based on the large error on the normalisation of the luminosity function (Gruppioni et al., 2013). We do acknowledge that the consistency is partly due to the large error bars in this normalisation.

We make an additional 10 mock catalogues where we modify the mean temperature in the relations of Symeonidis et al. (2013) to measure the effect of the average temperature of DSFGs on the observed number counts. In Table 2 we show the total number counts as function of (mean) temperature. It is clear that the number of observed galaxies is a strong function of temperature and it is therefore important to get a better understanding of the distribution of galaxy types at high redshift to fully understand the number counts.

In Figure 9 we show the resulting and number counts for our mock catalogues show in Table 2. Our mock catalogue is consistent at but over predicts the number of bright sources at , even when we raise the temperature of our mock catalogues with 5 K we keep over predicting the number of sources at 50 mJy.

We use our mock model as input for EAZY to predict the observed luminosity function using our method. On top of the 200 galaxies we have already drawn at every redshift and luminosity bin we draw an additional 100 galaxies for every very bright bin (), an additional 300 galaxies for the bins and an additional 500 galaxies for the bins. these extra galaxies lead to a total mock size to test the luminosity function of 630,500 galaxies. These extra galaxies give us extra statistics on the lower end of the luminosity function, where galaxies are intrinsically not bright enough to be detected with our detection method but might be very occasionally scattered up by noise. In Figure 10 we compare the predicted luminosity with the calculated luminosities for our galaxies.

From Figure 10 we can see that the simulated galaxies are scattered up in luminosity due to confusion and instrumental noise. This is a flux boosting effect, well-known in sub-mm surveys (e.g. Coppin et al., 2005, 2006). From our mock catalogue we derive that 61 per cent of the mock galaxies which observational properties fulfil our selection criteria are intrinsically not bright enough and are scattered up due to confusion and instrumental noise. We use the average boosting factor (difference between input and output Luminosity of our Mock) to correct our SRFD in Figure 8.

4.1.2 Comparison with simulations

The Simulated Infrared Dusty Extragalactic Sky (SIDES, Béthermin et al., 2017) includes a 274 deg simulation to match the size of the HeLMS field. The size of the model and its capability to simulate the observed FIR and submillimetre flux densities makes it ideal for comparison with our observations.

The main SIDES model predicts the FIR and submillimetre emission in a 2 deg light cone, which simulates clustering by using abundance matching to populate dark matter haloes with galaxies according to their star formation evolution model. This model is accurate in describing the number counts at 350 and 500 m. This 2 deg light cone is not a large enough volume to get accurate predictions for our rare sources. Béthermin et al. (2017) tackled this problem by producing the 274 deg simulation to predict number counts for much rarer (brighter) sources but this larger simulation does not contain any clustering estimates.

The number of sources in the 274 deg SIDES model which fulfil the Asboth et al. (2016) criteria is 22, and all are strongly lensed. This number goes down to 11 in the case we use our 63 mJy cut on top of the Asboth et al. (2016) criteria. These numbers are an order of magnitude lower that the bright red sources found in the HeLMS field.

Those results do not account for the effect of flux boosting by both instrumental and confusion noise. Béthermin et al. (2017) calculated this effect of flux boosting by adding random (Gaussian) instrumental and confusion noise to the fluxes. This increased the number count to 114 sources which fulfil the Asboth et al. (2016) criteria and 35 sources when we add 63 mJy constraint. The 2 deg SIDES model was used to calculate the effect of clustering on these number counts. They found that the confusion which arises from clustering increases the number of red sources by a factor of 1.7. This leaves them with an estimate of 229 sources which is within of the 477 sources found in Asboth et al. (2016). This boosting factor of 1.7 is however not high enough to boost the 35 sources in the 274 deg SIDES model to the 200 sources found in the HeLMS field.

Our method of drawing correlated confusion noise estimates provides us with a different way of using the 274 deg SIDES model to predict the number of sources in the HeLMS field. We do this by adding both random Gaussian instrumental noise, and our correlated confusion noise estimates to the SIDES 274 deg catalogue. This noise increases the number of sources from 11 to 172 (where the noise only accounts for different sets of random numbers and Poisson noise, and does not account for any other uncertainties in the SIDES model), which is very close to 200 sources which were detected with our selection criteria (see Table 2 and Figure 9).

17 per cent of these 172 sources are strongly lensed and the mean redshift is . Figure 11 shows the full redshift distribution of our data compared with the SIDES model and our mock catalogue.

From Figure 11 we can see that the redshift distribution of the mock has a larger tail to higher redshifts than our observations. We test if there is any significant net bias we calculate the mean of the observed mock and input mock redshifts, we calculate the error on this mean using jack-knife samples. We find a different value for the mean redshift (4.17 0.04 q.v. 3.69 0.08), which is smaller than the RMS of the refshifts of 0.6, but nevertheless statistically significant. Flux boosting can happen at every wavelength band but because of our 500 m riser selection we are biased towards selecting galaxies which are boosted at 500 m. These selected galaxies look therefore redder than they truly are, which results in a over estimate of the redshift. This argument mainly holds for galaxies which are intrinsically not red or bright enough to fulfil our selection criteria. For all galaxies we see the same trend as in Figure 6, where our redshifts are overestimated at high redshift and under estimated at low redshifts. As we stated in more detail in Section 3.4 this trend is partly due to selection effects and due to the prior pushing us towards mean and not “extreme” redshift estimates.

The 274 deg SIDES model model has a comparable high redshift tail, but this model peaks at lower redshift, causing the mean redshift to be lower (3.1 0.07 q.v. 3.6 0.04 from our observations).

4.2 SDSS and WISE quasars

We cross-matched the 188 galaxies with the SDSS quasar catalogue (Pâris et al., 2017) and found two matches within 20 arcsec. We test the change on a random alignment with a a SDSS quasar by taking 50,000 random positions in the HeLMS field and see how many of these random positions match with a SDSS quasar within a 20 arcsec radius. The number of matches is 127, leading to a probability of 0.25 per cent that there is a random alignment within 20 arcsec. Using this statistics we would expect that there is a 38 per cent chance that at least one of our object is randomly aligned with a SDSS quasar and there is a probability of 8 per cent for at least two alignments.

HELMS_RED_80 is located 3 arcsec from SDSS_J005036.93+014449.1 which has a redshift of 3.4351. Our estimated photometric redshift is 3.65, which is within 1 agreement with the quasars spectroscopic redshift. The quasar is furthermore detected in WISE-1, WISE-2 and WISE-3. We use the intrinsic quasar SED derived in Symeonidis et al. (2016) in combination with the WISE magnitudes to calculate the AGN contribution to the FIR luminosity. This contribution is estimated at and is a factor of 3 lower that our measured Luminosity. We thus conclude that it is likely that HELMS_RED_80 is associated with SDSS_J005036.93+014449.1 and that the quasar contaminates our SFR estimate.

HELMS_RED_421 is located 12 arcsec away from SDSS_J000127.11-010603.1 which has a redshift of 1.934. Our estimated photometric redshift is 2.95, which is in 1.3 tension with the quasars spectroscopic redshift. The separation of 12 arcsec is furthermore in 2 tension with our positions. The quasar is not detected in any WISE bands, but there is a nearby (z = 0.163) SDSS galaxy 9.0 arcsec away from our SPIRE detection which is detected in all 4 WISE bands (WISE_J000127.76-010607.5). Furthermore, WISE_J000126.74-010612.2 is located 9.6 arcsec away and is detected in WISE-1 and WISE-2 and WISE_J000127.44-010626.6 is located 12.4 arcsec away and has besides a WISE-1 and WISE-2 detection a 3.3 detection in WISE-3. The location of our SPIRE source lies in the middles between those 4 WISE/SDSS sources, indicating that this source is likely contaminated by several of those galaxies. We tested the probabilistic de-blender XID+ (Hurley et al., 2016)using the default flat uniform flux prior ( as used for the HELP database) to disentangle the SPIRE flux densities over the four sources. XID+ with a uniform flux prior, assigns the flux evenly among them as they are all located at roughly the same distance from the centre of the SPIRE emission. We note XID+ can be run with more sophisticated priors, using both SED and redshift information, however this requires thorough analysis and so we leave the nature of this SPIRE detection for future work.

HELMS_RED_421 may be associated with SDSS_J000127.11-010603.1 but would be consistent with a spurious coincidence. The percentage of the FIR luminosity which is caused by the (potential) quasars is a function of the AGN luminosity (Rosario et al., 2012; Symeonidis et al., 2016; Symeonidis, 2017), which we do not know. We therefore exclude the source from our final corrected SFRD.

4.3 Sub-mm interferometry

We use the high-resolution SMA data, the ALMA and the CARMA redshifts to more closely examine the properties of the subset of galaxies possessing this information. The images and SED fits of the 6 galaxies with interferometry data are shown in Figure 12. We now discuss the sources individually below:

• HELMS_RED_1: The photometric redshift of is consistent with the spectroscopic redshift of 4.163 which is obtained with CO(4-3) and CO(5-4) line detections (Riechers et al. in prep). The 500 m flux density of 192 mJy suggests that the object is lensed (e.g. Negrello et al., 2017). This source was also detected with ACT with flux densities of , mJy and 72.326.26 mJy at 148 (2.0), 218 (1.4) and 278 (1.1) GHz (mm), respectively. Our best fit SED predicts flux densities of 7, 24 and 47 mJy at those frequencies, which are considerably lower. The SMA flux density at 1.1 mm is mJy, which is less than half that of the ACT value at 278 GHz which is observed at a similar wavelength but with a much larger beam. The predicted 1.1 mm flux density form our best fit SED is 46.6 mJy. The nearest WISE-1 or SDSS source near to the SMA position is 15.6 arcsec away. The SMA position is 3.9 arcsec away from the SCUBA-2 position.

• HELMS_RED_2: The photometric redshift of is consistent with the spectroscopic redshift of 4.373, which is obtained with CO(4-3) and CO(5-4) line detections (Riechers et al. in prep). The 500 m flux density of 132 mJy means the object is likely to be lensed. The SMA flux density is and the predicted 1.1 mm flux density form our best fit SED is 53.9 mJy. The nearest WISE-1 or SDSS object near the SMA position is 2.0 arcsec away, the location of the source is J005258.53+061317.5 and has a WISE-1 AB magnitude of . The SMA position is 2.0 arcsec away from the SCUBA-2 position.

• HELMS_RED_4: The photometric redshift of is in 1.7 tension with the spectroscopic redshift of 5.162, which is obtained with CO(5-4) and CO(6-5) line detections (Asboth et al., 2016). The 500 m flux density of 116 mJy makes the object likely to be lensed. The SMA flux density is mJy and the predicted flux density at 1.1 mm from our best fit SED is 29.8 mJy. The nearest WISE-1 or SDSS object near the SMA position is 1.0 arcsec away, the location of the source is J002220.73-015520.2 and has a WISE-1 AB magnitude of . The SMA position is 1.5 arcsec away from the SCUBA-2 position.

• HELMS_RED_10: The photometric redshift is . The SMA flux density of and the predicted 1.1 mm flux density form our best fit SED is 24.5 mJy. The nearest WISE-1 or SDSS object near the SMA position is 8.7 arcsec away. The SMA observations are not centred on the SCUBA-2 position and the brightest peak is 4.7. The SMA position is 13.4 arcsec away from the SCUBA-2 position. It is unclear if the SMA sources is the same source as our SPIRE/SCUBA-2 detection more detail of this sources will be provided in Greenslade et al., in prep.

• HELMS_RED_13: Our photometric redshift of . The SMA flux density is mJy and the predicted flux density at 1.1 mm from our best fit SED is 19 mJy. The nearest WISE-1 or SDSS object near the SMA position is 3.6 arcsec away. The SMA position is 2.9 arcsec away from the SCUBA-2 position.

• HELMS_RED_31: This object has a single line detection which might be either the CO(5-4) or the CO(4-3) transition (Asboth et al., 2016) suggesting a redshift of 3.798 or 4.997. The photometric redshift of is consistent with the lower redshift from and in a small () tension with . The nearest WISE-1 or SDSS source is 4 arcsec away from the SCUBA-2 position.

4.4 Extreme sources

We isolate a subset of potentially high-redshift extremely bright galaxies. This subset consists of galaxies which have a clear detection with SCUBA-2 () as well as a redshift PDF which has 50 per cent of its probability at . In total we find 21 galaxies fulfilling those conditions, which includes HELMS_RED_2 , 4, 10 and 31. Figure 13 shows the WISE-1, SPIRE and SCUBA-2 cut-outs of these sources, excluding the ones we already discussed in Section 4.3.

These sources might contain some of the highest redshift DSFGs ever detected. Therefore this catalogue provides a high priority sample for spectroscopic follow-up with ALMA. High-resolution follow up observations are also required for accurately determining the blending fraction (see Section 5.1) for these types of sources.

Our candidate with the highest chance of being a z 6 galaxy is HELMS_RED_69. Its redshift is estimated to be 5.19 and 19 per cent of its redshift PDF lies above a redshift of 6. Another remarkable feature of HELMS_RED_69 is that its 500 m flux density is 1.5 times higher than that of HFLS3. There is a possibility that this source has been lensed by a foreground galaxy as we find an SDSS counterpart at a distance of 3.0 arcsec.

5 Discussion

5.1 Blending and Lensing

Due to the relatively large beam of the 500 m data there is a high probability that in many cases some parts of the measured flux density comes from randomly aligned galaxies or companion galaxies of the main source (confusion, Nguyen et al., 2010).

ALMA observations (Karim et al., 2013; Hodge et al., 2013) of bright LABOCA ( 12 mJy) sources in the 0.25 degree LESS survey (Weiß et al., 2009) showed that these sources contain emission from several fainter sources with an upper limit of 9 mJy per source, in later work this fraction of sources breaking up is found to be less significant (Simpson et al., 2015). This indicates that there might be a maximum SFR for DSFGs of M yr (Chabrier IMF). Bussmann et al. (2015) found that 20 out of 29 bright SPIRE sources ( = 52-134 mJy) break down into multiple ALMA sources, and of the 9 isolated sources 5 have a magnification factor larger than 5. Simpson et al. (2015) found that 61 per cent of their sample of bright galaxies (median 0.4 mJy) consist of a blend of 2 or more sources in the ALMA maps. Their sample was selected to be representative of the bright end of the 1 degree deep 850 m S2CLS field. The brightest detection with ALMA had a flux density of 12.9 0.6 mJy and is considerably brighter than the sources observed in Karim et al. (2013). Michałowski et al. (2017) found that bright DSFGs found in SCUBA-2 blank fields (around 10 mJy) typically have a second component of about 1-2 mJy. Furthermore, they found that the bright end of the source counts is hardly affected by replacing from SCUBA-2 flux densities with those from ALMA. The survey was taken over an area of 2 deg.

The bright end of the Karim et al. (2013), Simpson et al. (2015) and Michałowski et al. (2017) sources are fainter than 20 mJy, and are thus much fainter than our Group 1 galaxies. Hence it would be interesting to see if our brightest sources are also characterised by having a second component of about 1-2 mJy or a 61 per cent blending fraction.

Prior-based source extraction (XID+ Hurley et al., 2016) to investigate multiplicities of bright Herschel sources at m in the COSMOS field show that the brightest component contributes roughly 40 per cent of the source flux density (Scudder et al., 2016).

The multiplicity due to blending seen in these studies is a potential concern. Blending of objects at the same redshift will not seriously impact on the redshift determination, although we will determine the luminosity and star formation of the combined system, rather than a single object. Blending of two (or more) objects at different redshifts will produce composite SEDs which are likely to elicit an intermediate redshift estimate. We derive from our mock observations that 60 per cent of our detected galaxies are likely to be scattered up to our selection criteria due to flux boosting partly caused by blending with foreground objects. Some of those boosting factors are as large as 0.5 dex, but can be explained by instrumental and confusion noise. An example of such a large effect might be HELMS_RED_421 where the SPIRE position is in the middle of three WISE sources and a SDSS detect quasar.

The advantage with our sample is that we probed a much wider field, over 100 times wider than COSMOS and S2CLS and more than 1000 times bigger than the area targeted by the ALMA observations of the LESS field. Our sample is therefore expected to be comprised of a much rarer and more luminous and less confused population of sources. However, a proper investigation of the blending of these objects is deferred until we are able to obtain high-resolution data of a significant sub-set.

HFLS3 has an observed flux density of 35.4 mJy at 850 m (Robson et al., 2014), which is comparable to our Group 1 galaxies. Riechers et al. (2013) found that HFLS3 is only marginally magnified by a factor of 1.2-1.5 by a foreground lens. This magnification factor was updated by Cooray et al. (2014) to a factor of 2.2 0.3, which yields a SFR of 1320 yr.

Another explanation for the high SFRs in our sample is that the galaxies might be lensed. We do not expect to detect weak lensing from high-redshift lenses or even unambiguous confirmation of large magnifications from nearby lenses from our current data. We can however asses the likely incidence of lensing statistically by looking at the density of WISE-1 and SDSS sources near our SCUBA-2 detections. For this we use our Group 1 galaxies which have SCUBA-2 detections. For this Group we can use the SCUBA-2 positions with a statistical positional accuracy of (Ivison et al., 2007) which is of the order of 2 arcsec, which is comparable with the JCMT pointing accuracy of 2-3 arcsec. We combine these 2 uncertainties () and assume that either the optical/NIR counterpart or the lens should lie within a 7 arcsec annulus around our source position.

We use this 7 arcsec aperture to count the number of WISE-1 detected objects near our group 1 sources and find that 53 per cent have a nearby WISE galaxy. We calculate the significance of this number by using the same aperture at 64 (same number as Group 1 galaxies) random positions in the HeLMS field 1000 times. With these 1000 runs we can calculate both the expectation value and the 84.1, 97.8 and 99.9 percentiles.

In Figure 14 we show our results for our 7 arcsec aperture and several larger apertures. It is clear that there is a significant overdensity of WISE-1 sources near our Group 1 galaxies. The total space density of WISE-1 sources is a factor of 3 higher, which is a strong indication that part of this sample is lensed (Wang et al., 2011). We perform the same measurement on a 100 mJy subset and find that 7 of the 9 sources have a WISE-1 counterpart within the 7 arcsec annulus. This leads to an even higher WISE-1 space density compared to our Group 1 sample, but due to the low number of galaxies this is less significant than our Group 1 overdensity. González-Nuevo et al. (2014) found a similar result by cross-correlating SDSS and GAMA galaxies with mJy H-ATLAS sources. They found a spatial correlation, which for non-overlapping redshift distributions can be explained by weak gravitation lensing, where the weak regime has a lensing factor smaller than 2.

We do, however, wish to stress that this is a statistical measurement, and we lack accurate enough positions and morphologies for our DSFGs to do a proper lensing analysis to yield a magnification estimate. Furthermore, we lack sufficiently deep optival data to find potential lenses at z 0.5. For our estimates of the and SFR we do not take this significant foreground source detection into account.

An overlap between the redshift distribution of WISE-1 sources and 500 m risers could provide another explanation for the excess WISE-1 sources seen near our high-redshift DSFGs. Ménard et al. (2013) used the cross-correlation between SDSS quasars and WISE sources to recover the redshift distribution of WISE sources. They found that there is a potential sub-sample of WISE sources with a redshifts 2, indicating that it is possible that our high-redshift sample could be detected in WISE. However, their red WISE sources with z 2 are at least 1.2 magnitudes brighter in the WISE-2 band than in WISE-1 and several orders of magnitude brighter in the WISE-3 band.

We find that 34 of our 64 Group 1 galaxies have at least one WISE-1 source within the 7 arcsec aperture. Of these 33 closest WISE sources, only 2 have a WISE-2 magnitude 1.2 brighter than WISE-1. The remaining 31 are either undetected in WISE-2 (and we would have detected them if they were 1.2 magnitude brighter than in WISE-1) or they have a WISE-2 magnitude that is not bright enough to fall in the red sample. The non-red WISE sources have a mean redshift of 0.5 and are very unlikely to have a redshift above 1.5 Ménard et al. (2013).

High-resolution follow-up is required to properly assess the incidence of lensing and to resolve any blending issues and determine the merging, interacting or stable disk-like morphologies of the systems (e.g. Dye et al., 2017; Oteo et al., 2017, Leung et al. in prep).

5.2 Space density

Ivison et al. (2016) compared their space density of DSFGs at to the space density of UVJ selected massive galaxies at 3.4 4.2 (Straatman et al., 2014) which are predicted to form their stellar mass around redshift 5. Straatman et al. (2014) found a space density for massive quiescent galaxies of Mpc, which is an order of magnitude higher than Ivison et al. (2016) and 2 orders of magnitude higher than our results of Mpc.

We can therefore make a similar conclusion to that of Ivison et al. (2016) for H-ATLAS ( mJy), i.e. that the HeLMS ( mJy) red sample cannot account for the massive quiescent galaxies found at . This can be confirmed by our measured SFRs (5000 yr), which for a of 100 Myr generate a higher stellar mass of yr than the (Straatman et al., 2014) sample. This suggests that part of our sample might go through a short phase (100 Myr) of extremely high star formation, or is lensed, or is more massive (and thus rarer) than the population probed by Straatman et al. (2014).

6 Conclusions

We have observed 188 high-redshift, dusty, star-forming galaxy candidates with the SCUBA-2 camera at the JCMT. The sample had been selected to be very red and bright at Herschel SPIRE wavelengths and was taken from the 274 deg HerMES HeLMS field. We achieve a 1 rms depth of 4.3 mJy and detected 87 per cent of our candidates with S/N 3.

We developed a new method of incorporating correlated confusion noise into our SED fitting procedure.

We applied EAZY with a range of galaxy templates to determine the full redshift PDF. The addition of the longer wavelength 850m data improves our photometric redshifts which are systematically lower than with SPIRE data alone and reduces the estimated uncertainties from to 0.19. Our photometric redshifts are consistent with the four spectroscopic redshifts available in our sample.

With this final PDF we compute the redshift, FIR luminosity and SFRs of our sample. From these we computed the SFRD and showed that the population of 500 m risers with 63 mJy contribute less than 1 per cent of the total SFRD at any epoch. The number density of 500 m risers is consistent with a model extrapolated from the Gruppioni et al. (2013) FIR luminosity function and with the Béthermin et al. (2017) empirical model, contradicting previous tensions with physically motivated models. consistency with the models arises from our novel way in adding both confusion and instrumental noise were 60 per cent of the galaxies are predicted to be scattered up to our selection criteria due to flux boosting.

The excess number of WISE-1 sources near our DSFGs gives a strong indication that some of our sample may be lensed or that there is a surprising overlap with the redshift distribution of WISE-1 sources.

High resolution SMA observations of 5 of our sources reveal that two of the sources have a WISE-1 source at the same position. In all cases the SMA flux density at 1.1 mm is lower than predicted from our best fit SED.

We identify a subset of 21 excellent very high-redshift DSFGs candidates, of which two are already identified as z4 DSFGs. This group is clearly detected by SCUBA-2 with a high probability that they lie above a redshift of 4. These 21 galaxies would be ideal targets for interferometric imaging and spectroscopy to get a better understanding of these high-redshift objects with extreme ( Myr) star formation.

Observing the high-redshift, dust-obscured Universe remains an important challenge for current day astronomy. Interferometry has the resolution and sensitivity to answer our questions about the SFR and nature of these obscured galaxies. With telescopes like ALMA it is still impractical to cover a large enough area of the sky to find a representative population of extreme star-forming sources. With our new catalogue we now possess an ideal target list for high-resolution and spectroscopic follow up.

7 acknowledgements

This project has received funding from the European Unionâs Horizon 2020 research and innovation programme under grant agreement No. 607254; this publication reflects only the authorâs view and the European Union is not responsible for any use that may be made of the information contained therein. S.D acknowledges support from the Science and Technology Facilities Council (grant number ST/M503836/1). S. O. and JMS acknowledges support from the Science and Technology Facilities Council (grant number ST/L000652/1). JMS acknowledges the support from the PATT travel grant. D.R. acknowledges support from the National Science Foundation under grant number AST-1614213 to Cornell University. J.L.W is supported by an STFC Ernest Rutherford fellowship, and acknowledges additional support from a European Union COFUND/Durham Junior Research Fellowship under EU grant agreement number 609412, and from STFC (ST/L00075X/1). H.D. acknowledges financial support from the Spanish Ministry of Economy and Competitiveness (MINECO) under the 2014 RamÃ³n y Cajal program MINECO RYC-2014-15686. KC acknowledges support from the Science and Technology Facilities Council (grant number ST/M001008/1). MTS was supported by a Royal Society Leverhulme Trust Senior Research Fellowship (LT150041).

The James Clerk Maxwell Telescope is 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 Observatories 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 participating universities in the United Kingdom and Canada.

The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen) responsible for the payload module and for system testing at spacecraft level, Thales Alenia Space (Turin) responsible for the service module, and Astrium (Toulouse) responsible for the telescope, with in excess of a hundred subcontractors

SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University of Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, University of Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK); and Caltech, JPL, NHSC, University of Colorado (USA). This development has been supported by national funding agencies CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA).

References

• Ahn et al. (2012) Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2012, ApJS, 203, 21
• Aretxaga et al. (2003) Aretxaga, I., Hughes, D. H., Chapin, E. L., et al. 2003, MNRAS, 342, 759
• Asboth et al. (2016) Asboth, V., Conley, A., Sayers, J., et al. 2016, MNRAS, 462, 1989
• Aversa et al. (2015) Aversa, R., Lapi, A., de Zotti, G., Shankar, F., & Danese, L. 2015, ApJ, 810, 74
• Baugh et al. (2005) Baugh, C. M., Lacey, C. G., Frenk, C. S., et al. 2005, MNRAS, 356, 1191
• Becker et al. (1995) Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559
• Berta et al. (2013) Berta, S., Lutz, D., Santini, P., et al. 2013, A&A, 551, AA100
• Béthermin et al. (2011) Béthermin, M., Dole, H., Lagache, G., Le Borgne, D., & Penin, A. 2011, A&A, 529, A4
• Béthermin et al. (2012) Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23
• Béthermin et al. (2017) Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89
• Blain & Longair (1993) Blain, A. W., & Longair, M. S. 1993, MNRAS, 264, 509
• Bothwell et al. (2013) Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047
• Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503-1513
• Burgarella et al. (2013) Burgarella, D., Buat, V., Gruppioni, C., et al. 2013, A&A, 554, A70
• Bussmann et al. (2015) Bussmann, R. S., Riechers, D., Fialkov, A., et al. 2015, ApJ, 812, 43
• Carlstrom et al. (2011) Carlstrom, J. E., Ade, P. A. R., Aird, K. A., et al. 2011, PASP, 123, 568
• Casey et al. (2012a) Casey, C. M., Berta, S., Béthermin, M., et al. 2012, ApJ, 761, 139
• Casey et al. (2012b) Casey, C. M., Berta, S., Béthermin, M., et al. 2012, ApJ, 761, 140
• Casey (2012) Casey, C. M. 2012, MNRAS, 425, 3094
• Casey et al. (2014) Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45
• Cesarsky et al. (1996) Cesarsky, C. J., Abergel, A., Agnese, P., et al. 1996, A&A, 315, L32
• Chapin et al. (2011) Chapin, E. L., Chapman, S. C., Coppin, K. E., et al. 2011, MNRAS, 411, 505
• Chapin et al. (2013) Chapin, E. L., Berry, D. S., Gibb, A. G., et al. 2013, MNRAS, 430, 2545
• Conley et al. (2011) Conley, A., Cooray, A., Vieira, J. D., et al. 2011, ApJ, 732, L35
• Cooray et al. (2014) Cooray, A., Calanog, J., Wardlow, J. L., et al. 2014, ApJ, 790, 40
• Condon et al. (1998) Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693
• Coppin et al. (2005) Coppin, K., Halpern, M., Scott, D., Borys, C., & Chapman, S. 2005, MNRAS, 357, 1022
• Coppin et al. (2006) Coppin, K., Chapin, E. L., Mortier, A. M. J., et al. 2006, MNRAS, 372, 1621
• Coppin et al. (2008) Coppin, K., Halpern, M., Scott, D., et al. 2008, MNRAS, 384, 1597
• Cox et al. (2011) Cox, P., Krips, M., Neri, R., et al. 2011, ApJ, 740, 63
• Cutri et al. (2013) Cutri, R. M., & et al. 2013, VizieR Online Data Catalog, 2328
• Dempsey et al. (2013) Dempsey, J. T., Friberg, P., Jenness, T., et al. 2013, MNRAS, 430, 2534
• Dowell et al. (2014) Dowell, C. D., Conley, A., Glenn, J., et al. 2014, ApJ, 780, 75
• Dye et al. (2017) Dye, S., Furlanetto, C., Dunne, L., et al. 2017, arXiv:1705.05413
• Eales et al. (2010) Eales, S., Dunne, L., Clements, D., et al. 2010, PASP, 122, 499
• Franceschini et al. (1991) Franceschini, A., Toffolatti, L., Mazzei, P., Danese, L., & de Zotti, G. 1991, A&AS, 89, 285
• Fudamoto et al. (2017) Fudamoto, Y., Ivison, R. J., Oteo, I., et al. 2017, MNRAS, 472, 2028
• Geach et al. (2017) Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789
• González-Nuevo et al. (2014) González-Nuevo, J., Lapi, A., Negrello, M., et al. 2014, MNRAS, 442, 2680
• Griffin et al. (2010) Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3
• Greve et al. (2012) Greve, T. R., Vieira, J. D., Weiß, A., et al. 2012, ApJ, 756, 101
• Gruppioni et al. (2013) Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23
• Haas et al. (2014) Haas, M., Leipski, C., Barthel, P., et al. 2014, ApJ, 790, 46
• Hayward et al. (2013) Hayward, C. C., Behroozi, P. S., Somerville, R. S., et al. 2013, MNRAS, 434, 2572
• Hodge et al. (2013) Hodge, J. A., Karim, A., Smail, I., et al. 2013, ApJ, 768, 91
• Holland et al. (2013) Holland, W. S., Bintley, D., Chapin, E. L., et al. 2013, MNRAS, 430, 2513
• Holland et al. (2017) Holland, W. S., Matthews, B. C., Kennedy, G. M., et al. 2017, MNRAS, 470, 3606
• Hurley et al. (2016) Hurley, P. D., Oliver, S., Betancourt, M., et al. 2016, MNRAS,
• Ivison et al. (2007) Ivison, R. J., Greve, T. R., Dunlop, J. S., et al. 2007, MNRAS, 380, 199
• Ivison et al. (2011) Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913
• Ivison et al. (2016) Ivison, R. J., Lewis, A. J. R., Weiss, A., et al. 2016, ApJ, 832, 78
• Karim et al. (2013) Karim, A., Swinbank, A. M., Hodge, J. A., et al. 2013, MNRAS, 432, 2
• Kennicutt (1998) Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189
• Kirkpatrick et al. (2017) Kirkpatrick, A., Pope, A., Sajina, A., et al. 2017, ApJ, 843, 71
• Klaas et al. (1997) Klaas, U., Haas, M., Heinrichsen, I., & Schulz, B. 1997, A&A, 325, L21
• Lacey et al. (2010) Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2010, MNRAS, 405, 2
• Lapi et al. (2014) Lapi, A., Raimundo, S., Aversa, R., et al. 2014, ApJ, 782, 69
• Leung & Riechers (2016) Leung, T. K. D., & Riechers, D. A. 2016, ApJ, 818, 196
• Lonsdale et al. (1984) Lonsdale, C. J., Persson, S. E., & Matthews, K. 1984, ApJ, 287, 95
• Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
• Ménard et al. (2013) Ménard, B., Scranton, R., Schmidt, S., et al. 2013, arXiv:1303.4722
• Michałowski et al. (2017) Michałowski, M. J., Dunlop, J. S., Koprowski, M. P., et al. 2017, MNRAS, 469, 492
• Narayanan et al. (2010) Narayanan, D., Hayward, C. C., Cox, T. J., et al. 2010, MNRAS, 401, 1613
• Nayyeri et al. (2016) Nayyeri, H., Keele, M., Cooray, A., et al. 2016, ApJ, 823, 17
• Negrello et al. (2010) Negrello, M., Hopwood, R., De Zotti, G., et al. 2010, Science, 330, 800
• Negrello et al. (2017) Negrello, M., Amber, S., Amvrosiadis, A., et al. 2017, MNRAS, 465, 3558
• Nguyen et al. (2010) Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, A&A, 518, L5
• Oliver et al. (2012) Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614
• Oteo et al. (2017) Oteo, I., Zwaan, M. A., Ivison, R. J., Smail, I., & Biggs, A. D. 2017, ApJ, 837, 182
• Pâris et al. (2017) Pâris, I., Petitjean, P., Ross, N. P., et al. 2017, A&A, 597, A79
• Pilbratt et al. (2010) Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1
• Polletta et al. (2007) Polletta, M., Tajer, M., Maraschi, L., et al. 2007, ApJ, 663, 81
• Riechers et al. (2013) Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329
• Riechers et al. (2017) Riechers, D. A., Leung, T. K. D., Ivison, R. J., et al. 2017, ApJ, 850, 1
• Robson et al. (2014) Robson, E. I., Ivison, R. J., Smail, I., et al. 2014, ApJ, 793, 11
• Rosario et al. (2012) Rosario, D. J., Santini, P., Lutz, D., et al. 2012, A&A, 545, A45
• Rowan-Robinson et al. (1997) Rowan-Robinson, M., Mann, R. G., Oliver, S. J., et al. 1997, MNRAS, 289, 490
• Sayers et al. (2014) Sayers, J., Bockstiegel, C., Brugger, S., et al. 2014, Proc. SPIE, 9153, 915304
• Scoville et al. (2007) Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1
• Scudder et al. (2016) Scudder, J. M., Oliver, S., Hurley, P. D., et al. 2016, MNRAS, 460, 1119
• Simpson et al. (2015) Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 807, 128
• Smail et al. (1997) Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJ, 490, L5
• Smail et al. (2011) Smail, I., Swinbank, A. M., Ivison, R. J., & Ibar, E. 2011, MNRAS, 414, L95
• Soifer et al. (1984) Soifer, B. T., Neugebauer, G., Helou, G., et al. 1984, ApJ, 283, L1
• Straatman et al. (2014) Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2014, ApJ, 783, L14
• Strandet et al. (2016) Strandet, M. L., Weiss, A., Vieira, J. D., et al. 2016, ApJ, 822, 80
• Strandet et al. (2017) Strandet, M. L., Weiss, A., De Breuck, C., et al. 2017, ApJ, 842, L15
• Su et al. (2017) Su, T., Marriage, T. A., Asboth, V., et al. 2017, MNRAS, 464, 968
• Symeonidis et al. (2013) Symeonidis, M., Vaccari, M., Berta, S., et al. 2013, MNRAS, 431, 2317
• Symeonidis et al. (2016) Symeonidis, M., Giblin, B. M., Page, M. J., et al. 2016, MNRAS, 459, 257
• Symeonidis (2017) Symeonidis, M. 2017, MNRAS, 465, 1401
• Vieira et al. (2013) Vieira, J. D., Marrone, D. P., Chapman, S. C., et al. 2013, Nature, 495, 344
• Wang et al. (2011) Wang, L., Cooray, A., Farrah, D., et al. 2011, MNRAS, 414, 596
• Wardlow et al. (2013) Wardlow, J. L., Cooray, A., De Bernardis, F., et al. 2013, ApJ, 762, 59
• Weiß et al. (2009) Weiß, A., Kovács, A., Coppin, K., et al. 2009, ApJ, 707, 1201
• Weiß et al. (2013) Weiß, A., De Breuck, C., Marrone, D. P., et al. 2013, ApJ, 767, 88
• Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868-1881
• York et al. (2000) York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579