BLAST: THE REDSHIFT SURVEY

Excluding the cosmic microwave background, the main peaks in the extragalactic background radiation are in the optical and far-IR/submillimeter wavebands with roughly the same amount of energy in each (Dwek et al. 1998; Fixsen et al. 1998), implying that approximately half of the total energy emitted by galaxies since their formation has been absorbed by dust and then reradiated at longer wavelengths. This energy budget strongly suggests that to completely understand galaxies and their evolution it is crucial to understand the nature of the sources that make up the cosmic infrared background (henceforth the CIB). However, in the 13 years since the discovery of this background (Puget et al. 1996) it has proved difficult to answer this question, partly because of the technical challenges of working at these wavelengths and partly because the atmosphere is opaque over much of the wavelength range from 20 μm to 1 mm, with only the 850 μm atmospheric window having routine transmission of over 50%.

After the discovery of the CIB, much of the early progress in determining the nature of the sources that compose it came from the ground-based surveys with the SCUBA camera on the James Clerk Maxwell Telescope (Smail et al. 1997; Hughes et al. 1998; Barger et al. 1998; Eales et al. 1999). These surveys resolved about 30% of the background at 850 μm and close to 100% of the background if one includes the small numbers of sources detected in lensing surveys (Blain et al. 1999; Knudsen et al. 2008). However, the full potential of these surveys has been hard to achieve due to the poor angular resolution combined with the faintness of the optical counterparts, which has made it a challenge both to identify the correct optical counterparts and to measure their redshifts. The most extensive redshift survey of the SCUBA surveys (Chapman et al. 2005) found a median redshift of ≃2.2, and in general the SCUBA sources are luminous dusty galaxies seen in the early universe that are even more luminous than the Ultraluminous Infrared Galaxies (ULIRGs) found in the universe today (Coppin et al. 2008). Evidence from X-ray observations (Alexander et al. 2005) and from mid-infrared spectroscopy with Spitzer (Pope et al. 2008; Menendez et al. 2007; Menendez et al. 2009) suggests that while a large fraction of these sources appear to contain active nuclei, most of the energy emitted by these objects ultimately comes from young stars rather than an obscured active nucleus. The star formation rates implied by the luminosities of these objects are often as much as 1000 M_☉ yr⁻¹ (Alexander et al. 2005), enough to build a large galaxy in only 1% of the age of the universe. Many authors have argued that the space density of these sources and their implied star formation rates show that they are probably the ancestors of present-day elliptical galaxies (Lilly et al. 1999; Scott et al. 2002; Dunne et al. 2003).

The SCUBA surveys, however, had two major limitations. First, the energy in the background (I_νν) at 850 μm is only one-thirtieth of the energy in the background at its peak at ≃200 μm, and so the sources detected in the SCUBA surveys may not be representative of the CIB as a whole. Dye et al. (2007) used a stacking argument to show that the sources that constitute 30% of the background at 850 μm make up at most 18% of the background at 160 μm. Chapman et al. (2005) used indirect arguments to reach the stronger conclusion that the sources that make up 60% of the background at 850 μm contribute only 6% of the background at 200 μm. Second, since most of the sources detected in the SCUBA surveys are at very high redshifts, we actually know remarkably little about the submillimeter properties of the nearby universe. To produce a fair sample of the nearby universe that is not biased by the presence of a small number of clusters or unusually empty regions, it is necessary to survey a large area of sky, which was not possible with SCUBA because of its small field of view. Therefore, estimates of the local luminosity function at submillimeter wavelengths ($\rm 100\ \mu m < \lambda < 1\ mm$), which are crucial for investigating the cosmic evolution in this waveband, are based either on extrapolations from the survey with the IRAS satellite at shorter wavelengths or on submillimeter observations of samples of galaxies selected in other wavebands (Dunne et al. 2000; Vlahakis et al. 2005). Both of these approaches have obvious drawbacks.

There has recently been a major step forward in this field as the result of observations with the Balloon-borne Large Aperture Submillimeter Telescope (BLAST; Devlin et al. 2009). BLAST has carried out surveys of the extragalactic sky in two fields, one near the South Ecliptic Pole and one centered on the southern field of the Great Observatories Origins Deep Survey (GOODS-South). Each survey covered about 10 deg², and, for comparison, the largest SCUBA survey only covered ≃0.33 deg² (Coppin et al. 2006). The BLAST surveys were at three wavelengths—250, 350 and 500 μm—and since the shortest wavelength is close to the peak of the CIB, the sources detected in these surveys are likely to be more representative of the CIB than the sources detected in the SCUBA surveys.

The BLAST survey of GOODS-South has been particularly useful because of the wealth of data at other wavelengths that exists in this field. There have been several studies of the statistical properties in the BLAST bands of sources from catalogs defined from Spitzer 24 μm observations (Devlin et al. 2009; Marsden et al. 2009; Pascale et al. 2009). These have shown that the 24 μm sources may well contribute all of the CIB (see also Dole et al. 2006). Therefore, whereas the sources found in samples at one end (850 μm) of the far-IR/submillimeter waveband are not representative of the CIB, those at the other end (24 μm) do seem to be. By combining the BLAST results with Spitzer 70 μm data and a mixture of photometric and spectroscopic redshifts, Pascale et al. (2009) made the first direct measurements of the history of dust-obscured energy output, finding a gradual increase from z = 0 to z = 1. Finally, Viero et al. (2009) investigated the clustering of star-forming galaxies from the power spectra of the BLAST maps, and Patanchon et al. (2009) used the distribution of fluctuations in the maps to estimate the submillimeter number counts.

Whereas these studies looked at the statistical properties of the BLAST maps or the statistical properties in the BLAST bands of galaxies in catalogs selected in other wavebands, a sixth paper (Dye et al. 2009) looked for the counterparts in other wavebands of the individual sources detected in the BLAST survey. This is a challenge because the angular resolution of BLAST (FWHM of 36, 42, and 60 arcsec at 250, 350, and 500 μm, respectively) is larger than the angular resolution of SCUBA at 850 μm (FWHM of 14 arcsec). Nevertheless, using the standard frequentist technique (Section 2) that has been used for other submillimeter surveys, Dye et al. 2009 (henceforth D09) succeeded in finding radio and/or 24 μm counterparts for 227 out of 351 sources detected at 5σ in the BLAST survey centered on GOODS-South. The authors used the spectroscopic and photometric redshifts that exist for many of these counterparts to show that 75% of them lie at z < 1, while only a handful of SCUBA sources lie at such a low redshift. The luminosities of these counterparts are also lower than those of the SCUBA galaxies, being more typical of luminous infrared galaxies (LIRGs) than ULIRGs. An important point is that the catalog used by D09 and in this paper is a 5σ catalog in the sense that σ is the instrumental noise, whereas an additional source of noise is the fluctuations in the map produced by faint sources. The effect of both types of noise on the fluxes of the sources in the catalog is one of the issues we will address in this paper.

This paper represents a continuation of the work described in D09. We present the results of a redshift survey of the counterparts to the BLAST sources with the AAOmega multi-object spectrometer on the Anglo-Australian Telescope. We use the spectroscopic data for a number of different purposes. First, we use the spectroscopy to investigate the nature of the galaxies that are the counterparts to the BLAST sources. Second, we use the spectroscopic data and the imaging data that exist for this field to address the issue of multiple counterparts to the BLAST sources. This is a familiar problem from attempts to find counterparts to SCUBA sources (Ivison et al. 2007) and occurs when the frequentist approach finds multiple possible counterparts to the submillimeter source. The possible causes are either that the submillimeter source actually consists of two submillimeter sources confused together—an obvious strong possibility given the poor angular resolution–or that only one of the possible counterparts is a submillimeter source with the second galaxy being physically associated in some way (possibly in the same galaxy group) with the first. Third, we use the spectroscopic redshifts to investigate the accuracy of the photometric redshifts used in D09. Fourth, we use a combination of spectroscopic and photometric redshifts to make the first estimates of the luminosity function at these wavelengths, and also a first estimate of the dust-mass function (the space density of galaxies as a function of dust mass).

The layout of this paper is as follows. In Section 2, we revise the frequentist identification technique and give a list of secondary counterparts that complements the list of primary counterparts given in D09. Section 3 describes the redshift survey. Section 4 describes the results of the analysis based on the redshift survey, including the first estimates of the galaxy luminosity function in this waveband. Section 5 contains a brief discussion and our conclusions. We assume everywhere the standard concordance cosmology: Ω_M = 0.28,  Ω_Λ = 0.72,  H₀ = 72 km s⁻¹ Mpc⁻¹.

Full details of the multi-wavelength data sets are given in D09. Briefly, the BLAST survey of GOODS-South consisted of a wide-area map of 8.7 deg² with a deeper confusion-limited map of 0.8 deg². D09 lists a catalog of all the sources detected at >5σ in any of the three BLAST bands. From the point of view of the detection of counterparts, there are two distinct regions. The central 30 × 30 arcmin² of the deep BLAST survey was surveyed by the Far-Infrared Deep Extragalactic Legacy Survey (FIDEL; Magnelli et al. 2009), and is the same region that was surveyed in the 17-band optical survey COMBO-17 (Wolf et al. 2004). We call this region the "FIDEL area." Outside this region, the whole BLAST survey area has been covered by the Spitzer Wide-area InfraRed Extragalactic Survey(SWIRE; Lonsdale et al. 2004) in all the Spitzer bands, although only ≃4 deg² were surveyed by the Spitzer team in optical bands (u, g, r, i, z). The radio catalogs discussed in D09 also consist of a deeper central region covering 0.33 deg² and a wider shallower catalog covering ≃4 deg².

Full details of the identification procedure were given in D09. Here we revise the main points. We searched for 24 μm and radio counterparts to the BLAST sources using the frequentist approach of Lilly et al. (1999), which is based on the method of Downes et al. (1986). The method is to search for possible counterparts close to the submillimeter position and then use a Monte Carlo analysis to estimate the probability that the possible counterpart is there by chance and is actually not genuinely associated with the submillimeter source. The advantage in this situation of this frequentist approach over Bayesian approaches (Sutherland & Saunders 1992) is that it does not require much information about the positional errors, which in this case are poorly known because of the effects of source confusion. The details of the procedure are given here for a radio catalog but are the same for a 24 μm catalog.

• 1.  
  Select a random position within the area common to the BLAST and radio catalogs.
• 2.  
  Find the minimum of the quantity S = r²_sepn(>f), where r_sep is the separation between a radio source and the random position, f is the flux density of the radio source, and n(>f) is the surface density of radio sources brighter than this radio source. An important point is that only radio sources within a maximum separation radius of r_max are included (see below).
• 3.  
  Repeat steps 1 and 2 for N realizations to determine the distribution of S for the radio sources.

We determined this distribution separately for the FIDEL area and the area outside FIDEL. If we then have a real potential counterpart with a value for S of S_i, we can estimate the probability that it is simply there by chance from the distribution of S generated by the Monte Carlo simulation, D(S). The probability that the potential counterpart is simply there by chance is

$$\begin{equation} P(S<S_i) = {1 \over N} \int ^{S_i}_0 D(S) dS. \end{equation} \tag{ 1 }$$

The crucial point that was investigated in D09 was the choice of r_max. Even though the accuracy of the submillimeter positions is uncertain, we do know enough about the accuracy of the positions to insist on a maximum value for r_sep; otherwise a bright radio source such as Cygnus A would yield a low value of P even if it were many degrees away from the BLAST source. The choice of the value of r_max is a balance between not missing genuine counterparts and including too many false IDs. D09 describes a method for determining the value of r_max at which the expected number of excluded genuine counterparts equals the number of included false counterparts. A byproduct of this analysis was an estimate of the distribution of offsets between the positions of the BLAST sources and the counterparts:

$$\begin{equation} n(r) \propto r e^{-r^2 \over 2 \sigma ^2} \end{equation} \tag{ 2 }$$

with σ≃ 8 arcsec. This agrees well with a prediction (D09) based on the analytical formula for the positional errors of submillimeter sources derived by Ivison et al. (2007). We derived a value for r_max of 20 arcsec in the FIDEL area and a value of 25 arcsec outside this area.

In D09, we listed the counterparts with P < 0.05 for the BLAST sources in the FIDEL area and the counterparts with P < 0.1 for the BLAST sources outside this area. The different values of P_max in the two regions were chosen because of the different surface densities of 24 μm and radio sources in the two regions. By summing the values of P for our list of 227 possible counterparts, we predict that ≃5 are incorrect. The counterparts listed in D09 were the primary counterparts, the counterpart to each BLAST source that had the lowest value of P and satisfied the condition P < P_max. In Table 1, we list the counterpart with the next lowest value of P and P < P_max, if one exists. There are 69 of these "secondary counterparts" compared to 227 primary counterparts. Approximately one third of the BLAST sources have more than one possible counterpart. Figure 1 shows the optical or mid-infrared images of all BLAST sources with more than one counterpart. We will present images of all the counterparts in a later paper (L. Moncelsi et al. 2010, in preparation). For the sources in the FIDEL region, these images are taken from the optical R-band image from the COMBO-17 survey (Wolf et al. 2004). For the sources outside the FIDEL region, the image is taken by preference from an image taken in the r-band by the SWIRE team (Lonsdale et al. 2004) and, if that does not exist, from the Spitzer 3.5 μm image. We discuss the reason for multiple counterparts in Section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On 2008 November 24, we used the AAOmega Spectrometer on the Anglo-Australian Telescope to observe targets from a preliminary list of counterparts to the BLAST sources. AAOmega (Sharp et al. 2006) consists of 392 fibers that feed the light from targets within a field 2° in diameter to a blue and a red camera via a dichroic. We used the 580 V and the 365 R gratings for the blue and the red cameras respectively, which gave a wavelength coverage from 370 to 880 nm and a resolution (λ/δλ) of 1300.

We adopted the following scheme for placing the fibers. We only placed fibers on targets with sufficiently accurate positions (the fibers are 2.1 arcsec in diameter). If a BLAST source had a radio counterpart, we placed the fiber on the radio position. If a BLAST source had only a 24 μm counterpart, we searched for an optical or 3.6 μm counterpart within 3 arcsec of the 24 μm position; we placed the fiber on the optical position if an optical counterpart existed and, if not, on the 3.6 μm position. We were constrained in our placement of fibers by the geometry of the BLAST survey area. The wide-area survey is much larger than the field of view of AAOmega, whereas the deep central area in which there is the greatest density of counterparts is significantly smaller (0.87 deg²) than the AAOmega field of view. To observe our main target list, we used three configurations of fibers. In each configuration, we observed the same targets in the BLAST deep area but a different set of targets in the surrounding area, observing the central targets with 12 exposures of 1800 s and the surrounding targets with four exposures of 1800 s. There were far more fibers than counterparts to 5σ BLAST sources with accurate positions, and so this list consisted of counterparts to BLAST sources with signal to noise >3.5. Even then we had spare fibers, so we placed the remaining fibers on the counterparts to sources detected in the SWIRE survey (Lonsdale et al. 2004). The counterparts to the BLAST sources in this list were always the primary counterparts, the ones with the lowest probability of being there by chance. We also observed a second list of targets, which contained many of the secondary counterparts (Section 2), but we only succeeded in observing this list for two exposures of 1200 s (we had planned to come back to this field on succeeding nights but the rest of the run was lost because of the weather).

We reduced the data using the standard data-reduction pipeline for AAOmega, 2dfdr, which is described on the AAOmega Web site.⁹ We then extracted redshifts from the spectra using the RUNZ software package (S. Croom 2009, private communication). This package automatically extracts a redshift from each spectrum by fitting continuum templates to the spectrum and by looking for emission lines. A crucial aspect of the program is that the user is able to inspect the result and use his/her judgment to, if necessary, change the redshift, for example if one of the lines used in the fit looks like an artifact. Two of us (S.A.E. and L.M.) did this independently and agreed in all cases. We obtained spectroscopic redshifts for 399 galaxies, 82 of which are primary or secondary counterparts to 5σ BLAST sources. Out of the 669 targets in our main list, a mixture of BLAST and SWIRE galaxies, we obtained spectroscopic redshifts for 339, a success rate of 51%.

Table 2 contains the list of spectroscopic redshifts we have obtained for the counterparts to the 5σ BLAST sources. The counterparts are mostly the primary counterparts, but nine are the secondary counterparts. The other redshifts, which are for a mixture of SWIRE galaxies and counterparts to BLAST sources that now fall below the 5σ cutoff, are given in Table 3. Figure 2 includes spectra of a representative sample of the BLAST counterparts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is impossible to quantify accurately the probability of a redshift being correct, because the final redshift is a combination of the automatic continuum and emission-line fitting by RUNZ plus the subjective judgment of the user. We adopted the quality assessment system used by RUNZ, in which a redshift with a quality of 5 is defined as being a "definite redshift," one with a quality of 4 is defined as being "almost certain with roughly a 95% probability of being correct" and one with a quality of 3 as being "somewhat less certain but probably correct." We have listed in Tables 2 and 3 our estimates of the quality of each redshift using this system. In many cases, the value is simply the one produced by RUNZ; in others it is our modification of the RUNZ value based on an inspection of the spectrum. Of the 82 redshifts we measured for the BLAST counterparts, two have a quality flag of 3, 10 have a quality flag of 4, and the remainder have a quality flag of 5.

Taylor et al. (2009) have measured spectroscopic redshifts for 21 of the same galaxies for which we have measured redshifts. In 19 cases, our redshifts and those of Taylor et al. agree. In the case of the first discrepancy, our redshift is 0.672, whereas that of Taylor et al. is 0.553. We inspected our spectrum, and the redshift has a quality flag of 5 and is based on detection of an [O ii] 372.7 emission line and several absorption features. We are therefore confident that our redshift is correct. In the case of the second discrepancy, our redshift of 0.205 is of a poor quality (quality flag of 3), whereas the redshift of Taylor et al. is 0.620. Our redshift is based on a fit to the continuum and a single emission line, which we have assumed is [O iii] 500.7. Taylor et al. have clearly assumed that the line is [O ii] 372.7. On reinspection of our spectrum, we decided this is more likely to be correct, and so the redshift we have listed is based on this assumption.

The spectra in Figure 2 are a representative sample of the spectra we obtained for the BLAST counterparts, consisting of spectra of the 1st, 6th, 11th etc. galaxies listed in Table 2. They show quite clearly that we were undoubtedly more successful in obtaining spectroscopic redshifts for galaxies with low redshifts than for galaxies at high redshifts. It was fairly easy to measure the redshift of a galaxy at z < 0.3 because of the presence of the bright Hα and [N ii] 658.3 lines, the [SII] 671.6, 673.1 doublet and often many other absorption and emission lines. At 0.3 < z < 1.0, it was a little more difficult: the Hα and [N ii] 658.3 lines are redshifted out of the accessible waveband, and the only bright emission line is the [O ii] 372.7 line. In this range of redshifts, the redshift usually came from this line plus a fit to the continuum. At even higher redshifts, the [O ii] line is redshifted out of the waveband, and it was only possible to extract a redshift if the object is a quasar because quasars have several broad emission lines, such as $\rm C\,{{\mathsc {iii}}}] 190.9$, which appear in the accessible waveband at z>1. Thus, the success rate for obtaining redshifts for BLAST counterparts that do not have active nuclei probably falls from close to 100% at z = 0 to close to 0% at z ≃ 1. Therefore, for investigating the evolution of the luminosity function (Section 4.4), we still rely on photometric redshifts.

4.1. Inferences from the Spectra

The absolute flux scale of our spectra and its dependence on wavelength are uncertain, since we did not observe spectrophotometric standards. Nevertheless, we were able to measure two quantities from the spectra that do not require this: the equivalent width of the Hα line and the ratio of the flux in the [N ii] 658.3 line to the flux in the Hα line, two lines which are very close in wavelength. In Table 2, we have listed the Hα equivalent width of each galaxy, corrected to the galaxy's rest frame, and the value of this line ratio. If there are no values for these quantities, it is either because the galaxy is at too high a redshift for these to be measured or because there was a problem with the spectrum or because the Hα line is broad.

The Hα equivalent width is useful because it gives a measure of the star formation rate in the galaxy relative to the average star formation rate since the galaxy was formed. Figure 3 shows a histogram of the Hα equivalent width for the BLAST primary counterparts and for a sample of ≃25,000 galaxies with 0.05 < z < 0.095 drawn from the Sloan Digital Sky Survey (SDSS) and the 2dF Galaxy Redshift Survey (2dFGRS; Balogh et al. 2004). Note that because of the very large number of galaxies with equivalent widths less than 5 Å in the latter sample, the plot for this sample is on a logarithmic scale. These galaxies with low equivalent widths represent the "old, red, and dead" population in which the star formation rate was much higher in the past than it is today. It is now clear that there is a dichotomy between the properties of this population and those of actively star-forming galaxies (Kauffmann et al. 2003). The most striking difference between the two samples is that there are very few BLAST galaxies with equivalent widths less than 5 Å, and so the BLAST counterparts are almost exclusively drawn from the actively star-forming population. However, if we exclude the galaxies in both samples with equivalent widths <5 Å, the two samples appear quite similar, with the mean equivalent width being 20 Å for the SDSS/2dFGRS sample and 24 Å for the BLAST sample. Therefore, judged by the Hα equivalent width, the BLAST counterparts are star-forming galaxies but do not appear to be extreme ones. We will address the issue of whether the BLAST galaxies are exceptional or run-of-the-mill galaxies in more detail in a subsequent paper (L. Moncelsi et al. 2010, in preparation).

Zoom In Zoom Out Reset image size

Five of the BLAST counterparts clearly contain powerful active galactic nuclei (AGNs) because they have spectra typical of quasars. We looked for less powerful active nuclei by measuring line ratios. The classic way of determining whether the emission lines from a galaxy are dominated by emission from gas that is photoionized by an active nucleus or by gas that is heated by young stars is to look at the galaxy's position on a line ratio diagram. For example, the two classes fall in separate regions in a diagram of [N ii] 658.3/Hα versus [O iii] 500.7/Hβ (Baldwin et al. 1981). We almost never had measurements of all four of these lines, but we did have measurements of the first ratio for many of the galaxies. In their study of the galaxies detected in the SDSS, Miller et al. (2003) argued that if this line ratio is >0.63 the galaxy must lie in a region of the four-line diagram dominated by AGNs (the reverse is not true because if the line ratio is below this value, the galaxy may still lie in a region of the diagram that is dominated by AGNs). We were able to measure this ratio for 54 of our sample of 73 primary counterparts for which we have redshifts. If we include the objects with a value for this line ratio >0.63, the quasars and two other objects for which there is some evidence for an AGN (broad Hα in one case and a large value for the ratio of the [O iii] 500.7 and Hβ lines in the other), there is evidence for an AGN in 19 out of 73 primary counterparts. This percentage (26%) is rather greater than the percentage of SDSS galaxies that contain AGNs, which is approximately 18% (Miller et al. 2003), but given the differences between how these samples were selected we do not believe we can draw any profound conclusions from this result.

4.2. A Test of Photometric Redshifts

We have not measured enough spectroscopic redshifts to be able to investigate cosmic evolution without recourse to photometric redshifts. A potential problem with using photometric redshifts is that the BLAST galaxies probably contain large amounts of dust, and most photometric redshift techniques have only been tested on galaxies detected in optical surveys. The only method that has been tested on galaxies that may be similar to those detected by BLAST is that of Rowan-Robinson et al. (2008), who estimated redshifts for galaxies detected in the Spitzer SWIRE Legacy Survey.

We used the spectroscopic redshifts to test the accuracy of the two sets of photometric redshifts that we used in our investigation of the evolution of the luminosity function (Section 4.4). The first set was obtained by Wolf et al. (2004) from the COMBO-17 survey, a survey of a field $\rm 30 \times 30$ arcmin² in size through 17 optical filters. The second set were obtained by Rowan-Robinson et al. (2008) for the much larger area (≃4 deg²) covered by the optical images taken as part of the SWIRE survey. The quality of the imaging data used to estimate the second set of redshifts (three broadband optical images plus the Spitzer images at 3.6 and 4.5 μm) is more typical of the imaging data available over large areas of sky. Since there are many new large-area surveys, such as those with Herschel and VISTA (Sutherland 2009), which will be reliant on photometric redshifts, a test of the latter set is particularly interesting. Our test is unbiased because our new spectroscopic redshifts are not part of the sets that were used to tune either of the original methods.

Figure 4 shows a comparison of the photometric and spectroscopic redshifts for the two data sets. We have included all the galaxies for which we have spectroscopic redshifts (i.e., Tables 2 and 3). We have only included spectroscopic redshifts with a quality of four or greater (Section 3) to ensure that any discrepancies are caused by errors in the photometric redshifts rather than errors in the spectroscopic ones. We only included SWIRE photometric redshifts that were based on photometry in at least three optical bands (Rowan-Robinson et al. 2008). In our error analysis, we used $\delta = { z_{\rm phot} - z_{\rm spec} \over 1 + z_{\rm spec} }$ as our measure of the discrepancy between the photometric and spectroscopic redshift, and we treated any photometric redshift with $\rm |\delta | > 0.15$ as a catastrophic error.

Zoom In Zoom Out Reset image size

We found that the percentage of catastrophic errors is ≃8% for the COMBO-17 photometric redshifts and ≃15% for the SWIRE redshifts. However, many of the discrepant SWIRE photometric redshifts are for objects flagged by Rowan-Robinson et al. (2008) as probable quasars (confirmed in 9 out of 10 cases by our spectroscopy and probably also true in the tenth case), and once these objects are removed the percentage of catastrophic errors falls to 9%. The very discrepant object in Figure 4(a) was also shown by our spectroscopy to be a quasar but was not flagged as such in the COMBO-17 survey. Excluding all objects for which there are catastrophic errors, we estimate that the errors for the two methods ($\sqrt{\langle \delta ^2 \rangle }$) are 0.031 for the COMBO-17 data set and 0.056 for the SWIRE data set. In calculating the luminosity function (Section 4.4), we have used redshift bins with a width of 0.2 in redshift, and we therefore conclude that the errors in the photometric redshifts are not likely to be the limiting factor in the accuracy of our analysis.

4.3. Multiple Counterparts

A common problem with submillimeter surveys is that there are sometimes two possible counterparts to the submillimeter source, each of which has a low probability of being there by chance. Ivison et al. (2007) found that approximately 10% of sources detected in the SHADES 850 μm survey have more than one counterpart. For the BLAST survey of GOODS-South, the percentage of sources with multiple counterparts is ≃30% (Section 2). Of the 69 sources with multiple counterparts, 13 have more than two counterparts, although, as we show below, in three cases the 24 μm/radio sources are all associated with a single large galaxy.

This high percentage might, of course, be due to the clustering of 24 μm or radio sources in the catalogs used for the identification analysis; one of the counterparts might be genuine with the second simply being there because of the clustering within the radio/24 μm catalog. We have investigated this by adapting the Monte Carlo simulation described in D09. We lay down points randomly in the area covered by the 24 μm and radio catalogs, carrying out separate simulations for the radio and 24 μm catalogs and for the area covered by the FIDEL survey and for the area outside this survey (Section 2). The points represent artificial submillimeter sources. We then apply the frequentist identification technique to look for counterparts at this list of positions, using the same search radii that we used for the real data. For all the counterparts that satisfy the condition that the probability is <P_max, we determine the percentage of cases in which there is also a second counterpart that satisfies this condition. We find that the percentages in the FIDEL region are 10.7% for the radio catalog and 3.3% for the 24 μm catalog, with the percentages outside the FIDEL area being 7.2% for the radio catalog and 5.8% for the 24 μm catalog. These percentages should also be good estimates of the percentages of the real counterparts that have secondary counterparts because of clustering in the radio/24 μm catalogs. Applying these percentages to the real list of counterparts, we estimate that ≃15 of the 69 sources with multiple counterparts are caused by this effect. Therefore, the true number of sources with multiple counterparts is ≃54, ≃24% of the total number of sources with counterparts. This is still larger than the percentage found by Ivison et al. (2007) for the SHADES survey, which is not surprising because of the poorer angular resolution of BLAST (36 arcsec at 250 μm verses 14 arcsec for SHADES), although the difference in effective linear resolution is somewhat less because the BLAST sources tend to be at lower redshift (D09); for example, the linear resolution for a BLAST source at z = 0.3 is only 33% greater than the linear resolution for a SHADES source at z = 2.

There are a number of possible explanations of multiple counterparts. One of these is the possibility that there is a single genuine counterpart, which has been gravitationally lensed by a nearby galaxy (Blain 1998). The arguments against this are (1) the angular distances between the counterparts are often much greater than the typical distances between a lensed source and a lense and (2) that the BLAST galaxies are typically spirals or interacting galaxies, whereas most lenses are predicted to be ellipticals (Blain 1998). There are two explanations that cannot be ruled out.

• 1.  
  The two counterparts are physically associated in some way. They might either be in the same cluster or they might be two galaxies that are gravitationally interacting. If this explanation is correct, the counterparts need not both be submillimeter sources.
• 2.  
  Both of the counterparts are submillimeter sources which are not linked physically in any way.

We can investigate which of these is correct by first inspecting the images of all the BLAST sources with multiple counterparts, and second by comparing the redshifts of the primary and secondary counterparts. The images (Figure 1) immediately reveal a few interesting facts. First, there are several sources where there is clearly only a single galaxy. In these cases, it seems almost certain that the two apparent counterparts are actually radio or 24 μm sources within the galaxy. For example, BLAST 4 (see Table 1 for full name) has two radio counterparts which lie close to the center of a large galaxy. Inspection of the radio image shows that these two sources are actually two peaks in a single source that is extended in the same direction as the optical structure. Not surprisingly, our spectroscopic redshifts of these two counterparts are virtually the same. BLAST 2 and BLAST 53 also seem to be cases where the two apparent counterparts are actually sources within a single galaxy. Second, there are several cases where the counterparts seem to be clearly interacting. In four cases—BLAST 6, 9, 103, and 127—there is clear morphological evidence for a gravitational interaction between the two counterparts.

Unfortunately, this still leaves a large number of multiple counterparts for which there are no clear morphological clues. We can make more progress by comparing the redshifts, either spectroscopic or photometric, of the primary and secondary counterparts. We have redshifts for both counterparts for 27 systems, although unfortunately spectroscopic redshifts for both counterparts in only three cases. Figure 5 shows there is a clear correlation between the two redshifts (Spearman's ρ = 0.52; probability of the two variables being uncorrelated is ≃0.4%) with 14 out of the 27 systems lying sufficiently close (given the errors on the photometric redshifts) to the line on which the redshifts of the two counterparts are the same. Therefore, since there is no reason why there should be a correlation if the second explanation is correct, we conclude that in at least half the cases where there are multiple counterparts, the counterparts appear to be physically associated.

Zoom In Zoom Out Reset image size

4.4. Luminosity Functions

In this section, we make a first attempt to estimate luminosity functions at the three BLAST wavelengths. There are two major obstacles to overcome. First, there are still many BLAST 5σ sources that do not have counterparts, and some of the counterparts do not have redshifts. Second, the fluxes of the BLAST sources are systematically biased upward by the effect of noise, both instrumental noise and the fluctuations in the map caused by other submillimeter sources. This is the well-known Eddington bias (Eddington 1940), in which the effect of the steep differential source counts is that more sources in a flux-limited sample have had their fluxes increased by noise than decreased by noise. The large effects of Eddington bias in the BLAST maps are why our earlier papers concentrated on the statistical properties of the BLAST maps rather than the properties of individual sources (e.g., Patanchon et al. 2009). Here, if we wish to use the information about cosmic evolution provided by our redshift survey, we are forced to confront its effects. The confusion of sources discussed in Section 4.3 is often treated as a separate problem, but this is really just a form of Eddington bias in which the noise comes from discrete sources. This upward bias in the fluxes of sources is often called "flux-boosting." We have developed a number of techniques for overcoming these obstacles and believe that our conclusions at the end of this section are not invalidated by any of these effects.

For each wavelength, the sample from which we start is the sample of sources detected at 5σ at this wavelength and which fall in a 4.2 deg² area covered by either the SWIRE optical images or the COMBO-17 survey. The point of restricting the investigation in this way is that outside this area there are no photometric redshifts. The three samples are thus subsamples of the list in Table 3 of D09. Table 4 in this paper lists the statistics of these samples: the number of sources; the number with either a radio or a 24 μm counterpart or both; the number with either a photometric redshift or a spectroscopic redshift, with the number of spectroscopic redshifts in the brackets. The table illustrates one of the problems mentioned above, that not all the sources have counterparts and some of the counterparts do not have redshifts, either spectroscopic or photometric.

Table 4. The 5σ Samples at Each Wavelength

Wavelength	Sources	Counterparts	Redshifts
250 μm	115	94	82 (49)
350 μm	89	62	48 (27)
500 μm	107	52	39 (12)

Notes. Reading from the left, the columns are: the wavelength of the sample; the number of sources detected at >5σ at this wavelength and that are in an area covered by deep optical images; the number of these sources with radio and 24 μm counterparts or both; the number of counterparts with either photometric or spectroscopic redshifts (the number of spectroscopic redshifts is in brackets).

Download table as:  ASCIITypeset image

We will address the problem of the lack of redshifts first because it is easiest to deal with. It seems likely that counterparts without redshifts are galaxies that are at very high redshift, and are thus too faint at optical wavelengths for photometric redshift methods to work. We can test this by comparing the mid-infrared colors of the counterparts with and without photometric redshifts. Figure 6 shows all the counterparts on a plot of S_3.6/S_4.5 versus S_5.8/S_8.0, colors which Pascale et al. (2009) show depend on redshift. The counterparts without photometric redshifts cluster in the bottom right of the figure, which Figure 3 of Pascale et al. shows corresponds to z>1. Therefore, we can assume that as long as we restrict our estimates of the luminosity function to z < 1, our estimates should not be affected if we omit these objects.

Zoom In Zoom Out Reset image size

The problem of the missing counterparts is more complicated because there are several possible causes. The first is that the BLAST sources without counterparts are at such high redshifts that their 24 μm and radio fluxes fall below the limits of the 24 μm and radio catalogs. Evidence for this is the fact that the percentage of BLAST sources with counterparts falls with increasing wavelength, in line with the predictions of models that the fraction of sources at very high redshift should increase with increasing wavelength (D09). If this is the cause, then we should again have no problems if we restrict our estimates of the luminosity function to low redshifts. The second possibility, suggested by the simulations in Appendix B, is that some of these sources are not genuine sources or are a confused combination of instrumental noise and many very dim sources. If these sources are not real sources, of course, omitting them from our estimates of the luminosity function is the right thing to do. However, a third possible explanation which must be correct for at least some of the BLAST sources is that the 24 μm or radio counterpart is in the catalog but our frequentist method of finding the counterparts (Section 2) has failed to find it. It is possible to estimate the number of counterparts that are present in the catalogs but are missed by the selection procedure from the predicted distribution of distances between the BLAST sources and their counterparts that was given in Equation (2). Each counterpart found by the method should be multiplied by a correction factor to compensate for the counterparts that could not have been found by our frequentist method. This is given by

$$\begin{equation} c_i = { \int ^{\infty }_0 r e^{-r^2 \over 2 \sigma ^2} dr \over \int ^{r_{\rm cut,i}}_0 r e^{-r^2 \over 2 \sigma ^2} dr}, \end{equation} \tag{ 3 }$$

in which σ≃ 8 arcsec (Section 2), and r_cut,i is the maximum radius at which this counterpart could have been found by this method. For counterparts that are bright 24 μm or radio sources, r_cut,i is the same as the search radius, r_max. However, a faint 24 μm source would have been dismissed as a possible counterpart at a smaller distance from the BLAST source than r_maxbecause its probability of being there by chance would have exceeded P_max. Therefore, r_cut, i is given by the smaller of the search radius, r_max, and

$$\begin{equation} r^{\prime } = \left({P_{\rm max} \over P_i}\right)^{1/2} r_i, \end{equation} \tag{ 4 }$$

in which r_i is the distance of the counterpart from the BLAST position, and P_i is the probability that the potential counterpart is not genuinely associated with the BLAST source. The number of counterparts missed by the selection procedure is given by

$$\begin{equation} N_{\rm miss} = \sum _{i=1}^n c_i - n, \end{equation} \tag{ 5 }$$

in which n is the number of counterparts.

Figure 7 shows a plot of 1/c_i versus redshift for all counterparts of 250 μm sources with either photometric or spectroscopic redshifts. There is a small correction factor for all counterparts, simply because we have used a maximum search radius (r_max) in the procedure, and so some true counterparts will have been missed no matter how bright they are at radio wavelengths or at 24 μm. However, for some counterparts the correction factor is much larger because they are faint 24 μm/radio sources and therefore only have sufficiently low values of P_i if they are very close to the BLAST position. The fraction of counterparts with high values of c_i increases with redshift, showing this incompleteness effect is worse at higher redshifts. The number of missing counterparts, according to Equation (5), is 388, which is clearly too high because the number of missing counterparts in Table 4 is only 21. However, the number of missing counterparts is dominated by the three counterparts with extremely high values of c_i. If we omit these as statistical fluctuations, we obtain a value of 24. Therefore, although the number of counterparts missed by the method is clearly uncertain, partly because of the effect of huge statistical fluctuations caused by small offsets and partly because the value of σ we have used in Equation (3) has a large error, it is clear that some of the missing counterparts to the 250 μm sources must be in the radio/24 μm catalogs. We will return to how we correct for these missing counterparts after discussing the best way of estimating the luminosity function.

Zoom In Zoom Out Reset image size

The standard method of estimating the luminosity function is fairly simple. Suppose one wishes to estimate the value of the luminosity function (the space density of galaxies as a function of luminosity) in a particular range of redshift and luminosity. If there are n galaxies in this luminosity–redshift bin, the standard estimate of the value of the luminosity function in this bin is

$$\begin{equation} \phi (L_1<L<L_2,z_1<z<z_2) \Delta {\rm log}_{10} L = \sum _{i=1}^n {1 \over V_i}, \end{equation} \tag{ 6 }$$

in which V_i is the comoving volume in which the ith galaxy could both have been detected by the survey and still have been found within the range of redshifts for this bin. The error on this estimate is usually given as ${1 \over \sqrt{n}}$, although this is only an approximation because the values of V_i for the sources are usually different.

However, we have preferred to use a different estimate of the luminosity function which has several major advantages for deriving the luminosity function from submillimeter surveys. In this method, which was first suggested by Page & Carrera (2000), the luminosity function is given by

$$\begin{equation} \phi (L_1<L<L_2,z_1<z<z_2) \Delta \mbox{log}_{10} L = {n \over V}, \end{equation} \tag{ 7 }$$

in which n is again the number of galaxies in the bin, and V is the accessible comoving volume averaged over the luminosity range in this bin. The important difference is that V is now not calculated from the measured luminosities of the sources. Page and Carrera have shown that this estimator is always better than the one given in Equation (6), and the error on the estimate is now truly ${1 \over \sqrt{n}}$. A major advantage for submillimeter surveys is that whereas V_i in Equation (6) depends on the luminosities of the sources, which are often uncertain because of flux-boosting, V in Equation (7) does not depend on the measured luminosities of the sources. This is not a complete solution for the problem of flux-boosting, which we will discuss more below, because the number of sources, n, in a bin will only be correct if their luminosities have been measured correctly; but both the methods suffer from this problem, while only the standard method suffers from the problem that an error in the luminosity of a source also produces an uncertainty in V_i. V in Equation (7) is given by

$$\begin{eqnarray} V &=& {1 \over \Delta {\rm log}_{10}(L)} \int ^{L_u}_{L_d} \int _{\rm survey} \int ^{{\rm min}\lbrace z_u,z(L,S_{\rm min}(A))\rbrace }_{z_d} {c \over H_0} \nonumber\\ &&\times {D^2 \over \sqrt{\Omega _M (1+z)^3 + \Omega _{\Lambda }}} dz dA d{\rm log}_{10}(L) \end{eqnarray} \tag{ 8 }$$

in which Δlog₁₀L is the width of the bin in luminosity; dA is an element of the BLAST survey area; S_min is the minimum flux density a galaxy could have and still be detected in this area element¹⁰; L_d, L_u, z_d, and z_u are the limits of the bin in luminosity and redshift; and all the cosmological terms have their usual meanings.

We chose to estimate the luminosity function in five redshift slices: $\rm 0 < z < 0.2$, 0.2 < z < 0.4, 0.4 < z < 0.6, 0.6 < z < 0.8, 0.8 < z < 1.0. To partially overcome the problem of the Eddington bias, rather than using the flux densities in the original catalog (D09), we remeasured the fluxes from the BLAST images at the positions of the counterparts. To tackle the problem of the missing counterparts, we made a minimal and maximal estimate of the luminosity function. We made the minimal estimate by making no correction at all for the missing counterparts and the maximal estimate by making the assumption that all the missing counterparts in Table 4 are actually in the 24 μm/radio catalogs but have just been missed by our identification technique. In calculating the average comoving accessible volume we assumed the average spectral energy distribution found in D09. In calculating the luminosities of the individual BLAST sources, we either used the SED of the individual source given in D09 or, if that was not possible, the average SED.

We made our maximal estimate of the luminosity function by correcting for the missing counterparts in the following way. We replace Equation (7) by

$$\begin{equation} \phi \Delta {\rm log}_{10} L = {\sum _{i=1}^n c_i \over V}, \end{equation} \tag{ 9 }$$

in which c_i is the correction factor for the ith counterpart that falls in that particular luminosity–redshift bin, and the sum is over all the counterparts that fall in that bin. We adopted the following scheme for estimating the correction factors. We started by deriving the values using the same method that was used to construct Figure 7, except for the three counterparts with very high value of c_i. For these we used the average values of c_i for the rest of the counterparts at z>0.2. We then scaled all the values of c_i by a constant factor so that the number of missing counterparts predicted by Equation (5) matched the number of missing counterparts in Table 4. In doing this, we are implicitly assuming that all the missing counterparts in Table 4 are actually in the 24 μm and radio catalogs but have just been missed by our identification analysis. While this is quite possibly true at 250 μm, it seems unlikely it is true at 500 μm because of the much greater percentage of missing counterparts. This assumption is why this method yields a maximal estimate of the luminosity function.

The final question is how to deal with flux-boosting. Although we have remeasured the fluxes at the positions of the counterparts, our simulations (Appendix B) show that the fluxes are still too high. We dealt with this issue by making estimates of the luminosity function both with and without making a correction for flux-boosting. Appendix B describes how we have estimated the effect of flux-boosting in the BLAST images.

Figures 8 and 9 show our estimates of the luminosity functions at the three wavelengths in the five redshift slices. Figure 8 shows our estimates of the luminosity function when no correction is made for flux-boosting. Figure 9 shows the effect of including the correction for flux-boosting. On the left-hand side of each figure is the minimal estimate, without any correction for the missing counterparts; on the right-hand side is the maximal estimate. These are the first measurements of the galaxy luminosity function at these wavelengths. We have also plotted estimates of the low-redshift luminosity function at these wavelengths, which have been obtained from the IRAS PSCZ survey and the results of the only large submillimeter survey of nearby galaxies (Appendix A).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Inspection of the figures shows that corrections for flux-boosting and the missing counterparts make very little obvious difference to the estimates of the luminosity function. The reason for this is that the plots are on logarithmic axes and the luminosity bins are very wide (0.33 in dex). This agreement gives us confidence that the two obvious features of the luminosity functions are correct. The first is that the agreement between the measurements of the luminosity function in the low-redshift slice (0 < z < 0.2) and the extrapolation of the local luminosity function from shorter wavelengths is surprisingly good. The second is that there appears to be cosmic evolution, in the sense that the space density of the most luminous sources increases steadily with redshift. The fact that there appears to be evidence for cosmic evolution in all the sub-panels of Figures 8 and 9 suggests that this is a robust result.

As one additional check, we have calculated the 250 μm luminosity function using only the BLAST sources outside the area covered by the FIDEL survey (Section 2). This tests whether the evolution could be the result of some peculiarity associated with our deepest optical/IR data set, which also covers part of the BLAST survey where the confusion is worst (Section 2). Figure 10, which contains no correction for flux-boosting or missing counterparts, shows that the evolution is still present even if we do not use the FIDEL data set.

Zoom In Zoom Out Reset image size

This evolution in the space density of the most luminous sources has been seen before in the Spitzer bands (Huynh et al. 2007; Le Floc'h et al. 2005), but we might suspect that our results are adding something new because whereas the monochromatic luminosity in the Spitzer bands is extremely sensitive to the temperature of the dust, in the BLAST bands the monochromatic luminosity is equally sensitive to the mass of dust that is present. Thus Figure 8 suggests that there may be strong cosmic evolution not only in the luminosities of galaxies but in the masses of dust in the galaxies. We can test whether this is so by calculating the space density of galaxies as a function of dust mass. We can do this using a straightforward adaptation of the formalism above. The monochromatic luminosity at a frequency ν is connected to the mass of dust in a galaxy by the relation

$$\begin{equation} L_{\nu } = B_{\nu }(T_d) \kappa _{\nu } M_d, \end{equation} \tag{ 10 }$$

in which B_ν is the Planck function, T_d and M_d are the dust temperature and dust mass, and κ_ν is the dust-mass opacity coefficient. Although galaxies clearly contain dust with a range of dust temperatures, Dunne & Eales (2001) have shown that most of the dust, even for a ULIRG like Arp 220, has a temperature of only ≃ 20 K. In using Equation (10) to make the connection between dust mass and luminosity, we have assumed a dust temperature of 20 K and the value of the dust-mass opacity coefficient from James et al. (2002), extrapolating this to the BLAST frequencies assuming that it scales as ν². We have used the sample of galaxies at 250 μm (Table 4) because the percentage of sources with counterparts is highest at this wavelength. The dust-mass function (the space density of galaxies as a function of dust mass) in a bin in the mass–redshift plane is then given by

$$\begin{equation} \phi (M_1 < M_d < M_2,z_1<z<z_2) \Delta {\rm log}_{10} M_d = {n \over V}, \end{equation} \tag{ 11 }$$

in which n is the number of galaxies with dust masses and redshifts that fall within this bin and V is calculated using Equation (8). Figure 11 shows the results for the five redshift slices without making any correction for missing counterparts. There is clearly strong evolution, in the sense that the space density of the galaxies with the highest dust masses increases steadily with redshift. Pascale et al. (2009) concluded that there was no evolution in the comoving density of dust in the universe. However, Pascale et al. effectively measured ϕ〈M_d〉 in each redshift slice, and Figure 11 shows that this does not change very much: the average dust mass of the galaxies detected at low redshift is lower than at high redshift but their space density is higher. It is only by comparing the space density at different redshifts but at the same dust mass that it is possible to see the evolution.

Zoom In Zoom Out Reset image size

Note that the result in Figure 11 is very insensitive to our assumptions about temperature, because on the Rayleigh–Jeans tail the Planck function in Equation (10) only depends on the first power of dust temperature. The strength of the evolution would be less if the temperature of the bulk of the dust at high redshift were higher than at low redshift. But even if the temperature were a factor of two greater at high redshift, the effect on the high-redshift points in Figure 11 would be to move them a factor of two to the left, which is not enough to remove the result. It is possible to think of scenarios in which the evolution was caused by temperature. Suppose that as one moves to higher redshift, the fraction of BLAST galaxies that contain a luminous but obscured quasar gradually increases, and by a redshift of ∼1 the temperature of the dust is a factor of 10 greater than at z = 0. This would explain the evolution seen in Figure 11. However, because of the strong dependence of bolometric luminosity on temperature, this increase in temperature would correspond to an increase in bolometric luminosity of at least a factor of 10⁵. Therefore, it is much harder to explain the evolution visible in Figure 11 as a temperature effect than as an increase in the number of galaxies with high dust masses.

We have carried out a redshift survey of the sources found in the BLAST survey of GOODS-South. Our basic results are as follows.

• 1.  
  The equivalent widths of the Hα line show that the counterparts to the BLAST sources are mostly star-forming galaxies with a mean equivalent width similar to that for the star-forming galaxies found in the SDSS and the 2dF Galaxy Redshift Survey. Therefore, the BLAST counterparts at low redshift appear to be star-forming galaxies but not particularly extreme ones.
• 2.  
  Approximately one quarter of the BLAST counterparts contain an active nucleus, judged either by the line ratios or the presence of broad emission lines.
• 3.  
  We have made an unbiased estimate of the errors in the redshifts produced by the photometric redshift methods developed from the COMBO-17 survey (Wolf et al. 2004) and the SWIRE survey (Rowan-Robinson et al. 2008). Using $\delta = { z_{\rm phot} - z_{\rm spec} \over 1 + z_{\rm spec} }$as our measure of the discrepancy between the photometric and spectroscopic redshift, we found that 8% of COMBO-17 photometric redshifts and 9% of SWIRE redshifts had catastrophic errors in the sense that |δ|>0.15. Excluding these catastrophic errors, we found that the errors ($\sqrt{<\delta ^2>}$) were 0.031 for the COMBO-17 redshifts and 0.056 for the SWIRE redshifts.
• 4.  
  We have used the redshifts and the images to investigate the 30% of BLAST sources that have two or more counterparts. We conclude that there is evidence in at least half the cases that the two counterparts are physically associated, either because they are interacting or because they are in the same large-scale structure.
• 5.  
  We have made the first estimates of the luminosity function at the three BLAST wavelengths and in five redshift slices. We find strong evolution, in the sense that the space density of the most luminous sources increases steadily with redshift out to z = 1.
• 6.  
  We have also investigated the evolution of the dust-mass function with redshift, finding gradual evolution in the space density of the galaxies with the highest dust masses out to z = 1.

The most interesting result is probably the last one. It is well known that the luminosity density of the universe evolves strongly with redshift, whether observed in the optical waveband or the far-IR/submillimeter wavebands, and that the space density of the most luminous sources evolves strongly with redshift (Lilly et al. 1996; Huynh et al. 2007; Le Floc'h et al. 2005; Magnelli et al. 2009). But this increased luminosity–density need not necessarily be associated with an increase in the amount of interstellar matter in galaxies. If the increased luminosity density is caused by an increase in the global star formation rate, it is possible, for example, that this is caused by a larger number of galaxy interactions at high redshift, which trigger starbursts, and not necessarily by the larger amount of interstellar material in galaxies. However, Figure 11 shows that the space density of galaxies with high dust masses, and thus presumably large reservoirs of interstellar material, is also evolving strongly with redshift.

We are grateful to Heath Jones for his help with the observations and Rob Sharp for his help with the 2dfdr data-reduction pipeline. This work makes use of the Runz redshifting code developed by Will Sutherland, Will Saunders, Russell Cannon and Scott Croom, and we are grateful to Scott Croom for making this available to us. We thank Seb Oliver for making available the SWIRE optical images and for a useful conversation about the SWIRE photometric redshift method and Luca Cortese for emergency help with IRAF. We acknowledge the support of NASA through grant numbers NAG5-12785, NAG5-13301 and NNGO-6GI11G, the NSF Office of Polar Programs, the Canadian Space Agency, the National Sciences and Engineering Research Council of Canada and the UK Science and Technology Facilities Council. This paper relies on observations made with the AAOmega spectrograph on the Anglo-Australian Telescope, and we thank the staff of the telescope and especially those involved in the development of the spectrograph.

The local luminosity function is well known in the wavelength range 10 ⩽ λ < 100 μm because of the all-sky IRAS survey. However, the problem in simply estimating the local luminosity function at the BLAST wavelengths by extrapolating an IRAS luminosity function is that it is not obvious what spectral energy distribution to use; the galaxies we know most about are those that were detected by IRAS but these are likely to contain warmer dust than the galaxies detected at longer wavelengths in the BLAST survey. Serjeant & Harrison (2005) suggested a way of overcoming this problem by using the results of the largest submillimeter survey of nearby galaxies: the SCUBA Local Universe Galaxy Survey (SLUGS; Dunne et al. 2000; Dunne & Eales 2001; Vlahakis et al. 2005). This is the method we have adopted here.

We started with the ∼10,000 galaxies detected in the IRAS PSCZ survey (Saunders et al. 2000), which was a survey of 82% of the sky down to a flux density of S_60 μm⪰0.6 Jy. We only included in our analysis galaxies that had both 60 and 100 μm detections and, to avoid the effects of peculiar motions and evolution, velocities between 300 and 30,000 km s⁻¹. We calculated the accessible volume for each galaxy using both the flux and the velocity limits. We then used the tight relation that exists between the ratio of 60 and 100 μm flux and the ratio of 60 and 850 μm flux which was discovered in SLUGS (Dunne & Eales 2001) to estimate the flux of each galaxy at 850 μm. The precise form of the relationship we used is the one given by Vlahakis et al. (2005):

$$\begin{equation} {\rm log}_{10} {S_{60} \over S_{100}} = 0.365 {\rm log}_{10} {S_{60} \over S_{850}} - 0.881. \end{equation} \tag{ A1 }$$

We fitted a simple two-component dust model to the 60, 100, and 850 μm values, which allowed us to estimate the luminosity of each galaxy at the BLAST wavelengths. We then used Equation (6) to estimate the local luminosity function at the three BLAST wavelengths. Finally, we fitted the modified Schechter function that Saunders et al. (1990; Equation (6.1)) found was a good fit to the 60-μm luminosity function to the three BLAST luminosity functions. This is the function plotted in Figures 8–10.

Dye et al. (2009) presented a 5σ catalog of BLAST sources, which we have also used in this paper. This catalog is affected by Eddington bias (Section 4.4), and we have used a Monte Carlo simulation of the BLAST fields to investigate the effects of Eddington bias on the catalog, and in particular to investigate the influence of these effects on our estimates of the luminosity function described in Section 4.

We used the source counts from Patachon et al. (2009) and our maps of instrumental noise to generate Monte Carlo realizations of the deep BLAST image and of the shallower wide image. We have not incorporated any clustering of the BLAST sources, since Patanchon et al. (2009) concluded that clustering has only a small effect in models of the effect of Eddington bias. After generating the maps, we used the source-finding software used on the real BLAST images to generate a catalog of 5σ sources. We call these sources the "output sources." For each of these sources we then found all the sources that were used to create the realization within 0.5×FWHM (full-width half maximum of the telescope beam) of the position of the output source; we call these sources the "input sources."

In estimating the luminosity function (Section 4.4), we used fluxes measured from the BLAST maps at the positions of the counterparts. We modeled this procedure by making the assumption that the brightest of the input sources represents the true submillimeter emission associated with the counterpart. We remeasured the submillimeter fluxes from the artificial submillimeter maps at the position of the brightest input source. These fluxes, which we call F_out, thus represent the fluxes we used to estimate the luminosity functions shown in Figure 8. We call the true fluxes of these sources F_in. Figure A1 shows the ratio of F_in to F_out plotted against the signal-to-noise of the output source for both the wide and the deep images.

Zoom In Zoom Out Reset image size

These two figures show two effects of Eddington bias. First, the panels for the wide and deep surveys suggest that some of the 5σ sources in both regions are either instrumental noise, promoted by Eddington bias to appear as actual sources, or a confused combination of instrumental noise and many faint sources. This is most apparent in the figure for the wide survey, which is dominated by instrumental rather than confusion noise. This figure shows a cluster of sources with ${F_{\rm in} \over F_{\rm out}} \simeq 0$, which can clearly not be associated with a single luminous source. This feature can also be seen in the figure for the deep survey (see the inset to the right-hand figure). An alternative empirical way to determine the fraction of spurious sources is, of course, to look for optical counterparts: the number that does not have counterparts gives an upper limit to the percentage of sources that are effectively instrumental noise (Table 4). The results of this simulation are why we suggest in Section 4.4 that some of the sources with missing counterparts may actually not be genuine sources. If this is true, the luminosity functions on the left-hand sides of Figures 8 and 9 are likely to be the best estimates of the genuine luminosity functions.

The second effect is the bias on the fluxes measured at the counterpart. Although we adopted this procedure to mitigate the effect of flux-boosting, flux-boosting can still clearly be a big effect. To quantify this effect, we have measured the average value of $F_{\rm in} \over F_{\rm out}$ as a function of the signal to noise of the output source. We excluded all sources with ${F_{\rm in} \over F_{\rm out}} < 0.2$, which we argue are either instrumental noise or that do not represent single luminous sources. Figure A2 shows the results at 250 μm for the deep survey. We have used tables constructed from figures like this to "deboost" the fluxes used to estimate the luminosity functions shown in Figure 9.

Zoom In Zoom Out Reset image size

There is also a third effect, which is that an output source is composed of more than one input source. Figure A3 shows a histogram of the ratio of the true fluxes of the second brightest and brightest input source. The figure suggests that most output sources are dominated by a single input source, although 21% of the sources in the artificial catalog made for the deep survey have a second input source that is over 50% of the brightness of the brightest input source. We have made no correction for this effect, although a simple thought experiment suggests that this effect effectively operates in the opposite direction to the flux-boosting effect. Suppose an output source is composed of three input sources of equal brightness. If we make the assumption that these sources also all have the same redshift, the correction we should make to a point on the luminosity function is to move it to a luminosity that is three times lower and to a number-density that is three times higher. Since this correction is roughly parallel to a typical luminosity function, the net effect is relatively small. Therefore, by only correcting for the flux-boosting effect in Figure 9, we are essentially putting an upper limit on the effect of Eddington bias on the luminosity function.

Zoom In Zoom Out Reset image size