THE HAWAII SCUBA-2 LENSING CLUSTER SURVEY: NUMBER COUNTS AND SUBMILLIMETER FLUX RATIOS

In order to estimate the number of fake sources contaminating the number counts that we measured from our science maps, we need to generate source-free maps with only pure noise for each of our fields. These maps are sometimes referred to as jackknife maps in the literature. Following Chen et al. (2013a, 2013b), we subtracted two maps that were each produced by coadding roughly half of the flux-calibrated data. In doing so, the real sources are subtracted off, and the residual maps are source-free maps. We then rescaled the value of each pixel by a factor of $\sqrt{{t}_{1}\times {t}_{2}}/({t}_{1}+{t}_{2})$, with t₁ and t₂ representing the integration time of each pixel from the two maps. Finally, we applied the matched-filter with the same procedure for the science maps.

We have shown in previous work that sources detected above a 4σ level have a low contamination rate (Chen et al. 2013a, 2013b). However, for computing number counts, we can use a lower detection threshold where there are still significantly more true sources than false detections. In Chen et al. (2013b), we extracted sources down to ∼2σ. However, here we use 3σ as our detection threshold.⁸ We experimented with different detection thresholds and binning in order to extend the faint end of the counts while keeping good signal-to-noise ratio (S/N) throughout all the flux bins at the same time. We found that using 3σ leads to better S/N in the number counts than using lower thresholds, especially at the faint end.

We first generated the PSFs by averaging all the primary calibrators, the ones we used for deriving the FCFs. Following the methodology of source extraction in Chen et al. (2013a, 2013b), we identified the pixel with the maximum S/N, subtracted this pixel and its surroundings using the PSF centered and scaled at the position and value of this pixel, and then searched for the next maximum S/N. We iterated this process until the 3σ threshold was hit. We ran the source extraction on both the science maps and and the pure noise maps. The ratio of the total number of sources from the pure noise maps and from the science maps (=N_false/N_total) is ∼55% (22%) at 450 (850)μm. The effect of false sources is subtracted in the computation of number counts, as we will describe in Sections 3.4 and 3.5.

In order to obtain the intrinsic flux of a lensed source, both the lens model of the cluster and the redshift of the source are required. However, since we do not have redshift measurements for individual sources, we simply adopted estimated median redshifts of the lensed 450 μm and 850 μm sources to compute our de-lensed number counts. Note that what we need to estimate are the "observed median redshifts" of the lensed populations, which would be higher than the real median redshifts of the distributions (e.g., Weiß et al. 2013). This is because sources at higher redshifts have higher probability of being lensed and generally have higher lensing magnifications, causing a selection bias (e.g., Hezaveh & Holder 2011). We leave the discussion of the blank-field sources and their redshift distributions to Section 4.2.

We estimated the two median redshifts by exploring the flux ratios between 450 μm and 850 μm for all the 4σ detected sources in the cluster fields. At 450 (850)μm, we took the flux and position of a detected source and then measured the flux value at the same position on the 850 (450)μm map. In Figure 1, we plot the 850 μm-to-450 μm and 450 μm-to-850 μm flux ratios against the product of S/N at 450 μm and at 850 μm. In each bin of (S/N)_{450 μm}  × (S/N)_{850 μm}, we took the median of the flux ratios and calculated the error using bootstrapping. We can see that the flux ratio increases with increasing S/N product and then flattens. The lower (negative) measured flux ratios at lower (negative) (S/N)_{450 μm} ×(S/N)_{850 μm} are a result of the mismatch between the positions of the 450 μm and 850 μm flux peaks due to lower S/N. We compared the measured flux ratios with what we measured from simulated maps, which were produced by populating the pure noise maps with sources with constant flux ratios. A detailed description of how we performed such simulations is left to the Appendix. In Figure 1, red dashed lines are the input constant flux ratios of the simulations, and the red solid lines are what we measured from the simulated maps, which show good agreement with our measurements from the science maps. We therefore conclude that the values of the two dashed lines correspond to the median flux ratios of the 450 μm and 850 μm selected populations in the cluster fields.

Zoom In Zoom Out Reset image size

To convert flux ratios to redshifts, we assumed a modified blackbody spectral energy distribution (SED) of the form S_ν ∝ (1 − e^−τ(ν))B_ν(T), where τ(ν) = (ν/ν₀)^β and ν₀ = 3000 GHz. Assuming β = 1.5, we determined the redshifts from our estimated median flux ratios for a dust temperature of 30 K, 40 K, or 50 K. The corresponding redshifts for the 450 μm sources are z ∼ 1.5, 2.2, or 2.8. At 850 μm, we obtained z ∼ 2.0, 2.8, or 3.5. The final number counts shown in this paper are based on source plane redshifts of 2.2 and 2.8 for 450 μm and 850 μm, respectively. We chose these values because they are the central values of the different SED models used. However, we will show in Section 3.6 that using z = 1.5, 2.8 at 450 μm and z = 2.0, 3.5 at 850 μm does not change our results significantly, and that the computation of the number counts is not sensitive to the adopted source plane redshifts.

To compute the de-lensed, differential number counts, we corrected all the source fluxes in the cluster fields using the publicly available software LENSTOOL (Kneib et al. 1996), which allows us to generate magnification maps with the angular sizes of our SCUBA-2 maps. We therefore used the lens models from the LENSTOOL developers (CATS team) for A1689 (Limousin et al. 2007), A2390 (Richard et al. 2010), MACS J1423.8+2404 (Limousin et al. 2010), and the three Frontier Fields (Hubble Frontier Field archive⁹ ). For each source from a science or pure noise map, we calculated its number density by inverting the detectable area, which is the area in which this source can be detected above the 3σ threshold. For a source in a cluster field, the detectable area is defined on the source plane. We then computed the number counts by summing up the number densities of the sources in each flux bin with errors based on Poisson statistics (Gehrels 1986). Finally, we subtracted the counts of the pure noise maps from the counts of the science maps to produce the pure source counts.

While the discrepancy in the magnifications between different lens models can be a factor of a few at the cluster center, the effect on the measured number counts is not significant. This is the same as the effect caused by the different source plane redshifts, as we discussed in Section 3.3. Although there are uncertainties in the lens models, the source plane redshifts, and the positions of the submillimeter sources, the de-lensed flux and detectable area of a source are directly related, causing little change in the slope and normalization of the measured number counts.

The pure source counts we computed above, however, still do not represent the true underlying submillmeter populations because the fluxes of the sources are boosted and there is incompleteness. The cause of the flux boost is the statistical fluctuations of the flux measurements for flux-limited observations, known as the Eddington bias (Eddington 1913). Following Chen et al. (2013a, 2013b), we ran Monte Carlo simulations to estimate the underlying count model at each wavelength. We first generated a simulated map by randomly populating sources in the pure noise map, drawn from an assumed model and convolved with the PSF. The count model we used is in the form of a broken power law

$$\begin{eqnarray}&&\tfrac{{dN}}{{dS}}=\left\{\begin{array}{l}{N}_{0}{\left(\tfrac{S}{{S}_{0}}\right)}^{-\alpha }\,{\rm{if}}\,S\leqslant {S}_{0}\\ {N}_{0}{\left(\tfrac{S}{{S}_{0}}\right)}^{-\beta }\,{\rm{if}}\,S\gt {S}_{0}.\end{array}\right.\end{eqnarray} \tag{ 1 }$$

For the cluster fields, we populated the sources in the source plane and projected them onto the image plane using LENSTOOL.

For each simulated map we reran our source extraction and computed the recovered counts using the same method and flux bins used for the science map. We repeated the simulation 50 times for each input model and then averaged the recovered counts from these simulations. In order to measure the actual counts we adopted an iterative procedure. Using the ratios between the averaged recovered counts and the input counts, we renormalized the observed raw counts in each bin from the science map. We then did a χ² fit to the corrected observed counts using a broken power law. This fit was then used as the input model for the next iteration. We continued until the corrected counts matched the corrected counts of the previous iteration within 1σ throughout all the flux bins. It took only two or three iterations to converge for each field.

We show the corrected number counts for all the fields together at both 450 μm and 850 μm in Figure 2. Thanks to the lensing magnification, we are able to detect counts down to fluxes fainter than 1 mJy and 0.2 mJy in several fields at 450 μm and 850 μm, respectively. The solid lines represent the best-fit broken power law models for the counts. In each panel, we also show the best-fit model for the counts computed with a lower source plane redshift (z = 1.5 at 450 μm and z = 2.0 at 850 μm) with the dotted line and the best-fit model to the counts computed with a higher source plane redshift (z = 2.8 at 450 μm and z = 3.5 at 850 μm) with the dashed line. We can see that the results are not very sensitive to the assumed redshifts.

Zoom In Zoom Out Reset image size

In order to better constrain our count model at each wavelength, we combined the counts from all the fields. The results are shown in Figure 3, which are weighted averages of the corrected counts from each field (black circles). We assigned a weight for each flux bin of a field in the following way. For each field, we used the final count model we obtained (Section 3.5) to run the same simulation 50 times and obtained the 1σ scatter of the recovered counts in each flux bin. We then normalized this scatter by the average of the recovered counts. The inverse square of the scatter is adopted as the weight. We also show various results from the literature. The combined number counts and the best-fit parameters of the models are summarized in Tables 2 and 3.

Zoom In Zoom Out Reset image size

Table 2.  Combined Differential Number Counts at 450 μm and 850 μm

S450 μm	dN/dS	S850 μm	dN/dS
(mJy)	(mJy−1 deg−2)	(mJy)	(mJy−1 deg−2)
0.94	${40579}_{-13496}^{+17834}$	0.16	${468599}_{-147442}^{+215679}$
4.63	${1438}_{-510}^{+741}$	0.23	${217910}_{-70715}^{+105487}$
8.77	${263.9}_{-79.2}^{+105.7}$	0.55	${33138}_{-9975}^{+15129}$
14.45	${76.58}_{-15.96}^{+18.52}$	0.85	${9650}_{-3444}^{+4984}$
19.76	${30.53}_{-9.39}^{+10.02}$	1.27	${4576}_{-826}^{+1114}$
24.05	${13.55}_{-6.09}^{+6.17}$	1.87	${2646}_{-309}^{+345}$
34.53	${2.17}_{-1.10}^{+1.72}$	2.56	${1209}_{-129}^{+140}$
		3.63	${552.1}_{-44.5}^{+47.6}$
		4.96	${238.6}_{-24.8}^{+31.3}$
		6.06	${155.4}_{-16.9}^{+18.8}$
		7.85	${54.52}_{-13.99}^{+21.06}$
		8.93	${24.16}_{-4.59}^{+6.12}$
		11.14	${12.88}_{-4.84}^{+7.79}$
		15.52	${5.31}_{-1.79}^{+3.02}$

Download table as:  ASCIITypeset image

Table 3.  Best-fit Broken Power Laws for the Combined Differential Number Counts at 450 μm and 850 μm

Wavelengths	N0	S0	α	β
(μm)	(mJy−1 deg−2)	(mJy)
450	33.3${}_{-18.1}^{+44.0}$	20.1${}_{-5.1}^{+6.8}$	2.34${}_{-0.17}^{+0.14}$	5.06${}_{-1.52}^{+5.03}$
850	342${}_{-80}^{+56}$	4.59${}_{-0.38}^{+0.26}$	2.12${}_{-0.07}^{+0.06}$	3.73${}_{-0.47}^{+0.59}$

Download table as:  ASCIITypeset image

Semi-analytical simulations have shown that source blending could impact the number counts obtained from single-dish observations (Hayward et al. 2013a; Cowley et al. 2015). Some recent studies with ALMA observations have also discussed the effect of multiplicity on the number counts (Hodge et al. 2013; Karim et al. 2013; Simpson et al. 2015). In Chen et al. (2013b), we used the SMA detected sample in CDF-N (Barger et al. 2014) to obtain the multiple fraction as a function of flux at S_{850 μm} > 3.5 mJy, and we computed the multiplicity-corrected CDF-N 850 μm number counts above 3.5 mJy, assuming that all the blends split into two equal components. There is one incorrect coefficient in Equation (4) of Chen et al. (2013b), which should instead be written as

$$\begin{eqnarray}\displaystyle \frac{{{dN}}_{{\rm{corr}}}(S)}{{dS}} & = & \displaystyle \frac{{{dN}}_{{\rm{orig}}}(S)}{{dS}}\times (1-{f}_{{\rm{mul}}}(S))\\ & & +{f}_{{\rm{mul}}}(2S)\times 2\times \displaystyle \frac{{{dN}}_{{\rm{orig}}}(2S)}{{dS}},\end{eqnarray} \tag{ 2 }$$

where f_mul is the multiple fraction of the SMA detected SCUBA-2 sources as a function of flux, and dN_corr/dS and dN_orig/dS are the multiplicity-corrected and the original SCUBA-2 counts, respectively.¹⁰ However, this correction does not significantly change the result of Chen et al. (2013b). The systematic changes introduced by multiplicity are still smaller than the statistical errors of the counts.

Computing multiplicity corrections is difficult because the multiple fractions at different fluxes are still not well determined. For SCUBA-2 selected sources, Simpson et al. (2015) found a multiple fraction of ${61}_{-15}^{+19}$% (17 out of 28) at S_{850 μm} > 4 mJy using ALMA, while Barger et al. (2014) found a multiple fraction of only ${12.5}_{-6.8}^{+12.1}$% (3 out of 24) at S_{850 μm} > 3.5 mJy using SMA. The much lower multiple fraction from Barger et al. (2014) can be caused by the sensitivities of their SMA maps, which only allow detections of secondary SMGs brighter than 3 mJy. The multiple fraction is simply sensitive to the depth of the follow-up interferometric observations. However, Chen et al. (2013b) showed that most of the sources with a single SMA detection in CDF-N have flux measurements that statistically agree with those made by SCUBA-2. Using a larger SMA detected sample in CDF-N (31 4σ detected sources; L. L. Cowie et al. 2016, in preparation), we again compare the fluxes measured by SCUBA-2 and by SMA for the sources with a single SMA detection. We show the comparison in Figure 4. The two fluxes statistically agree with each other for most of the sources. The median flux ratio of SMA to SCUBA-2 is 0.96 ± 0.06. This suggests that most secondary sources that are missed by SMA would be faint and unlikely to affect the bright-end counts. These sources would be unlikely to increase the faint-end counts significantly, either, since they contribute a small fraction of the faint sources based on the broken power law model.

Zoom In Zoom Out Reset image size

Instead of computing the multiplicity corrections by assuming multiple fractions, here we use another approach to examine the effect of multiplicity on the number counts. We took the SMA detected sample in CDF-N (L. L. Cowie et al. 2016, in preparation) and their corresponding SCUBA-2 sources to compute two sets of number counts. We corrected the SCUBA-2 counts for Eddington bias using simulations. For the SMA sources, we simply took their fluxes and computed the detectable areas and number counts as if they were detected in our SCUBA-2 map. Note that these sources comprise a incomplete sample at S_{850 μm} > 3 mJy because only 29 out of 81 SCUBA-2 sources above this flux (19 out of 29 at S_{850 μm} > 6 mJy) have SMA observations. We did not apply any completeness correction since we are only attempting to see whether there is a significant difference between these two sets of counts. The comparison is shown in Figure 5. We can see small deviations at both the faint and bright ends, but the two sets of counts are essentially consistent within their uncertainties.

Zoom In Zoom Out Reset image size

Although Simpson et al. (2015) found a multiple fraction of 61% at S_{850 μm} > 4 mJy, their differential counts constructed from ALMA and SCUBA-2 still statistically agree with each other at S_{870 μm} ≲ 15 mJy (see their Figure 6). The effect of multiplicity is more obvious in the cumulative counts at S_{870 μm} ≳ 10 mJy. If we plot the cumulative EBL, the ALMA and SCUBA-2 results of Simpson et al. (2015) would deviate at S_{870 μm} > 5 mJy, which is well above the confusion limit of the JCMT. As we will discuss in Section 4.1, the majority of the EBL comes from sources fainter than our detection limits at both wavelengths, and this conclusion is not affected by multiplicity.

Zoom In Zoom Out Reset image size

We plot the cumulative EBL as a function of flux at both wavelengths in Figure 6. Without gravitational lensing, our surveys would be limited to sources with S_{450 μm} ≳ 10 mJy or S_{850 μm} > 2 mJy. However, sources fainter than these limits contribute the majority of the EBL at both wavelengths. If we use the measurement by Fixsen et al. (1998) as the total EBL, ∼90% of the 450 μm background comes from sources fainter than 10 mJy and ∼80% of the 850 μm background comes from sources fainter than 2 mJy. These numbers would not change significantly even if we consider the effect of multiplicity. Our result also suggests there is at least ∼50% of the EBL with S_{450 μm} < 1.0 mJy. If we integrate the 450 μm differential count down to the lower flux limit in the upper panel of Figure 2, there is still at least ∼40% of the EBL with S_{450 μm} < 0.7 mJy that is not resolved by our SCUBA-2 maps. Most of these faint SMGs would have L_IR < 10¹² L_⊙, corresponding to LIRGs or normal galaxies. However, note that all these fractions of the EBL we calculate here depend on which measurement of the total EBL we assume, as well as its uncertainty.

We note that the faint-end slopes (α) of the number counts should become less than one at fluxes fainter than 1 mJy at 450 μm and 0.1 mJy at 850 μm. In Figure 7, we show the combined differential number counts (from Figure 3) multiplied by the flux. Because the total EBL equals the integral of dN/dS × S over S, dN/dS × S must turn over at some point such that the derived cumulative EBL does not significantly exceed the total EBL measured by the COBE satellite. Using second order polynomial fits to the log–log plots in Figure 7, the estimated turnovers (where α becomes one) are at ∼0.8 mJy at 450 μm and ∼0.06 mJy at 850 μm. Sources with these fluxes would contribute the most to the EBL. If we again assume a modified blackbody SED with β = 1.5 and a dust temperature of 40 K, S_{450 μm} ∼ 0.8 mJy and S_{850 μm} ∼ 0.06 mJy correspond to L_IR ∼ 1.3 × 10¹¹ L_⊙ at z = 2.2 and ∼3.4 × 10¹⁰ L_⊙ at z = 2.8, respectively.

Zoom In Zoom Out Reset image size

We also show other measurements and model predictions of the EBL from the literature in Figure 6. All of the observational results are consistent with ours within 1σ (note that the 1σ spread of all the other cumulative EBL curves are not shown in Figure 6). There is, however, significant disagreement between the 450 μm model from Béthermin et al. (2012a) and our result. This difference is mainly caused by the discrepancy in the differential number counts at S_{450 μm} ∼ 2–15 mJy (see Figure 3), where the model slightly overproduces sources.

Viero et al. (2013) quantified the fraction of the EBL that originates from galaxies identified in the UV/optical/near-infrared by stacking K-selected sources on various Spitzer and Herschel maps at different wavelengths. They were able to resolve 65% ± 12% of the EBL at 500 μm (2.60 nW m⁻² sr⁻¹ or 0.434 MJy sr⁻¹) based on the measurement by Lagache et al. (2000). If we correct their result using the EBL measured by Fixsen et al. (1998), their sample contributes ∼70% of the EBL at 500 μm (2.37 nW m⁻² sr⁻¹ or 0.395 MJy sr⁻¹), which includes ∼10%, 40%, and 20% coming from normal galaxies, LIRGs, and ULIRGs, respectively. For comparison, we can assume an extreme case, where all of the 450 μm sources lie at z = 2.8 and have a modified blackbody SED with β = 1.5 and T = 50 K (see Section 3.3). In such a case, galaxies with L_IR < 10¹² L_⊙ would have S_{450 μm} ≲ 3.4 mJy, which still contribute ∼75% of the EBL and cannot be fully accounted for by the normal galaxies and LIRGs in Viero et al. (2013). If we assume a lower dust temperature or a lower redshift, galaxies with L_IR < 10¹² L_⊙ would contribute even more to the EBL. This is consistent with recent SMA (Chen et al. 2014) and ALMA (Kohno et al. 2016; Fujimoto et al. 2016) observations, which suggest that many faint SMGs may not be included in the UV star formation history.

Because the majority of sources that contribute the submillimeter EBL have L_IR < 10¹² L_⊙, a full accounting of the cosmic star formation history requires a thorough understanding of the galaxies with FIR luminosities corresponding to LIRGs and normal star-forming galaxies. It is therefore critical to determine how much the submillimeter- and UV-selected samples overlap at this luminosity range. Future work using interferometry is needed to determine the fraction of faint SMGs that is already included in the UV population as a function of submillimeter flux.

As described in Section 3.3, we used a statistical approach to explore the submillimeter flux ratios for both 450 μm and 850 μm selected samples from the cluster fields. If we do the same exercise on the blank-field data and again use a modified blackbody SED with β = 1.5, T = 40 K, the median redshifts of the 450 μm and 850 μm populations would be 2.0 and 2.6, respectively. These are in rough agreement with the median redshifts (2.06 ± 0.10 and 2.43 ± 0.12) from the simulations of Zavala et al. (2014). Béthermin et al. (2015) also presented the median redshift of dusty galaxies as a function of wavelength, flux limit, and lensing selection bias based on their empirical model (Béthermin et al. 2012a). The 4σ detection limits for our blank-field 450 μm, cluster 450 μm, blank-field 850 μm, and cluster 850 μm images are ∼18, 10, 2, and 2 mJy, respectively. According to Figure 3 of Béthermin et al. (2015), these flux cuts correspond to median redshifts of z ∼ 1.9, 1.8 ≲ z ≲ 2.4, z ∼ 2.4 and 2.4 ≲ z ≲ 2.8, respectively. These also roughly agree with our estimated median redshifts. Note that for our cluster fields, the corresponding redshifts are shown as intervals, because the relations in Figure 3 of Béthermin et al. (2015) are for all galaxies and "strongly lensed" galaxies, while our SCUBA-2 maps extend out to radii of ∼6' and therefore detect both strongly and weakly lensed sources.

In Figure 8, we show the 850 μm-to-450 μm and 450 μm-to-850 μm flux ratios versus the observed 450 μm and 850 μm fluxes, respectively, for both our cluster and blank-field data. The main difference between Figures 1 and 8 is that we correct the data points in Figure 8 for the effect of image noise on the flux ratio measurements. This is again done by running simulations in which we generated sources with constant flux ratios, and a detailed description is left to the Appendix. In Figure 8, we also show some predicted flux ratios of a modified blackbody SED with β = 1.5 and a dust temperature of 40 K at several redshifts. At 850 μm, a clear relation between the flux ratio and the observed flux can be seen in both the cluster and blank fields. This relation can be explained by a redshift evolution if the SEDs of these galaxies do not change significantly with redshift. Since the observed 850 μm flux of an SMG remains almost invariant over z = 1–8 due to the strong negative K-correction, the variation of the flux ratio we see here might be a result of "cosmic downsizing" (Cowie et al. 1996), where SMGs at higher redshifts have higher gas fractions and therefore higher luminosities and star formation rates (e.g., Heavens et al. 2004; Bundy et al. 2006; Franceschini et al. 2006; Dye et al. 2008; Mobasher et al. 2009; Magliocchetti et al. 2011). Note that although the variation can be explained by an evolution in dust temperature, it would be interpreted as brighter sources having lower temperatures, which conflicts with many studies of dusty star-forming galaxies (e.g., Casey et al. 2012; U et al. 2012; Lee et al. 2013; Symeonidis et al. 2013).

Zoom In Zoom Out Reset image size

Another possible factor that can cause the redshift variation we see here would be lensing bias, in which brighter sources contain a higher fraction of high-redshift, lensed galaxies. For the sources in our cluster fields, although their redshifts are needed to compute the precise lensing magnifications, the changes in their magnifications are much more sensitive to the source positions than to the redshifts. As a consequence, we can use the magnification maps for z = 2.2 and z = 2.8 that we generated using LENSTOOL (see Section 3.4) to quantify how strong the lensing effect is for each source. We do not see any correlation between the observed flux and the magnification for our lensed sources. This suggests that these brighter sources in the cluster fields are generally not more strongly lensed and are simply brighter intrinsically.

We also cannot rule out the possibility of galaxy–galaxy strong lensing events (which are not included in the cluster lens models) that cause lensing bias in both the blank and cluster fields. Theoretical predictions (e.g., Blain 1996; Perrotta et al. 2002, 2003; Negrello et al. 2007; Paciga et al. 2009; Béthermin et al. 2012a; Wardlow et al. 2013) showed that the fraction of strongly lensed sources increases with the observed submillimeter flux. Wide-area, flux-density limited surveys with Herschel have successfully discovered many bright, strongly lensed SMGs (e.g., Negrello et al. 2010; Conley et al. 2011; Wardlow et al. 2013). However, at the flux range of our SCUBA-2 sources, these galaxy–galaxy strong lensing events should be rare and should have little effect on the observed redshift distribution. If we take the count models from Béthermin et al. (2012a), 850 μm sources with fluxes of 3, 5, 7, 9, 11, and 13 mJy (corresponding to the flux bins in Figure 8) have lensing fractions of ∼1%, 2%, 3%, 5%, 7%, and 10%, respectively. Although a fraction of 10% might cause a significant effect, the lensing fractions in the four faintest 850 μm flux bins in Figure 8 are too small to produce the variation of redshift we see in both cluster and blank fields. Therefore, we conclude that lensing bias only has a minor effect on the observed redshift distribution and a downsizing scenario is the most likely cause.

On the other hand, we do not see a clear relation between the submillimeter flux ratio and the 450 μm flux, mainly because of the large uncertainties due to small number statistics. Deeper 450 μm maps obtained in the future should improve our measurements. However, a nearly flat distribution for the lensed sources shown in Figure 8 is in agreement with Roseboom et al. (2013). This result might still be consistent with a downsizing scenario, given that the observed 450 μm flux of an SMG does decrease with the redshift. The observed distribution of flux ratios might be flattened due to a mixture of high-redshift bright and low-redshift faint objects in the same flux bin. A similar trend is also seen in Figure 3 of Béthermin et al. (2015), where the median redshift of 450 μm sources between flux-density cuts of 10 and 50 mJy (the flux range shown in Figure 8) changes less than that of 850 μm sources between flux-density cuts of 2 and 14 mJy (the flux range shown in Figure 8) in both full samples and strongly lensed samples. Again, if we take the count models from Béthermin et al. (2012a), 450 μm sources with fluxes of 15, 25, 35, and 45 mJy (corresponding to the flux bins in Figure 8) have galaxy–galaxy lensing fractions of ∼1%, 2%, 4%, and 7%, respectively, which have little effect on the observed redshift distribution.