We present a study of the far-infrared (IR) properties of a stellar mass selected sample of 1.5 < z < 3 galaxies with log (M*/M⊙) > 9.5 drawn from the Great Observatories Origins Deep Survey (GOODS) Near Infrared Camera and Multi-Object Spectrometer (NICMOS) Survey (GNS), the deepest H-band Hubble Space Telescope survey of its type prior to the installation of Wide Field Camera 3 (WFC3). We use far-IR and submm data from the Photoconductor Array Camera and Spectrometer (PACS) and Spectral and Photometric Imaging Receiver (SPIRE) instruments on-board Herschel, taken from the PACS Evolutionary Probe (PEP) and Herschel Multi-Tiered Extragalactic Survey (HerMES) key projects, respectively. We find a total of 22 GNS galaxies, with median log (M*/M⊙) = 10.8 and z = 2.0, associated with 250 μm sources detected with signal-to-noise ratio (SNR) > 3. We derive mean total IR luminosity log LIR(L⊙) = 12.36 ± 0.05 and corresponding star formation rate (SFR)IR + UV = (280 ± 40) M⊙ yr−1 for these objects, and find them to have mean dust temperature Tdust ≈ 35 K. We find that the SFR derived from the far-IR photometry combined with ultraviolet (UV)-based estimates of unobscured SFR for these galaxies is on average more than a factor of 2 higher than the SFR derived from extinction-corrected UV emission alone, although we note that the IR-based estimate is subject to substantial Malmquist bias. To mitigate the effect of this bias and extend our study to fainter fluxes, we perform a stacking analysis to measure the mean SFR in bins of stellar mass. We obtain detections at the 2–4σ level at SPIRE wavelengths for samples with log (M*/M⊙) > 10. In contrast to the Herschel detected GNS galaxies, we find that estimates of SFRIR + UV for the stacked samples are comparable to those derived from extinction-corrected UV emission, although the uncertainties are large. We find evidence for an increasing fraction of dust obscured star formation with stellar mass, finding , which is likely a consequence of the mass–metallicity relation.
Star formation rates (SFRs) in galaxies can be measured using many different methods (see e.g. Kennicutt ). The most easily accessible tracer at high redshift (z > 1) is rest-frame ultraviolet (UV) emission, which correlates with the number of young, massive stars and hence the global SFR of a galaxy. However, in dusty galaxies, this requires a significant correction due to absorption of UV photons by dust, which can be estimated using the correlation between the UV and far-infrared (IR) luminosity ratio (LIR/LUV, where LIR is conventionally defined over the wavelength range ) and the UV slope (β; typically determined from a power-law fit of the form fλ ∝ λβ between 1500 and 2800 Å), which has been measured from local starburst galaxies (e.g. Meurer, Heckman & Calzetti ; Calzetti et al. ). Observations at far-IR wavelengths are generally thought to quantify the amount of obscured star formation more directly, as UV radiation associated with young stellar populations is absorbed by interstellar dust and re-emitted at far-IR wavelengths, and have revealed that much of the star formation activity that occurred at z > 1 is obscured (e.g. Le Floc'h et al. ; Pérez-González et al. ; Caputi et al. ; Magnelli et al. , ).
Observations over the last decade spanning a wide range in redshift and galaxy environments have shown that stellar mass is a key parameter for predicting the properties of a given galaxy. At low redshift (z < 0.1), the most massive galaxies tend to be red and located in denser environments than bluer, lower mass galaxies (e.g. Baldry et al. ). Although the colour–density relation weakens as redshift increases, a strong colour–mass relation is still seen at z ∼ 2 (e.g. Grützbauch et al. ). For galaxies which are actively forming stars, SFR is seen to be correlated with stellar mass up to z ∼ 3 (e.g. Daddi et al. ; Magdis et al. ; Oliver et al. ; Bauer et al. ; Karim et al. ; Rodighiero et al. ). Environment, while certainly important (as seen by the dominance of early type, passively evolving galaxies in clusters), seems to be more weakly correlated with other galaxy properties in comparison to stellar mass, particularly at high redshift (e.g. Peng et al. ; Grützbauch et al. ). This suggests that studies of the assembly of stellar mass, much of which occurs in obscured bursts of star formation, are crucial for developing our understanding of the galaxy formation process.
In this paper we use far-IR photometry from the Herschel Space Observatory (Pilbratt et al. ): Herschel Multi-Tiered Extragalactic Survey (HerMES; Oliver et al. ) and Photoconductor Array Camera and Spectrometer (PACS) Evolutionary Probe (PEP; Lutz et al. ) key projects to investigate obscured star formation in a stellar mass selected galaxy sample: the Great Observatories Origins Deep Survey (GOODS) Near Infrared Camera and Multi-Object Spectrometer (NICMOS) Survey (GNS; Conselice et al. ). The GNS sample is selected in the H band and is estimated to be complete for galaxies with stellar masses down to log (M*/M⊙) = 9.5 at z < 3 (Conselice et al. ; Grützbauch et al. ; Mortlock et al. ). Bauer et al. () carried out a study of star formation activity in the GNS sample over the redshift range 1.5 < z < 3. This coincides with the peak of cosmic star formation activity as measured in the UV (e.g. Bouwens et al. ); note, however, that in the IR a flat plateau in the SFR density is seen from 1 < z < 2 (e.g. Béthermin et al. ; Magnelli et al. ). The Bauer et al. () study primarily used rest-frame UV luminosity (corrected for extinction according to the UV slope) to estimate SFRs. In addition, they estimated obscured SFRs for the ≈20 per cent of their sample that were detected at 24 μ m using the Multiband Imaging Photometer for Spitzer (MIPS) instrument on board Spitzer, finding that the inferred total SFR (SFRIR + UV) is on average 3.5 times larger than the SFR derived from the UV-slope extinction-corrected UV flux (SFRUV, corr). This factor of 3.5 may be overestimated, as several previous studies have shown that while 24 μ m flux densities can be reasonably extrapolated to measure LIR (and hence SFRIR) for galaxies at z < 1.5, this is not the case at higher redshift (e.g. Papovich et al. ; Murphy et al. , ), where LIR as estimated from 24 μ m photometry alone can be a factor of ∼5 higher than LIR measured for the same sources when additional longer wavelength photometry is available to constrain the spectral energy distribution (SED) fits. The discrepancy is greater for ultraluminous infrared galaxies (ULIRGs, which have LIR > 1012 L⊙). Similar results have been reported in studies using Herschel data (e.g. Elbaz et al. , ; Nordon et al. , ).
Star formation in the massive (M* > 1011 M⊙) galaxies on which most of the GNS fields are centred (see Section ) has been investigated using far-IR data from the Balloon-borne Large-Aperture Submillimeter Telescope (BLAST; Viero et al. ) and Herschel (Cava et al. ), who found that disc-like galaxies (selected by the use of the Sérsic index) have significantly higher SFRs than spheroidal-like galaxies. In this work we aim to improve the characterization of obscured star formation as a function of stellar mass at 1.5 < z < 3, using the combination of Herschel photometry and the wide stellar mass range spanned by the full GNS sample (log (M*/M⊙) > 9.5).
The structure of this paper is as follows. In Section we give a brief overview of the GNS and the Herschel data used in this work. We investigate the properties of the GNS galaxies detected at 250 μm using Herschel in Section . We extend the study to lower luminosity galaxies through a stacking analysis which is presented in Section . We present our conclusions in Section .
The galaxy sample used in this work is taken from the GNS (Conselice et al. ), which consists of 60 F160W (H band) pointed observations in the GOODS fields (Giavalisco et al. ) using the NICMOS instrument on-board the Hubble Space Telescope (HST). The footprint of the GNS overlaid on the Spectral and Photometric Imaging Receiver (SPIRE) 250 μm maps is shown in Fig. 1. Each GNS field is ≈50 arcsec on a side, and covers the region around one or more massive galaxies (M* > 1011 M⊙) at 1.7 < z < 2.9, initially selected using a variety of colour selection techniques: distant red galaxies (DRGs; Papovich et al. ), Infrared Array Camera (IRAC) extremely red objects (IEROs; Yan et al. ) and BzK galaxies (Daddi et al. ). While this selection is not homogeneous, Conselice et al. () show that this combination of colour selection techniques leads to an almost complete sample of massive (M* > 1011 M⊙) galaxies: no single one of these colour selection methods selects more than 70 per cent of the massive galaxy population that would be selected in a photometric redshift survey, while a subsequent stellar mass selection in these fields based on photometric redshifts found an almost identical massive galaxy sample to the initial colour-based selection (Conselice et al. ).
In addition to providing high-resolution near-IR photometry of the massive galaxies targeted in each GNS pointing, the depth of the survey allows galaxies with much lower stellar masses to be detected: GNS is complete for galaxies with stellar masses down to log (M*/M⊙) = 9.5 at z < 3 (Grützbauch et al. ; Mortlock et al. ). The stellar mass measurements are described in detail in Conselice et al. (); briefly, a grid of Bruzual & Charlot () stellar population models, with exponentially declining star formation histories (τ-models, with 0.01 < τ(Gyr) < 10), spanning a wide range in metallicity (−2.25 < [Fe/H] < +0.56), were fitted to the BVizH photometry for each galaxy.
In this paper we use a sample of 860 1.5 < z < 3 galaxies with log (M*/M⊙) > 9.5 drawn from the GNS (the redshift range is chosen to match previous analyses of this catalogue presented in e.g. Bauer et al. ; Grützbauch et al. ; Mortlock et al. ). We include galaxies with both spectroscopic and photometric redshifts, using the former where possible. We do not cut galaxies with low photometric redshift probability (P, the χ2 probability outputted by hyperz, the code used to compute the GNS photometric redshifts; Bolzonella, Miralles & Pelló ), because a comparison of the spectroscopic and photometric redshifts showed that the scatter of the residuals is similar regardless of the cut in P (σz = 0.0451 when using only galaxies with P > 95 per cent, compared to σz = 0.06 using the full sample; see Bauer et al. ; Grützbauch et al. ). Note that 450 galaxies in this sample have P > 95 per cent.
To reduce contamination of the sample by active galactic nucleus (AGN), we remove galaxies found within a 2 arcsec matching radius of X-ray sources listed in the 2 Ms Chandra catalogues of Alexander et al. (, GOODS-N) and Luo et al. (, GOODS-S). These catalogues have flux limits of ≈1.4 × 10−16 erg cm−2 s−1 in the 2–8 keV band, and are therefore deep enough to allow sources brighter than erg s−1 to be detected at z ∼ 2 (assuming a power-law spectrum with Γ = 2).
Later in this paper, we measure SFRs for GNS galaxies from the Herschel IR data and compare these with UV-based SFR measurements from Bauer et al. () for the same galaxy sample. Here, we briefly summarize the method used to estimate these UV-based SFRs.
Bauer et al. () estimated unobscured UV SFRs from K-corrected Advanced Camera for Surveys (ACS) z850-band flux measurements, applying the SFRUV–L2800 relation of Kennicutt (), where L2800 is the UV luminosity at 2800 Å. These were corrected for obscuration by dust using the UV slope (β) to estimate the amount of extinction, where β was measured from the 1600 to 2800 Å luminosities of the best-fitting model SED for each galaxy. A similar methodology to Calzetti et al. () was used to convert β values into extinction estimates at 2800 Å. The typical uncertainty on the UV-slope extinction-corrected SFR estimates (SFRUV, corr) is ∼30 per cent (Bauer et al. ).
The Herschel photometry used in this work is taken from two key projects. The PACS (Poglitsch et al. ) Evolutionary Probe (PEP; Lutz et al. ) provides 100 and 160 μm data covering both GOODS fields, as well as 70 μm coverage of GOODS-S. Simulations show that in GOODS-N, the flux limits at 80 per cent completeness are 4.5 and 7.0 mJy at 100 and 160 μm, respectively, while in GOODS-S the corresponding limits are 1.5, 2.0 and 4.8 mJy at 70, 100 and 160 μm. We also use 250, 350 and 500 μm SPIRE (Griffin et al. ) imaging data which were obtained as part of the Herschel Multi-Tiered Extragalactic Survey (HerMES;2 Oliver et al. , ). Unlike the PACS data, the SPIRE data are dominated by confusion noise from unresolved background sources. The calibration of the SPIRE instrument is described in Swinyard et al. ().
Photometry was performed on all the Herschel maps, using prior positions derived from the MIPS 24 μm catalogue of Magnelli et al. () for source extraction. This 24 μm catalogue is extracted from the GOODS-Legacy program observations (PI: M. Dickinson), and reaches a 5σ depth of about 30 μJy. Note that by requiring a 24 μm detection for source extraction in the Herschel maps, a small fraction of sources will be missed at the GOODS depth (<10 per cent; e.g. Roseboom et al. ; Magdis et al. ; Béthermin et al. ). A blind extraction might be able to find such sources, at the expense of significantly noisier photometry due to source blending. Fluxes in the PACS maps were measured by fitting scaled point spread functions (PSFs) at each object position, as in Magnelli et al. (). In the case of the longer wavelength HerMES data, photometry was performed on all sources simultaneously, with the 24 μm catalogue being used to provide reliable deblending, using a slightly modified version of the method described in Roseboom et al. (). The changes to the method are described in Roseboom et al. (); briefly, a global (rather than local) background estimate was used in producing the catalogues used in this work, and a different (and faster) model selection algorithm was used in the fitting procedure. Using this deblending method, reliable fluxes can be extracted close to the formal ≈4–5 mJy SPIRE confusion noise (measured after a 3σconf source cut, where σconf is the confusion noise measured without this cut; Nguyen et al. ). The 24 μm prior positional information reduces the impact of confusion noise, and so the approximate 3σ limit for the SPIRE catalogue at 250 μm used in this work is ≈9 mJy in both fields. We use this catalogue to investigate the properties of GNS galaxies detected at 250 μm in Section .
In Section we present a stacking analysis of GNS galaxies in bins of stellar mass, and we use data from other infrared surveys to broaden the wavelength coverage outside of the Herschel bands. In both the GOODS-N and GOODS-S fields we use Spitzer MIPS 24 μm maps, taken from the Far Infrared Deep Extragalactic Legacy Survey (FIDEL DR2; PI: Mark Dickinson; for GOODS-S) and the GOODS-Spitzer survey (for GOODS-N). In addition, in GOODS-N we make use of the combined AzTEC/MAMBO 1160 μm map of Penner et al. (), while in GOODS-S we use the 870 μm Large Apex Bolometer Camera (LABOCA) map from LABOCA Extended Chandra Deep Field-South (ECDFS) Submillimetre Survey (LESS; Weiß et al. ). To simplify the stacking analysis, the MIPS and PACS maps (in surface brightness units) are cross-correlated with the appropriate area normalized PSF such that each pixel in the resulting map represents the maximum likelihood flux density (in Jy) of an isolated point source at that position. For the publicly available AzTEC/MAMBO and LESS maps, this operation has already been performed.
PROPERTIES OF SPIRE DETECTED GNS GALAXIES
We cross-match the GNS catalogue with the HerMES/PEP catalogue using a simple 2 arcsec matching radius. Since the HerMES/PEP catalogue was extracted using MIPS 24 μm prior positions, a small matching radius, appropriate to the astrometric accuracy achievable with MIPS at 24 μm, can be used (e.g. Bai et al. ). We select robust detections at 250 μm from the catalogue using cuts of S250 > 3ΔS250, where ΔS250 is the flux uncertainty (including confusion noise), i.e. S250 > 8–9 mJy (see Section ), and χ2 < 5 (i.e. the goodness of fit of the source solution within the neighbourhood of the source; see Roseboom et al. ). We find that a total of 22 GNS galaxies with 1.5 < z < 3 and log (M*/M⊙) > 9.5 are matched across both the GOODS-N and GOODS-S fields; this corresponds to ≈2.5 per cent of the GNS sample within these stellar mass and redshift cuts. We note that if we repeat the selection at 350 μm, we obtain a sample of 14 objects, only one of which is not in common with the 250 μm selected sample. This additional source is ID 283 in the GNS catalogue, and has photometric redshift zp = 1.55 ± 0.15 and stellar mass log (M*/M⊙) ≈ 10.6.
Fig. 2 shows 10 × 10 arcsec2 NICMOS F160W postage stamp images centred on each detected GNS galaxy, with the position of the HerMES source and the 2 arcsec matching radius indicated. In almost all cases each GNS galaxy is unambiguously identified with the HerMES source; there are only two cases (IDs 4180 and 5310) where two galaxies of similar brightness are located within the matching circle. We estimated the fraction of potentially spurious matches by randomizing the positions of the submm sources and repeating the cross-matching procedure 1000 times. We found a mean number of 3 ± 2 of the 250 μm sources were randomly associated with GNS galaxies in this test (where the uncertainty is the standard deviation). This can be treated as an upper limit, as it assumes no correlation between objects detected in the submm and near-IR – and so the real fraction of spurious matches is likely to be lower.
Table 1 lists the properties (redshift, stellar mass, rest-frame colour) and flux densities of the individual detected sources. The median redshift of the detected objects is z = 2.02, and the median stellar mass of the detections is log (M*/M⊙) = 10.8. We note that in comparison to the bulk of the GNS sample (Section ), these objects typically have lower photometric redshift probabilities, with median P = 61.
Fig. 3 shows the location of the detected objects in the (U − B) colour–stellar mass plane. Clearly, relatively more massive galaxies with red rest-frame (U − B) colours are detected, as shown in Fig. 4. We find that roughly 13 per cent of the sample with log (M*/M⊙) > 11 and (U − B)rest > 0.85 (the fiducial colour criterion adopted for dividing quiescent and star-forming galaxies in Kriek et al. ) are detected at 250 μm. Given their far-IR flux densities, these objects are clearly not quiescent, and we expect them to have high dust masses and high SFRs, with their red colours being as a result of dust extinction. However, it is possible that the dominant origin of the IR emission is hot dust associated with AGN, rather than star formation, although this is not likely: e.g. Symeonidis et al. () found that all of their 70 μm selected galaxy sample were primarily powered by star formation. Although X-ray AGN were removed from the sample at the outset (Section ), we checked for additional AGN using colours in the Spitzer IRAC bands (Stern et al. ), using data from the GOODS Spitzer Legacy program (Dickinson et al. ). Fig. 5 shows the [3.6]–[4.5], [5.8]–[8.0] colour–colour plot of the 250 μm detected GNS galaxies. We find that six objects fall within the region typically occupied by AGN. We do not remove these objects from the sample, as some studies have shown that AGN mainly contribute to the IR flux at wavelengths (Netzer et al. ; Mullaney et al. , see also Hatziminaoglou et al. ); we will instead note these objects in the following analysis (see also Section ).
We note that it is possible that the presence of either an AGN or starburst may lead to the stellar masses of some of the detected sources being overestimated. Other studies, which explicitly correct for the effect of power-law emission from AGN, find that neglecting such corrections can lead to differences of 10–25 per cent in stellar mass estimates of submm galaxies (e.g. Hainline et al. ). We show in Section that more sophisticated SED modelling, using rather different assumptions to those used in deriving the GNS stellar masses, verifies that the 250 μm detected GNS galaxies are genuinely massive systems (see also the discussion concerning stellar mass estimates of AGN hosting GNS galaxies in Bluck et al. ).
SED fitting). We also fit the SEDs using the templates of Chary & Elbaz (, hereafter CE01), as a consistency check on our results.
We fit the SEDs using χ2 minimization, allowing the dust temperature to vary in the range 10–70 K. We ignore the 24 μm flux densities when fitting the SEDs using models of the form of equation , since at z > 1.5 we do not expect the modified blackbody model to be a reasonable description of the SED at this wavelength in the observed frame. However, we do include the 24 μm fluxes when fitting to the CE01 templates, as these include the contribution from polyaromatic hydrocarbon (PAH) features. Note that we include SED points with SNR <3 in the fitting – given the requirement of a 24 μm detection and prior position, so long as the uncertainties on these points are accurately estimated, then the additional information they provide should help to better constrain the SED than either neglecting these points, or replacing them with 3σ upper limits. We comment on the effect of this on our results in Section .) IMF. We therefore apply a correction of −0.23 dex to SFRs estimated using equation to account for the Chabrier () IMF assumed in this work (see e.g. Kriek et al. ). ). There are many caveats for the dust mass estimates obtained in this way, such as the uncertainty in the value of κ250; the fact that equation can underestimate the true dust mass due to the presence of warm dust in galaxies being neglected in the modified blackbody model (equation ) and the large K-correction to the redshift range of our study. Although the absolute values of Mdust are highly uncertain, we use the relative values obtained by this method to give an indication of the relation of Mdust with M*, assuming that the dust properties are similar in galaxies of different stellar mass in our redshift range of interest (see Section ).
We estimate errors on the parameters derived from the SED fits using Monte Carlo simulations. For each observed SED we generate 1000 random realizations, assuming that the errors on the fluxes are Gaussian. For objects with only photometric redshifts, we simultaneously randomize the redshift of the fitted model SED according to the scatter of σz = 0.06 measured by Grützbauch et al. (). We adopt the 68.3 percentile range from the distribution of parameter values obtained from the random realizations as the corresponding ±1σ uncertainty.
The SED fitting shows that the GNS galaxies individually detected in HerMES are ULIRGs, spanning the range 11.9 < log LIR(L⊙) < 12.9, with mean log LIR = 12.36 ± 0.05 L⊙, where the quoted uncertainty is the standard error on the mean. We estimate total SFRs for these galaxies under the assumption that this corresponds to the sum of the SFR derived from the far-IR SEDs and the UV-based unobscured SFR measurements from Bauer et al. (). We find that the mean total SFR for these galaxies is SFRIR + UV = 280 ± 40 M⊙ yr−1. Removing the six galaxies with IRAC colours consistent with AGN has no significant effect: with these objects excluded, we find SFRIR + UV = 260 ± 50 M⊙ yr−1. This is a factor of >2 larger than the mean UV-slope extinction-corrected SFR estimates from Bauer et al. () for these same galaxies, i.e. SFRUV, corr = 120 ± 30 M⊙ yr−1. We obtain results within <1σ of these values for all of these properties if we take into account the fraction of potential spurious matches (Section ) in a Monte Carlo fashion.
We checked the sensitivity of these estimates to the adopted submm selection criteria. We find consistent results for the smaller sample of eight galaxies detected with SNR > 5 at 250 μm (mean log LIR = 12.39 ± 0.09 L⊙, mean SFRIR + UV = 290 ± 60 M⊙ yr−1), and for the sample of 14 galaxies detected at SNR > 3 at 350 μm (mean log LIR = 12.34 ± 0.07 L⊙, mean SFRIR + UV = 260 ± 40 M⊙ yr−1). We also checked the effect of including SED points with SNR < 3 in the fits (see Section ) – replacing them with 3σ upper limits, we obtain mean log LIR = 12.40 ± 0.05 L⊙, with corresponding mean SFRIR + UV = 300 ± 40 M⊙ yr−1, for the whole sample of 22 galaxies.
Dividing the sample by rest-frame colour, we see no evidence for different IR properties for galaxies detected at 250 μm with red or blue colours, although of course the sample is very small. We find mean log LIR = 12.33 ± 0.09 (SFRIR + UV = 270 ± 60 M⊙ yr−1) for the 11 galaxies with (U − B)rest > 0.85, and mean log LIR = 12.34 ± 0.06 (SFRIR + UV = 260 ± 40 M⊙ yr−1) for the 11 galaxies with (U − B)rest < 0.85.
We conclude that SFRIR + UV is significantly higher than SFRUV, corr for our sample. Wuyts et al. () also found that SFRUV, corr is underestimated compared to SFRIR + UV for galaxies with similar total SFRs and redshifts to our sample. However, several other recent studies find the reverse situation. For example, Murphy et al. () observed a sample of 0.66 < z < 2.6 24 μm selected sources with additional 70 μm photometry, and found that their measurements of SFRUV, corr are a factor of >2 higher than SFRIR + UV. They concluded that the dust corrections applied to their sample (from the Meurer et al. relation) were overestimated for many objects. Nordon et al. () found similar results from a study using PACS observations of massive galaxies at 1.5 < z < 2.5 in GOODS-N, finding SFRUV, corr is overestimated by a factor of about 2 for galaxies with SFRUV > 40 M⊙ yr−1, assuming a Calzetti UV attenuation law (note however that Wuyts et al. showed that this result may in part be driven by the relatively bright Ks < 22 limit adopted in Nordon et al. ). Buat et al. () reached similar conclusions from a study of 250 μm selected z < 1 galaxies from HerMES with UV photometry from the Galaxy Evolution Explorer (GALEX) satellite. At lower redshift (z < 0.35), Wijesinghe et al. () found only a weak correlation with large scatter between the UV slope (β) and LIR/LUV, which would also lead to overestimated SFRUV, corr. However, for UV selected samples (e.g. Lyman break galaxies; LBGs) which are not ULIRGs, dust corrections from the local Meurer et al. () relation appear to be valid at z ∼ 2 (e.g. Overzier et al. ; Reddy et al. ). Reasonable agreement between SFRUV, corr and SFR derived from stacked radio and 24 μm observations is also seen up to z ∼ 3 for LBGs (Magdis et al. ).
We expect large IR-derived SFRs for the galaxies we detect at 250 μm given their redshift and the 3σ flux limit, which is ≈9 mJy at 250 μm in the GOODS-N field. This leads to a large Malmquist bias (with some flux boosting due to the low SNR) in comparison to the UV-derived SFRs, which reach to ∼1 M⊙ yr−1 (Bauer et al. ). Fig. 6 shows the SFRIR limit as a function of redshift for a modified blackbody model SED (equation ) with Tdust = 35 K, normalized to a 250 μm flux density of 9 mJy. Highlighted in this plot are the SFRIR and SFRUV, corr values for the SPIRE-detected galaxies; and clearly in most cases SFRUV, corr is much lower than the fiducial SFRIR corresponding to the 250 μm flux limit. This makes the comparison between these two SFR measures for our sample difficult to interpret. There is one clear exception, where SFRUV, corr is roughly a factor of 3 larger than SFRIR – this is ID 5918, which, from inspection of the ACS imaging, seems to be a multiple component merger system, with regions of significant unobscured star formation (see Fig. 7). It may be that only one component of this system is the source of the FIR emission, but it is not possible to determine which using the current data.
Fig. 8 shows the comparison of SFRIR + UV and M* for the SPIRE-detected GNS galaxies with the wider GNS sample, where for the latter SFRUV, corr is used as the estimate of the total SFR. We see that almost all of the SPIRE-detected galaxies scatter above the SFR–M* relation measured by Daddi et al. (), which is as expected given the approximate SFRIR limit shown in Fig. 6.
For the 16 galaxies with flux measurements in all SPIRE bands, we find dust temperatures in the range 23–48 K, with mean 35 ± 6 K (where the quoted uncertainty is the standard deviation). Note however that only four of these galaxies have SNR > 3 in all SPIRE bands, and so the individual temperature estimates are poorly constrained, with typical statistical uncertainty ≈5 K. We find that replacing the SNR < 3 SED points in the fits with 3σ upper limits (see Section ) gives Tdust values for individual galaxies in this subsample that agree within <1σ of the values obtained when the low SNR SED points are included. For a sample selected with SNR > 3 at 350 μm, we find mean Tdust = 33 ± 7 K, while for a sample with SNR > 5 at 250 μm, we find mean Tdust = 34 ± 7 K. The single GNS galaxy which is detected at SNR > 3 at 350 μm but is not in our 250 μm selected sample (ID 283; see Section ) has a slightly lower dust temperature (Tdust = 20 ± 5 K).
The mean temperature we find is somewhat lower than the typical temperature of ULIRGs at z < 1 (Tdust ≈ 42 K; e.g. Yang et al. ; Clements, Dunne & Eales ); although note that β is fitted for in the former work, whereas in the latter it is fixed at β = 1.5, as we assume here. This is not unexpected given the high redshift of the sample and the selection at SPIRE wavelengths (Symeonidis, Page & Seymour ). The dust temperatures we find are similar to those found for other samples at z > 1 (Chapin et al. ; Amblard et al. ; Chapman et al. ; Hwang et al. ). Adopting β = 2.0 in the modified blackbody model (equation ) gives mean Tdust about 4 K lower.
We find a mean dust mass for these galaxies of Mdust ∼ 3 × 108 M⊙, which is comparable to the characteristic mass in the dust mass function of measured at z ∼ 2.5 by Dunne, Eales & Edmunds (, note the value quoted here is taken from table 3 of Dunne et al. ). However, the range in Mdust spans more than an order of magnitude, and the individual values are highly uncertain. The median Mdust/M* ratio for these galaxies is ∼5 × 10−3 and spans the range 4 × 10−4–3 × 10−2. This is similar to the range found by Rowan-Robinson et al. (), with a sample reaching to z ∼ 2 and using a different method to estimate Mdust. Fixing the value of β = 2.0 in the modified blackbody model (equation ) would increase the mean dust mass that we find by ≈60 per cent.
Joint optical–IR SED fitting
We tested the sensitivity of the results described above to the simple modified blackbody model used in the SED fitting (Section ) by jointly fitting the full optical–IR SEDs (BVizH from HST, IRAC channels 1–4, MIPS 24 μm, plus the Herschel photometry) using cigale (Noll et al. ), a code which fits the attenuated optical light from stars and dust emission associated with star formation and AGN simultaneously. The available models for use within cigale differ from those assumed for deriving the GNS stellar masses (see Section ) and the SFRs estimated in this work (see Section ). We used the Maraston () stellar population models to fit the optical part of the spectrum, the Dale & Helou () templates to fit the dust emission, and a Kroupa () IMF. Some example SED fits are shown in Fig. 9.
We find that cigale gives stellar masses that span the range 10.0 < log (M*/M⊙) < 11.5, with median log (M*/M⊙) = 10.9, confirming that these systems have high stellar masses, as measured in the GNS using a different SED fitting code (Conselice et al. ). A two-sample Kolmogorov–Smirnov (KS) test reveals that the stellar mass distributions are not significantly different (p = 0.33), although there is a scatter of 0.23 dex in the residuals between the two stellar mass estimates for each galaxy.
The SFRs estimated by cigale are systematically lower than the results obtained using the modified blackbody model, presumably as a result of the different stellar population model, IMF, and dust emission spectral templates used, but the mean SFR (210 ± 30 M⊙ yr−1) is similar to that found from the modified blackbody SED fits, despite this. The fraction of the IR luminosity due to warm dust associated with AGN estimated by cigale spans the range 3–30 per cent, with median ≈5 per cent. This suggests that star formation is the primary source of the IR emission in these objects, as found in other studies (e.g. Netzer et al. ; Mullaney et al. ).
As shown in Section , only 2.5 per cent of the 1.5 < z < 3, galaxy sample is detected in the 250 μm maps used in this work, and the detected galaxies are ULIRGs with large stellar masses (∼1011 M⊙). We therefore performed a stacking analysis to extend our study to galaxies with lower stellar masses and fainter far-IR luminosities. An additional advantage of the stacking analysis is that the results are less biased than those obtained from a small number of sources detected at low SNR. The stacking was performed on maps from which sources were not subtracted. Note that in contrast to the analysis in Section , additional maps at longer wavelengths than SPIRE were used in the stacking analysis (see Section ).
We divide the 1.5 < z < 3 GNS galaxy sample into four bins of stellar mass, reaching to the log (M*/M⊙) > 9.5 limit to which the survey is complete (Grützbauch et al. ; Mortlock et al. ). Fig. 10 shows the location of the mass-limited subsamples in the (M*, z) plane, compared to the full GNS catalogue covering both GOODS fields. Because of the low SNR of the resulting stacked detections (see Section ), we are not able to divide the sample into redshift bins, nor examine subsamples of passive versus actively star-forming galaxies (although note that the latter is investigated using the GNS galaxy sample by Bauer et al. , using UV-based SFR measurements). Table 2 lists the properties of the mass-limited subsamples we stack.
The far-IR data used in this work have low angular resolution, particularly in the SPIRE bands where the beam sizes are 18, 25 and 36 arcsec at 250, 350 and 500 μm, respectively, resulting in relatively large confusion noise. The source densities of GNS galaxies per beam are also large (median nine sources per beam at 250 μm), and if the effect of clustered confused sources is not accounted for, the resulting stacked fluxes will be biased.
We use the global stacking and deblending algorithm of Kurczynski & Gawiser (, hereafter KG2010) to mitigate the effect of this bias (for other approaches to this problem see Béthermin et al. ; Bourne et al. ). We generalized the method to simultaneously stack and deblend all of the mass-limited samples (see Table 2), in addition to two ‘non-target’ samples of objects. The first of these non-target galaxy samples is drawn from the 24 μm catalogue of Magnelli et al. (, which is also used to provide prior positions for source extraction in the Herschel maps used in this paper). This catalogue provides coverage outside of the GNS footprint, and allows infrared bright galaxies beyond the edges of the GNS fields to be deblended. Since 24 μm bright sources are correlated with sources detected at PACS and SPIRE wavelengths, these objects are the most likely to contaminate stacked flux measurements of the mass-limited samples at far-IR wavelengths. The second non-target galaxy sample consists of all GNS galaxies which are not 24 μm sources and not included in the stellar mass selected samples (i.e. with z < 1.5 or >3, and/or log (M*/M⊙) < 9.5). X-ray detected objects that are not included in the stellar mass selected samples (Section ) were also included in this second non-target sample.
We estimate errors on the stacked fluxes by bootstrapping: we run the stacking and deblending algorithm 1000 times, assigning the flux at each object position uniformly at random (with replacement) from the observed fluxes in each sample. During this process, the positions of all sources in the samples are kept fixed, and so the attenuation factors used in deblending sources (αkj in KG2010) remain constant (i.e. it is only the flux values that are bootstrap resampled). We adopt the 68.3 percentile as the uncertainty in the stacked flux. We also estimated errors by jackknifing (i.e. from the distributions of stacked fluxes obtained after removing a single source from each stacking sample in turn), finding slightly smaller error bars – the detection significances inferred using the jackknife error estimates are 0.1–0.2σ higher than those obtained using the bootstrap error estimates.
We test the robustness of the mean stacked flux measurements by randomizing the object positions in each of the stacking samples (both target and non-target samples) and running the stacking algorithm, repeating this process 1000 times. For simplicity, we perform this test using the GOODS-N sample only. We show the results for each of the stellar mass samples in the SPIRE bands (since these are the most likely to suffer from the effects of confusion as they have the largest beams) in Fig. 11. With the exception of the lowest stellar mass bin, we find that the probability of a chance spurious stacked detection is higher for the lower resolution channels. The detection probabilities inferred from this null test are consistent with those obtained from stacking on real object positions and assuming the bootstrap error estimate; the maximum difference is 0.3σ, with detection significances inferred from the random stack tests being higher.
We perform simple simulations to check that we can recover SED parameters such as LIR and Tdust without significant bias. We create simulated maps with the same pixel scales as the real GOODS-N maps and insert Gaussian sources with the appropriate full width at half-maximum (FWHM) for each channel at the positions of real objects in the GNS catalogue. The simulated sources are modelled using the modified blackbody SED (equation ). We note that this is somewhat idealized, as we do not include different SEDs from those used in the fitting procedure.
For the stellar mass selected samples, anticipating the LIR measurements obtained for the real maps (shown in Section ), we set each model SED to have log LIR(L⊙) = 11.0, 11.5, 11.7 and 11.9 for galaxies in stellar mass bins log (M*/M⊙) 9.5–10.0, 10.0–10.5, 10.5–11.0 and >11, respectively. We draw Tdust for each galaxy in each stellar mass subsample from a uniform distribution, with a slightly different (–) range used for each bin: (15–45 K), (20–50 K), (25–55 K), (30–60 K), in ascending order of stellar mass. This ensures that each bin has different mean Tdust, for clarity in the right-hand panel of Fig. 12.
Models for galaxies in the non-target sample of 24 μm bright sources have log LIR(L⊙) = 11, which is the median value we find for these sources when estimating their LIR from their 24 μm flux densities alone (where we estimate LIR for each source as the median value over the full range of CE01 templates). We do not include the non-target galaxies that were not detected at 24 μm in the simulated maps. Each model source is redshifted to its corresponding z in the GNS catalogue. We apply a Gaussian random scatter of (1 + zp) 0.06 in redshift to galaxies with only photometric redshift estimates to simulate the effect of incorrect redshifts, where the amount of scatter is as found by Grützbauch et al. () from a comparison of a subset of GNS galaxies with spectroscopic redshifts (see Section ). For sources in the 24 μm detected non-target sample without redshift information, we assign their model SED a redshift selected at random from the redshift distribution of GNS galaxies detected at 24 μm.
Fig. 12 shows the results of running our stacking and SED fitting code (Section ) on the simulated maps. We find that we recover LIR to within ±30 per cent down to the lowest stellar mass bin. We see that there is a small positive bias in Tdust, with the recovered value being at most about 7 K lower than the mean input Tdust. This bias is absent if we set Tdust to a fixed value for all galaxies in each bin, and is likely to be a consequence of the smearing of the stacked SED shape due to the different redshifts and dust temperatures of the model SEDs that go into each stack.
Table 3 lists the mean stacked flux densities for each stellar mass selected subsample in each field. We find consistent results between the northern and southern fields given the large uncertainties, although the stacked fluxes in the south are typically fainter than in the north for most stellar mass samples (see Fig. 13). The stacked SNR values are low: in the north, we obtain ≈2–3σ detections across almost all SPIRE and PACS bands for only the two most massive stellar mass bins. However, the detection significance increases to ≈4σ in some channels for the second highest log M* bin when the combined sample is used. The SNR in the lowest mass bin is only ≈1σ across the PACS and SPIRE bands when using the combined sample.
Despite the low SNR for each individual SED point, we proceed to fit the SEDs, in order to derive rough estimates of LIR and SFRIR for each stellar mass bin. We include the low SNR points in the fits, rather than excluding them, or treating them as upper limits. Under the assumption that the estimated error bars are reasonable (note that here they are obtained in a consistent way across all wavelengths), this should not bias the fit. We fit the SEDs for each stack using nearly the same method that was used for the SPIRE-detected galaxies (Section ). We make one change to the fitting procedure in order to account for the wide redshift range covered by the galaxy sample: during the Monte Carlo procedure used to estimate error bars on the fitted parameters (i.e. LIR, Tdust, the uncertainties of which feed through to SFRIR and Mdust), we bootstrap sample the redshift applied to the model SEDs from the distribution of redshift values in each stellar mass bin. This approximately doubles the size of the uncertainties on LIR and SFR in comparison to those obtained when the redshift is held fixed at the mean redshift of the galaxy sample. Fig. 14 presents the stacked SEDs and best-fitting results using the modified blackbody templates for the northern, southern and combined samples.
We obtain estimates of SFRIR for each sample with typically a factor of 2 uncertainty, despite the low SNR measurements in each individual band. We find that the difference between the stacked flux densities measured for the GOODS-N and GOODS-S fields (Fig. 13) leads to lower SFRs for most stellar mass bins in the GOODS-S sample. However, there is little tension between the SFRs measured in each field: the largest discrepancy is between the highest stellar mass bins, but even in this case the difference in the SFR estimates is significant only at the <2σ level. We find that the SFRIR estimates obtained using the modified blackbody model and the CE01 templates are consistent.
We estimate mean total SFRIR + UV for the stacked samples by adding to each sample the mean UV-based estimate of unobscured SFR from Bauer et al. () for the same galaxies in each stellar mass bin. Fig. 15 shows the resulting comparison with the mean UV-slope extinction-corrected estimates (SFRUV, corr) from Bauer et al. () for the same galaxies. We see a rough agreement between the two measurements given the large uncertainties, although while in GOODS-N SFRIR + UV is higher than SFRUV, corr, the opposite is true in GOODS-S. Much of this difference comes from a factor of ∼2 difference in SFRUV, corr between the two fields, with SFRUV, corr being higher in GOODS-S than GOODS-N. For all stellar mass bins apart from log (M*/M⊙) > 11, the difference in SFRUV, corr between the fields is significant at the ≈3σ level. The difference in SFRIR + UV between the fields is less significant, at most 1.6σ. Also, in GOODS-S, the highest SFR is seen for the second most massive log M* bin, in both SFRIR + UV and SFRUV, corr, although neither of these SFR estimates are significantly different from those measured for the most massive log M* bin.
We checked for differences between the GOODS-N and GOODS-S samples that could lead to these effects. It is not likely that they arise from different redshift distributions: a two sample KS test gives p = 0.21, i.e. the distributions are not significantly different. Another possibility is environmental effects: the GOODS-S field contains a galaxy overdensity at z = 1.6 (Kurk et al. ) which lies within our redshift range. This structure is thought to be a forming cluster of galaxies, and so the denser environment on average relative to the GOODS-N field may lead to a higher fraction of quiescent galaxies in GOODS-S, and therefore lower average SFR. However, we find that excising the region within 2 Mpc projected radius of this structure makes no significant difference to the derived SFRs. It seems likely that the difference between the results for each field can be ascribed to the small area covered by the GNS.
Fig. 16 shows the SFRIR + UV–M* relation obtained for the stacks in each field. Similarly to Bauer et al. (), we see a shallower relation compared to the SFR–M* relation of Daddi et al. (), who measured for star-forming galaxies at z ∼ 2. This is not surprising, because the GNS sample is selected differently, purely by stellar mass, and therefore includes quiescent in addition to star-forming galaxies (see also Bauer et al. ). Furthermore, the different methods used to measure SFR are also subject to different selection effects. Using weighted least-squares regression, we find the relation16 and are consistent within the errors.
A straightforward comparison of our results with other works is not possible due to differences in sample selection, stellar mass completeness and the wide redshift range used here. Karim et al. () performed a stacking analysis in 1.4-GHz data using a 3.6-μm-selected sample of galaxies in the Cosmic Evolution Survey (COSMOS) field (Scoville et al. ); however, this survey suffers from incompleteness for log (M*/M⊙) < 10.4 at z > 1.5. Attempting a rough comparison of our measurements with this work, we find that our SFRIR + UV estimates are systematically lower, for similar stellar mass and redshift ranges. However, the discrepancy is only significant at the 2–3σ level for our most massive bin (log (M*/M⊙) > 11, where we find SFR a factor of ∼2 less than Karim et al. ), and there is reasonable agreement for the 10.5 < log (M*/M⊙) < 11 bin. Kurczynski et al. () studied a sample of star-forming BzK galaxies (sBzKs) in the ECDFS, comparing several methods of measuring SFR using essentially the same stacking algorithm we used in this work. Their sample was not stellar mass selected, but we find that our SFRIR + UV estimates for log (M*/M⊙) > 10.5 galaxies are in good agreement with their measurements (obtained using IR data from MIPS, BLAST and LESS) at the same redshift as our study, after accounting for the Salpeter () IMF assumed in Kurczynski et al. ().
Ratio of obscured to unobscured star formation and relation to stellar mass
We plot the ratio of obscured to unobscured star formation (SFRIR/SFRUV) as a function of stellar mass for the stacked samples in Fig. 17. Since the uncertainties are large, this is not well constrained from our data. For both GOODS fields combined, we find the relation), who suggest that the mass–metallicity relation is responsible, with higher mass (metallicity) galaxies having larger dust column densities and correspondingly larger SFRIR/SFRUV ratios (see also Pannella et al. ). For galaxies with log (M*/M⊙) > 11, we find the range spanned across the GOODS-N and GOODS-S fields is SFRIR/SFRUV ∼ 6–20. For comparison, Reddy et al. () find SFRIR/SFRUV = 4.2 ± 0.6 for a sample of galaxies observed as part of the GOODS–Herschel project.
Although we derive estimates for Tdust in each stellar mass bin from the SED fits (Section ), they are not well constrained, with uncertainties ∼10 K. All of the stellar mass samples in each field have Tdust in the 20–40 K range, consistent within errors across the stellar mass range, and consistent with the mean value found for the individually detected sources (Section ). We note that simulations suggest that the Tdust estimates from the stacked SEDs may be biased low, perhaps by roughly 7 K (Section ).
The estimates of Mdust we obtain are fairly low in comparison to , the characteristic mass in the dust mass function, as measured by Dunne et al. () for 0 < z < 0.5 and at z ∼ 2.5 by Dunne et al. (). The largest value of Mdust that we measure (≈1.3 × 108 M⊙), corresponding to the log (M*/M⊙) > 11 bin, is a factor of >3 lower than measured by Dunne et al. () at similar z, and also lower than measured at 0.4 < z < 0.5 (Dunne et al. ). This may be as a result of the purely stellar-mass-based sample selection used here, which contains both passive and actively star-forming galaxies; naturally, the samples used in Dunne et al. (, ) consist of galaxies selected in the submm, and are therefore dominated by dusty star-forming galaxies.
The relation we see between Mdust and M* is very poorly constrained (), owing to the large uncertainties in the dust masses, but suggests a mildly decreasing Mdust/M* ratio with increasing stellar mass, with Mdust/M* falling from ∼5 × 10−3 to ∼7 × 10−4 over the stellar mass range 9.5 < log (M*/M⊙) < 11.
We have investigated the far-IR properties of a stellar mass selected sample of 1.5 < z < 3 galaxies drawn from the GOODS NICMOS Survey – the deepest H-band HST survey of its type prior to the installation of the WFC3 instrument – using deep Herschel 70–500 μm photometry from the HerMES and PEP key projects. We found the following.
Only 22 galaxies from the sample are detected at SNR > 3 at 250 μm. They are ULIRGs (median log LIR(L⊙) = 12.4), have high stellar masses (median log (M*/M⊙) = 10.8) and are located at z ≈ 2.
From fitting the SEDs of the SPIRE detected galaxies, we find they have mean SFRIR + UV a factor of >2 higher than the UV-slope extinction-corrected estimates of Bauer et al. (). However, we note that the IR-based SFR estimate suffers from a significant Malmquist bias, making the interpretation difficult. The mean dust temperature of the 16 objects with flux estimates in all HerMES and PEP bands (Tdust = 35 ± 6 K) is slightly lower than found for ULIRGs at z < 1.
Using a stacking algorithm which attempts to deblend sources, we find marginal detections (2–4σ) at SPIRE wavelengths when stacking the galaxy sample in bins of stellar mass, even for the highest stellar mass bins (log (M*/M⊙) > 10.5).
Despite the low SNR of the stacked flux measurements in each band, we obtain estimates of SFRIR for the stacked samples with factor of ∼2 uncertainties. We find that SFRIR + UV measured for the stacked samples is in reasonable agreement with measurements of SFRUV, corr for the same galaxy sample by Bauer et al. ().
We find a relatively shallow slope for the SFR–M* relation () compared to previous studies (e.g. Daddi et al. ), which is likely due to selection effects, as our purely stellar mass selected sample contains a mixture of passive and actively star-forming galaxies.
We find evidence for an increase in the ratio of obscured to unobscured star formation with increasing stellar mass (). This is most likely a consequence of the mass–metallicity relation, with higher mass and metallicity galaxies being more obscured.
Since the far-IR and submm data used in this paper are amongst the deepest that will be obtained by Herschel, it is clear that to make further progress in characterizing the far-IR properties of low stellar mass (log (M*/M⊙) < 10) galaxies at z ∼ 2 using Herschel, a much larger galaxy sample is needed, as will be provided by the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. ; Koekemoer et al. ).
We thank the referee for many helpful comments which have improved this paper. We thank Amanda Bauer for providing the UV-based SFR measurements of GNS galaxies and useful discussions. MH and CJC acknowledge financial support from the Leverhulme Trust and STFC. Support for the GNS was also provided by NASA/STScI grant HST-GO11082.
SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK) and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK) and NASA (USA).
PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy) and IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy) and CICYT/MCYT (Spain).