© ESO 2019

1. Introduction

Over the last few decades, advances in the sensitivity of observations (specially in the IR part of the spectrum) and progress in data processing have allowed us to probe the high redshift universe in greater detail. The direct observation of galaxies in the redshift range z ∼ 1 − 10 gives us the opportunity to study the cosmic history of star and galaxy formation at different cosmic epochs (see for example de Zotti et al. 2010; Eales 2015). However, despite the constant increase in diameter of the telescopes and increase of sensitivy of the detectors, observations of the distant universe are still flux-limited, rendering only those objects that are above the detection threshold. In a universe in which the inverse-square law prevails, a flux limit implies that the highest redshift galaxies accessible to any observatory will be among its faintest detectable objects. This situation is alleviated for sources selected in the submillimetric range of the electromagnetic spectrum thanks to the strong, negative K correction, which leads to high-redshift galaxies being relatively easy to detect at submm wavelengths as compared with their low-redshift counterparts (Blain & Longair 1993). In addition, lucky alignments of background objects with foreground lenses can push the limits further by enhancing the flux of objects that could not be detected otherwise. But even with the aid of the negative K correction and gravitational lensing, signal processing techniques are a fundamental tool to reach the faintest and most distant galaxies. This is particularly true in the microwave and far infrared (IR) parts of the electromagnetic spectrum, where the fluctuations from the cosmic infrared background (CIB) create a confusion noise whose level is comparable to the flux density of the typical high redshift galaxies.

The standard single-frequency detection methods for point sources in the CMB and far IR are based on wavelet techniques (Vielva et al. 2003; Barnard et al. 2004; González-Nuevo et al. 2006) or on the matched filter (or MF hereafter, Tegmark & de Oliveira-Costa 1998; Herranz et al. 2002; Barreiro et al. 2003; López-Caniego et al. 2006, see also Herranz & Vielva 2010 for a review.). Wavelets are well suited for the detection of compact sources due to their good position-scale determination properties, whereas the MF is the optimal linear detector-estimator because it provides the maximum signal-to-noise (S/N) amplification for a source with a known shape (usually the point-spread function, or PSF hereafter, of the telescope) embedded in statistically homogeneous and spatially correlated noise. By default, these techniques are applicable only to single-frequency sky images: even for multiwavelength observatories such as the Herschel Space Observatory (Pilbratt et al. 2010) or Planck (Tauber et al. 2010), the standard detection pipelines have produced individual source catalogs for each frequency band (see e.g., Planck Collaboration VII 2011; Planck Collaboration XXVIII 2014; Planck Collaboration XXVI 2016; Maddox et al. 2018). The next logical step is to boost the signal of faint sources by combining the different bands into a single detection, that is, “multifrequency detection”. Most of the blind component separation algorithms that are used for diffuse components in microwave and far IR astronomy can not deal with the high diversity of spectral behaviors associated to the different populations of extragalactic compact sources (see for example Leach et al. 2008). However, over the last few years a number of multifrequency compact source detection techniques have been proposed in the literature (Herranz & Sanz 2008; Herranz et al. 2009; Lanz et al. 2010, 2013; Planck Collaboration Int. LIV 2018). A review on the topic can be found in Herranz et al. (2012). In particular, if the spatial profile and the spectral energy distribution (SED) of the sources are known, and if the cross-power spectrum is known, or can be estimated from the data, the optimal linear detection method is the matched multifilter (or MMF hereafter, Herranz et al. 2002). Lanz et al. (2010) also showed that the MMF can be generalized for the case where the SED of the sources is not known. This generalization outperforms the single-frequency MF in terms of S/N and can be used to infer the spectral index of synchrotron-dominated radio sources, as shown in Lanz et al. (2013). However, in this paper we will incorporate a specific SED to the MMF in order to derive a photometric redshift estimation of dusty galaxies and high-redshift star forming galaxies detected in the IR part of the spectrum¹. We will do so by applying the multifrequency MMF filter to the first and second data releases of the Herschel Astrophysical Terahertz Large Area Survey (the Herschel-ATLAS or H-ATLAS, Eales et al. 2010), the largest single key project carried out in open time with the Herschel Space Observatory. We restrict our multifrequency analysis to the three wavelength bands covered by the SPIRE instrument aboard Herschel (Griffin et al. 2010), centered around 250, 350 and 500 μm. As discussed in Hopwood et al. (2010), Lapi et al. (2011), González-Nuevo et al. (2012), Pearson et al. (2013) and Donevski et al. (2018), the SPIRE bands are ideal for capturing the peak in the SED corresponding to dust emission of star-forming galaxies at z ∼ 2, that is redshifted from its rest-frame wavelength around 70–100 μm to the SPIRE wavelengths: This is the redshift range where galaxies have formed most of their stars. At higher redshifts, dusty star-forming galaxies (DSFGs) occupy the most massive halos and are among the most luminous objects found at z ≳ 4 (Michałowski et al. 2014; Oteo et al. 2016; Ikarashi et al. 2017). These high-redshift DSFGs have markedly red colors as seen by SPIRE, with rising flux densities from 250 to 500 μm (the so-called “500 μm-risers”), and have received a great deal of attention in the recent years (see for example Ivison et al. 2016; Negrello et al. 2017; Strandet et al. 2017). The DSFGs, and particularly the 500 μm risers uncovered by Herschel, are providing much insight into the early star forming history of the universe. However, sensitivity and limited angular resolution severely constrain the power of this type of objects as astrophysical probes. The sensitivity of SPIRE allows for the direct detection of only the brightest, and thus rarest objects, at the bright end of the luminosity function. By means of our multifrequency MMF technique, we intend to enhance the detectability and statistical significance of very faint red objects in the H-ATLAS source catalog and so expand the list of reliable 500 μm-riser candidates.

Although a non-negligible part of the faint H-ATLAS sources at z >  1 could be detected thanks to having been amplified by weak lensing (González-Nuevo et al. 2014, 2017), most of the faint high-z candidates in the H-ATLAS catalog have not been strongly lensed (with magnification factors larger than a few) by foreground halos (Negrello et al. 2017). In the other end of the flux density distribution, gravitational lensing plays an important role by magnifying distant galaxies that could be otherwise below the detection threshold or, at the very least, be observed with a significantly smaller flux (Negrello et al. 2007, 2014, 2017; Hopwood et al. 2010; Cox et al. 2011; Conley et al. 2011; Lapi et al. 2011; González-Nuevo et al. 2012; Bussmann et al. 2012, 2013; Vieira et al. 2013; Wardlow et al. 2013; Canalog et al. 2014; Messias et al. 2014; Dye et al. 2015; Nayyeri et al. 2016; Spilker et al. 2016). Gravitational lensing is a powerful astrophysical and cosmological probe particularly rewarding at submillimeter wavelengths. As mentioned before, submillimeter telescopes such as Herschel have limited spatial resolution and consequently high source confusion, which makes it difficult to probe the dusty star forming galaxies. However, due to the relatively low probability of lensing (with typical magnification factors of a few), the identification of gravitational lenses is difficult and usually results in a few candidates. At high fluxes, wide-area submillimeter surveys can simply, and easily detect strong gravitational lensing events, with close to 100% efficiency, as was proved by Hopwood et al. (2010). These are often strongly lensed galaxies (SLGs) with magnification factors of order ten that can be more easily detected owing to their magnified flux. The identification of these lenses is of great interest for multiple reasons. They offer the possibility to study in greater detail distant galaxies and resolve some of their features. Also, the background galaxies can be used to reveal the internal structure of the lenses. Having a large catalog of SLGs will be important in future studies. For example, caustic crossing events on these galaxies can be used to study, not only distant luminous stars, but also the constituents of the lens itself. If a sizable fraction of dark matter is made of compact objects, caustic crossing events can be used to set limits in their fraction on a range of masses from subsolar mass to tens of solar masses (through microlensing). This mass range can be difficult to probe otherwise (see for instance Diego et al. 2018).

In this work we aim at producing a catalog of distant and faint IR sources. Also, we select from those the ones that are more likely to be gravitationally lensed. Our identification of lensed candidates takes advantage of our newly inferred redshifts. As it will be described in Sect. 2 of this paper, the use of a parametrized SED template in the MMF technique allows us to provide a photometric estimation of the redshift of sources, that facilitates the identification of possible lensing systems (lens plus background source). Some of these systems will be confirmed (either as lens systems, or random alignments) in the near future with ground observations. Also, we should notice that our search for lensed systems was restricted to galaxies with known spectroscopic redshift from the SDSS. Extending the number of galaxies with redshift by using other data sets, like for instance with the data from the GAMA fields, would extend the number of candidates to lens systems not only contained in our various samples of high-z Herschel candidates.

The structure of this paper is as follows. Section 2 describes the theoretical framework of the MMF technique and the frequency dependence model used to estimate the redshifts and fluxes of the SPIRE sources. In Sect. 3 we describe the images and data of the H-ATLAS survey on which we have applied the method. In Sect. 4 we expose the results of testing the method with simulations. In Sect. 5 we compare our MMF-derived flux densities and photometric redshifts with the H-ATLAS fluxes and spectroscopic redshifts of previously studied objects. Selection of high-z Herschel sources is explained in Sect. 6 and the search for possible lensed Herschel sources is reported in Sect. 7. Finally, conclusions are detailed in Sect. 8.

2. Method

The MMF is the optimal linear detection method when the frequency dependence and the spatial profile of the sources are known, and the cross-power spectrum of the noise is known or can be estimated from the data. In the Fourier space the MMF can be written as follows:

$$\begin{array}{cc}& \mathbf{\Psi }(q)=\alpha {\mathbf{P}}^{-1}\mathbf{F},\hfill \\ \hfill & {\alpha }^{-1}\equiv {\displaystyle \int \mathrm{d}q{\mathbf{F}}^t{\mathbf{P}}^{-1}\mathbf{F},}\hfill \\ \hfill & {\sigma }^2=\int \mathrm{d}q{\mathbf{\Psi }}^t\mathbf{P}\mathbf{\Psi }=\alpha \hfill \end{array}$$(1)

where Ψ(q) is the column vector of the filters Ψ(q)=[ψ_ν(q)], F is the column vector F = [f_ντ_ν], being f_ν the frequency dependence and τ_ν the source profile at each frequency ν, P⁻¹ is the inverse matrix of the cross-power spectrum P and σ² is the variance of the output-filtered image. In Eq. (1) and in the following discussion, q≡|q| is the modulus of the Fourier wave vector; since we are assuming circularly symmetric source profiles, and since the cross-power spectrum only depends on the modulus q, all the formulas can be expressed in terms of q instead of the full vector. However, it would be easy to generalize our formulas for non symmetric profiles just by replacing q by q in the equations. Finally, α in Eq. (1) can be interpreted as the normalization that is requested in order to guarantee that the filters Ψ are unbiased estimators of the flux density of the sources under study. Further details can be found in Herranz et al. (2002); Lanz et al. (2010, 2013).

Rewriting the vector F = [f_ντ_ν] in the matrix form F = T(q)f(ν), with diagonal matrix T(q) = diag [τ₁(q), ..., τ_N(q)] and f = [f_ν] the vector of frequency dependence, we are able to include all the dependence in q in the matrix T pulling it completely apart from the dependence in ν. This way Eq. (1) can be rewritten as:

$$\begin{array}{c}\hfill \mathbf{\Psi }(q)=\alpha {\mathbf{P}}^{-1}\mathbf{T}\mathit{f},{\alpha }^{-1}={\displaystyle \int \mathrm{d}q{\mathit{f}}^t\mathbf{T}{\mathbf{P}}^{-1}\mathbf{T}\mathit{f}={\mathit{f}}^t\mathbf{H}\mathit{f}}\end{array}$$(2)

where matrix H = ∫dqTP⁻¹T and we used the facts that T^t = T and that vector f does not depend on q.

This reformulation of the Eq. (1) is very convenient for implementation of the MMF. The most time-consuming part of the filtering is the calculation of the matrices P and T since they must be calculated for all values of q. In the case we are considering in this paper the only quantity that varies during the maximization process is the redshift of the source we want to estimate. This allows us to compute the integrals of matrix H only once for each set of images of the source considered.

The MMF takes as argument a set of N images corresponding to the same area of the sky observed simultaneously at N different frequencies and returns a single filtered image where the source is optimally enhanced with respect to the noise. For N images, the frequency dependence f_ν has N degrees of freedom. Choosing one of the frequencies under consideration as fiducial frequency of reference allows to reduce to N − 1 the number of independent degrees of freedom.

The total filtered map is the result of a two-phase process. The first phase is the slowest one but, having separated the dependence in q from the dependence in ν in Eq. (2), it only needs to be done once for each set of images of the source considered. It consists on the calculation of a prefiltered map without any frequency dependence information, and for what is necessary to have previously calculated the Fourier transforms of the N images and the filters without frequency dependence. The second phase is faster and requires only the calculations of the normalization α and the linear combination of prefiltered maps using a given frequency dependence f_ν. The two necessary requirements to guarantee that the filtered field is optimal for the detection of point sources are that the filtered map is an unbiased estimator, on average, of the amplitude of the source (unbiased filter) and that the variance of the filtered map around that value is as small as possible, that is, that it is an efficient estimator of the amplitude of the source (maximum efficiency filter).

Summing up, in the first step, each individual frequency image is filtered with a linear filter, and in the second step all the resulting filtered maps are combined so that the signal is boosted and the noise tends to cancel out.

The frequency dependence f_ν of the sources is not known a priori just with the information of the images. A template model from Pearson et al. (2013) developed to estimate redshifts using only the SPIRE fluxes from Herschel has been used as frequency dependence for all the sources considered:

$$\begin{array}{c}\hfill S_{\nu }=A_n[B_{\nu (1+z)}(T_{\mathrm{h}})[\nu ·(1+z){]}^{\beta }+a{B}_{\nu (1+z)}(T_{\mathrm{c}})[\nu ·(1+z){]}^{\beta }]\end{array}$$(3)

where S_ν is the flux at a redshift frequency ν(1 + z), z is the unknown redshift of the source, A_n is a normalization factor, B_{ν(1 + z)} is the Planck function, β = 2 is the emissivity index, T_h = 46.9 K and T_c = 23.9 K are the temperatures of the hot and cold dust components, and a = 30.1 is the ratio of the mass of cold dust to the mass of hot dust.

This template has emerged from a subset of 40 bright Herschel-ATLAS sources with very well known redshifts in the range 0.5 <  z <  4.3. The redshifts of 25 of them, with z <  1, were obtained through optical spectroscopy. The redshifts of the other 15 objects, in the range 0.8 <  z <  4.3, were estimated from CO observations. This SED has also already been used and studied in several previous works (Eales 2015; Ivison et al. 2016; Bianchini et al. 2016, 2018; Negrello et al. 2017; Fudamoto et al. 2017; Bakx et al. 2018; Donevski et al. 2018).

Given that all the sources used to build this model are among the most luminous H-ATLAS sources at their respective redshifts, a bias may arise from the fact that the model may not be representative of the less luminous sources. For instance, low-z H-ATLAS sources have cooler SEDs than the template derived in Pearson et al. (2013) from their high-z spectroscopic sample. It is important to bear in mind that the many different types of sources distributed in the sky constitute a very heterogeneous set of objects that do not have a common spectral behavior. This is the reason why the detection and estimation of the flux of point sources is a difficult task. In this sense, it should be noted that this template model is not expected to be a physically real SED but simply a representative model that can be used as a statistical tool for estimating redshifts from SPIRE fluxes.

We use source positions given by the H-ATLAS catalog and follow the procedure described above. For each source we maximize the S/N of the filtered map, defined as

$$\begin{array}{c}\hfill S/N=\frac{A}{\sigma },\end{array}$$(4)

with respect to the frequency dependence f_ν. In the previous equation A is the amplitude and σ the standard deviation of the point source in the image. Since for the frequency dependence we use the SED template (Eq. (3)) with fixed a, β, T_k and T_c parameters, the only free parameters in the optimization are the source amplitude A and its redshift z. In fact, the amplitude is not really a variable, because for any given set of images it is determined by z for any iteration of the filter through Eq. (1). By construction, the resulting A coincides with the source’s flux density when the optimization is completed. Therefore, the only variable in the optimization is z and the maximization of the filtered S/N of a given source is tantamount to finding its redshift, provided Eq. (3) is a valid description of its SED. In the end, we have a maximized filtered image of the source, with an amplitude A that corresponds to the flux density of the source in the chosen fiducial frequency. The fluxes of the source at the other frequencies can be obtained by multiplying this amplitude by the frequency dependence vector f_ν, which is normalized to the fiducial frequency. This method is robust only in the case of point sources, that is, those whose spatial profile in each frequency agrees with the beam profile in that frequency.

3. Data

Herschel-ATLAS is the extragalactic survey covering the widest area undertaken with Herschel Space Observatory (Pilbratt et al. 2010), imaging 659.25 deg² of the sky distributed in five fields: three (GAMA9 with 53.43 deg², GAMA12 with 53.56 deg² and GAMA15 with 54.5 deg²) on the celestial equator, a large field (180.1 deg²) centered on the north Galactic pole (NGP) and an even larger field (317.6 deg²) centered on the south Galactic pole (SGP). Images have been taken in five far-infrared (far-IR) to submm photometric bands, 100, 160, 250, 350 and 500 μm, using the PACS and SPIRE instruments in parallel mode. PACS measurements have not been used in the main analysis of this work. The main reason is that the SED model (Eq. (3)) from Pearson et al. (2013) exploited to estimate the redshifts has been developed to use only the SPIRE fluxes. This was owing to that not all H-ATLAS sources have flux measurements at PACS wavelengths and only a few per cent of them were detected at greater that 5σ in these bands. However, SPIRE bands themselves are ideal for capturing the emission peak belonging to the high-redshift sources aimed in this work.

Both Data Release 1 (DR1) and the recently released Data Release 2 (DR2) have been used in this analysis. Herschel-ATLAS DR1 includes the three equatorial fields covered by the Galaxy And Mass Assembly (GAMA; Driver et al. 2009, 2016) spectroscopic survey. The three fields are ∼162 deg² combined, and are approximatively located around 9h, 12h and 15h in α. The associated catalog, described in Valiante et al. (2016) and Bourne et al. (2016), covers all three regions and includes 120 230 SPIRE sources, which have at least a S/N = 4σ (including confusion noise) in any of the 250, 350 or 500 μm maps. Herschel-ATLAS DR2 covers the two fields centered in the NGP and SGP, which are about 450 deg² combined. The maps are described in Smith et al. (2017) while the submillimeter catalog is described in Maddox et al. (2018)² and include 118 980 (NGP) and 193 527 (SGP) sources, respectively. These sources have also at least a S/N = 4σ detection in all of the SPIRE bands. The complete H-ATLAS catalog contains a total of 432 737 sources, most of them being point sources. After removing extended sources and stars, the catalog contains 410 997 sources.

As explained in greater detail in Valiante et al. (2016), sources were detected using the MADX algorithm (Multiband Algorithm for source Detection and eXtraction) applied to the SPIRE maps. The first step of this method is to use Nebuliser to remove the diffuse Galactic dust emission from all maps in the three bands, resulting in raw images with the local large-scale background subtracted (“backsub” maps). Then the images are convolved with a proper matched-filter for each band (Chapin et al. 2011). Maps of the variance in each of these convolved maps are also created. During the convolution, the contribution of each pixel of the input image is weighted by the inverse of the square of the instrumental noise in that pixel. The resulting maps are background subtracted maps and noise-weighted maps filtered with a customized matched filter (“fbacksub” maps). In the next step (in MADX), the maps at 350 and 500 μm are interpolated onto images with the same pixel scale as the 250 μm one, and the three images and their corresponding variance maps are then combined together to form a single S/N or “detection” image. In practice, images at 350 and 500 μm are given a zero weighting regarding source detection, that is, the detection image in MADX is simply the 250 μm image. The position of the source in this image will be used to estimate the fluxes of the source in the 350 and 500 μm maps. A list of potential sources is produced by finding all peaks in the detection image with S/N >  2.5σ. A Gaussian fit is carried out in each of these peaks to provide an estimate of the source position and their flux densities are measured at the positions of those peaks in all the SPIRE bands. Finally, only sources with S/N >  4σ in at least one of the three SPIRE bands are kept in the final catalog.

We have worked with the backsub maps instead of the fbacksub ones in order to test our own multifrequency matched-filter’s performance without any other alteration but the subtraction of the large-scale background emission. The method used to subtract this large-scale emission does not affect the flux density of point sources. The units of the maps are Jy/beam. We converted these fluxes to Jy/pixel by dividing the values in the maps by the ratio between the beam area and the pixel area in arcsec² (469/36, 831/64 and 1804/144 at 250, 350 and 500 μm, respectively). The maps have pixel sizes of 6, 8 and 12 arcsec at 250, 350 and 500 μm, respectively. All maps must have the same pixel size so as to be able to combine the three-channel images of a source into one single filtered image. Thus we re-binned 350 and 500 μm maps to a pixel size of 6 arcsec, the same pixel size than the 250 μm map. This repixelization may cause small alignment errors between the pixel positions of the center of the source in the different channels. These pixel misalignments have already been considered and monitored in the method. We have achieved a perfect alignment for ∼90% of the H-ATLAS sources leaving the rest with deviations not greater than 2 pixels in one or some of the channels.

Once all maps have the same pixel size and are in units of Jy/pixel we apply our new algorithm on the positions of the 4σ detections produced by the MADX algorithm in order to test our method. Positions (α and δ) of all point sources identified in the maps are known and taken from the H-ATLAS catalog, converted from degrees into pixels and used to select the objects in the backsub maps. We extract square patches of 256 × 256 pixels centered on the position of the source for the three bands (250, 350 and 500 μm) and align them to run our MMF algorithm. When a H-ATLAS source is close to the edge of the H-ATLAS footprint, the zeros in the map are replaced by white noise generated with the same mean and standard deviation than the rest of the map (within the footprint).

Prior to the filtering step, a prefiltered map without any frequency dependence is built using the Fourier transforms of the three image patches and the matrix H. Cross-power spectrum of the images and matrices T are needed to get matrix H. T matrices are calculated using the point spread functions (PSF) at 250, 350 and 500 μm as source profiles at each frequency. And the inverse of the cross-power spectrum matrix is calculated for each position of a pixel in the images. In the final phase of the filtering process, we introduce the frequency dependence shown in Eq. (3). The 250 μm channel is chosen as the frequency of reference. Since the redshift z of the source is unknown, this last step is repeated for equally distributed redshift values in the range 0 ≤ z ≤ 7, with a step of 0.007 until the redshift which maximizes the S/N, that is, the optimal redshift, for the considered source is found. The result is an estimation of the redshift of the source, its frequency dependence vector f_ν derived from Eq. (3) and a maximized filtered image of the source. The flux at 250 μm is directly the amplitude of the filtered image and the fluxes at 350 and 500 μm can be estimated multiplying this amplitude by the corresponding components of f_ν.

4. Simulations

In order to test our method, we used simulated data with a well known SED. Simulations are useful for both, identifying possible biases and also to estimate the errors in the reconstructed redshift. Simulations were done using only GAMA’s backsub maps from H-ATLAS DR1. The recently released NGP and SGP fields from H-ATLAS DR2 were not used for the simulations but this should have no impact on our results.

We started each simulation with a randomly chosen square patch of the desired size (256 × 256 pixels in our case) from any of the three equatorial fields surveyed (GAMA 9, GAMA 12 or GAMA 15). The same patch region was selected for the three submm photometric bands. Since the three SPIRE channels have different pixel sizes, and the MMF needs to work with a common pixel size, we re-bin the 350 and 500 μm maps to have the same pixel size than the 250 μm map. Alignment errors between the pixel positions of the source in the different channels (which may harm the MMF filtering result), can take place due to this repixelization but, as we already explained in Sect. 3, they have already been considered for the H-ATLAS sources, as they are for simulated sources. Thus, all maps used in simulations have a pixel size of 6 arcsec and are in units of Jy/pixel. Then a source with the corresponding beam profile (according to the PSF of the channel), an adequate amplitude (in order to obtain fluxes like the H-ATLAS ones), and a fixed redshift and SED (Eq. (3)), is placed in the middle of each one of the three patches. From this moment we followed the same procedure, described in Sect. 3, as with any H-ATLAS source. As it is done with the H-ATLAS sources, if the selected map patch contains zeros (i.e., it is near the edge of one of the GAMA fields), these are replaced by white noise with dispersion given by the map background. We performed 5000 simulations, as described before, for each one of the redshifts considered within the ranges 1 ≤ z ≤ 4.5, with a 0.1 step, and 4.5 ≤ z ≤ 7, with a 0.5 step. For each input redshift value considered, z_in, we compute the mean value of the 5000 output redshift values, z_out, and the standard deviation. The difference between the redshifts estimated after applying our MMF method to the simulated sources and the input redshift is shown in Fig. 1 as a function of z_in.

Fig. 1. Difference between the redshift recovered with our MMF method (zout) and the input redshift (zin) as a function of zin. 5000 simulations are run for each zin value in the range 1 ≤ zin ≤ 7. The mean value (zout) is computed and shown as blue dots. The solid line shows a polynomial fit to these mean values. This fit is later used to correct for this bias. zout − zin errors (1σ) are plotted as a shaded region.

The bias observed above z ≈ 4 could be due (in part) to the fact that the Pearson model (Eq. (3)) is built based on some of the most luminous H-ATLAS sources and a restricted range in redshift (0.5 ≤ z ≤ 4.3). However, more importantly, photometric redshifts derived from a SED have problems when the peak of the IR emission is not bracketed by the three SPIRE bands. This peak falls in the SPIRE bands between z ∼ 1 and z ∼ 4, and it is precisely in this redshift range where our method seems more robust returning unbiased redshift estimates. Beyond z ≈ 4 a positive bias can be appreciated which can be as high as ≈0.6 at z ≈ 7. Using a polynomial fit, we find that our estimations of the redshift after applying the MMF can be corrected through:

$$\begin{array}{c}\hfill z_{\mathrm{true}}=0.0013·z_{\mathrm{out}}^4-0.018·z_{\mathrm{out}}^3+0.056·z_{\mathrm{out}}^2+0.943·z_{\mathrm{out}}+0.057\end{array}$$(5)

where z_true is the unbiased redshift estimation of the corresponding H-ATLAS source.

5. Comparison with known-redshift H-ATLAS sources

We compare the redshifts obtained by the MMF method with a set of 32 Herschel-ATLAS sources with known spectroscopic redshifts from Pearson et al. (2013), Negrello et al. (2017) and Bakx et al. (2018). Several of the sources selected are ubiquitous in all these references. Ten of these 32 sources are chosen from Negrello et al. (2017), 17 are sources with z_spec >  0.8 used in Pearson et al. (2013) to build their template and five are taken from Bakx et al. (2018). Redshifts and flux densities estimated with the MMF for these sources are shown in Table 1.

The differences between photometric redshifts estimated with the MMF and the measured spectroscopic redshifts for these 32 objects are shown in Fig. 2. The top plot shows Δz/(1 + z_spec)=(z_photo − z_spec)/(1 + z_spec) before the bias correction. The mean and median are μ = 0.004 and μ_1/2 = −0.017 respectively, with an rms scatter of σ = 0.143. The bottom plot shows the same quantity after the bias correction. Since at z <  4 the bias correction is small, the improvement is small in this redshift range. Nevertheless, the mean and median (μ = 0.009, μ_1/2 = −0.003), and scatter (σ = 0.138) are slightly better than in the sample without bias correction. These statistical parameters are also included in the corresponding redshift column of Table 1. If we take the definition for outliers (those with |Δz/(1 + z_spec)|  >  0.3) used in Ivison et al. (2016), only one of the objects considered is identified as an outlier both for the analysis with biased and unbiased redshifts. Error bars in the top panel are calculated from using Eqs. (1) and (3) while error bars in the bottom panel are derived from simulations described in Sect. 4. We note how the error bars in the bottom panel are more representative of the dispersion around the zero value than the error bars in the top plot. This result indirectly confirms that the error bars derived from the simulation are the most meaningful ones for our estimated redshifts. For the meaning of the error bars in the top panel, see the following note:

Fig. 2. Difference, Δz/(1 + zspec), as a function of zspec between the biased (top) or unbiased (bottom) photometric redshifts estimated with our matched multifilter (MMF) and the spectroscopic redshifts from sources in Table 1. The statistical parameters noted illustrate the systematic overestimates or underestimates, mean μ and median μ1/2, and the degree of scatter, σ, of the photometric redshifts ($z_{\mathrm{MMF}}^{\mathrm{biased}}$ and $z_{\mathrm{MMF}}^{\mathrm{unbiased}}$) found using the MMF and the Pearson et al. (2013) SED. Error bars in the top panel are calculated from using Eqs. (1) and 3) while error bars in the bottom panel are derived from simulations described in Sect. 4. The gray band depicts the 1σ region centered in the median μ1/2 value for each distribution.

A note on error bars: from Eq. (1), and using a parametric SED such as Eq. (3) it is possible, under some general (but not necessarily true) assumptions, such as the statistical independence of the background noise at the three different SPIRE channels, to estimate the degree of uncertainty of our photometric redshift z_MMF. The error bars of all redshift estimations from Tables 1–3, except for $z_{\mathrm{MMF}}^{\mathrm{unbiased}}$, have been obtained this way. However the statistical uncertainty of an estimator, and its actual error with respect to the groundtruth are not necessarily the same thing. The uncertainty given to an estimator can be under, or overestimated depending on the validity of the statistical assumptions made. On the other hand, the estimator may be biased and this bias may not be accounted for in the calculation of the uncertainty. When possible, it is preferable to calculate the error of the estimator using real, already known values of z_spec or, if few spectroscopic redshifts are available, by means of realistic simulations. This is the approach followed in this section to obtain the errors of the unbiased redshifts $z_{\mathrm{MMF}}^{\mathrm{unbiased}}$ from Table 1, showing that for SPIRE, and in the redshift range 1 ≤ z ≤ 7, the actual average error of the estimation of z is under control and typically smaller than the uncertainty calculated from Eqs. (1) and (3).

In order to test the robustness of our results, we have repeated the comparison with the spectroscopic sources, but changing the maps from which we extract these sources and the convolution functions used. We run our method, but using several combinations of the backsub and fbacksub maps as well as SPIRE’s PSFs and MFs. The results obtained with the different configurations are shown in Table 2. We also obtain redshift estimates for these configurations by applying the Pearson’s χ² test statistic, but without using the MMF, and taking into account only the flux density measurements from different maps or the tabulated fluxes from the H-ATLAS catalog, and comparing them with the fluxes predicted by the Pearson et al. (2013) SED. These last results are shown in Table 3.

As can be seen by comparing the redshift estimations for the 32 known-spectroscopic sources from Tables 1–3, both our unbiased and biased estimates obtained with the MMF method on the backsub maps and using the SPIRE PSFs outperform all the other redshift estimations derived with the alternative combinations of maps and convolution functions. Only redshifts estimated by applying a Pearson’s χ² test statistic with flux measurements from fbacksub maps ($z_{{\chi }^2}^{\text{fb}}$) and with tabulated fluxes from the H-ATLAS catalog ($z_{{\chi }^2}^{\text{t}}$) get comparable results. Nevertheless, our unbiased MMF redshifts are the ones which get lower offsets and scatter, and agree with spectroscopic redshifts with the greatest accuracy.

Focusing now our attention on flux densities, practically all recovered fluxes with the MMF (S_MMF) are consistent with the corresponding tabulated fluxes from H-ATLAS catalog (S_SPIRE), in the sense that the IR peak is recovered at the right corresponding wavelength for 29 out of the 32 sources considered from Table 1. On the other hand, and as expected, all IR peaks are recovered in the right band for the flux densities $S_{{\chi }^2}^{\text{fb}}$ taken from the fbacksub maps (see Table 3). The comparison between our MMF estimates of the flux densities and those from the H-ATLAS catalog is shown in Fig. 3. It can be seen how our flux estimations seem to be systematically below the values from the H-ATLAS catalog. This slight underestimate is expected since the noise reduction carried out by the MMF must lead to flux densities lower than the H-ATLAS ones. The average flux underestimates between the flux densities estimated from the MMF method and the H-ATLAS fluxes are 17 ± 13 mJy at 250 μm, 18 ± 9 mJy at 350 μm and 14 ± 14 mJy at 500 μm.

Fig. 3. Flux measurements in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels derived with the MMF versus the corresponding flux measurements from the H-ATLAS catalog in logarithmic scale for the 32 known-spectroscopic sources from Table 1. The dashed line marks perfect correlation. All fluxes are in units of mJy/beam according to the beam profile of the respective channel.

6. High-z candidates in H-ATLAS

In this section we describe our strategy to find high-z candidates in the H-ATLAS data. In the two subsections below, we explore two different strategies. A first subsample (Sect. 6.1) is defined, where the candidates have to be visible in all three bands in SPIRE (this will define our bright subsample of high-z candidates) and for which a more reliable estimate of the redshift can be obtained. In the second subsample (Sect. 6.2), we focus on the 500 μm risers where the highest flux is found in the band with the longest wavelength. Although this does not guarantee that the source is at high redshift, all of the most distant objects in H-ATLAS will be 500 μm risers, as the peak of the IR emission will be at wavelengths longer than 500 μm.

6.1. The bright subsample: High-z candidates with photometric redshift estimates

In order to define a sample of reliable high-z candidates in the H-ATLAS, and motivated by the work exposed in Negrello et al. (2017), we follow a strategy based on applying different cuts to the official photometric catalog and the results of running the MMF. This strategy is based on flux cuts, color cuts, and agreement of fit between the photometric measurements and an assumed SED. By applying a series of cuts to the full sample, we reduce the number of candidates until we arrive to a small subsample of objects which meet all our criteria. Since we focus on high-z candidates, some of the cuts are designed to remove low-z sources. We describe these cuts in detail below:

(i) Since our method only works for the case of strictly point sources and high-z galaxies will appear as unresolved, first we discard all sources we know for sure are non-pointlike, that is, sources which have aperture fluxes at the three SPIRE wavelength different from the point source fluxes. We retain only those sources for which their aperture radius has “−99” value in the H-ATLAS catalog, which means that the aperture flux and the point source flux are the same. We also remove sources identified as stars and those with null or negative fluxes in any of the channels. This results in a sample of 410997 objects from H-ATLAS (see Fig. 4), on all of which our MMF method is applied.

Fig. 4. Evolution of the color–color diagram of the H-ATLAS sources studied as cuts are applied in order to get a sample of robust high-z candidates. The dashed pink line is the track of the Pearson et al. (2013) SED for redshifts in the range [0.5,4.5] (in increasing order from the top-right to the bottom-left corner). The vertical lateral coloured bar present in all plots is a scale of the S/N of the sources exhibited, achieved with the MMF.

The redshift distribution found for these H-ATLAS objects is shown in Fig. 5. An important peak close to zero-redshift can be seen. This figure shows clearly that there are a lot of sources that could be either, (a) low redshift sources or (b) sources for which their frequency dependence does not resemble the Pearson model considered in Eq. (3), and hence are not adequate to be used with our method (resulting in erroneous redshifts and fluxes).

Fig. 5. Redshift distribution, according to the redshift estimates obtained with the MMF method, of the 410 997 sources from H-ATLAS selected after removing the non-point objects and the ones identified as stars.

(ii) After removing stars and non-point sources, we proceed to make a preselection using the photometric information of the H-ATLAS catalog. We want our sources to be at high redshift and bright enough so our redshift estimations are robust. Thus we select those which have a S/N greater than 5σ in all three SPIRE channels, which leaves us with just 9159 sources, 2.2% of the total (see Fig. 4). This is the most stringent cut.

(iii) Good candidates are required to have a flux ratio between 250 and 350 μm bands less than or equal to 1.5 (i.e S₂₅₀/S₃₅₀ ≤ 1.5). This cut has the effect of excluding local galaxies at low redshifts. (see Fig. 4).

(iv) One important requirement for our preselection is to ensure that the chosen sources have a photometric behavior close to the response offered by the Pearson model (Eq. (3)) used to estimate their redshifts, since, as we discussed earlier in Sect. 2, the method does not work equally well for all H-ATLAS sources (see Fig. 6). This can clearly be seen reflected in the large number of sources far away from the Pearson et al. (2013) SED model in Fig. 4. We discard sources that are at a distance larger than 0.3 from the Pearson model in the color–color diagram (see Fig. 4), according to:

Fig. 6. Normalized SED from Pearson et al. (2013), as defined in Eq. (3), at z = 2 in contrast with the normalized tabulated fluxes at 250, 350 and 500 μm (vertical dotted) of two H-ATLAS sources at zMMF = 2: one (J144102.9+012805) that fits well to the model and other one (J233138.6−354345), whose points have been slightly displaced in x-axis to get better clarity, that does not fit properly to the model according to our criteria.

$$\begin{array}{c}\hfill \sqrt{{(\frac{S_{250}^{\mathrm{SPIRE}}}{S_{350}^{\mathrm{SPIRE}}}-\frac{S_{250}^{\mathrm{MMF}}}{S_{350}^{\mathrm{MMF}}})}^2+{(\frac{S_{350}^{\mathrm{SPIRE}}}{S_{500}^{\mathrm{SPIRE}}}-\frac{S_{350}^{\mathrm{MMF}}}{S_{500}^{\mathrm{MMF}}})}^2}<0.3.\end{array}$$(6)

The number 0.3 is a compromise between a more stringent requirement that would result in a smaller number of candidates and a more relaxed requirement that would increase the number of candidates but at the expense of increasing the number of sources with unreliable redshift estimations.

(v) The last requirement in our preselection is to exclude the presence of possible blazars. As showed in Negrello et al. (2017), the leaking of blazars into a catalog of high-z candidates can be reduced by demanding our sources to have S₃₅₀/S₅₀₀ >  1, unless this ratio is already above the Pearson model in the color-color diagram (see Fig. 4).

At this point, after the cuts (i) through (v), 5079 sources remain in the sample. Since this sample will be used later, we denote it the “full high-z” sample. These cuts are not perfect at removing low-z objects but the sample should be dominated by high-z candidates.

(vi) Finally, to reduce this level of contamination, and select high-S/N sources for which photometry is expected to be robust (and consequently the photometric redshift as well), we impose the condition that the S/N defined in Eq. (4) must be greater than or equal to 15σ in the filtered image after our MMF has been applied, except if the H-ATLAS SPIRE position has an association with a galaxy of known spectroscopic redshift within a separation of 5 arcsec. This cut leaves 370 objects (see Fig. 4). This selection constitutes our “robust high-z” sample of high-z candidates, with a redshift distribution of $\overline{z}=2.13$ and σ = 0.65, and will be used later to identify possible lensed systems.

This sample is partially shown in Appendix A including the estimated redshifts and flux densities. The entire catalog is available at the CDS. We have performed an additional study for the objects within this sample by comparing the Pearson’s χ² value obtained considering only the SPIRE flux densities with the one obtained taking into account also the PACS flux densities. For this we have used the flux densities from the H-ATLAS catalog as the observed data and the frequency dependence provided by the Pearson SED (Eq. (3)), at the photometric redshift estimated by our MMF, as the theoretical data. The result of this study is shown within the online catalog through the flag “Reliability”. Those sources for which the χ² improves or remains the same when using PACS fluxes are flagged with a “0”, whether the χ² worse slightly but it is still acceptable they are flagged with a “1”, if the χ² is much worse they are flagged with a “2” and if the source does not have PACS flux densities we can use, it is flagged with a “−1”.

Within this sample, 201 candidates are in the GAMA fields (60 in the GAMA9, 58 in the GAMA12 and 83 in the GAMA15), 82 in the NGP and 87 in the SGP. The number density of sources in the GAMA fields is higher than in NGP and SGP after this last cut because the number of associations with objects having spectroscopic redshifts is higher in the GAMA fields. Among all the objects of this sample we find 35 QSOs. Figure 7 shows the redshift distribution of the robust high-z sample in order to compare it with the redshift distribution of the initial sample shown in Fig. 5.

Fig. 7. Redshift distribution, according to the redshift estimates obtained with the MMF method, of the 370 high-z H-ATLAS sources from the robust high-z sample selected after imposing all cuts enumerated in Sect. 6.1.

A direct comparison between our estimates of the flux densities in all channels and the tabulated fluxes from the H-ATLAS catalog is shown in Fig. 8. A clear linear trend is observed and the accordance is pretty good. As happens with the spectroscopic redshift sources from Fig. 3, an overall underestimation of our MMF flux densities, greater for fainter sources, can be observed. As the number of sources is greater, here the effect is most remarkable. The average flux underestimates between the flux densities estimated from the MMF method and the H-ATLAS fluxes are 10 ± 9 mJy at 250 μm, 12 ± 9 mJy at 350 μm and 9 ± 8 mJy at 500 μm. Since the MMF combines information from all three wavelengths, which allows to reduce the background and boost the signal, instrumental, foreground and confusion noises are better removed so flux density estimates are less affected by Eddington bias than H-ATLAS flux densities. This underestimation with respect to H-ATLAS flux densities is stronger toward low flux densities, which supports the Eddington bias hypothesis, but is also observed to a lesser extent for high flux densities, suggesting a possible degradation of the MMF photometry that could be related to the way we re-pixelize the 350 and 500 μm images and combine them during the multifiltering step (see Sect. 3).

Fig. 8. Flux measurements in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels after our MMF has been applied versus the corresponding tabulated flux measurements from the H-ATLAS catalog in logarithmic scale for the 370 high-z H-ATLAS from the robust high-z sample (see Appendix A). The linear behavior with zero-intercept is drawn with a black dashed line. All fluxes are in units of mJy/beam according to the beam profile of the respective channel.

Figure 9 shows the improvement in S/N achieved with our MMF method for the robust high-z sample in contrast with the S/N of the three μm SPIRE bands. An average improvement of 76% in the S/N has been achieved for this sample with our MMF technique compared to the 500 μm band. Besides, an average improvement of 16% and a slight improve of 0.2% have been obtained for the 350 μm and 250 μm, respectively.

Fig. 9. S/N in the filtered image after our MMF has been applied versus the S/Ns in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels in logarithmic scale for the 370 high-z H-ATLAS sources from the robust high-z sample. The linear behavior with zero-intercept is drawn with a black dashed line.

In the end, we have a selected sample that includes several hundreds of interesting objects from H-ATLAS which both agree with the Pearson et al. (2013) SED used to estimate the redshift and have high redshifts and S/Ns.

6.2. Faint subsample: “500 μm-risers”

Apart from the robust high-z sample explained above in Sect. 6.1, we also looked for faint sources at 250 and 350 μm but bright at 500 μm in the H-ATLAS data, the so-called “500 μ-risers”. Our selection criterion looks for sources whose detection is at least barely significant at 500 μm and that are not clearly detected at 250 and 350 μm in the H-ATLAS catalog. We select objects with S/N₅₀₀ ≥ 4σ, S/N₂₅₀ ≤ 4σ and S/N₃₅₀ ≤ 4σ in the H-ATLAS catalog and apply our multifrequency MMF filter to them, in order to enhance the statistical signification of those detection candidates. Those sources with S/N ≥ 5σ after the MMF filtering and that satisfy the condition S₂₅₀ ≤ S₃₅₀ ≤ S₅₀₀ are considered to be statistically significant enough to be firm candidates to be 500 μm-risers. This way, we get a sample of 695 reddened SPIRE objects. We must not forget the limitations of the Pearson et al. (2013) SED (Eq. (3)) used to estimate the redshifts so by selecting again the sources which fit better to the model in the color-color diagram (Eq. (6)) we are left out with 237 objects. This selection constitutes our 500 μm-riser sample of robust high-z candidates, with a redshift distribution of $\overline{z}=4.62$ and σ = 0.71. This sample is partially shown in Appendix A including redshift and flux density estimates. The entire catalog is available at the CDS. The same additional χ² study, considering PACS flux densities, performed for the “robust high-z” sample has been applied to this sample, and the result is shown within the online catalog through the same flags explained in Sect. 6.1. Within this sample, 97 objects are from the GAMA fields (27 in the GAMA9, 37 in the GAMA12 and 33 in the GAMA15), 68 from the NGP and 135 from the SGP.

The comparison between our estimates of the flux densities in all channels and the tabulated fluxes from the H-ATLAS catalog for the 500 μm-riser sample is shown in Fig. 10. A much larger scattering than the one seen in Fig. 8 for the robust high-z sample can be observed. But this behavior was expected as we are aiming to sources which have a barely significant detection at 500 μm and are not detected at 250 and 350 μm. The average flux underestimates between the flux densities estimated from the MMF method and the H-ATLAS fluxes are 4 ± 4 mJy at 250 μm, 0.4 ± 4 mJy at 350 μm and 3 ± 5 mJy at 500 μm.

Fig. 10. Flux measurements in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels after our MMF has been applied versus the corresponding tabulated flux measurements from the H-ATLAS catalog in logarithmic scale for the 237 high-z H-ATLAS from the 500 μm-riser sample (see Appendix A). The linear behavior with zero-intercept is drawn with a black dashed line. All fluxes are in units of mJy/beam according to the beam profile of the respective channel.

Figure 11 shows the comparison between the S/N reached with our MMF method and the S/N in all three SPIRE channels. It seems logical that the improvement achieved with our method in S/N for these 500 μm-riser objects (Fig. 11) should be better than for the objects from the robust high-z sample (Fig. 9), as they are near the H-ATLAS detection limit. This is confirmed since we have achieved average improvements of 25%, 55% and 76% in the S/N for the 500, 350 and 250 μm, respectively. This clearly reflects that it is in this kind of faint objects where our MMF method accomplishes bigger impact in terms of signal significance.

Fig. 11. S/N in the filtered image after our MMF has been applied versus the S/Ns in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels in logarithmic scale for the 237 high-z H-ATLAS from the 500 μm-riser sample. The linear behavior with zero-intercept is drawn with a black dashed line.

Unlike the previous robust high-z sample that sought to select bright objects in all bands, among which the probability of finding lensed systems is relatively high, now we pursue faint high-z objects. Most of them are not expected to be lensed by foreground sources but to be intrinsically luminous. These robust high-z and 500 μm-riser samples are, in fact, built from starting requirements mutually excluding. Sources from the robust high-z sample are initially required to have S/N₂₅₀, S/N₃₅₀ and S/N₅₀₀ greater than 4σ while sources from the 500 μm-riser sample are demanded to have S/N₅₀₀ ≥ 4σ, S/N₂₅₀ ≤ 4σ, and S/N₃₅₀ ≤ 4σ. But this does not mean, for example, that there are not candidates to lensed sources among H-ATLAS 500 μm-riser galaxies. It should be pointed out, for instance, the case of J090045.4+004125 (α = 135.191, δ = 0.6897), a dusty star-forming galaxy at z = 6 revealed by strong gravitational lensing and detected in the GAMA field being part of a subsample of the H-ATLAS 500 μm-riser galaxies (Zavala et al. 2018). This object does not appear in the robust high-z sample but is part of our 500 μm-riser sample under the identifying name J090045.5+004131 with a redshift estimation of z = 6.35 via our MMF technique.

An important effort was made recently by Ivison et al. (2016) in order to take advantage of the 250, 350 and 500 μm Herschel-ATLAS imaging survey and select extremely red objects. That work focused on studying the space density of luminous dusty star-forming galaxies (DSFGs) at z >  4 by selecting galaxies from the H-ATLAS survey with extremely red far-infrared colors and faint 350 and 500 μm flux densities, called ultra-red galaxies. It is important to bear in mind that they used a modified version of the MADX algorithm to identify their sources, so some of their sources are not in the official H-ATLAS catalog. This fact explains why we have been able to locate in the H-ATLAS catalog only 78 of the 109 sources shown in their sample.

None of these 78 red sources is, of course, in our robust high-z sample since they were detected at a S/N ≥ 3.5σ at 500 μm, being mostly brighter in this band than in the others, and our robust high-z sample was built with relatively bright sources in all bands. Instead, and as expected, there is an overlap between our 500 μm-riser sample and these red sources from Ivison et al. (2016). These 78 objects were required to be above 3.5σ in any of the three SPIRE bands and our 500 μm-riser sample was built demanding S/N₅₀₀ ≥ 4σ, S/N₂₅₀ ≤ 4σ and S/N₃₅₀ ≤ 4σ. To begin with, 54 out of these 78 objects are not included in our 500 μm-riser sample because all they have S/N₃₅₀ >  4 and thus are excluded by our criterion, which leaves us with 24 possible objects. Only nine of these remaining 24 objects from Ivison et al. (2016) (J090045.5+004131, J090304.5−004616, J114038.8−022804, J114350.3−005210, J114353.5+001250, J114412.1+001812, J115614.0+013900, J142710.6+013806 and J004615.0−321825) are included in our 500 μm-riser sample. If sources were not demanded to behave like the Pearson SED, there would be 16 objects.

In the next section we will study the correlations between the high-redshift sources and their possible lenses. We will focus on the robust high-z sample because we find no significant correlation for the case of the 500 μm-riser sample. This suggests us that this last sample is mostly not lensed as expected due to the lower flux criterion used to select its sources.

7. Possible lensed galaxies

7.1. Preliminary comparison with previous works

In the previous section we presented our samples of high-redshift candidates (full, robust high-z, and 500 μm-riser samples). Since the optical depth of strong lensing grows with the redshift of the background source, these samples of high-z candidates may contain some lensed galaxies. In fact, this is known to be particularly true for Herschel sources, where the brightest high-z sources correspond to SLGs (Negrello et al. 2007, 2014, 2017; Hopwood et al. 2010; Cox et al. 2011; Conley et al. 2011; Lapi et al. 2011; González-Nuevo et al. 2012; Messias et al. 2014; Dye et al. 2015).

Here we compare the robust high-z sample (370 candidates) with similar catalogs found in the literature. Our robust high-z sample contains 62 of the 80 candidate SLGs with flux density above 100 mJy at 500 μm presented in Negrello et al. (2017). 17 of the candidates in the robust high-z sample are part of the sample of 20 confirmed SLGs (Negrello et al. 2017). The only three confirmed strong lens systems that are not included in our sample are J085358.9+015537 (flagged as a star), J142935.3−002836 (which is a major merger system at z = 1.027 Messias et al. 2014 and is excluded by our cut iii) and J125135.3+261457 (excluded by our cut v). Among the confirmed lensed galaxies, J114637.9−001132 at z = 3.26 is interesting since it is associated to a candidate high-z proto-cluster³(Fu et al. 2012; Herranz et al. 2013; Clements et al. 2016; Greenslade et al. 2018).

In addition, six of the eight objects labeled in Negrello et al. (2017) as likely to be lensed and 39 of the 51 objects defined as unclear are included in our robust high-z sample of 370 candidates. The two missed objects labeled as likely to be lensed were excluded in our cut iv. The only one object from Negrello et al. (2017) confirmed to not be a strongly lensed galaxy (J084933.4+021442) is nor part of our sample because it is flagged as a star in the H-ATLAS catalog. It is indeed a binary system of Hyper Luminous Infrared Galaxies (HyLIRGs) at z  =  2.410 (Ivison et al. 2013). Our sample also contains five sources from the SGP field (J004736.0−272951, J011424.0−333614, J235623.1−354119, J001010.5−360237 and J014849.3−331820) which meet the flux criterion demanded by Negrello et al. (2017) but are not in their proposal.

In González-Nuevo et al. (2012), the authors applied to the H-ATLAS Science Demonstration Phase field (≃14.4 deg²), which covers part of the GAMA9 field, a method for efficiently selecting faint candidate SLGs. This method was called HALOS (Herschel-ATLAS Lensed Objects Selection). They found 31 candidate SLGs, whose respective candidate lenses were identified in the VIKING near-infrared catalog and proposed that the application of HALOS over the full H-ATLAS surveyed area would increase the size of the sample up to ∼1000 SLGs. Eight of these sources are included in our robust high-z sample of 370 sources: J090302.9−014127, J090311.6+003906, J090740.0−004200, J091043.0−000321, J091304.9−005343 (all of them confirmed as strongly lensed in Negrello et al. 2017), J085855.3+013728, J090957.6−003619 and J091331.4−003644.

The H-ATLAS catalog can be used to find potential lens systems (lens plus lensed galaxy) using the already available optical associations with SDSS (Blanton et al. 2017) for each SPIRE source (Bourne et al. 2016). These associations are sought via a Likelihood-Ratio analysis of optical candidates within 10 arcsec of all SPIRE sources with S/N ≥ 4 at 250 μm. Bourne et al. (2014) studied the fact that redder and brighter submm sources have optical associations with greater positional offsets than would expected if they were due to random positional errors. They concluded that lensing is the most plausible cause for increased offsets of red submm sources and that the problem of misidentifying a galaxy in a lensing structure as the counterpart to a higher redshift submillimeter galaxy may be more common than previously thought. Most of these optical associations do not have spectroscopic information (i.e. secure redshift), however, there are 180 objects in our robust high-z sample for which this condition is fulfilled (mostly because of the cut vi). Spectroscopic redshifts are obtained from many different surveys, like SDSS DR7, SDSS DR10, 6dFGS, 2SLAQ or GAMA. 138 sources out of these 180 have a reliable spectroscopic redshift (Z_QUAL ≥ 3) in the range 0.1 ≤ z≤ 1.1, that is significantly smaller than the photometric redshift estimated by our MMF method. Hence, these associations may correspond to possible lens systems since the redshifts of the alleged lens and the high-z candidate are so different. In those cases where the optical association is not the same object as the SPIRE source, it will be an object at a smaller redshift and close (in angular separation) to the SPIRE source. The conditions would be given for the lens effect to take place and these cases should be studied in detail to verify it. However, since these associations are already given in the catalogs themselves and their spectroscopic redshifts come from many different sources, we are going to proceed to look for our own associations.

The above discussion shows how our robust high-z sample has the potential to host many unknown lensed galaxies. Most of the previously known Herschel lensed galaxies were unveiled by the 500 μm flux density criterion (S₅₀₀ >  100 mJy), which has proven to be a simple (but powerful) method of selecting strongly lensed candidates. Here we rely on a cross-correlation study based on matching distant IR sources with foreground potential lenses located at distances that make them consistent with being a lens system.

7.2. Statistical lensing analysis. Correlation analysis with SDSS

Additional evidence for significant lensing in our two samples (full and robust high-z) can be obtained through a simple correlation analysis with a catalog of foreground galaxies. If the Herschel sources are tracing the magnification pattern produced by a population of lenses at z <  1, one would expect an excess of IR sources detected around regions of magnification larger than one. Alternatively, the alleged high-z source could be instead a lower redshift associated with the lens. In this case, the excess found in the correlation would be produced by contamination of our sample (i.e low-z sources being misinterpreted as high-z sources).

For the catalog of potential foreground lenses (z <  1), we use lenses extracted from the SDSS. By lenses, we mean here either individual galaxies or groups of galaxies (see below). Since SDSS does not cover the SGP field, we consider only the IR sources which come from the GAMA and NGP fields. After removing IR high-z candidates from the SGP field, the full high-z sample is reduced to 2828 candidates while the robust high-z sample is left with 283 candidates. For a simple estimation of the correlation, we compare the number of matches found within an aperture and for different aperture radii, N_m(R), with the expected number from a random distribution (N_r(R), see Eq. (7)). This random number is obtained by the following equation:

$$\begin{array}{c}\hfill N_r\left(R\right)=N_{\mathrm{H}}\frac{A_{\mathrm{c}}\left(R\right)}{A_{\mathrm{H}}}\end{array}$$(7)

where N_H is the number of H-ATLAS high-z candidates, A_c(R) is the total area covered (within the footprint of H-ATLAS) by the disks of radius R around the SDSS sources, and A_H = 341.65 square degrees is the total area of H-ATLAS survey excluding the SGP field. By construction, N_r(R)≤N_H.

On the other hand, the number of matches (N_m(R)) between the H-ATLAS sources and the SDSS lenses is obtained by computing the number of associations between both catalogs as a function of radius by centering disks of radius R on the SDSS lenses and counting the number of H-ATLAS sources which fall within the disk. Any significant excess over the expected value in the random case is evidence for either lensing or contamination. The uncertainty, or significance, with respect to the random case is given by the Poissonian error (i.e., the uncertainty is given by ${\sigma }_r(R)=\sqrt{N_r(R)}$). If the excess is due to contamination, this hypothesis can be tested, since one would expect the separation between the positions of the Herschel sources and the SDSS lenses to be comparable to the positional error in Herschel (which is significantly larger than the corresponding error in SDSS), that is 2–3 arcsec. In these cases, the Herschel source may actually be the SDSS lens. If, on the contrary, a high-zHerschel candidate is found at more than 3 arcsec from the SDSS source, lensing is possibly responsible for that association. Some of the associations should be due to pure random alignments but this number can be estimated by the Poissonian expectation number discussed above.

For the SDSS lenses, we use two catalogs of potential lenses derived from SDSS. The first catalog focuses on rare but massive potential lenses at z ≤ 0.6 while the second catalog focuses on less massive, but more abundant, potential lenses with z ≤ 1.1. We set a lower limit to the redshift of the potential lenses since below certain redshift, strong lensing becomes inefficient due to the increase in the critical surface mass density (z_min ∼ 0.1).

For the association with massive lenses, we use the SDSS DR8 redMaPPer cluster catalog with 26 111 objects (Rykoff et al. 2014). This catalog is the result of applying the Red-sequence Matched-filter Probabilistic Percolation (redMaPPer) cluster finding algorithm to the 10 400 deg² of photometric data from the Eighth Data Release (DR8, Aihara et al. 2011) of the SDSS. The redMaPPer algorithm has been designed to handle an arbitrary photometric galaxy catalog, with an arbitrary number of photometric bands (≥3), and performs well provided the photometric bands span the 4000 Å break over the redshift range of interest. It adapts therefore well to a survey such as the Sloan Digital Sky Survey. Because the number of objects in this catalog is not very large, we use all of them in the cross-correlation which cover a range of 0.08 ≤ z ≤ 0.6 in redshift and 19.85 ≤ λ ≤ 299.46 in cluster richness. N_H = 881 of the 26 111 halos fall in the footprint of H-ATLAS. We find no significant excess when cross-correlating redMaPPer with our catalog of high-z H-ATLAS sources. Given the fact that 17 of the 20 strongly lensed candidates from Negrello et al. (2017; confirmed as such) are in our selected sample, and none of them has a match with redMaPPer, this confirms that the lenses in Negrello et al. (2017) are not massive halos, but rather relatively small halos (like elliptical galaxies for instance).

Our second search for potential lensed galaxies focuses on the low-mass regime of the lenses. From SDSS DR14 (Abolfathi et al. 2018) we select a larger catalog of galaxies with known spectroscopic redshifts. We focus on galaxies with known redshift in order to minimize possible contamination from galaxies that are misinterpreted as having z >  0.15 and also to reduce the computation time. The sample contains 1 776 242 galaxies from the Sloan Digital Sky Survey Data Release 14 with 0.15 ≤ z ≤ 1.1. As mentioned above, we limit the minimum redshift to 0.15 since below this redshift most galaxies are expected to be subcritical (and not produce strong lensing effects). Among all of them, N_H = 50 175 are the galaxies that fall in the footprint of H-ATLAS. We cross-correlate our full high-z subsample of 2828 H-ATLAS sources with the SDSS catalog of galaxies and compare it with the expected number in the case of no correlation (i.e, the random case described above). The ratio of the observed (N_m(R)) and random matches (N_r(R)) between this catalog and our full high-z selection sample is shown in Fig. 12 for different radii. In Fig. 13 we exhibit the same but for our robust high-z subsample of 283 candidates. Both Figs. 12 and 13 show a non-one signal for aperture radii of several arcminutes which is unexpected and an example of the lensing-induced cross-correlations between high-z submillimeter galaxies and low-z galaxy population (Wang et al. 2011; González-Nuevo et al. 2014, 2017; Bourne et al. 2014).

Fig. 12. Ratio between the number of matches found (Nm(R)), after cross-matching the full high-z subsample of 2828 H-ATLAS sources with a sample of 1 776 242 known-redshift galaxies from SDSS DR14, and the number of matches expected (Nr(R)) from a random distribution (Eq. (7)) for different aperture radii R around the SDSS sources.

Fig. 13. Ratio between the number of matches found (Nm(R)), after cross-matching the robust high-z subsample of 283 candidates with a sample of 1 776 242 known-redshift galaxies from SDSS DR14, and the number of matches expected (Nr(R)) from a random distribution (Eq. (7)) for different aperture radii R around the SDSS sources.

There is a clear increase in significance when considering the robust high-z subsample of 283 candidates. A sharp increase in the excess of matches is found at distances below 60 arcsec. The smaller amplitude of the excess in the full high-z sample with 2828 sources suggests that this sample may be more contaminated by low-z candidates.

Focusing on the smaller radii, we find 40 associations at a separation lower than 20 arcsec, between the 50 175 known-spectroscopic-redshift galaxies from SDSS DR14 that fall in H-ATLAS footprint and our robust high-z subsample of 283 H-ATLAS sources. One of them (J145420.6−005203) is identified as a QSO in the H-ATLAS catalog. We chose to consider this separation since most of the lenses would have an Einstein radius less than 20 arcsec, which is the radius around which the strongest magnifications are expected. From among these matches, 28 have a separation greater than the positional error in Herschel (> 3 arcsec) so lensing is possibly responsible for that association. ∼4 associations should be caused due to pure random alignments so it is expected that a considerable number of these associations are lensed. These 40 matches are shown in Table 4. And snapshots of them, centered on the SDSS DR14 sources, are shown in Figs. B.1–B.3.

We have used the SDSS DR14 asinh magnitudes in the r-band of these 40 low redshift optical sources shown in Table 4 to get a rough estimation of the Einstein radius of each possible lensed system. Firstly, we have estimated from van Uitert et al. (2015) the corresponding corrections for the redshift of their spectra (i.e., the k-correction) and for the intrinsic evolution of their luminosity (i.e., the e-correction) in order to correct the r-band magnitudes. Once the magnitudes of the optical sources are corrected, we have calculated their fluxes and then their luminosities (through their luminosity distances D_L). At this point, we used the luminosity-to-halo mass relation M_eff = M_0, L(L/L₀)^β_L parametrized in van Uitert et al. (2015) to estimate the mass of each SDSS galaxy for the corresponding luminosity previously obtained. The pivot luminosity L₀ is the same for every object while the M_0, L and β_L parameters depend on the spectroscopic redshift of the galaxy. Finally, we supposed the galaxy behaves as a singular isothermal sphere to estimate the Einstein radius (see Narayan & Bartelmann 1999). We assumed the virial radius of the galaxy to be r = 1.3(M/10¹⁵ M_⊙)^1/3 Mpc in order to estimate its velocity dispersion ${\sigma }_{\mathrm{v}}=\sqrt{GM/2r}$. The Einstein radius can be then estimated through ${\theta }_{\mathrm{E}}=(4\pi {\sigma }_{\mathrm{v}}^2{D}_{\mathrm{ds}})/(c^2{D}_{\mathrm{s}}$), where D_ds and D_s are the angular diameter distances between the lens and the source, and observer and source, respectively. These distances are calculated with the spectroscopic redshift of the SDSS galaxy acting as lens ($z_{\mathrm{spec}}^{\mathrm{SDSS}}$) and the photometric redshift of the source estimated with our MMF ($z_{\mathrm{MMF}}^{\mathrm{unbiased}}$). These Einstein radius estimates are included in Table 4 and endorse the choice of 20 arcsec as upper limit for the separation between H-ATLAS sources and their partners from SDSS DR14. 15 of the 40 matches have a separation radius smaller than the corresponding estimation of the Einstein radius. For these 15 we will expect to find counter-images outside the Einstein radius. Possible reasons why we do not find those could be: (i) the foreground galaxy is not a lens, (ii) there is substructure around the lens that can create fluctuations in the flux or (iii) the Einstein radius is overestimated. Einstein rings corresponding to this radii have been drawn in the snapshots from Appendix B.

All H-ATLAS objects from Table 4, except four of them, had a previous association with a galaxy with spectroscopic redshift at less than 10 arcsec. In our search for possible lens systems, we extend our association radius to 20 arcsec. We find that 2 of the previously unmatched H-ATLAS sources (J132111.2+265009 and J134429.5+303034) can now be matched with a galaxy having spectroscopic redhsift within 20 arcsec. For the other two remaining unmatched sources (J131533.9+233254 and J141351.9−000026), we find associations within 10 arcsec with SDSS galaxies having spectroscopic redshift. These last two associations may have not been included in the H-ATLAS catalog because SDSS DR7 and SDSS DR10 were used and we cross-match our robust high-z sample with known-low-redshift galaxies from SDSS DR14.

There are five sources for which there is a mismatch between the redshift published in the H-ATLAS catalog and the redshift of the match found in this work. These sources are J090038.0−003522, J091331.4−003644 (SLG candidate from González-Nuevo et al. 2012), J114512.9−002446, J115101.7−020024 (labeled as unclear in Negrello et al. 2017) and J145420.6−005203. In all cases, our matches are beyond the 10 arcsec radius used in the H-ATLAS association. All these five sources lie in the GAMA fields for which more redshift information is available other than the SDSS redshifts (and that was not used in our association). In fact, we can confirm that spectroscopic redshifts for the first four sources come from the GAMA project (Driver et al. 2009, 2016) but the last one comes from the WiggleZ Dark Energy Survey.

After visual inspection, among the 40 matches we identify interesting cases that are likely to be lensed. Three of them (J090311.6+003907, J134429.5+303034 and J141351.9−000026) are already confirmed as strongly lensed systems in Negrello et al. (2017). Some other SDSS stamps reveal overdensities of galaxies at z <  1 that could indicate the presence of a group of galaxies acting as a lens for the corresponding candidate. These are mainly J083726.1+015641, J085229.0+010217, J085250.9−010000, J115819.5−000127, J130138.7+302808, J131407.0+271208, J131533.9+233254, J133231.6+350843, J133806.8+351530, J141832.9+010212 and J143845.8+013503.

The properties and complex nature of the strong gavitationally lensed system J090311.6+003907 (proposed as candidate SLGs in González-Nuevo et al. 2012) have been studied in detail in Dye et al. (2015) by modeling ALMA long baseline imaging. This source has an optical association in the H-ATLAS catalog at 1.87 arcsec with z = 0.3 which is compatible with both the assocciation we found after the cross-matching and the association identified in Negrello et al. (2017).

The source J134429.5+303034 has no association in the H-ATLAS catalog (within 10 arcsec) with an optical counterpart with known-redshift. Although we find an SDSS galaxy (with unknown redshift) at 0.43 arcsec distance. At larger radii, we find an association with a SDSS galaxy having z_spec = 0.67 at 16.96 arcsec. In Negrello et al. (2017) they find a potential lens with z_spec = 0.6721 but much closer to the H-ATLAS source (0.43 arcsec). A closer look at the SDSS images reveals how the potential lens can be barely resolved in SDSS images. We propose that the small lens indentified in Negrello et al. (2017) is part of a larger halo (or lens) at z = 0.67.

The source J141351.9−000026 resembles the situation of the previous source. It has no optical association with any galaxy having known spectroscopic redshift in the H-ATLAS catalog but we find a galaxy (with unknown redshift) in the SDSS images at 1.12 arcsec distance. For this source we find an association with a SDSS galaxy at 7.61 arcsec at z = 0.547. This association differs from the one in Negrello et al. (2017) which finds a much fainter galaxy at 1.12 arcsec (also at z = 0.547). A closer inspection of the optical high resolution image in Negrello et al. (2017) reveals a giant arc at the position of the IR source. The morphology of the arc (with an arc radius of order 10 arcsec and aligned in a tangential direction with respect to our lens candidate) suggests that our association is the correct one. However, the galaxy found in Negrello et al. (2017) could still be contributing to the magnification since it seems to lie at, or near, the critical curve (for the redshifts of the lens and IR source).

The object J115101.7−020024 is also present in Negrello et al. (2017) proposal. It is labeled as an unclear case because no optical association was found. But we have found an association for this source with z_spec = 0.4543 at 13.63 arcsec which differs from the optical association assigned in H-ATLAS catalog.

J131533.9+233254 is an interesting example. The H-ATLAS catalog contains no optical association. We find an association in SDSS at 8.61 arcsec and with z  =  0.551. The snapshot of this association in Fig. B.2 reveals a number of galaxies in the vicinity of the SDSS galaxy, suggesting that the lens may be a group or a cluster at z  =  0.55. The stamp for J141351.9−000026 (discussed above) resembles this situation although with fewer SDSS galaxies. Other interesting candidates are J083726.1+015641, J085229.0+010227 and J115819.5−000127 which show a concentration of SDSS galaxies near the IR source. The already mentioned object J091331.4−003644 and J085855.4+013729 are the only ones candidates from González-Nuevo et al. (2012) presented in the robust high-z sample for which we found associations with spectroscopic redshifts.

The selection shown in Appendix B contains only a small sample of potential candidates to be lensing systems. Our search for this type of systems is limited by the fact that we restrict ourselves to SDSS galaxies with spectroscopic redshifts. Future surveys with many more galaxies having spectroscopic redshifts (like EUCLID), or even photometric redshifts with accurate photo-z (like J-PAS) can increase significantly the number of lensed systems.

8. Conclusions

In this paper we have extended previous work on the multifrequency matched multifilter (MMF) technique to submm astronomy, in particular the Herschel-ATLAS fields. Our multifrequency technique allows us to, on one hand, boost the statistical significance and S/N of the detections and, on the other hand, provide an estimation of the photometric redshift. This approach relies on semi-analytic modeling of the thermal peak of the SED of dusty galaxies (see e.g., Pearson et al. 2013). We remark that the use of SED templates is not strictly necessary for the MMF technique⁴, but it becomes very useful for the purpose of redshift estimation. We have tested our technique with both realistic simulations and also with H-ATLAS sources for which their redshifts are known, concluding that the MMF indeed leads to an improved S/N with respect to single-frequency detection and that the photometric redshifts are relatively accurate, specially for those H-ATLAS sources with a photometric behavior close the SED used in the range 0.8 <  z <  4.3. Outside of this range, we have studied the bias in the MMF photo-z estimation and provided correction factors up to z ∼ 7.0.

410 997 point-like sources from the DR1 and DR2 H-ATLAS fields (Valiante et al. 2016; Maddox et al. 2018) have been re-analyzed with this multifrequency filter. Two different sets of selection criteria based on flux, color, S/N, and agreement between observed photometry and the SED model have been carried out in order to select the H-ATLAS sources for which the MMF method performs the best and illustrate its potential usefulness. A total of 607 H-ATLAS sources, for which new photometry and photometric redshifts are provided, are found. They are split in two samples of interest:

• A sample of 370 bright high-z objects with $\overline{z}=2.13$ and σ_z = 0.65. This sample contains 62 of the 80 strongly lensed candidates identified in Negrello et al. (2017), and including 17 out of the 20 confirmed strongly lensed galaxies. An average improvement of 76% in the S/N has been achieved for this sample with our MMF method compared to the 500 μm band. Besides, an average improvement of 16% and a slight improve of 0.2% have been obtained for the 350 μm and 250 μm, respectively. By cross-correlating this high-z sample with a catalog of low-known-redshift galaxies from SDSS DR14, we find 40 potential candidates to being gravitational lens systems which have been shown and discussed. Some SDSS stamps around these candidates reveal overdensities of galaxies at z <  1 that could indicate the presence of a group of galaxies acting as a lens.

• A second sample of 237 faint high significance candidates to 500 μm risers (with $\overline{z}=4.62$ and σ_z  =  0.71) which were previously near the H-ATLAS detection limit but are now confirmed with the MMF as high significance detections. We have achieved average improvements of 25%, 55% and 76% in the S/N for the 500 μm, 350 μm and 250 μm, respectively. This clearly reflects that is in this kind of faint objects where our MMF method reaches its full potential in terms of S/N improvement.

The MMF has proven to be able to enhance the sensitivity of fainter objects and improve its photometry. According to the comparative analysis of the photometric redshift estimates with known-redshift H-ATLAS sources, it may also return more robust photometric redshifts than the ones that would be obtained using the fluxes of the H-ATLAS catalog, as proven by the mean μ and median μ_1/2 values shown in Tables 1, 3.

Running a non-blind search, using the already-known source positions in the public H-ATLAS catalog, allows us to directly compare our method with the detection pipeline of the public H-ATLAS catalog. However, given that, in the end, H-ATLAS detection image simply reduces to the 250 μm channel (Valiante et al. 2016), in the future we plan to apply this method in a blind search mode over the H-ATLAS maps. Thereby we will be able to exploit the information from the other wavelengths and be able to detect faint sources that we could be missing now.

The work done here could also be extended in the future focusing only in the sources already confirmed. One of the possible deeper analysis or improvements would be to repeat the MMF execution using other SEDs different from the Pearson et al. (2013) model. For instance, SEDs derived from ultraluminous dusty galaxies Arp200; SMM J2135-0102, “The Cosmic Eyelash” at z  =  2.3 (Ivison et al. 2010) or H-ATLAS J142413.9+022304 alias G15.141 at z  =  4.23 (Cox et al. 2011), which are characteristic of the local submillimeter galaxies. This procedure would show us how much the redshifts and flux densities change depending on the SED used and would help us to select sources in a more robust way. Nevertheless, we think that the Pearson et al. (2013) model is the best model choice since it is not a certain SED from a single source but a empirical model based on 40 H-ATLAS sources and it has already been used in several previous works (Eales 2015; Ivison et al. 2016; Bianchini et al. 2016, 2018; Negrello et al. 2017; Fudamoto et al. 2017; Bakx et al. 2018; Donevski et al. 2018). This further analysis would also allow to check whether the methodology bias observed in simulations is because of the particular shape of the function employed.

Acknowledgments

We acknowledge partial financial support from the Spanish “Ministerio de Economía y Competitividad” (MINECO) project AYA2015-64508-P. DH also acknowledges support from the RADIOFOREGROUNDS project, funded by the European Comission’s H2020 Research Infrastructures under the Grant Agreement 687312. JGN and LB also acknowledge financial support from the I+D 2015 project AYA2015-65887-P (MINECO/FEDER). JGN also acknowledges financial support from the Spanish MINECO for a “Ramon y Cajal” fellowship (RYC-2013-13256). We acknowledge Santander Supercomputacion support group at the University of Cantabria who provided access to the supercomputer Altamira Supercomputer at the Institute of Physics of Cantabria (IFCA-CSIC), member of the Spanish Supercomputing Network, for performing simulations and analyses. The Herschel-ATLAS is a project with Herschel, which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. The H-ATLAS website is http://www.h-atlas.org/. GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalog is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programs including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website is http://www.gama-survey.org/. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the US Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

References

Appendix A: Selected high redshift sources

We present here two tables: the first one with the 20 first objects out of the 370 high-z H-ATLAS sources from the robust high-z sample and the second one containing the 20 first objects out of the 237 high-z H-ATLAS sources from the 500 μm-riser sample.

Appendix B: Postage stamps of candidate lensed galaxies

We show here the postage stamps of the 40 candidate lensed galaxies selected after cross-matching the robust high-z subsample of 283 high-z candidates with a sample of known-low redshift galaxies from SDSS DR14.

Fig. B.1. SDSS images (47 × 47 arcsec2) of part of the 40 associations found between the robust high-z subsample of 283 high-z candidates and sample of 1 776 242 known-low-redshift galaxies from SDSS DR14. H-ATLAS sources and SDSS associations are marked with a pink and a green square, respectively. Einstein rings according to the Einstein radii estimated in Table 4 are drawn as red circles.

Fig. B.2. SDSS images (47 × 47 arcsec2) of part of the 40 associations found between the robust high-z subsample of 283 high-z candidates and sample of 1 776 242 known-low-redshift galaxies from SDSS DR14. H-ATLAS sources and SDSS associations are marked with a pink and a green square, respectively. Einstein rings according to the Einstein radii estimated in Table 4 are drawn as red circles.

Fig. B.3. SDSS images (47 × 47 arcsec2) of part of the 40 associations found between the robust high-z subsample of 283 high-z candidates and sample of 1 776 242 known-low-redshift galaxies from SDSS DR14. H-ATLAS sources and SDSS associations are marked with a pink and a green square, respectively. The object J145420.6−005203 is identified as a QSO in the H-ATLAS catalog. Einstein rings according to the Einstein radii estimated in Table 4 are drawn as red circles.

All Tables

All Figures

Fig. 1. Difference between the redshift recovered with our MMF method (zout) and the input redshift (zin) as a function of zin. 5000 simulations are run for each zin value in the range 1 ≤ zin ≤ 7. The mean value (zout) is computed and shown as blue dots. The solid line shows a polynomial fit to these mean values. This fit is later used to correct for this bias. zout − zin errors (1σ) are plotted as a shaded region.

In the text

Fig. 2. Difference, Δz/(1 + zspec), as a function of zspec between the biased (top) or unbiased (bottom) photometric redshifts estimated with our matched multifilter (MMF) and the spectroscopic redshifts from sources in Table 1. The statistical parameters noted illustrate the systematic overestimates or underestimates, mean μ and median μ1/2, and the degree of scatter, σ, of the photometric redshifts ($z_{\mathrm{MMF}}^{\mathrm{biased}}$ and $z_{\mathrm{MMF}}^{\mathrm{unbiased}}$) found using the MMF and the Pearson et al. (2013) SED. Error bars in the top panel are calculated from using Eqs. (1) and 3) while error bars in the bottom panel are derived from simulations described in Sect. 4. The gray band depicts the 1σ region centered in the median μ1/2 value for each distribution.

In the text

	Fig. 3. Flux measurements in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels derived with the MMF versus the corresponding flux measurements from the H-ATLAS catalog in logarithmic scale for the 32 known-spectroscopic sources from Table 1. The dashed line marks perfect correlation. All fluxes are in units of mJy/beam according to the beam profile of the respective channel.
In the text

Fig. 4. Evolution of the color–color diagram of the H-ATLAS sources studied as cuts are applied in order to get a sample of robust high-z candidates. The dashed pink line is the track of the Pearson et al. (2013) SED for redshifts in the range [0.5,4.5] (in increasing order from the top-right to the bottom-left corner). The vertical lateral coloured bar present in all plots is a scale of the S/N of the sources exhibited, achieved with the MMF.

In the text

	Fig. 5. Redshift distribution, according to the redshift estimates obtained with the MMF method, of the 410 997 sources from H-ATLAS selected after removing the non-point objects and the ones identified as stars.
In the text

Fig. 6. Normalized SED from Pearson et al. (2013), as defined in Eq. (3), at z = 2 in contrast with the normalized tabulated fluxes at 250, 350 and 500 μm (vertical dotted) of two H-ATLAS sources at zMMF = 2: one (J144102.9+012805) that fits well to the model and other one (J233138.6−354345), whose points have been slightly displaced in x-axis to get better clarity, that does not fit properly to the model according to our criteria.

In the text

	Fig. 7. Redshift distribution, according to the redshift estimates obtained with the MMF method, of the 370 high-z H-ATLAS sources from the robust high-z sample selected after imposing all cuts enumerated in Sect. 6.1.
In the text

Fig. 8. Flux measurements in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels after our MMF has been applied versus the corresponding tabulated flux measurements from the H-ATLAS catalog in logarithmic scale for the 370 high-z H-ATLAS from the robust high-z sample (see Appendix A). The linear behavior with zero-intercept is drawn with a black dashed line. All fluxes are in units of mJy/beam according to the beam profile of the respective channel.

In the text

	Fig. 9. S/N in the filtered image after our MMF has been applied versus the S/Ns in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels in logarithmic scale for the 370 high-z H-ATLAS sources from the robust high-z sample. The linear behavior with zero-intercept is drawn with a black dashed line.
In the text

Fig. 10. Flux measurements in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels after our MMF has been applied versus the corresponding tabulated flux measurements from the H-ATLAS catalog in logarithmic scale for the 237 high-z H-ATLAS from the 500 μm-riser sample (see Appendix A). The linear behavior with zero-intercept is drawn with a black dashed line. All fluxes are in units of mJy/beam according to the beam profile of the respective channel.

In the text

	Fig. 11. S/N in the filtered image after our MMF has been applied versus the S/Ns in the 250 (left), 350 (middle) and 500 μm (right) SPIRE channels in logarithmic scale for the 237 high-z H-ATLAS from the 500 μm-riser sample. The linear behavior with zero-intercept is drawn with a black dashed line.
In the text

	Fig. 12. Ratio between the number of matches found (Nm(R)), after cross-matching the full high-z subsample of 2828 H-ATLAS sources with a sample of 1 776 242 known-redshift galaxies from SDSS DR14, and the number of matches expected (Nr(R)) from a random distribution (Eq. (7)) for different aperture radii R around the SDSS sources.
In the text

	Fig. 13. Ratio between the number of matches found (Nm(R)), after cross-matching the robust high-z subsample of 283 candidates with a sample of 1 776 242 known-redshift galaxies from SDSS DR14, and the number of matches expected (Nr(R)) from a random distribution (Eq. (7)) for different aperture radii R around the SDSS sources.
In the text

	Fig. B.1. SDSS images (47 × 47 arcsec2) of part of the 40 associations found between the robust high-z subsample of 283 high-z candidates and sample of 1 776 242 known-low-redshift galaxies from SDSS DR14. H-ATLAS sources and SDSS associations are marked with a pink and a green square, respectively. Einstein rings according to the Einstein radii estimated in Table 4 are drawn as red circles.
In the text

	Fig. B.2. SDSS images (47 × 47 arcsec2) of part of the 40 associations found between the robust high-z subsample of 283 high-z candidates and sample of 1 776 242 known-low-redshift galaxies from SDSS DR14. H-ATLAS sources and SDSS associations are marked with a pink and a green square, respectively. Einstein rings according to the Einstein radii estimated in Table 4 are drawn as red circles.
In the text

Fig. B.3. SDSS images (47 × 47 arcsec2) of part of the 40 associations found between the robust high-z subsample of 283 high-z candidates and sample of 1 776 242 known-low-redshift galaxies from SDSS DR14. H-ATLAS sources and SDSS associations are marked with a pink and a green square, respectively. The object J145420.6−005203 is identified as a QSO in the H-ATLAS catalog. Einstein rings according to the Einstein radii estimated in Table 4 are drawn as red circles.

In the text