The Herschel-ATLAS lensed galaxies sample

 This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 707601
Building the sample
I searched for strong lensing events over the entire H-ATLAS area (Negrello et al. 2017), using the simple selection in flux density at 500μm described here.
I first identified all the sources with F500 ≥ 100mJy in the H-ATLAS fields. A Herschel/SPIRE color-color plot of the sample, shown in Fig. 1, gives a crude indication of the redshift of the sources: galaxies with z > 1 have their dust peak redshifted above the 250μm Herschel/SPIRE band compared to lower  redshift galaxies, and, therefore, display lower values of the 250μm to 350μm flux density ratio. The figure shows that the distribution of the F250/F350 values is bimodal. This is because the sample comprises two populations of sources at very different redshifts:
Figure 1. Herschel/SPIRE colour–colour diagram of all the sources with F500 ≥ 100 mJy identified in the H-ATLAS fields. The colours are indicative of the redshift of the sources and allow us to differentiate between local galaxies (blue squares) and high redshift sources, i.e. candiate lensed galaxies (red dots) and blazars (green triangles). The black line is the track of the Pearson et al. (2013) empirical template for redshifts in the range [0.5,4.5] in steps of 0.5 (in increasing order from the top-right to the bottom-left corner), as marked by the black dots [from Negrello et al. (2017)]
  A population of nearby (z < 0.1) galaxies with F250/F350 > 1.5. A first visual inspection of available optical images showed that such sources are extended, with a disk and/or spiral arms clearly visible in most cases. Then, a quick search in the NASA Extragalactic Database confirmed that they are local galaxies with spectroscopic redshifts z < 0.1. 
•  A population of high-redshift (z > 1) galaxies with F250/F350 < 1.5. This comprises the candidate lensed galaxies I am looking for (red dost in the figure) and a small group of blazars (green triangles in the figure). The latter were easily singled out by checking for strong radio emission in shallow radio surveys, like NVSS and FIRST.
I identified 80 candidate lensed galaxies with flux density F500μm≥100mJy over the ~600 deg2 of the H-ATLAS. The are represented by the red dots in Fig. 1. Their spatial distribution within the H-ATLAS fields is shown in Fig. 2.
Currently 21 of them – which includes 80% of those with F500μm > 150mJy – are confirmed to be lensing systems based on available optical/near-infrared to sub-millimetre/millimetre spectroscopic and high-resolution imaging data, while follow-up observations are still ongoing to confirm the nature of the others.
Only one object in our sample, HATLASJ084933.4+021442, has been so far confirmed to not be a strongly lensed galaxy. It is indeed a binary system of Hyper Luminous Infrared Galaxies5 (HyLIRGs) at z = 2.410 reported by Ivison et al. (2013), which is believed to represent the early stage in the formation of the core of a massive galaxy cluster. 
The list of candidate lensed galaxies is available hereand can be exported in either excel or cvs format. 
Figure 2. Herschel/SPIRE maps of the H-ATLAS fields with the position of the 80 candidate lensed galaxies marked in yellow [from Negrello et al. (2017)].
Properties of the sample
•  Number counts
There is an apparently striking variation in the number density of candidate lensed galaxies between the different H-ATLAS fields, that can be noticed by looking at Fig. 2. In order to investigate this issue, I have computed the integral number counts of the candidate lensed galaxies for each of the H-ATLAS fields. The results are shown in Fig. 3: top panels for the three equatorial GAMA fields; bottom-left and bottom-middle panels for the NGP and SGP fields, respectively. Error bars and upper limits correspond to the 68 per cent confidence interval, assuming Poisson statistics.
The GAMA 9h field and the NGP field appear to be particularly rich in lensed galaxies with F500 ≥ 200 mJy. On the contrary, the SGP field is short of such bright sources. In fact, only one candidate lensed galaxy with F500 ≥ 200 mJy is found in the SGP field, while 13 have been identified in the other fields (most of which are already confirmed lensed galaxies).
In the bottom-right panel of the same figure, I show the comparison between the number counts extracted from the SGP field and those derived by combining the three equatorial fields with the NGP field. In this case, the error bars and the upper limits have been drawn to represent the 95 per cent confidence interval. It is clear that, despite the large difference in number density for F500 150 mJy, the two samples are still consistent with each other at the ∼2σ level. The observed fluctuations are also consistent with the results of simulations based on the modelled lensed number counts (red shaded region in the bottom rightmost panel).
Figure 3.  Integral number counts of the candidate lensed galaxies with F500 ≥ 100mJy identified in the five H-ATLAS fields: GAMA9h (top-left), GAMA12h (top-middle), GAMA 15h (top-right), NGP (bottom-left) and SGP (bottom-middle). The bottom-right panel shows the comparison between the number counts measured in the SGP field and those derived by combining the three GAMA fields with the NGP fields. Error bars and upper limits correspond to the 68 per cent confidence interval assuming Poisson statistic (Gehrels 1986), with the exception of the bottom-right panel where they are shown at the 95 per cent confidence level. In all panels, the solid black curve is the prediction for the abundance of unlensed DSFGs at z > 1.5, based on the model of Cai et al. (2013). The red curves show the predicted number counts of lensed galaxies for different values of the maximum magnification, μmax, experienced by the background sources, from 10 to 30. The mass profile of the lenses is assumed to be the superposition of an NFW profile (for the dark-matter component) and a de Vaucouleurs profile (accounting for baryons), as described in Lapi et al. (2012). The shaded region in the bottom-right panel illustrates the variation in the number density of lensed galaxies based on a set of 100 simulations performed over an area of 300 deg2 and assuming μmax = 12.5 [from Negrello et al. (2017)].
Figure 4. Integral number counts of the candidate lensed galaxies with F500 ≥ 100 mJy identified in all the H-ATLAS fields (purple dots). Error bars correspond to the 95 per cent confidence interval. The meaning of the curves is the same as in Fig. 3. The squares represent measurements of the number density of unlensed DSFGs derived from P(D) analysis (Glenn et al. 2010) [from Negrello et al. (2017)].
Once I combine the samples from all the five H-ATLAS different fields I obtain the number counts shown in Fig. 4. The number density is in agreement with expectations based on a physical model for the unlensed DSFG population (Cai et al. 2013) coupled with a NFW plus de Vaucouleurs model for the mass profile of the lenses (Lapi et al. 2012), assuming a maximum magnification in the range 10 to 15. These values of are in agreement with the magnification factors estimated from the lens modelling and source reconstruction of individual lensed DSFGs (Enia et al. 2018). 
An excel version of these counts and of the corresponding differential counts can be downloaded from here together with model predictions.
•  Redshifts
Most of the (lensed) sources in the sample do not have a spectroscopic redshift yet. Follow-up observations are currently ongoing to measure their redshifts via the detection of carbon monoxide (CO) emission lines. For now I have used the empirical template of Pearson et al. (2013) to derive a photometric redshift. 
Pearson et al. (2013) have used a sub-set of 40 H-ATLAS sources with spectroscopic redshifts in the range 0.5 < z < 4.2 to construct an empirical SED for high-redshift Herschel-selected sources to be used to estimate photometric redshifts from Herschel/SPIRE data. The template is the sum of two modified blackbody spectra with temperatures Tcold = 23.9 K and Thot = 46.9 K, and dust emissivity index fixed to β = 2. The ratio between the normalisation of the two components is 30.1. The track shown in Fig. 1 was derived using this template. The errors on the photometric redshifts are calculated as 0.12 × (1 + zphot), 0.12 being the rms scatter in the (zphot − zspec)/(1 + zspec) values as measured by Pearson et al. for sources with zspec > 1.
Figure 5. Redshift distribution of the candidate lensed galaxies with F500 ≥ 100 mJy. Upper panel: distribution of the photometric redshift sample (orange) and of the spectroscopic redshift sample (blue). Lower panel: distribution of the photometric+spectroscopic redshift sample (dots with error bars). The shaded yellow region represents the 68 per cent confidence interval assuming Poisson statistics. The red histograms are predictions based on the galaxy formation and evolution model of Cai et al. (2013) coupled with the lensing formalism of Lapi et al. (2012), assuming different maximum magnifications: μmax = 10, 15 [from Negrello et al. (2017)].
In Fig. 5 I show the redshift distribution of the candidate lensed galaxies with F500 ≥ 100 mJy, with photometric redshifts replaced by spectroscopic ones where available. I haven’t included the Ivison et al. (2013) source in this plot.
The distribution is derived by performing 10000 simulations. Each time, the redshift of the sources is resampled at random, from a Gaussian probability distribution with a mean equal to the measured photometric/spectroscopic redshift and dispersion equal to the associated error. The simulated redshifts are then binned into a histogram and the mean value in each bin is taken as the estimate of the number of objects in that bin. Error bars, corresponding to the 68 per cent confidence interval, are derived following the prescriptions of Gehrels et al. (1986) for Poisson statistics.
The galaxies in the sample have redshifts ranging from ~1 to ~4.5 with a median value zmedian = 2.53. Similar median redshifts are obtained for the spectroscopic redshift sample and the photometric redshift sample alone, i.e. zmedian = 2.49 and 2.54, respectively (the two distributions are shown in the upper panel of Fig. 3).
•  Infrared Luminosities and Dust Temperatures
Figure 6. Dust temperature versus total IR luminosity(8−1000μm, not corrected for lensing) of the candidate lensed galaxies extracted from the H-ATLAS fields, excluding the Ivison et al. (2013) source. Squares represent objects with spectroscopic redshift while the dots indicate objects with photometric redshift. Error bars correspond to the 68 per cent confidence interval and take into account uncertainties on both the redshifts and the Herschel photometry. The IR luminosity and the dust temperature are derived assuming an optically thin greybody model with dust emissivity index fixed to β = 1.5 [from Negrello et al. (2017)].
I used the spectroscopic and photometric redshift information to derive the total IR luminosity, LIR (integrated over the rest-frame wavelength 8 to 1000μm), of the candidate lensed galaxies. I assumed a single temperature, optically thin, modified blackbody spectrum with dust emissivity index β = 1.5 and fit it to the Herschel photometry of each source. The IR luminosity and dust temperature, Tdust, are kept as free parameters.  In order to account for the uncertainty on the photometric/spectroscopic redshift, as well as that on the photometry, I performed for each source 1000 simulations by resampling at random the distributions of redshift values and measured flux densities, both assumed to be Gaussian. The final values of LIR and Tdust are shown in Fig. 6, and were obtained as the median of the derived distributions of 1000 best-fitting values.
All the galaxies in the sample have infrared luminosities in excess of 1013 Lsun, and even 1014 Lsun in few cases, thus making them the apparently brightest star-forming galaxies in the Universe. However, these estimates are not corrected for the effect of lensing: with an expected (or measured, where possible) magnifications μ~5-15, the galaxies in the sample are more likely to be intrinsically ultra luminous infrared galaxies (ULIRGs: 1012 ≤ LIR/L ≤ 1013).
The median value of the dust temperature is Tdust = 34.6 K, consistent with what was found for other sub-mm/mm-selected unlensed (Magnelli et al. 2012) and lensed (Weiß et al. 2013; Canameras et al. 2015; Nayyeri et al. 2016) galaxies.
•  Magnifications and Physical Sizes
In Enia et al. (2018) we performed lens modelling and source reconstruction of Sub-millimetre Array (SMA) data for 13 of the 500μm-selected lensed galaxies in the H-ATLAS sample. We adopted a pixellated source brightness distribution, following the semi-linear inversion method of Warren & Dye (2003) with the adaptive pixel scale scheme by Nightingale & Dye (2015) [the method is described here], and we extended it to work in the Fourier space of interferometry. Some examples of the source reconstruction are shown in Fig. 7.
Figure 7. Examples of source reconstruction of H-ATLAS lensed galaxies performed on SMA data. From left to right: input SMA image (created using a natural weighting scheme); reconstructed image; reconstructed background source with contours at 3σ (black curve) and 5σ (white curve). The caustics and the critical lines are shown in brown and in red, respectively. The white hatched ellipse in the bottom left corner of the leftmost panels represents the SMA synthesized beam [from Enia et al. (2018)]. 
The modelling allowed us to constraint the size of dust emitting region in the source plane and to estimate the magnification factor, μ. 
We defined an effective radius of the reconstructed sources based on the area in the source plane where emission is detected above 5σ. We also fit the reconstructed source surface brightness with an elliptical Gaussian model. We derive a median value reff ∼ 1.77 kpc and a median Gaussian FWHM ∼ 1.47 kpc. The magnification factors are found to be in the range ~3 to ~10.
After correction for magnification, our sources have intrinsic star formation rates SFR ∼ 900–3500 Msun yr−1 , resulting in a median SFR surface density SFR ∼ 132 Msun yr−1 kpc−2 (or SFR ∼ 218 Msun yr−1 kpc−2 for the Gaussian fit). This is consistent with that observed for other star-forming galaxies at similar redshifts, and is significantly below the Eddington limit for a radiation pressure regulated starburst.
Final remarks
The boost in brightness due to gravitational lensing is what makes strongly lensed sources so valuable in probing the properties and physical conditions of the interstellar medium in high redshift star-forming galaxies. In fact, spectroscopic follow-up observations of such bright targets are less time demanding, allowing the swift detection of emission lines from e.g. carbon monoxide (Lupu et al. 2012; Harris et al. 2012; Yang et al. 2017), water vapor (Omont et al. 2011; Omont et al. 2013; Yang et al. 2016) and dense molecular gas tracers HCN, HCO+ and HNC (Oteo et al. 2017).
Thanks to the increase in spatial resolution offered by gravitational lensing, our sample is also ideal for studying the morphological and dynamical properties of z ∼ 2 DSFGs down to sub-kpc scales, as demonstrated, for example, by the analysis of the high-resolution observations of HATLAJ090311.6+003907 (aka SDP.81) obtained with the Atacama Large Millimeter Array (e.g. Dye et al. 2015 and references therein) as part of its long baseline campaign (ALMA partnership 2015).