A Near-infrared, x-ray, & radio matching and k-corrections

Hidden starbursts and active galactic nuclei at from the Herschel1-VVDS-CFHTLS-D1 field: Inferences on coevolution and feedback2

Key Words.:
Galaxies: evolution - Galaxies: high-redshift - Galaxies: starburst - Galaxies: active - Techniques: spectroscopic - Techniques: photometric - Submillimeter: galaxies

We investigate of the properties of 2000 Herschel/SPIRE far-infrared-selected galaxies from in the CFHTLS-D1 field. Using a combination of extensive spectroscopy from the VVDS and ORELSE surveys, deep multiwavelength imaging from CFHT, VLA, Spitzer, XMM-Newton, and Herschel, and well-calibrated spectral energy distribution fitting, Herschel-bright galaxies are compared to optically-selected galaxies at a variety of redshifts. Herschel-selected galaxies are observed to span a range of stellar masses, colors, and absolute magnitudes equivalent to galaxies undetected in SPIRE. Though many Herschel galaxies appear to be in transition, such galaxies are largely consistent with normal star-forming galaxies when rest-frame colors are utilized. The nature of the star-forming “main sequence” is studied and we warn against adopting this framework unless the main sequence is determined precisely. Herschel galaxies at different total infrared luminosities () are compared. Bluer optical colors, larger nebular extinctions, and larger contributions from younger stellar populations are observed for galaxies with larger , suggesting that low- galaxies are undergoing rejuvenated starbursts while galaxies with higher are forming a larger percentage of their stellar mass. A variety of methods are used to select powerful active galactic nuclei (AGN). Galaxies hosting all types of AGN are observed to be undergoing starbursts more commonly and vigorously than a matched sample of galaxies without powerful AGN and, additionally, the fraction of galaxies with an AGN increases with increasing star formation rate at all redshifts. At all redshifts () the most prodigious star-forming galaxies are found to contain the highest fraction of powerful AGN. For redshift bins that allow a comparison (), the highest galaxies in a given redshift bin are unobserved by SPIRE at subsequently lower redshifts, a trend linked to downsizing. In conjunction with other results, this evidence is used to argue for prevalent AGN-driven quenching in starburst galaxies across cosmic time.

1 Introduction

Questions related to the processes responsible for transforming galaxies from the large populations of low-mass, gas-dominated, blue star-forming galaxies found at high redshifts to the massive, quiescent, bulge-dominated galaxies observed in the local universe remain numerous. Gas depletion resulting from normal, continuous, or quasi-continuous star formation extended over a large fraction of the Hubble time, perhaps coupled with an initial phase of vigorous star formation, is a seemingly viable way of explaining the cessation of star-formation activity in galaxies and is invoked by several studies as the primary method of building up stellar mass in galaxies (e.g., Legrand et al. 2001; Noeske et al. 2007a, 2007b; James et al. 2008; Rodighiero et al. 2011; Gladders et al. 2013; Oesch et al. 2013). This scenario is not, however, the only method of transforming the gas-rich galaxies observed in the high-redshift universe (e.g., Daddi et al. 2010a) into those whose baryonic component is dominated by their stellar matter. In some cases, star formation is thought to proceed in a sporadic manner, particularly at high redshifts, with undisturbed, low-level, continuous star formation being the exception in galaxies rather than the primary mechanism of building stellar mass (e.g., Kolatt et al. 1999; Bell et al. 2005; Papovich et al. 2005, 2006; Nichols et al. 2012; Pacifici et al. 2013). Violent episodes of star formation, known as starbursts, can dramatically alter the stellar mass of galaxies with sufficient gas reservoirs during the relatively short period that constitutes the lifetime of a starburst (typically 100-500 Myr; Swinbank et al. 2006; Hopkins et al. 2008; McQuinn et al. 2009, 2010b; Wild et al. 2010; though perhaps up to 700 Myr for particular populations, e.g., Lapi et al. 2011; Gruppioni et al. 2013). Such events contribute significantly to the global star formation rate () of galaxies both in the high-redshift universe (e.g., Somerville et al. 2001; Chapman et al. 2005; Erb et al. 2006; Magnelli et al. 2009) and locally (e.g., Brinchmann et al. 2004; Kauffmann et al. 2004; Lee et al. 2006, 2009) and are furthermore thought to be a typical phase of evolution undergone by massive quiescent galaxies observed at low redshifts early in their formation history (e.g., Juneau et al. 2005; Hickox et al. 2012).

While the operational definition of the term starburst changes significantly throughout the literature (see, e.g., Dressler et al. 1999 vs. Daddi et al. 2010b and the discussion in Lee et al. 2009 and references therein), typically the phenomenon is attributed to galaxies that are prodigious star formers, i.e., forming stars at rates of 100-1000 yr, though dwarf galaxies with substantially lower star formation rates have also fallen under the dominion of this term (see, e.g., McQuinn et al. 2010a). Beyond this common thread, weaving a coherent picture of how starbursting galaxies are formed and how they evolve is difficult due to the large range of properties observed in the population. Though galaxies undergoing a starburst event are commonly associated with a late-type morphology (i.e., spiral, irregular, or chaotic), a non-negligible fraction of starbursts are, surprisingly, associated with early-types (Mobasher et al. 2004; Poggianti et al. 2009, Abramson et al. 2013). Additionally, starbursts are observed in hosts that span a wide range of stellar masses (e.g., Lehnert & Heckman 1996; Elbaz et al. 2011, Hilton et al. 2012; Ibar et al. 2013), color properties (e.g., Kartaltepe et al. 2010b; Kocevski et al. 2011a), and environments (e.g., Marcillac et al. 2008; Galazzi et al. 2009; Hwang et al. 2010a; Lemaux et al. 2012; Dressler et al. 2013). And though some starbursting populations show morphologies indicative of a recent merging or tidal interactions, others appear primarily as undisturbed, grand-design spirals, with little evidence of recent interactions with their surroundings (see, e.g., Ravindranath et al. 2006; Kocevski et al. 2011a; Tadhunter et al. 2011). Thus, it is likely that a variety of different processes are responsible for inducing starburst events in galaxies, though which mechanisms are responsible and how prevalent their effects are, remain open issues.

Just as the relative importance of the processes that serve to induce starbursts in galaxies remain ambiguous, so too are the mechanisms responsible for the cessation of star-formation activity in galaxies. There are many processes capable of negatively feeding back on a galaxy once a starburst is underway, but what role each of these plays in frustrating or quenching star formation is not fully understood. These processes include photoionization by hot stars (Terlevich & Melnick 1985; Filippenko & Terlevich 1992; Shields 1992), galactic shocks (Heckman 1980; Dopita & Sutherland 1995; Veilleux et al. 1995), post-asymptotic giant branch stars (Binette et al. 1994; Taniguchi et al. 2000), and emission from an active galactic nucleus (AGN; Ferland & Netzer 1983; Filippenko & Halpern 1984; Ho et al. 1993; Filippenko et al. 2003; Kewley et al. 2006). The latter is of particular importance, as AGN feedback is often invoked as the primary culprit in transforming the morphologies of galaxies (e.g., Fan et al. 2008; Ascaso et al. 2011; Dubois et al. 2013), in color space (e.g., Springel et al. 2005; Georgakakis et al. 2008; Khalatyan et al. 2008; though see Aird et al. 2012 for an opposing view), and in either stimulating (i.e., positive feedback) or truncating (i.e., negative feedback) star-formation activity in their host galaxies.

However, the observational connection between starbursts or galaxies that have recently undergone a starburst (i.e., post-starburst or K+A galaxies; Dressler & Gunn 1983; Dressler & Gunn 1992) and AGN activity is, to date, largely circumstantial. This is not due to to a lack of studies investigating this connection, nor to the quality of the data in those studies, but is rather due to the difficulty in proving a causal relation between the two phenomena. While AGNs, especially those of the X-ray variety, have been observed in hosts that exhibit transitory properties either in their broadband colors (e.g., Silverman et al. 2008, Hickox et al. 2009, Kocevski et al. 2009b, though see Silverman et al. 2009 for an alternative view) or in their spectral properties (e.g., Cid Fernandes et al. 2004; Yan et al. 2006; Kocevski et al. 2009b; Lemaux et al. 2010; Rumbaugh et al. 2012), determining timescales related to each phenomenon with precision in any particular galaxy remains extremely challenging. And it is this lack of precision which is fatal to directly proving causality (or lack thereof) between AGN activity and the cessation of star formation. Thus, many questions remain about the nature of the interaction between the AGN and its host in such galaxies. Does the AGN play a role in causing this transitional phase or does the nuclei become active after the host has already shutdown its star formation? If feedback originating from an AGN is primarily responsible for quenching star formation in galaxies, why are many transition galaxies observed without a luminous AGN? And how efficient is the quenching provided by an AGN, i.e., on what timescale can the two phenomena be observed coevally?

The difficulty of answering these questions is compounded by the confusion of separating starburst and AGN phenomena observationally. For those studies which have the luxury of integrated spectra and the requisite wavelength coverage, the latter condition becoming increasingly difficult to satisfy at high redshifts, starbursts can be separated from AGN, typically utilizing a BPT (Baldwin, Phillips, & Terlevich 1981) emission line ratio diagnostic or some variation thereof. However, even with high-quality spectral information, attributing the fractional contribution of processes related to star formation and those related to AGN to the strength of the emission lines is challenging. While attempts at quantifying relative contributions have been made using photoionization models (e.g., Kewley et al. 2001), much ambiguity still remains. In the absence of spectra, separating, or quantifying the relative strengths of AGN and starbursts have been attempted through spectral energy distribution (SED) fitting (e.g., Marshall et al. 2007; Lonsdale et al. 2009; Pozzi et al. 2010; Johnson et al. 2013). This method, however, rapidly loses its effectiveness for galaxies whose SEDs are not dominated by one of the two phenomena, a difficulty which limits its usefulness for investigations on the connection between AGN and starbursts. Another approach is to use passbands that are (spectrally) far removed from typical wavelengths that AGN, when present, dominate. Over the past decade this has typically been achieved using Spitzer Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) 24m imaging to observe large populations of starbursting galaxies up to and beyond (e.g., Le Floc’h et al. 2005). Again, however, the difficulty of separating light originating from an AGN and that originating from a starburst persists. The SEDs of AGN obscured by dust, a population that is, not surprisingly, preferentially associated with the dusty starbursts selected by such imaging, peaks at 10-30m, roughly the center of the MIPS passband over a broad redshift range (i.e., ). Given that the reprocessed emission from HII regions in dusty starbursting galaxies peaks at 50-150m (e.g., Chary & Elbaz 2001), in order to minimize contamination from emission originating from an AGN, redder observations are necessary.



Figure 1: Overview of the sky coverage of the observations available on the CFHTLS-D1 field. More details on each observation are given in Table 1. Imaging data is represented by open regions bounded by the various lines. Spectroscopic observations are represented by shaded regions. The overlap region that is used to select our final Herschel sample is defined by the intersection of the CFHTLS optical imaging (dashed blue line), SWIRE NIR/MIR imaging (solid magenta line), and the Herschel far-infrared imaging (solid dark red line). In many cases throughout the paper the sample is further restricted to those galaxies that are covered by WIRDS NIR imaging (solid dark green line). The full extent of the SWIRE coverage is truncated here for clarity, though the coverage over the overlap region is accurately represented.

The recently launched Herschel Space Observatory provides, for the first time, the means to make these observations for statistical samples of starbursting galaxies in the local and distant universe. Already, observations from Herschel have been used to great effect to investigate the relationship between various types of AGN and their starbursting host galaxies. In galaxies with X-ray-selected AGN, numerous studies have found correlations between the strength of the starburst in galaxies, as determined by Herschel, and the strength of the AGN they host (Shao et al. 2010; Harrison et al. 2012b; Rovilos et al. 2012; Santini et al. 2012; though for an alternative viewpoint see Page et al. 2012). The host galaxies of radio- and IR-selected AGN have also been the subjects of investigation with Herschel (e.g., Hardcastle et al. 2012; Del Moro et al. 2013), though the connection between these AGNs and their host galaxies is less clear than in the case of X-ray-selected AGN. Additionally, large populations of starbursting galaxies without an active nucleus have been studied with Herschel, resulting in insights into the relationship between stellar mass and the in star-forming galaxies (e.g., Elbaz et al. 2011; Rodighiero et al. 2011; Hilton et al. 2012), specific star formation rates (s) and dust mass (Smith et al. 2012a), various indicators (Wuyts et al. 2011; Domínguez Sánchez et al. et al. 2012), and a variety of other properties (e.g., Casey 2012a; Casey et al. 2012b; Roseboom et al. 2012).

In this study we utilize wide-field Herschel/SPIRE observations taken of the CFHTLS-D1 field. The deep multiband optical/NIR imaging, extensive coverage of both radio and X-ray imaging, and exhaustive spectroscopic campaigns by the VVDS and ORELSE surveys on this field allow a unique view on both the properties of Herschel-selected galaxies and of the interplay between the star-forming properties of such galaxies and all flavors of AGN. In the first part of this paper we use all imaging and spectroscopic data on the field, along with the results of well-calibrated spectral energy distribution (SED) fitting, to broadly investigate the properties of galaxies detected by Herschel. Comparisons are made both between SPIRE-detected and equivalent samples of SPIRE-undetected galaxies and between various subsets of SPIRE-detected galaxies. In the second part of the paper, we use every method at our disposal (IR, radio, X-ray, and spectra) to select AGN and investigate the properties of starbursting galaxies who play host to such phenomena. The structure of the paper is as follows. In §2 the properties of the vast datasets available on the CFHTLS-D1 field are described, including the Herschel/SPIRE imaging, and, to varying degrees, the reduction of these data are also described. In §3 we describe various analyses which setup the framework of the sample and the subsequent interpretation of its properties, including the match of SPIRE sources, the selection of AGN, and the SED fitting process. In §4 the color, magnitude, stellar mass, and spectral properties of the full sample of 2000 SPIRE-matched galaxies is discussed. The contents of §5 is specifically devoted to those SPIRE-detected galaxies hosting AGN and to investigating both their prevalence and their properties relative to the full SPIRE population. Finally, §6 presents a summary of the results. Throughout this paper all magnitudes, including those in the IR, are presented in the AB system (Oke & Gunn 1983; Fukugita et al. 1996). We adopt a standard, pre-Planck, concordance CDM cosmology with = 70 km s, = 0.73, and = 0.27.

2 Observations and reduction

Over the past decade and a half, the Canada-France-Hawai’i Telescope Legacy Survey (CFHTLS)13-D1 field (also known as the VVDS-0226-04, the 02h XMM-LSS, or the Herschel-VVDS field) has been the subject of exhaustive photometric and spectroscopic campaigns, both by telescopes spanning the surface of the Earth and, more recently, those far above it. The observational advantages of this field are numerous. It is an equatorial field, allowing observations from telescopes located in both the southern and northern hemispheres. In addition, the field lies at a high (absolute) galactic latitude , resulting in low galactic extinction, and is reasonably far removed from the ecliptic, resulting in acceptable levels of zodiacal light. The field first became the subject of study in 1999 when it was targeted with the Canada-France-Hawai’i 12K wide-field imager (CFH12K; Cuillandre et al. 2000) to provide photometry (McCracken et al. 2003; Le Fèvre et al. 2004) for the initial phase of the VIMOS VLT Deep Survey (VVDS; Le Fèvre et al. 2005)14. Shortly thereafter, a subsection of the field was observed in multiple optical bands as part of the deep portion of the CFHTLS providing the basis for numerous future studies. Since the inception of the CFHTLS imaging campaign in this field, CFHTLS-D1 has been targeted by a large number imaging and spectroscopic surveys at wavelengths that span nearly the entirety of the electromagnetic spectrum. In this section we describe the current data available in the CFHTLS-D1 field. For those readers interested only in a brief overview of the data, Table 1 and Fig. 1 outline the basic properties and coverage, respectively, of each dataset used in the paper.

Telescope/Instrument Band Survey Limit15 16
XMM/EPIC CD (2-10 keV) XMM-LSS 2m 3.710 17
XMM/EPIC B (0.5-2 keV) XMM-LSS 1m 1.210 18
CFHT/MegaCam CFHTLS 0.37m 26.8
CFHT/MegaCam CFHTLS 0.49m 27.4
CFHT/MegaCam CFHTLS 0.65m 27.1
CFHT/MegaCam CFHTLS 0.77m 26.9
CFHT/MegaCam CFHTLS 0.91m 25.7
CFHT/WIRCam WIRDS 1.25m 24.2
CFHT/WIRCam WIRDS 1.63m 24.1
CFHT/WIRCam WIRDS 2.15m 24.0
Spitzer/IRAC SWIRE 3.6m 22.8 (3 Jy)
Spitzer/IRAC SWIRE 4.5m 22.1 (5 Jy)
Spitzer/IRAC SWIRE 5.8m 20.2 (29 Jy)
Spitzer/IRAC SWIRE 8.0m 20.1 (31 Jy)
Spitzer/MIPS SWIRE 24m 18.5 (141 Jy)
Herschel/SPIRE HerMES 250m 13.719 (12 mJy)
Herschel/SPIRE HerMES 350m 13.620 (13 mJy)
Herschel/SPIRE HerMES 500m 13.521 (13 mJy)
GMRT Band-4 VLA-VIRMOS 49 cm (610 MHz) 17.8 (275 Jy)
VLA/B-Array VLA-VIRMOS 21 cm (1.4 GHz) 19.0 (94 Jy)
Keck II/DEIMOS 1200 l mm ORELSE 0.79m 12322
VLT UT3/VIMOS LR-red VVDS Deep 0.65m 688923
VLT UT3/VIMOS LR-blue/LR-red VVDS Ultra-Deep 0.65m 69224
Table 1: CFTLS-D1 imaging and spectral data

2.1 Optical and near-infrared imaging

In this paper the deep CFHTLS imaging was chosen instead of the original VVDS CFH12K imaging due to the increased depth and number of bands. The deep CFHTLS imaging covers 1 deg, centered on =[02:26:00,-04:30:00], and consists of five bands, , with 5 point source completeness limits (i.e., ) of 26.8/27.4/27.1/26.1/25.7, respectively27. Magnitudes were taken from the penultimate release of the CFHTLS data (T0006, Goranova et al. 2009) and corrected for Galactic extinction and reduction artifacts using the method described in Ilbert et al. (2006). The pixel scale of the CFHTLS imaging is 0.186/pixel and the average seeing for the five bands ranges between 0.7-0.85 over the entire area of the imaging used for this paper, with progressively better seeing in progressively redder bands. For further details on properties of the CFHTLS-D1 imaging and the reduction process see the CFHTLS TERAPIX website28, Ilbert et al. (2006), and Bielby et al. (2012).

Imaging in the near infrared (NIR) of a subsection the CFHTLS-D1 field was taken from WIRCam (Puget et al. 2004) mounted on the prime focus of the CFHT. This imaging was taken as part of the WIRCam Deep Survey (WIRDS; Bielby et al. 2012), a survey designed provide deep NIR imaging for a large portion of all four CFHTLS fields. In the CFHTLS-D1 field, the one relevant to this paper, the WIRDS imaging covers roughly % of the area covered by the CFHTLS imaging, reaching 5 point source completeness limits of 24.2, 24.1, & 24.0 in the , , & bands, respectively. While the native pixel of WIRCam is 0.3/pixel, imaging was performed with half-pixel dithers (also known as micro-dithering due to the dither pattern resulting in moves of 9m at the scale of the detector) resulting in an effective pixel scale of 0.186/pixel for the WIRDS imaging. This pixel scale matches that of the CFHTLS imaging and, more importantly, is sufficient to Nyquist sample the 0.6 average seeing achieved by WIRDS at Mauna Kea in the NIR. For a full characterization of the WIRDS survey, as well as full details of the WIRDS imaging properties and reduction see Bielby et al. (2012).

Additional NIR imaging spanning a portion of both the CFHTLS and WIRDS coverage was drawn from the Spitzer Wide-Area InfraRed Extragalactic survey (SWIRE; Lonsdale et al. 2003). The SWIRE survey has imaged 49 deg, of which the CFHTLS-D1 field is a part, with Spitzer at 3.6/4.5/5.8/8.0 m from the Spitzer InfraRed Array Camera (IRAC; Fazio et al. 2004) and at 24m from the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004). All data were taken from the band-merged catalog of this field, released as part of SWIRE Data Release 2. While this release is deprecated, no changes were made to the CFHTLS-D1 SWIRE data (referred to as SWIRE-XMM-LSS by the SWIRE collaboration) in future data releases. The rough full-width half-maximum (FWHM) point spread functions (PSFs) of the IRAC images are 1.9, 2.0, 2.1, & 2.8 at 3.6, 4.5, 5.8, & 8.0 m, respectively, while the FWHM PSF of the MIPS image is . The data were reduced by the standard SWIRE pipeline (Surace et al. 2005)29. Magnitudes were derived using a combination of aperture-corrected fixed-aperture measurements for fainter sources () and variable diameter Kron apertures (Kron 1980; Bertin & Arnouts 1996) for brighter sources. The SWIRE imaging reaches 5 point source limiting magnitudes of 22.8, 22.1, 20.2, 20.1, & 18.5 (3, 5, 21, 29, & 141 Jy) in the 3.6, 4.5, 5.8, 8.0, & 24m bands, respectively30. For further details on the SWIRE data used in this study see Arnouts et al. (2007) and the footnote below. A schematic diagram showing the CFHTLS, WIRDS, & SWIRE coverage for the CFHTLS-D1 field is shown in Fig. 1. In this study we will be focusing exclusively on the 0.8 deg region where the available SPIRE, CFHTLS, and SWIRE observations overlap (i.e., the part of the blue CFHTLS coverage square in Fig. 1 that is to the right of intersecting SWIRE magenta line). This area is referred to in the remainder of the paper as the “overlap region”.

2.2 Spectroscopic data

While the deep ground-based broadband optical/NIR data described above are used to estimate redshifts of all detected objects solely from the photometry itself (hereafter photo-s, see §3.3), it is important to complement these photo-s with spectroscopic redshifts of a large subsample of these objects in order to check for biases, consistency, precision, and accuracy. While we did not, in this study, perform a dedicated spectroscopic survey of Herschel counterparts (as was done in, e.g., the study of Casey et al. 2012b, 2012c), the CFHTLS-D1 field is the subject of extensive spectroscopic campaigns from two different surveys, both of which are used in this study.

The first survey, and the survey from which a large majority of our spectroscopy are drawn, is the VVDS (see Le Fèvre et al. 2005, 2013a for details on the survey design and goals). The VVDS utilizes the Multi-Object Spectroscopy (MOS) mode of the VIsible MultiObject Spectrograph (VIMOS; Le Fèvre et al. 2003) mounted on the Nasmyth Focus of the 8.2-m VLT/UT3 at Cerro Paranal. In the CFHTLS-D1 field, the VVDS survey was broken into two distinct phases: the VVDS Deep survey (hereafter VVDS-Deep), a spectroscopic survey of % of the CFHTLS-D1 imaging coverage, magnitude limited to , and the VVDS Ultra-Deep survey (hereafter VVDS-Ultra-Deep), which covers a smaller portion of the field but targets galaxies to a magnitude limit of .

The two phases of the VVDS in this field used distinctly different spectroscopic setups and observing strategies. Briefly, VVDS-Deep observations were made with the LR-red grism, resulting in a spectral resolution of for a slit width of 1 and a typical wavelength coverage of Å. Integration times were 16000s per pointing. VVDS-Ultradeep observations were taken for both the LR-blue and LR-red grisms again with 1 slit widths, resulting in a spectral resolution of for both grisms. Due to the increased faintness of the target population the integration times of the VVDS-Ultra-Deep observations were significantly longer, typically 65000s per pointing per grism. The observation of both grisms resulted in a typical wavelength coverage of Å for the VVDS-Ultra-Deep spectra. For further details on the observation and reduction of the full VVDS dataset see Le Fèvre et al. (2005), Cassata et al. (2011), and Le Fèvre et al. (2013a).

In total, spectra of unique objects have been obtained in the CFHTLS-D1 field. In this study, only those VVDS spectra for which the redshift reliability has been determined to be in excess of 87% were used (see Le Fèvre et al. 2013a), which, in VVDS nomenclature, corresponds to flags 2, 3, & 431 (see Le Fèvre et al. 2005 and Ilbert et al. 2005 for more details on the VVDS flag system). Applying this reliability threshold results in 7352 and 722 spectra for the VVDS-Deep and VVDS-Ultra-Deep datasets, respectively, in the overlap region. Galaxies surrounding bright stars, for which the photometry was insufficient to perform accurate SED fitting (see §3.3), were removed from this sample, resulting in a sample of 6889 and 692 high-quality spectra in the overlap region for the VVDS-Deep and VVDS-Ultra-Deep samples, respectively.

Additional spectra were taken from the Observations of Redshift Evolution in Large Scale Environments (ORELSE; Lubin et al. 2009) survey. The ORELSE survey is an ongoing multi-wavelength campaign mapping out the environmental effects on galaxy evolution in the large scale structures surrounding 20 known clusters at moderate redshift (). One of these clusters, XLSS005 (), lies in a portion of the CFHTLS-D1 field. ORELSE spectra were obtained with the DEep Imaging Multi-Object Spectrograph (DEIMOS; Faber et al. 2003) on the 10-m Keck II telescope. Four slitmasks were observed of XLSS005 with an average integration time of s and an average seeing of . Observations were taken with the 1200 l mm grating tilted to a central wavelength of Å, yielding a typical wavelength coverage of Å. All targets were observed with slits, which results in a spectral resolution of . Data were reduced with a modified version of the Interactive Data Language (IDL) based spec2d package and processed through visual inspection using the IDL-based zspec tool (Davis et al. 2003; Newman et al. 2013). Only those spectra considered “high quality” (i.e., Q=-1,3,4) are used in this study (for an explanation of the quality codes see Gal et al. 2008), resulting in 317 high-quality spectroscopic redshifts (hereafter spec-s) of which 123 lie in the survey area and have sufficient photometric accuracy necessary for SED fitting. For further details on the reduction of ORELSE data see Lemaux et al. (2009).

Accounting for duplicate observations, in total the VVDS and ORELSE datasets combine for high-quality spectra of unique objects in the overlap region. Spectroscopic coverage of both phases of the VVDS as well as ORELSE is represented as shaded regions in Fig. 1.

2.3 Far-infrared imaging

Imaging of the CFHTLS-D1 field at 250, 350, & 500m was provided by the Spectral and Photometric Imaging REceiver (SPIRE; Griffin et al. 2010) aboard the Herschel space observatory (Pilbratt et al. 2010), taken as part of the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012). HerMES is a massive imaging campaign with Herschel designed to cover 340 deg of the sky to varying depths with both SPIRE and the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010). HerMES observations of the CFHTLS-D1 field (referred to as “VVDS” by the HerMES collaboration) were taken between the UTC 14th-15th of July 2010 (Obs. IDs 1342201438-1342201444) and designed as “level” (i.e., tier) 4 observations, where the level system runs from 1-7 and higher levels generally result in observation with increased area at the expense of depth. SPIRE observations were performed using Large Map mode, a mode which is described in §3.2.1 of the SPIRE Observer’s Manual32. The CFHTLS-D1 field was imaged with SPIRE at a relatively deep exposure time/area; the roughly 2 deg VVDS area that is considered to be within the good part of the VVDS mosaic (i.e., , the area in which the sampling is sufficiently deep relative to the mean sampling rate, see Oliver et al. 2012) was imaged with a total integration time of 10.4 hr (essentially equivalent depth as that of the Groth Strip and COSMOS HerMES, but shallower than the coverage in, e.g., the ECDFS and GOODS-N/S). The subsection of the 2 deg area imaged by SPIRE that intersects the CFHTLS-D1 field is inclined slightly () to the CFHTLS coverage, but, despite this, spans the entirety of the CFHTLS imaging (see Fig. 1). In addition, the entire deg 02h XMM-LSS field was imaged by both SPIRE and PACS to moderate depth (level 6) and released as part of Data Release 1 by the HerMES collaboration. These images were, however, too shallow to consider for the sample presented in this paper. For more details on HerMES and the Herschel observations relevant to this study see Oliver et al. (2012) and the SPIRE Observer’s Manual.

SPIRE observations of the CFHTLS-D1 field were obtained from the Herschel Data Archive (HDA). These observations consist of 14 scans which utilized the Herschel Large Map mode (described previously) and the nominal SPIRE scan rate, with an equal number of scans in quasi-orthogonal directions so as to enable the removal of low-frequency drifts. Level 1 data products processed by the Herschel Interactive Processing Environment were retrieved from the HDA and combined using the package Scanamorphos (Roussel 2013), which is an IDL-based set of routines specifically designed to process observations obtained from bolometer arrays. These routines, by using the single scan overlap regions, allow for the subtraction of the thermal and pixel drifts as well as the removal of low-frequency noise. Subsequently, Scanamorphos combines the different scans to produce a single mosaic at each available wavelength. The resulting mosaics, which span deg each, have a pixel scale of 4.5, 6.25, and 9 pix for the 250, 350, and 500 bands, respectively.

Far-Infrared (FIR) source catalogs in the three SPIRE bands were obtained with standard PSF fitting using the Image Reduction and Analysis Facility (IRAF; Tody 1993) daophot package (Stetson 1987). First, the 24-detected sources taken from the MIPS-SWIRE coverage of the CFHTLS-D1 field were used to create a list of priors for positions of possible FIR sources and a PSF fitting was performed at the positions of the priors for each of the three SPIRE mosaics. However, because of the relatively shallow depth of the SWIRE MIPS CFHTLS-D1 observations, many sources detected with Herschel lack counterparts at 24. We thus ran Sextractor on the residual image obtained after the first PSF fitting to search for any FIR sources that were not extracted with the 24 priors. This list of additional detections, comprising % of our total SPIRE catalog, was merged with the list of FIR sources detected from the initial 24 prior input catalog to create a full list of FIR detections. A final round of PSF fitting was then performed over each of the initial science images, this time using the final list of positions as position priors. With the exception of a few rather bright and nearby sources, the residual images obtained after this procedure did not reveal any obvious excess of signal related to faint sources remaining after the PSF fitting, nor negative signal that could be attributed to PSF over-subtraction. The bright objects where the PSF fitting did not perform well were removed from the input list of prior sources and their flux density measurements were treated separately with aperture photometry.

The Herschel data are calibrated in Jy/beam, so the flux density of each point source extracted by daophot was derived from the central pixel value of the scaled PSF which fit the object. The flux density uncertainties provided by daophot were verified against simulations that were also performed to estimate our source extraction completeness. In this simulation, modeled PSFs with flux densities ranging from 1 to 50 were randomly added to each mosaic and extracted with the real sources using the same PSF fitting technique as was used for building the Herschel source catalogs. The distribution of the differences between the input and measured flux densities for these artificially added point sources yielded 1 uncertainties of 4.0, 4.7 and 4.8 mJy in the 250, 350 and 500 bands, respectively. In a given bin of flux density, the completeness of the source extraction was inferred as the fraction of recovered sources for which the measured flux density agrees within 50% of that of the input PSF. This criterion led to 80% completeness limits of 18.5, 16 and 15 mJy in the 250, 350 and 500 bands, respectively33. These limits are slightly brighter than the 3 limits in each of the three Herschel bands (see below and Table 1), meaning that, with these cuts we are not performing a complete census of these sources, but are rather detecting 70% of such sources. The difference between these two numbers is largely driven by confusion, meaning that Herschel sources in areas of high SPIRE density will be underrepresented in our sample. Despite the high level of completeness for sources detected at in at least one of the three SPIRE bands, this confusion could induce some small level of bias to our sample if the SPIRE sources in areas of high real or projected density have appreciably different properties, either in the far infrared or in the other galaxy parameters used in this study. However, there is only one known cluster in the CFHTLS-D1 field (Valtchanov et al. 2004), which lies in an area of the sky that does not have a particularly high density of SPIRE sources. Thus, any sources missed due to this confusion will be missed due to projection effects, and as such, there is no a priori reason to believe such sources would not share the properties of the 70% of field galaxies recovered by imposing a significance threshold of .

The PSFs of the SPIRE 250, 350, & 500m band images are 18, 25, & 37, respectively. This is a factor of 10 more than our complementary images in shorter wavelengths, a subject that is dealt with extensively in §3.2. The SPIRE images reach measured 3 flux density limits of 12.0, 12.9, & 13.2 mJy at 3 in the 250, 350, & 500m bands, respectively, (or, alternatively, point source completeness limits of 18.3, 21.8, & 22.6 mJy, respectively). A total of 10431 SPIRE sources were extracted at all significances in the overlap region, of which 3843 were detected at in at least one of the three SPIRE bands. This latter definition will be adopted widely throughout the paper to define the sample detected “significantly” in SPIRE.

The above definition of significantly detected is more liberal than studies of Herschel-selected samples (e.g., Bétherman et al. 2012). As such, there is some concern that adopting such a limit would result in a higher number of spurious SPIRE sources. However, there are several reasons to believe this is not the case. Approximately half of the SPIRE sample is detected significantly in at least two SPIRE bands. Additionally, among the sample detected at 3 in at least one SPIRE band, the median S/N in the next most significant SPIRE band is . Roughly two-thirds of the sample detected significantly in only a single SPIRE band is detected at 5 in the MIPS-SWIRE imaging (see §3.2). Of the remaining objects, the median significance in the detection band is sufficiently high () that only a small number of objects in our full SPIRE sample are likely due to stochastic fluctuations (10 objects). The above arguments rely solely on the statistics of the MIR/FIR imaging. However, the sample we present here is not a Herschel-selected sample, but rather an optical/NIR sample that is matched to Herschel/SPIRE imaging (see §3.2). Thus, we do not require only a fluctuation in the SPIRE imaging, but a fluctuation at specific points in the sky. If we instead consider the probability of a spurious 3 Herschel/SPIRE source matching to an optical counterpart using either of the two matching methods detailed in §3.2, the number of spurious sources drops essentially to zero. As such, the effect of spuriously-detected SPIRE sources is ignored for the remainder of the paper.

2.4 Other ancillary data

Integral to our study of coeval starbursts and AGN in this paper were the X-ray and radio imaging provided by two different surveys. In this section we briefly describe the properties of this imaging.

X-ray observations of the CFHTLS-D1 field were taken from the X-ray Multi-mirror Mission space telescope (XMM-Newton; Jansen et al. 2001) as part of the XMM Medium Deep Survey (XMDS; Chiappetti et al. 2005). This survey, which lies at the periphery of the much larger XMM - Large Scale Structure survey (XMM-LSS; see Fig. 2 of Pierre et al. 2004 for rough survey geometries and details on survey motivation), utilized the MOS (Turner et al. 2001) and pn (Strder et al. 2001) cameras on the European Photon Imaging Camera (EPIC) aboard XMM to map out an area deg that is nearly coincident with the CFHTLS-D1 field. A majority of the survey area, and the entirety of the area used in this study, was imaged with a moderate integration time of 20-25 ks. This results in images with a point source flux limits of 3.710 and 1.210 ergs s cm in the soft (0.5-2 keV) and hard (2-10 keV) bands, respectively. This limit is equivalent to ergs s at , corresponding to a moderately powerful Seyfert.

X-ray fluxes were drawn from the catalog of Chiappetti et al. (2005), which provides fluxes and errors in the B (soft) and CD (hard) bands for 286 X-ray objects. The sources in this catalog are limited to those objects that satisfied a significance cut in one of the EPIC bands. Additionally, the area coverage of the XMDS sources are limited to the XMDS overlap with the original VVDS imaging, resulting in the pattern shown in Fig. 1. Duplicate entries from this catalog were removed and sources were re-matched to our own band-merged optical/NIR photometric catalog (see §3.1 and Appendix A). For further details on the observations and reduction of the XMM data we refer to Chiappetti et al. (2005).

Radio observations of the CFHTLS-D1 field were taken at 1.4 GHz and 610 MHz from the Very Large Array (VLA, B-Array) and the Giant Metrewave Radio Telescope (GMRT), respectively, as part of the VLT-VIRMOS survey (Bondi et al. 2003). Details of the observations and reduction of the VLA data are given in Ciliegi et al. (2005) and those of the GMRT data are given in Bondi et al. (2007). An updated band-merged version of the optical-radio catalog provided in Ciliegi et al. (2005) was generated. This updated catalog contains fluxes and errors in both 1.4 GHz and 610 MHz for all 1054 objects that meet a significance threshold of in the VLA imaging. A maximum liklihood matching procedure was used to match the radio objects with optical counterparts (see Ciliegi et al. 2005), which we have in turn matched to our updated band-merged optical/NIR photometric catalog (see §3.1 and Appendix A). The observations reach point source completeness limits of 275 and 94 Jy beam in the GMRT and VLA images, respectively. The PSF in both images is roughly . For further details on the observation and processing of the radio data see Ciliegi et al. (2005) and Bondi et al. (2007).

3 Analysis

In this study we combine datasets from across the full electromagnetic spectrum, combining data from telescopes that have PSFs that span nearly two orders of magnitude, ranging from for our ground-based NIR imaging to in the reddest Herschel band. In this section we discuss the process of matching and otherwise combining these data into a coherent, self-consistent dataset.

3.1 Near-infrared, x-ray, and radio imaging

The near-infrared, X-ray, and Radio imaging catalogs used in this study were drawn from a variety of other studies and adapted to our own purposes. A combination of methods was used to match these catalogs to optical sources in the CFHTLS images and to correct observed flux densities (or power densities in the case of radio) to rest-frame values for each source. The literature from which these catalogs were drawn as well as the various methods used to refine and adapt these datasets for use in this study is discussed in Appendix A.

3.2 Herschel counterparts

Two distinct methods were used to determine optical/NIR counterparts to those objects detected by Herschel. While these methods take dramatically different forms and are used to select significantly different populations, we will show later that the two methods are complementary and consistent. Our initial method of matching utilized SWIRE/MIPS 24m data. Because there is a strong physical motivation for expecting a positive correlation between 24m flux and flux in the FIR, those objects detected in 24m were used as a “parent” population to extract SPIRE sources at the position of the MIPS detection. This step was done during the extraction process described in section §2.3. Thus, any objects extracted through this methodology, by definition, had a counterpart in the MIPS data. Because the MIPS data was matched to the SWIRE/IRAC data, which was in turn matched to our optical CFHTLS photometric catalog, all Herschel objects extracted in this manner, except those MIPS detections that had no IRAC counterpart (comprising only % of the MIPS sample) for which no optical match was made or objects for which a sufficiently robust match between the IRAC and CFHTLS optical catalog could not be made, necessarily had a optical counterpart. Of the 2146 Herschel detections in the overlap region that were detected at in at least one of the three SPIRE bands and matched to a counterpart, 1682 (78.4%) were determined through this method.

Because a second set of extractions was performed on the residual Herschel mosaic subsequent to extracting all sources with 24m parents, a large population of SPIRE detections were without optical counterparts after the extraction process. Due to the shallowness of the SWIRE imaging in 24m relative to other similar fields (e.g., COSMOS), this 24m-“orphaned” Herschel population comprised a large fraction of the SPIRE sample over the full SPIRE mosaic (% of significantly detected objects), though some are likely galaxies which have low ratios (Magdis et al. 2011). In order to prevent severely limiting our sample size, a second round of matching was performed on the orphaned SPIRE sample. While there is an obvious reason to expect galaxies that are intrinsically brighter at m to be, on average, intrinsically brighter in the FIR, there exists no other sufficiently deep dataset in the CFHTLS-D1 field where a similar correlation is expected34. In the absence of an astrophysical motivation, we chose to adopt a statistical approach relying on the SWIRE/IRAC photometry and a well-established methodology of matching sources between images of significantly different PSFs (Richter 1975; Sutherland & Sanders 1992). This method, described in detail in Kocevski et al. (2009a) and Rumbaugh et al. (2012), uses a likelihood radio statistic of the form


where the indices and correspond to the IRAC and SPIRE object, respectively, is the positional error of the SPIRE object (we assumed that the error in the NIR position is negligible), and is defined as , the square root of the inverse of the number density of objects in the IRAC image brighter than the IRAC source, following the methodology of Rumbaugh et al. (2012). Positional errors for SPIRE objects depended on flux density, a dependence determined from simulations to take the form


where is the flux density of the SPIRE source in mJy and is defined to be in units of . The behavior and magnitude resulting from the above formula as a function of flux density are similar to those found in similarly deep fields by the HerMES collaboration (Smith et al. 2012b). Monte Carlo simulations were performed in an identical manner to the one described in Rumbaugh et al. (2012) and identical probability cuts were applied. For those SPIRE objects which were matched to more than one IRAC counterpart, the counterpart with the higher probability was selected as the genuine match if the primary match met the conditions given in Rumbaugh et al. (2012). All SPIRE objects with multiple IRAC counterparts that did not meet these conditions were removed from our sample. In total, 464 24m-undetected counterparts of significantly detected SPIRE sources were determined using this method, increasing our sample sample size by nearly a factor of 1.3 over the sample containing only those SPIRE objects with 24m parents. While this is a sizeable increase over the original sample, the fraction of orphaned SPIRE objects matched to an optical counterpart was only %, much lower than that fraction for the SPIRE-detected objects with 24m parents. This low fraction is as a result of two effects: the enormous PSF of SPIRE, which makes it incredibly difficult to properly identify which object (galaxy or star) the SPIRE source belongs to, and the conservativeness of the probabilistic matching scheme. For the latter point, the parameters of the scheme and the associated cutoffs were set sufficiently high to ensure matches are almost certainly real at the expense of lower likelihood matches (see discussion in Rumbaugh et al. 2012).

These two methods differ wildly in their approach, which causes some concern for induced differential bias in our full matched SPIRE sample. Consistency between the two methods was checked by blindly running the probabilistic matching technique on the full SPIRE sample and applying the same probability threshold as was applied for the orphaned matching. This catalog was then compared to the catalog of SPIRE objects with m parents. Remarkable consistency was found between the two methods, as % of the objects determined by using m priors that had IRAC counterparts were recovered through the full probabilistic process. While this is perhaps to be expected given that the m position for a given galaxy with an IRAC counterpart was, by virtue of the method used to generate the SWIRE catalog, forced to be at the same position as centroid of the object in the IRAC 3.6m band, the observed consistency of the two methodologies lends credence both to the choice of the weighting scheme used in the probabilistic method and the adopted reliability thresholds. In addition, since we do not use any galaxies in this paper that do not have an optical counterpart, i.e., galaxies that have a 24m and a SPIRE detection, for which we cannot test the consistency of the two methods, this test shows that the sample used in this paper would have been largely unaffected if we had instead chosen to solely adopt the probabilistic method for determining the optical counterparts to SPIRE detected sources from the onset. Thus, for the remainder of the paper no distinction is made between orphaned SPIRE objects and those with 24m parents. Between the two methods, 2146 of the 3843 SPIRE sources detected at in at least one of the three SPIRE bands were matched to an optical counterpart in the overlap region. Summarized in Table 2 are the number of objects determined at each stage in the process using the two methods.

A final note is necessary. Very recent results with the Atacama Large Millimeter/submillimeter Array (ALMA) have found that a large fraction (%) of sub-mm sources are resolved into multiple components when the resolution of the instrument used is increased dramatically (Hodge et al. 2013). The resolution of SPIRE, even in the shortest wavelength bandpass, is prohibitively large such that we cannot determine the pervasiveness of this phenomenon for the current sample, but given that the selection methods used in Hodge et al. (2013) are similar to our own, it is likely that a non-negligible fraction of the SPIRE-detected sources in our sample are in fact the composite emission from multiple optical/NIR sources. Though such an effect confuses the interpretation of our sample, it is, by in large, only necessary for our study that this effect induces no differential bias, as all the comparisons made in this study are either internal or to other sub-mm-selected surveys. Over the dynamic range probed by Hodge et al. (2013), the number of sub-mm sources found to have multiple counterparts, sources which were originally obseved with a similar PSF to the SPIRE 250m channel, was broadly representative in flux density of the sample as a whole. Since many of the results in this paper are predicated upon comparisons either between the properties of galaxies selected by their brightness in SPIRE or comparisons of the SPIRE brightnesses of different populations, the observed trend in Hodge et al. (2013) is encouraging. However, because ALMA observations have just begun, the number of sources studied in Hodge et al. (2013) was relatively small, and future work will certainly be necessary to fully characterize this effect as it relates to Herschel studies. Without further ability to determine the true number of SPIRE sources in our sample that result from two or more optical/NIR sources, we take refuge in trend observed from preliminary ALMA results and ignore the effects of multiple SPIRE counterparts for the remainder of the paper.

Sample 24m Prior 24m Orphaned35
Total SPIRE Sources36 2366 1477
SPIRE Counterparts37 1682 464
Final SPIRE Sample38 1367 386
Table 2: Summary of Herschel/SPIRE detections in the overlap region

3.3 Spectral energy distribution fitting and redshift distribution

Many of the results in this study are not directly reliant on the observed properties of the SPIRE sample, but rather on secondary physical quantities derived from those properties. In order to link observed properties to these secondary physical quantities we performed spectral energy distribution (SED) fitting on the entire SPIRE sample using the code Le Phare40 (Arnouts et al. 1999; Ilbert et al. 2006) in two steps. In this section the initial step of this process, used to derive all quantities for our sample except total infrared luminosities (TIRs) and s, is described and discussed.

The initial SED-fitting was performed solely on our ground-based optical and NIR photometry (i.e., CFHTLS+WIRDS) using the trained methodology similar to that described in Ilbert et al. (2006). Spitzer IRAC and MIPS magnitudes were not used at this stage as those data were significantly shallower than the CFHTLS/WIRDS imaging. The method used for deriving photo-s and almost all other physical parameters associated with our sample is identical to the one described extensively in Bielby et al. (2012) for the CFHTLS-D1 field. As such, we describe it only briefly here.

For the photo- fitting process, a combination of templates from Polletta et al. (2007) and Ilbert et al. (2009) were used. The final photo-s were drawn from the full probability distribution functions (PDFs) using a simple median as in Ilbert et al. (2010). For the CFHTLS-D1 field the reliability and precision of the derived photo-s were checked against the full VVDS spectroscopic sample (flags 3 and 4 only). Since the photo-s in this field were, by in large, trained on the VVDS spectral data, it is not suprising that the two sets of redshifts show excellent agreement. The normalized absolute median deviation (NMAD; Hoaglin et al. 1983), defined in Ilbert et al. (2013), is (or, using a more traditional metric, a median ) for galaxies with both a reliable spec- and photo- (see Bielby et al. 2012). The catastrophic failure rate, defined as the percentage of galaxies with (see Ilbert et al. 2013), is also quite low, roughly accounting for only 2.2% of the sample. Stellar masses and other physical parameters were derived using a stellar population synthesis models from Bruzual & Charlot (2003) and the methodology given in Ilbert et al. (2010) and Bielby et al. (2012). A rough range of parameters used in this fitting, e.g., , , metallicity, is given in Table 1 of Ilbert et al. (2010), with exceptions given in the text. As in Ilbert et al. (2010), we adopt a universal Chabrier initial mass function (IMF; Chabrier 2003) along with a Calzetti et al. (2000) reddening law.

It is extremely important for this study to note that, unlike what is observed in Casey et al. (2012c; hereafter C12) and discussed extensively therein, the disagreement between the spec-s and photo-s of our SPIRE-detected sample is not vastly larger than that of the full photometric sample. Applying the same metrics to our SPIRE-detected sample with reliable photo-s and spec-s, we find only a marginal increase in the NMAD ( or ) and the catastrophic failure rate (%) of the SPIRE-selected sample to the magnitude limit imposed for reliable photo-s (, see Ilbert et al. 2010). The former metric is nearly an order of magnitude lower than that of C12. While it is tempting to speculate on the source of this discrepancy, the photo-s derived in C12 are drawn from a variety of different methodologies using fields that have wildly different photometric coverage and depths. This confuses our ability to discern the source of imprecision in their methodology and makes even a qualitative comparison between the two datasets a daunting task. Therefore, we simply state here what is important for our sample: that the precision or catastrophic failure rate of the photo-s of our SPIRE-selected sample does not differ appreciably from the full photometric sample. This is integral to our study, as it not only gives us confidence in the derived photo-s of those SPIRE-selected galaxies that do not have spec-s, but allows us to make internal comparisons between SPIRE-selected and SPIRE-undetected galaxies without fear of induced bias.

As is discussed extensively in the literature (e.g., Vanzella et al. 2001; Coe et al. 2006), self-consistent magnitude measurements from band to band are extremely important to the SED fitting process. While the rewards of analyzing secondary quantities derived from SED fitting can be significant – stellar mass and rest-frame colors allow one to break degeneracies not possible when considering only observed-frame colors (e.g., Bell et al. 2004; Bundy et al. 2006; Lemaux et al. 2012; Bradač et al. 2013; Hernán-Caballero et al. 2013; Ryan et al. 2013) – there are also many potential hazards and subtleties associated with this type of analysis. However, the primary comparisons made in this study are internal. Thus, biases that introduce themselves during the process of fitting for these secondary quantities which affect the entire population are of less concern than those that affect populations differentially. To mitigate the latter, we treat the SPIRE-selected sample in an identical manner during the SED-fitting process as the full photometric sample.


[angle=0,width=1.0]FINAL.matching.zhist.3sigma.eps \includegraphics[angle=0,width=1.0]FINAL.matching.zhist.3sigma.specfaintvsbright.eps

Figure 2: Left:Photometric (open magenta histogram) and spectroscopic (dashed blue histogram) redshift distribution of the 1783 Herschel/SPIRE-matched sources in our sample. Only those objects which have reliable redshifts (see text) and are detected at in at least one of the three SPIRE bands are plotted. The various peaks in the photometric redshift distribution are not real but are rather an artifact of the method used to calculate the photometric redshifts as can be seen in the photometric redshift distribution of the full CFHTLS-D1 sample () plotted as the gray filled histogram. Right: Spectroscopic redshift histogram of two samples of SPIRE-detected galaxies at different optical brightness (orange and green hashed histograms). Also plotted is the spectroscopic redshift distribution of the full VVDS sample in this field for the same optical magnitude bins (black hashed and gray filled histograms). Despite median magnitudes that differ by in the band, the redshift distributions of the two SPIRE-detected samples appear largely similar, with an excess of galaxies at higher redshift in the (optically) fainter sample that is dwarfed by the excess observed in the full optical sample.

For the initial SED-fitting process, a process which resulted in our derived photo-s, rest-frame colors, and stellar mass measurements, aperture magnitudes of 2 were used for all input photometric measurements. While this choice can result in increased global bias relative to measurements using MAG_AUTO (Kron 1980; Bertin & Arnouts 1996) especially for lower redshift galaxies, the PSF across the 8-band CFHTLS/WIRDS imaging was essentially constant, meaning that there is no induced bias in our fitting process from differential loss of light at a given redshift. In addition, at the conclusion of the fitting process, quantities which are heavily dependent on aperture corrections, such as stellar mass and absolute magnitudes, were scaled by the difference41 between the aperture and MAG_AUTO magnitudes using the appropriate band (e.g., the median of the bands for stellar mass). Using this methodology allows for the determination of the best-fit model through the SED fitting process with a set of self-consistent magnitude measurements that are minimally effected by blending, while still allowing for a recovery of quantities which require a full account of the light of a galaxy. This statement is given credence by the photo-s derived using aperture magnitudes for the full CFHTLS-D1, which show considerably better agreement with the full VVDS spectral training set than did those derived using, e.g., MAG_AUTO.

A subset of several hundred SED fits, both for SPIRE-detected objects and those undetected in SPIRE, were examined by eye. While this was by no means a definitive test of bias, it was comforting that we found no major differences between the best-fit templates in the two subsamples relative to their photometric data points. This is reflected in similarity of the distribution of the SPIRE-selected sample relative to the full photometric sample. Both distributions are to a close approximation log normal, characterized by a mean of of 0.07 and -0.24 for the SPIRE-detected and full CFHTLS samples, respectively. Both values translate to statistically acceptable fits to the data and, though the SPIRE-detected objects have, on average, higher values, the difference between the mean of the two samples deviates by . With no a priori knowledge of the physical characteristics of the SPIRE-selected galaxies, we cannot go further in testing for differential bias induced on the physical parameters during the SED-fitting process. However, the consistency of the analysis used and the similarity of the primary statistic used to determine the quality of the SED-fitting process strongly suggests that no differential bias was induced during this process. Of the 2146 galaxies in the overlap region matched to a counterpart and detected at in at least one of the three SPIRE bands, 1800 had sufficiently precise photometry to perform this portion of the SED fitting. The spectral and photometric redshift distribution of these galaxies is shown in Fig. 2.

There are a few interesting notes to be made about the redshift distribution of SPIRE-matched galaxies. The first regards the photometric redshift distribution, which is not a smooth function, but rather varies appears with various peaks and troughs. While these peaks and troughs are mimicked to some degree by the spectroscopic redshift distribution, they are not necessarily a reflection of the underlying population. The method and statistical analysis used to calculate photometric redshifts, a method which is widely used, tend to siphon galaxies that are close in redshift into relatively narrow photometric redshift bins. The effect can also be seen in full CFHTLS-D1 photo- sample for galaxies brighter than plotted as the gray shaded histogram in the background of Fig. 2. The cause of this effect has to do with the narrowness of the various continuum breaks which are used as the primary discriminators of photometric redshift and the wideness of the broadband filters used to constrain the photometric redshift. This effect is mitigated to some extent when medium- or narrow-band filters are observed in addition to the broadband filters (see, e.g., Ilbert et al. 2009), a luxury that is not available for the CFHTLS-D1 field. However, the mean error on the photometric redshift of a given galaxy remains small (see §3.3) and, further, this effect is expiated by the use of large redshift bins for our analysis.

The second note of interest regards the spectroscopic redshift distribution, which appears qualitatively similar to the redshift distribution found in C12. At first glance this is perhaps not surprising; both samples are probing the redshift distribution of SPIRE-detected sources, albeit in different fields. However, the nature of the two surveys, and specifically the spectroscopy of the two surveys, is considerably different such that this similarity did not necessarily need to exist. As mentioned in §2.2, the study of C12 presented a systematic spectroscopic campaign of SPIRE-selected sources. In our paper, the spectroscopy was performed blindly with respect to SPIRE, in that no knowledge of the SPIRE-detected sources were used to target objects. While there are some additional criteria imposed on ORELSE spectroscopic targets, the vast majority of the spectroscopy in this paper is magnitude limited in the band, as opposed to the study of C12 in which galaxies were magnitude limited in the 250m band42. The similarity of the two spectroscopic redshift distributions then has two main consequences. The first consequence is comforting for this study. The spectroscopy of the CFHTLS-D1 field and the method of matching between optical and SPIRE counterparts reveals no bias in redshift with respect to the underlying population presented in C12, or only insomuch as their spectroscopic sample was biased with respect to the underlying population revealed by their photometric redshifts (though see the discussion earlier in the section regarding these measurements). The second consequence is of more interest astrophysically. While the (observed frame) wavelength coverage of the two samples is roughly similar, the spectroscopic data in the CFHTLS-D1 field probe much deeper than those of C12. The number of spectroscopically confirmed objects in C12 begins to drop dramatically at , whereas such objects are routinely observed in both phases of VVDS (see Le Fèvre et al. 2013a) and ORELSE (see Lubin et al. 2009). The broad similarity between the two distributions strongly suggest that the redshift distribution of SPIRE-detected sources does not vary rapidly with decreasing optical brightness. Plotted in the right panel of Fig. 2 are the spectroscopic redshift distributions of SPIRE counterparts in the CFHTLS-D1 field separated in two bins of optical brightness. The spectroscopic redshift distribution of galaxies with 2122.5 appears remarkably similar to galaxies with given that the two samples differ in their median magnitude by magnitudes. Only a slight increase in the number of galaxies at higher redshifts is observed in the histogram of the fainter sample relative to the brighter sample, a trend not shared by galaxies in general (see, e.g., Le Fèvre et al. 2013a; Newman et al. 2013). This lack of evolution to higher redshifts is further quantified by inspection of the spectroscopic redshift distribution of the full VVDS sample in the CFHTLS-D1 plotted for the same two optical magnitude bins in the right panel of Fig. 2. While the redshift distribution of SPIRE-detected counterparts is observed to peak at higher redshifts relative to the full VVDS sample for both the magnitude bins, the number of spectroscopically confirmed galaxies at increases only marginally, a factor of from the optically bright to the optically faint SPIRE sample. For the full VVDS sample, a sample which serves (largely) as the parent sample from which the SPIRE-counterparts are drawn and, as such, has the same observational biases, a factor of increase is observed in the number of galaxies over the same change in optical magnitudes (see also Le Fèvre et al. 2013b). In other words, it appears that large statistical samples of higher redshift () SPIRE-detected galaxies cannot simply be obtained by spectroscopically targeting galaxies at fainter optical/UV magnitudes. This phenomenon is perhaps due to the complex relationship that SPIRE-detected galaxies have with their optical/UV magnitudes, as dust can significantly modulate the latter at a given redshift, an effect which we will show later to be related, in turn, to the infrared luminosity of a given SPIRE-detected galaxy.

3.4 Total infrared luminosities and s

The second stage of the SED fitting process involved the full photometric dataset available for the CFHTLS-D1 field and was performed only on those galaxies detected with SPIRE. The goal for this stage was to derive total infrared (TIR) luminosities (), which allows us to calculate infrared-derived s for the SPIRE-selected sample. In principle, we would be better served by calculating the true global , a quantity which combines the ultraviolet (UV) light produced by young stars unobscured by dust and the component which is re-radiated at IR wavelengths by interaction with dust particles. It has been shown, however, that even in broad surveys of galaxies selected by (observed-frame) optical wavelengths, calculations of s from UV indicators lead to deficiencies of a factor of relative to the true, dust-corrected (e.g., Cucciati et al. 2012). As the SPIRE sample represents, by construction, the extreme of such galaxies in terms of their dust content, these values represent lower limits for the sample presented here. Empirical evidence of this using a combination of SPIRE imaging and spectral properties is given later in this paper (see §4.4.1). In addition, the SEDs and s of FIR-selected samples of galaxies appear dominated by the FIR portion of their spectrum (e.g., Elbaz et al. 2011; Rodighiero et al. 2010b, 2011; Smith et al. 2012a; Burgarella et al. 2013; Heinis et al. 2013). Also, as we will show later (see §4.4.1), the s of the SPIRE-detected galaxies in our sample determined from TIR luminosities are, on average, vastly higher (i.e., a factor of ) than those derived from the 3727Å [OII] UV emission feature. Thus, the contribution of UV light to the is ignored and for the remainder of the paper the IR-derived of the SPIRE-selected sample is equated with the global .


[angle=0,width=1.0]FINAL.SPIRE.02h.matched.lumhist.3sigmacut.eps \includegraphics[angle=0,width=1.0]FINAL.SPIRE.02h.matched.lumvsz.3sigmacut.eps

Figure 3: Left: Full total infrared (TIR) luminosity distribution for the 1753 galaxies in the final Herschel sample. Galaxies are binned as a function of (photometric or spectroscopic) redshift using the nominal redshift bins that are used thoughout most of the analysis in this study. Though the dynamic range of SPIRE observations is small, the large redshift baseline of this study allows us to study galaxies that span a factor of 10000 in TIR luminosity. Right: Similar to the left panel except the dynamic range of the Herschel sample is shown at each redshift. Herschel-detected galaxies are color coded in an identical manner to that of the left panel: galaxies with are represented by cyan Xs, galaxies with are represented by filled blue circles, galaxies with are represented by dark green open triangles, and galaxies with are represented by open red diamonds. The TIR limits used to contruct a volume-limited sample for each redshift bin are denoted by the dashed horizontal lines.

An additional complication to the calculation of the seemingly arises when discussing those SPIRE-detected galaxies that host a variety of active nuclei. Traditionally, separating emission coming from star-formation processes from that originating from AGNs has been extremely challenging. Even with high-quality spectroscopic information, which contains traditionally robust indicators, one has to rely on modeling (e.g., Kewley et al. 2001) or statistical methods (e.g., Kauffmann et al. 2003; Vogt et al. 2014) to determine the fractional contribution of each process to the strength of a given recombination or forbidden emission line. One of the greatest virtues of Herschel/SPIRE observations, however, is precisely the absence of this ambiguity. As has been shown for galaxies hosting X-ray-bright AGN (Mullaney et al. 2012; Johnson et al. 2013), spectroscopically confirmed samples of type-1 and type-2 AGN (Hatziminaoglou et al. 2010), those galaxies with IR-selected Compton-thick AGN (Sajina et al. 2012), and a variety of other types of AGN (Rowan-Robinson 2000) the SEDs of galaxies at the rest-frame wavelengths probed by the SPIRE observations, even for the highest redshifts considered in this paper (), are dominated by re-emitted IR photons originating from the cold dust component commonly associated with circumnuclear bursts in IR-bright galaxies rather than the warm dust component traditionally associated with the AGN torus. We can thus safely ignore the contributions of these types of AGN in the SPIRE bands.

The remaining population to be considered is then those galaxies hosting radio AGN. If we consider the simplest possible scenario, i.e., that the SED generated purely by the mechanism powering the radio AGN is governed by the average spectral slope observed for galaxies hosting a radio AGN (), the power density at THz frequencies (the frequency of the SPIRE bands) is 1% of that at 1.4 GHz. This translates to % of the TIR luminosities for all but the least luminous SPIRE sources in our full sample and is essentially negligible at redshifts . However, as is briefly discussed in the next section, radio emission originating from an AGN can exhibit a wide variety of spectral slopes, ranging from an extremely steep slope () to an inverted spectrum () in the case of optically thick synchrotron emission. In cases of the latter, it is possible that a large fraction of the SPIRE luminosity of the radio AGN hosts in our sample is dominated by self-absorbed synchrotron emission (see Gao et al. 2013 for a explanation of this phenomenon and Corbel et al. 2013 and López-Caniego et al. 2013 for astrophysical manifestations of it). On a global scale, however, this phenomenon seems to be an exception; the source density of such objects is extremely low (e.g., López-Caniego et al. 2013; Mocanu et al. 2013). In our sample, the median FIR colors (- and -) of the SPIRE-detected radio AGN hosts in our sample, or indeed the hosts of any type of AGN, are statistically indistinguishable from those hosting dormant nuclei. This is not likely to be the case if the mechanism driving the FIR emission of AGN hosts was disparate from that of the SPIRE population with dormant nuclei (though see discussion in Hill & Shanks 2012 and references therein for a possible alternative viewpoint). Another line of evidence comes from a more complex treatment of local ULIRGS involving multi-component SED fitting by Vega et al. (2008). In that study it was found that the bandpass spanning the rest-frame wavelengths 40-500 was the least sensitive to the presence of a radio AGN, such that only % of their sample needed AGN contributions to the TIR luminosity in excess of 10% to properly reproduce the SEDs of their sample. These wavelengths are essentially identical to the rest-frame wavelengths probed by SPIRE in this study. While it is unclear whether the properties of these lower redshift ULIRG radio AGN hosts are directly comparable to those in the redshift range of the galaxies presented in this study, the consistent picture painted by the various lines of evidence explored here allows us to confidently ignore the presence of any radio AGN, or indeed the presence of any other AGN, when deriving the of our sample from the SPIRE bands.

Because we now incorporate the SWIRE and SPIRE photometry in this stage of the SED fitting process, we move from using aperture magnitudes (as were used in the previous section) to MAG_AUTO for our CFHTLS/WIRDS photometry. At this stage of the fitting, we are attempting only to determine the TIR of our sample, the other physical parameters being set by the previous SED-fitting process. In order to determine this quantity, the inclusion of the SWIRE and SPIRE photometry is of paramount importance. This photometry is, however, aperture corrected, and thus it is necessary to use CFHTLS/WIRDS magnitudes which are corrected in a similar manner. Though in practice we find that the choice of MAG_AUTO over aperture magnitudes has only a small global impact on the derived (i.e., 5% mean offset) as the s are mostly dependent on the observed SWIRE/SPIRE magnitudes, we prefer to adopt the s derived from magnitudes which are self-consistent across the full range of wavelengths used by the SED-fitting process. The SED fitting was performed in the following manner. The photo- derived in the previous section or, where available, the spec-, was used as a prior and simultaneous fitting of the rest-frame UV/optical and IR portions of the SED of each SPIRE-selected galaxy was performed in a method similar to the one described in the previous section. The cutoff between the two portions of the SED was set to 7m in the rest-frame. This cutoff was chosen to place most of the polycyclic aromatic hydrocarbon (PAH) features observed in the typical spectrum of dusty starburst galaxies (see, e.g., Chary & Elbaz 2001) redward of the cutoff while minimizing contamination by stellar continua. The optical templates used were identical to those used for fitting the full photometric sample.

To these optical templates, added at this stage of the SED-fitting process were 64 different templates provided by Dale & Helou (2002) to simultaneously constrain the IR part of the spectrum. We chose these templates because the models incorporate a variety of different dust components at different temperatures for each model and because the templates are calibrated at SPIRE wavelengths using the Submillimetre Common-User Bolometer Array (SCUBA; Holland et al. 1999). These TIR luminosities were translated into s using the conversion of Kennicutt (1998), which is based on the models of Leitherer & Heckman (1995) and assumes a Salpeter (1955) initial mass faction (IMF)


where is the solar luminosity and is equal to ergs s. Employing instead a Chabrier (2003) IMF results in a % reduction in the TIR-derived (e.g., Magnelli et al. 2013), though the absolute value of the is not as important for this study as relative bias. While relative bias may enter through a varying IMF (e.g., Swinbank et al. 2008), without an ability to constrain the IMF for individual galaxies we adopt a Salpeter IMF for all galaxies to consistently compare with other surveys. Other suites of IR templates were experimented with to calculate TIR luminosities from our MIPS/SPIRE photometry. However, we found that the templates of Dale & Helou (2002) yielded the best results both in visual inspection of the fits and in comparing the template-derived TIR luminosities with luminosities derived from directly translating SPIRE flux densities in individual bands. In addition, the models of Dale & Helou (2002) showed the most internal consistency when comparing the 24m-derived TIR luminosity with the TIR luminosity derived from MIPS+SPIRE for lower redshift () galaxies.

It has been suggested that, given the range of model parameter space and the small number of photometric constraints, a functional fit to the FIR photometry is perhaps a more robust methodology of deriving physical parameters from FIR observations of galaxies (for a detailed discussion see Casey 2012a). This may be especially true for complex quantities like the dust temperature or dust mass, where standard models from the literature may fail to capture the subtleties of the physical conditions of sub-mm bright galaxies. Unfortunately, galaxies in our sample do not always have the number of photometric datapoints in the IR necessary to perform a functional fit to the data. However, we are concerned here only with of the SPIRE-detected galaxies, a quantity which is typically measured to within 10% (at most) between the two methodologies (Casey 2012a). Because we have no way of verifying the differences between the two methods for each of the galaxies in the full SPIRE sample, and, more importantly, because the main results presented in this study are insensitive to this level of precision on the overall normalization of , we chose to adopt the SED-fitting process for the measurements of the full SPIRE sample.

At the conclusion of this fitting process, a final SPIRE sample of 1753 of the 1783 galaxies selected in the previous section were selected by imposing that each galaxy, be detected at in at least one of the three SPIRE bands, have a well-defined redshift (a high-quality spectral redshift or a photo- with ) between , and have sufficiently accurate photometry to perform both the initial SED-fitting (from which we derive stellar masses, etc.) and to perform the fitting discussed here. The median redshift of the finaal SPIRE sample is , identical to that of the 1783 galaxies selected in the previous section. In Fig. 3 we show the SED-derived TIR luminosities as a function of redshift for our final SPIRE sample.

3.5 Active galactic nuclei selection & classification

In this section we briefly discuss the various AGN populations used for this study. The general philosophy adopted is to select AGN samples where purity is maximized without regard to completeness. This choice was made so that we need to make no allowance for ambiguity when discussing these populations relative to full SPIRE-selected sample. While it is true that we will miss, by virtue of these selections, those AGN which are less powerful (i.e., primarily lower radio power density populations), hybrid galaxies (e.g., Kewley et al. 2001; Davies et al. 2014) where the existence of the AGN cannot be proven definitively, and other exotic AGN or AGN-like populations which are discussed later in this section, it is enough for this study that the galaxies we do select here are definitively the hosts of powerful AGN. Four different selection techniques were used in this study to search for the presence of AGNs amongst the SPIRE sample. These were radio, NIR, X-ray, and spectral selection.

For readers interested only in the basic properties of each AGN population selected by the various methods used in this paper, these are summarized in Table 3. There are many subtleties associated with the selection of each type of AGN, subtleties which induce a variety of consequences during the selection process. As a large portion of this study is predicated upon properly quantifying and contextualizing these consequences, a detailed discussion of each AGN selection method used for this study as well as possible biases introduced by these selections follows in this section.


[angle=0,width=1.0]FINAL.radiovsTIR.SPIRE.AGNvsnon.eps \includegraphics[angle=0,width=1.0]FINAL.radioqhist.SPIRE.AGNvsnon.eps

Figure 4: Left: Total infrared (TIR) luminosity plotted against the rest-frame power density at 1.4 GHz of all galaxies with reliable redshifts detected significantly in both the VLA and Herschel imaging. Galaxies which we have defined as hosting radio AGN are represented by open red diamonds and galaxies whose radio emission likely originates either solely from star-forming processes or some combination of star formation and AGN activity are shown as blue Xs. The dashed line denotes our formal best fit relation for SPIRE radio starbursts. The solid line is the best-fit relation from Del Moro et al. (2013) for “radio normal” galaxies. Right: Histogram of for the SPIRE Radio AGN and SPIRE Radio starburst samples. The two histograms are normalized to each other. Though there is significant overlap between the two samples, in large part due to SPIRE radio starbursts that likely host lower luminosity AGN, the two distributions are inconsistent with being drawn from the same parent sample.

Initial investigations of the non-thermal emission from a radio AGN hinted that such emission exhibited a steeper spectral slope () than that originating from star-forming regions (Heckman et al. 1983). However, more recent work has found that radio emission originating from an AGN shows a wide variety of spectral indices (e.g., Ho & Ulvestad 2001; Clemens et al. 2008; Ibar et al. 2010). Because there is imaging in the radio of the CFHTLS-D1 field in two bands, we initially attempted to select radio AGNs based on a spectral slope () criterion alone. However, the large uncertainties in , primarily resulting from the shallowness of our 610 MHz imaging, made a selection for individual galaxies impossible. This selection was also not possible using a statistical treatment, as there was no observed correlation between radio power density and in our data: the average for sources with was identical within the uncertainty to that of lower power density sources. As is also noted in Condon (1992), it is possible to use the strength of optical recombination lines (or their proxies) in concert with radio power densities to determine the dominant source of radiation. Again, however, this selection was not possible with our data as not all of our SPIRE-detected sample was targeted for spectroscopy. In the absence of a more evolved selection method, we resigned ourselves to a brute force power density cut, where any source with a rest-frame was considered an AGN. This is consistent with maximal radio output of normal galaxies (Condon 1992; Kauffmann et al. 2008; Del Moro et al. 2013) and consistent with the selection used by other studies (e.g., Hickox et al. 2009) including other studies of SPIRE-selected galaxies (Hardcastle et al. 2012).

To test the robustness of this cut we calculated for all radio sources matched to a high significance () SPIRE detection, a ratio which broadly characterizes the amount of light in the TIR band () relative to that of the rest-frame 1.4 GHz radio band. This metric is used in some studies to discriminate between star-forming galaxies and those hosting AGN through the excess of radio emission in AGN-dominated galaxies relative to what is expected solely from processes related to star formation. Because a wide variety of observations are used especially to proxy the amount of light in the TIR band, a wide variety of definitions of this ratio are used, and each has its own subtleties with respect to selecting AGN phenomena (see discussion in Del Moro et al. 2012). Here, we define as


Because there is a strong correlation between the luminosity in the FIR bands and the luminosity in the TIR band, both in the models we use and in general in star-forming galaxies (see, e.g., Kennicutt 1998), we chose instead to calculate using the latter. In addition, because of the wide redshift baseline of our sample and the wide variety of corrections (or effective corrections applied from the SED-fitting process) employed in this study, was calculated with the observed and modeled rest-frame and values for all galaxies rather than observed-frame fluxes and flux densities. While these choices will result in some absolute offset in relative to measurements made in other studies (i.e., Del Moro et al. 2013), we are only interested here in the relative difference between the value of in galaxies we have selected as hosting an AGN and those we consider as normal SPIRE-selected starbursts. We broke up the SPIRE-radio matched sample into radio AGN () and those galaxies likely either to be dominated by star formation or have some combination of star formation and lower level AGN activity (). The resulting histograms are presented in Fig. 4. Visually, the two distributions appear to be quite different, with radio AGNs being biased, as expected, to lower values of (i.e., higher radio flux density for a given TIR flux density). To confirm this difference, we performed a Kolmogorov-Smirnov (KS) test on the two distributions, which resulted in a rejection of the null hypothesis that the two populations were drawn from the same parent population at a % confidence level (CL). The significance of this difference remains essentially unchanged if we instead chose for our star-forming criterion (i.e., the power density limit at which galaxies are likely dominated in the radio by their star formation; Condon 1992).

However, the sample of galaxies studied here is special relative to many other studies of radio AGN host galaxies. By initially selecting galaxies on their over a large redshift range, we are probing galaxies whose s well exceed even the most vigorous ULIRGs in the local universe on which this radio selection was initially calibrated. There is, therefore, some concern that a radio AGN host sample selected on a power density cut based on the observed properties of local galaxies may not be sufficiently pure for this study. There are several reasons, however, to believe this is not the case. The statistical difference between the distributions of the values of the galaxies selected as radio AGN hosts and those galaxies observed at lower power densities strongly suggested the two populations had different mechanisms powering their radio emission. This will be further demonstrated in the companion paper to this study (Lemaux et al. 2015) where we will remove the expected , as derived from and , from the observed (-corrected) of galaxies with , resulting in a residual power density of for % of the sample. This excess can only be attributed to either an evolution in the parameter as a function of redshift or (see Chapman et al. 2010 and Bourne et al. 2011 for a discussion on the former) or an AGN component. Conversely, for sources significantly detected in the radio but with , such an excess is not observed in % of the sample. In addition, if we instead select radio AGN hosts by a radio-excess criteria using the methodologies of Del Moro et al. (2013) none of the main results of this paper are substantially changed. Because of these reasons and because the measurements of individual sources are considerably less precise than measurements and, furthermore, because adopting a cut based on the former limits both internal comparisons and comparisons to the literature which can be made without fear of differential bias, we choose to retain the definition of radio AGN hosts as those galaxies fulfilling the criterion .

Radio starburst galaxies are then defined as any galaxy detected at radio wavelengths at , of which 251 lie within in the overlap region. This population will almost certainly contain a considerable number of lower luminosity AGN and hybrid galaxies. Similarly, the method of selecting X-ray AGN and IR AGN described later in this section will also leave open the possibility of lower level AGN activity in the non-AGN sample (mostly due to the shallowness of the XMM and SWIRE imaging). However, as stated earlier, the methods of selecting AGN presented here are designed to maximize purity of galaxies hosting powerful AGN. When we compare, in the latter sections of this paper, galaxies hosting powerful AGN to those that, to the limit of our data and selection methods, do not host a powerful AGN, any contamination in the latter sample from low-level AGN activity will likely only serve to decrease the observed differences between the two samples. As a result of the previous arguments regarding the purity of the radio AGN host sample, despite the large redshift baseline of this sample (), the redshift evolution of the radio power density (e.g., Condon 1989, 2002; Dunlop et al. 2003; Smolčić et al. 2009a) has little effect on our results. Or, rather, this evolution will simply result in more radio AGN being selected at high redshifts as the redshift distribution of the 10% of the radio AGN selected hosts whose subtracted is less than (see above) is indistinguishable from that of the remaining 90% of the sample. Of the 567 radio sources in the HerMES-SWIRE-CFHTLS overlap region that were matched to optical counterparts with well-measured redshifts, 316 satisfied the radio AGN criterion.

The selection of AGN obscured by dust has been the topic of much discussion since the advent of the Spitzer space telescope. The selection of such AGN are particularly important to this study, as the host galaxies that we are selecting with SPIRE are necessarily dusty. While pervasive dust in star-forming regions of host galaxies does not necessarily have a bearing on the dust properties of its AGN, it is reasonable to assume, given that the most violent starbursting events are circumnuclear (e.g., Lehnert & Heckman 1996; Kennicutt 1998; Kormendy & Kennicutt 2004; Haan et al. 2013), that obscured AGN may be more prominent than other AGNs in a FIR-selected sample of galaxies. Indeed, this assumption has been recently given observational support (Juneau et al. 2013). A variety of different methodology have been used to select obscured AGN, mostly relying on various features in the NIR originating from re-radiated emission of the dust-enshrouded nucleus (Lacy et al. 2004; Stern et al. 2005; Donley et al. 2008; Kirkpatrick et al. 2012). In this paper we adopt the hybrid criteria of Donley et al. (2012) because of the high purity and completeness demonstrated for samples obtained using their methodology. In Fig. 5 we illustrate this selection applied to our sample, plotting the IRAC flux density ratios for all galaxies detected in the overlap region at in all four IRAC bands along with the criteria given by Donley et al. (2012). Of the 1107 galaxies in the overlap region detected at in all four IRAC bands and matched to a optical counterpart with a well-defined redshift, 172 were selected as IR-AGN using this method.



Figure 5: Spitzer/IRAC NIR color-color diagram for all galaxies in the overlap region (see §2.1) with detections in all four IRAC bands. Both galaxies which are undetected in Herschel/SPIRE and galaxies which are marginally detected in Herschel/SPIRE () are plotted as blue Xs. Galaxies detected at in at least one of the three SPIRE bands are shown as open red squares. The solid black lines are adopted from Donley et al. (2012) and denote the region used to select IR-AGN. While a large fraction of galaxies with detections in all four IRAC bands are also detected in Herschel (because of our selection method) many galaxies, including many in the AGN region, are not. The depth of the SWIRE imaging appears sufficient to probe AGNs with NIR color properties spanning the entirety of the selection box.

Selection of AGNs in the X-ray was considerably easier than selection in the radio or NIR. Because of the shallowness of the XMM coverage in the CFHTLS-D1 field and the relatively low volume probed at low redshifts (), nearly all (%) of the objects detected in our X-ray catalog and matched to an optical counterpart with a well-measured redshift had X-ray luminosities in excess of what can be produced by normal star formation ( ergs s; Ranalli et al. 2003). Thus, any object detected in our X-ray catalog and matched to a SPIRE source, with the exception of only a few sources, was considered a X-ray bright AGN. It should be noted that this selection process is highly incomplete as the X-ray data are only complete to until redshift . In addition, it has been shown that X-ray weak or X-ray obscured AGNs selected by a combination of rest-frame optical line ratio diagnostics and rest-frame colors (Yan et al. 2011) or stellar mass (Juneau et al. 2013) constitute a non-negligible fraction, perhaps as much as %, of the total AGN population in exactly the type of infrared bright starbursts studied in this paper. Our selection is also incomplete with respect to those AGN, or AGN-like phenomena, whose presence is inferred by other rest-frame optical indicators (e.g., Yan et al. 2006; Lemaux et al. 2010). However, such objects should be less prevalent in our sample, as such phenomena are typically associated with quiescent red-sequence galaxies, a galaxy type that SPIRE observations are blind to. The lack of a full spectral sample and the large redshift baseline of our sample precludes the possibility of making similar selections for our sample. While we mitigate the effects of this completeness somewhat by also selecting AGN in the infrared and radio, all percentages quoted in this paper for the X-ray AGN population are hard lower limits which likely severely underestimate the true number of high-excitation, high-accretion-efficiency AGNs in our sample. Because of this, and because of the shallowness of the XMM data in the CFHTLS-D1 field relative to studies dedicated to the interplay between X-ray AGN and their Herschel-detected hosts (e.g., Harrison et al. 2012b), we do not attempt later in this paper to study the X-ray AGN population alone, but rather as part of a larger ensemble comprised of all AGN types.

Unobscured type-1 AGN, a population intimately related to the X-ray bright AGN population selected by XMM, were drawn from those SPIRE-matched galaxies which were targeted for spectroscopy either by VVDS or ORELSE. These AGN are characterized by their broadline emission features (typically with km s) and strong emission of those features which can result only from ionization by extremely high energy photons. In such objects, the AGN continuum emission also tended to dominate the stellar continuum in these objects, suggesting significant levels of activity. For both the ORELSE and VVDS surveys these galaxies were noted in the visual inspection process and were selected based on these notes. The spectrum for each AGN selected in this manner was re-inspected visually to confirm both the presence of an AGN and the quality of the redshift measurement. In total, 130 X-ray and 77 type-1 AGN were identified in the overlap region. In the left panel of Fig. 6 we plot a false color image of the Herschel/SPIRE mosaic in the overlap region. The mosaic was generated with linear scaling and made use of the python-based APLpy package (http://aplpy.github.com). The false color mosaic is contrasted by the right panel of Fig. 6 in which we plot the spatial distribution of all radio starbursts and all AGN selected in this section detected at a significance of in at least one of the three SPIRE bands against the backdrop of a density map composed from the 1753 SPIRE-detected galaxies fulfilling the criteria of the previous section. The density map was created using the methodology of Gutermuth et al. (2005), with the parameter dictating the number of nearest neighbors, , set to . Curiously, the AGN hosts appear to largely avoid the regions of highest SPIRE-detected galaxy density, especially in the case of radio AGN, an effect which cannot be caused by observational or selection issues. The environments of radio AGN in particular will be discussed in detail in a companion paper (Lemaux et al. 2015). Summarized in Table 3 are the total number of objects in the overlap region detected in each of the observations relevant for AGN selection and matched to an optical counterpart along with the total number of these objects which were selected as AGN hosts.


[angle=0,width=1.0]VVDS.2h.Herschel.colormosaic.fromAPLpy.edit2.eps \includegraphics[angle=0,height=0.98]XLSS005.Herscheldensity.AGNHerscheldet.wtypes.3sigma.eps

Figure 6: Left: False color image of the Herschel/SPIRE mosaic in the CFHTLS-D1 overlap region (see §2.1). The SPIRE 500m, 350m, and 250m images served as the RGB channels, respectively. Right: Sky distribution of multiwavelength galaxies also detected with SPIRE plotted against the backdrop of a density map of the 1753 galaxies in our final SPIRE sample. The spatial density of the full SPIRE sample is a convolution of the true spatial distribution (left panel) and a number of metholodological and observational effects. Bright stars, which prevent both spectroscopy and high-quality photometric redshifts, account for the presence of most of the voids seen in the background density. This includes voids around areas of otherwise high densities of SPIRE-detected sources, such as the grouping seen in the left panel at .
Sample 43
Radio 567 316
IRAC 28394 172
4-Channel IRAC 1107 172
VVDS Type-1 AGN 77 77
X-ray 137 130
SPIRE44 1753 18145
Table 3: Multiwavelength objects in the HerMES/SWIRE/CHFTLS overlap region

It is necessary at this point to make a final mention of subtleties regarding the SED-fitting process, this time for the AGN populations discussed in this section. Depending on the strength of the AGN and the dust properties of both the host and the region surrounding the nucleus, an active nucleus can dominate the stellar and nebular components of a host galaxy at a variety of wavelengths. This is particularly a problem for unobscured X-ray emitting AGN, quasars, and type-1 AGN, where the UV portion of the SED can almost entirely originate from the AGN (Vanden Berk et al. 2001; Gavignaud et al. 2006; Polletta et al. 2007). Obscured AGN can also present a problem, as the very features that allow their selection in the NIR can confuse ones ability to measure quantities such as stellar mass. We tested the SED fits of the parent samples of all the AGN (i.e., not only those detected by SPIRE) by comparing the photo-s and spec-s of those galaxies with reliable measurements for both and by comparing the value of the fits of the host galaxies of the various AGN with the full photometric sample. For radio- and IR-selected AGN the results were encouraging, no increase in the NMAD or the catastrophic outlier rate was observed for the host galaxies of those populations relative to the full photometric sample, nor was the statistic appreciably higher. However, this was not the case for the X-ray AGN and the type-1 AGN where both statistics were significantly higher. To alleviate this, we fit hybrid galaxy/QSO templates following the methodology of Salvato et al. (2009) using the templates of Polletta et al. (2007). This method resulted in a reduction of the catastrophic outlier rate of the full X-ray AGN host population with high-quality spectra from % to 6% and a reduction of the NMAD to . While both of these values are considerably higher than that of the full photometric sample, our large redshift bins and the fact that galaxies hosting X-ray AGN make up a small fraction of our overall SPIRE-detected AGN host sample (%), mean that this lack of precision is of limited consequence for our results. Because decomposing the stellar light from the light originating from the AGN is inexact, physical parameters other than redshift, such as stellar mass, are estimated for such galaxies with questionable accuracy. This, however, is of less consequence, as the stellar mass of such galaxies is largely considered only as part of the full SPIRE population which is comprised of only % of galaxies hosting X-ray AGN. For type-1 hosts, marked improvement in photo- reliability or accuracy was not observed upon the inclusion of the hybrid galaxy/QSO templates, as five out of the ten type-1 hosts detected at in one of the three SPIRE bands remained as catastrophic outliers. The type-1 hosts which are catastrophic outliers, however, make up an even smaller fraction of our full SPIRE and full AGN samples than the X-ray AGN hosts (0.2% and 2.7%, respectively). Furthermore, since the of these galaxies, which is the quantity most relevant to the analyses involving AGN host galaxies in this study, was calculated using spectroscopic redshifts when available, redshifts which the type-1 AGN hosts have by definition, these galaxies are retained as part of the full sample.

4 The properties of star-forming galaxies detected by SPIRE

We break the remainder of the paper into two parts. In this section we discuss the properties of the 1753 Herschel bright galaxies selected in the previous sections making no differentiation between those galaxies hosting AGN and those with dormant nuclei. Because the former is a small fraction of our sample (%), and because the fraction of AGN which may have their optical continua dominated by an AGN is an even smaller fraction of our sample (i.e., type-1 AGN and perhaps some X-ray AGN, 2%), the results in the first part of this section will be broadly applicable to normal dusty star-forming or starbursting galaxies, with a minimal contamination to our sample of stellar masses, rest-frame colors, or magnitudes by galaxies hosting an AGN. In the subsequent section we will divide our sample into two distinct sub-groups, those SPIRE-detected galaxies with and without an AGN, and compare the properties of each subsample.

4.1 Color, magnitudes, and stellar masses of the full SPIRE sample

We begin the investigation of the full SPIRE sample47 by plotting in Fig. 7 rest-frame color-magnitude diagram (CMD) and color-stellar-mass diagram (CSMD) of the full SPIRE sample in four different redshift bins. Rest-frame colors and magnitudes are derived from our SED-fitting process described in detail in §3.3. Also plotted in Fig. 7 are those galaxies in each redshift bin that went undetected in SPIRE but which fulfil all the other criteria of the full SPIRE sample (i.e., having a well-measured redshift between and sufficient photometry to perform SED fitting). For galaxies plotted in both the CMDs and CSMDs we additionally require galaxies be detected in the band at a magnitude brighter than the completeness limit of the CFHT-WIRDS imaging () and have a SED-derived stellar mass that is higher than the stellar mass completeness limit for each bin (see Fig. 7).

What is immediately apparent from the CMD is that the SPIRE-detected galaxies span a full range of colors and magnitudes, ranging from the brightest blue galaxies to the faintest red galaxies in each redshift bin. This is also true on inspection of the CSMD: the full SPIRE sample essentially spans the entire phase space encompassed by the SPIRE-undetected galaxies. The fraction of SPIRE-detected galaxies in the blue cloud (defined as roughly )48 is quite high in every redshift bin (70% on average), a consequence of the high rate of star formation in these galaxies. However, a non-negligible fraction (12.5% on average) of galaxies in the full SPIRE sample have colors consistent with the red sequence (defined as ) for redshifts where the red sequence is clearly separable from the locus of star-forming galaxies (i.e., ). The remaining galaxies in the full SPIRE sample (17.5% on average) have intermediate colors, lying in a region of phase space commonly referred to as the green valley. This portion of phase space is typically thought to house a transitory galaxy population, comprised both of galaxies that were previously on the red sequence but are undergoing low-level rejuvenated star formation and those galaxies that have finished their star formation and are slowly transitioning onto the red sequence (Silverman et al. 2008; Kocevski et al. 2009b; Vergani et al. 2010; Mendez et al. 2011; Fang et al. 2012; Gonçalves et al. 2012). While the number of galaxies located in the green valley is a sizeable fraction of the full SPIRE sample, this fraction only slightly exceeds the same quantity for all optical galaxies (%). Still, it is surprising that galaxies selected to be vigorously forming stars would (seemingly) be equally likely to be in a transition phase as galaxies which are either quiescent or undergoing low-level star formation.


[angle=90,width=1.0]SPIREvsnon.nolumbins.NUV-r.Mr.restframe.eps \includegraphics[angle=90,width=1.0]SPIREvsnon.nolumbins.NUV-r.stellarmass.restframe.eps

Figure 7: Rest-frame color magnitude (left) and color stellar mass (right) diagrams of the full SPIRE sample (filled red circles and dashed red histograms) and galaxies undetected in SPIRE (small black points and open black histograms). Galaxies are binned into four different redshift bins and histograms of each sample are normalized (such that the maximal value is unity) for each redshift bin. Only those galaxies detected at high significance in CFHT-WIRDS -band imaging are shown in both panels. In each panel the sample is restricted to those galaxies with stellar masses above the stellar mass limit for each redshift range (given in each panel). The arrow denotes the change in color and (where appropriate) absolute magnitude for an extinction value of using a Calzetti et al. (2000) extinction law. The arrow is vertical in the right panel because the measured extinction values are already incorporated into the SED-fitting process. The SPIRE-detected sample spans the full range of absolute magnitudes, stellar masses, and, surprisingly, colors in each redshift bin.


Figure 8: Rest-frame vs. color-color diagram for both the full SPIRE sample and galaxies undetected in the SPIRE imaging. The gray lines in each panel divide quiescent galaxies (above the lines) from star-forming galaxies (below the line) and is determined using a modified version of the methodology employed by Ilbert et al. (2013). The meanings of colors and symbols are nearly identical to Fig. 7. The one exception is the histogram plotted along the ordinate in each panel, in which galaxies are projected into an “effective color” defined as the offset from the dividing line. This representation of the sample serves to largely mitigate extinction effects, and as such, a large fraction (95%) of the full SPIRE sample are now in the phase space coincident with normal star-forming galaxies.

The prevalence of infrared-selected starburst galaxies in the green valley has been observed for a Spitzer/MIPS survey of 70 sources in the COSMOS field (Kartaltepe et al. 2010b). However, as is discussed extensively in that study, the high number of dusty starbursts in the green valley can be due to either properties intrinsic to the stellar populations of these galaxies or to dust extinction effects. In the case of the latter, the galaxies themselves are not inherently in a transitory phase, but rather appear so due to heavy extinction of the stellar continuum coming from the dusty component in the galaxy (see, e.g., Cardamone et al. 2010). While we could, in principle, correct the color and magnitude of each galaxy for extinction effects, the extinction values coming from the SED fitting process are not robust enough to do this in practice. In principle, we could also use the FIR-derived extinction values, values that we derive for a subsets of the full SPIRE sample later in this paper (see §4.4.1) to correct the colors and magnitudes of each SPIRE-detected galaxy. However, since we are comparing to an optically selected sample of galaxies undetected in SPIRE in each redshift bin, the colors and magnitudes of the comparison sample cannot be corrected in this manner, which limits the usefulness of the method.

Instead, as is done in Kartaltepe et al. (2010b), we employ a color-color diagram (CCD) to investigate which of these two possibilities is more applicable to the full SPIRE sample. Plotted in Fig. 8 is the vs. CCD of the full SPIRE sample against the backdrop of the SPIRE-undetected sample. As was done for both panels of Fig. 7, we require here that any galaxy plotted be detected in the band at a magnitude brighter than the completeness limit in order to ensure accurate rest-frame NIR magnitudes. Overplotted on the CCD is a dividing line between galaxies with quiescent colors and those consistent with young stellar populations (see Ilbert et al. 2013, Muzzin et al. 2013, and Arnouts et al. 2013, and discussions and references therein). The dividing lines given in Fig. 8 are slightly different than those presented in Ilbert et al. (2013) due to a different methodology used to determine the minima in color-color space between the two populations. In this paper, the CCD was rotated to search for the optimal angle of approach for the diagonal line in each redshift bin, requiring that this angle must remain the same in each bin. The optimal angle of approach was defined as the angle which maximized the disparity between the local minimum and the local maxima of the colors projected on the new axes (i.e., the axis defined by the diagonal line and an axis perpendicular to this line). The diagonal line was truncated for each redshift bin at the point where it began to intersect the star-forming locus, i.e., where continuing the line diagonally began to decrease the disparity between the minimum and local maxima of the rotated color-color histograms. The histograms along the ordinate of each panel in Fig. 8 no longer simply represent the distributions of colors. Rather an effective color is plotted, defined as the effective vertical offset of each galaxy from its projection onto the dividing line. The analysis of galaxies in this color-color space has the distinct advantage of almost completely eliminating extinction effects, as galaxies with higher extinction move parallel to the main discriminator nearly excluding the possibility of dust causing star-forming galaxies to appear quiescent.

In this context, the interpretation of full SPIRE sample changes dramatically. While in color-magnitude and color-stellar-mass space a large fraction of SPIRE-detected galaxies had colors consistent with the red sequence or green valley (%), this effect is largely gone in the CCD. The galaxies of the full SPIRE sample have colors that place them on the star-forming sequence (i.e., below the dividing line), with only 8.0%, 5.5%, 6.5%, and 0% of SPIRE-detected galaxies having colors consistent with the quiescent part of this diagram in the , , , and bins, respectively. This is in contrast to the SPIRE-undetected sample, in which 41.3%, 34.0%, 26.8%, and 3.8% of the galaxies have colors consistent with quiescent galaxies in the same bins, respectively. Transition galaxies are also more apparent in this diagram, as galaxies between the quiescent region and the star-forming sequence are there as a result of either age or metallicity effects. While there still remains a small fraction of galaxies that appear to be undergoing a transition, especially in the redshift range , the main conclusion from this diagnostic is that % of the galaxies in the full SPIRE sample have colors consistent with star-forming galaxies once extinction effects are minimized.

4.2 The star forming main sequence


[angle=90,width=1.0]Herschel.2h.SFRvsmass.wAGNtypes.zbins.3sigmacut.eps \includegraphics[angle=90,height=0.71]Herschel.2h.SSFRvsmass.wAGNtypes.zbins.3sigmacut.eps

Figure 9: Left: Total infrared (TIR)-derived plotted against stellar mass for volume-limited samples of all optically matched galaxies detected in at least one SPIRE band at for four different redshift bins. is shown on the right ordinate for reference. Because the conversion employs a Salpeter IMF, the stellar masses shown here are converted from a Chabrier to a Salpeter IMF. Black points represent galaxies detected in SPIRE, while galaxies detected in radio imaging and those galaxies hosting AGN are overplotted in larger symbols. The main sequence relationship for star-forming galaxies is adopted from E11 for galaxies at , an average of the relation in E11 and that relation of the lowest stellar mass bin of R11 for galaxies at , and the relation from the lowest stellar mass bin in R11 for galaxies at and is plotted as a solid gray line. The dashed line denotes the transition point between galaxies classified as normal star forming to those classified as starbursts (4 the main sequence threshold). There is no evidence in three of the four redshift bins of a correlation between and to the depth of our SPIRE imaging. Right: Specific plotted against stellar mass for the same samples as the left panel. Symbols and lines are identical to the left panel. The dynamic range of Herschel is small in any given redshift bin while the stellar mass range is not, which leads to the observed trend.

In this section we consider limited subsamples of the full SPIRE sample in the same four redshift bins as were used in the previous section. These cuts, which are shown in the right panel of Fig. 3, allow us to consider subsamples which are volume limited, such that the effects of Malmquist bias are minimized. This is essential for the study that follows. Plotted in the left panel of Fig. 9 are the stellar masses and TIR s of the volume-limited full SPIRE sample in four redshift bins. Due to the method in which s are calculated, SED-fit stellar masses are converted from a Chabrier IMF to a Salpeter IMF for the remainder of this section using the relation of Ilbert et al. (2010). One of the most striking features of the left panel of Fig. 9 plot is the high rate of star formation observed in galaxies hosting AGN relative to SPIRE-detected galaxies where no AGN is observed, a subject which will be covered extensively in §5.2. In this section, we make no distinction between galaxies hosting a powerful AGN and those without an active nuclei, and instead consider the question: is there a relationship between stellar mass and in the volume-limited SPIRE sample, and, if so, does that relationship evolve as a function of redshift?

In seminal Herschel study, Elbaz et al. (2011; hereafter E11) used a combination of Spitzer, AKARI, and Herschel observations to suggest the existence of an infrared main sequence for star-forming galaxies detected in the mid- and far-infrared. In this picture, stellar mass () is intimately linked to the amount of star formation in such galaxies, such that, if a galaxy is star forming, the higher the stellar mass the more vigorous the star formation. Multiple studies over a large range of redshifts using a variety of different selection methods have found that, in star-forming galaxies, increasing generally leads to an increase in (e.g., Brinchmann et al. 2004; Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007a; Santini et al. 2009; Gonzlez et al. 2011; Muzzin et al. 2012; Koyama et al. 2013), but this phenomenon was observed in E11 for the first time with a purely IR-selected sample of galaxies. These observations have in turn lead to significant tension with semi-analytic models of galaxy formation, specifically in the evolution of this relationship as a function of redshift and its overall normalization (for an excellent review of this topic, see Dav 2008). The uniqueness of the data presented in E11 was, however, that the sample presented was strictly IR-selected, probed a large range of redshifts, and probed a dynamic range of stellar masses and TIR-derived s of greater than two orders of magnitude. Using the observed - relationship and arguments based on the measured compactness of galaxies as measured in their IR photometry, E11 created a new method for defining starbursts by relating the current of a galaxy to its past averaged (Kennicutt 1983). The past averaged was proxied by the - infrared main sequence (MS), and therefore a galaxy significantly departing from the MS at a given redshift is considered a starburst galaxy (more specifically, when the ratio of the of a galaxy and a typical value of the for a MS galaxy at the same redshift exceeds two). This definition of starburst, which has been heavily used in the literature (see, e.g., Magnelli et al. 2012 and references therein), relies heavily both on the premise that the infrared MS exists and that the typical value of the for an MS galaxy can be measured precisely as a function of redshift. The latter relies on both a precise quantification of the relationship between and , the precision with which each parameter can be estimated, and the assumption that this relationship does not change based on the selection method of the sample nor the dynamic range of s or stellar masses probed in a given sample.

While our observations do not reach the depth achieved by E11 in the GOODS-N/S fields, we probe an area which is more than order of magnitude greater than the observations of E11. At the expense of fainter sources, our observations are more sensitive to rare galaxies which are undergoing extremely vigorous star-formation episodes. Galaxies with s in excess of 100 yr, a population which is largely absent in the E11 sample, are found in abundance in our observations for all redshifts greater than . To our knowledge, a statistical sample of this population has been investigated by only one other study (Rodighiero et al. 2011; hereafter R11, see also Rodighiero et al. 2014 for a followup study).

In our sample, we see little evidence for an infrared MS for the / range probed by our SPIRE observations. In all redshift bins the correlation between and is weak: the Spearman rank correlation coefficient is 0.17 on average and only one bin () shows a deviation from the null hypothesis that the two quantities are uncorrelated. The results of these tests remain essentially unchanged if we apply the stellar mass completeness limits from the previous section. This lack of correlation argues against the existence of the - MS for IR-selected galaxies at all redshifts that have s in excess of the limits used to form the volume-limited samples here49. While this lack of correlation observed in the full SPIRE sample cannot be due to differential effects as a function of stellar mass, e.g., imposing a stellar mass and limits have little effect on the lack of correlation, our data do not probe deep enough at any redshift to claim that such a correlation does not exist for galaxies anything but the most prodigious star formers. However, as mentioned previously, the large area coverage of our SPIRE reduces cosmic variance and allows for the characterization of starburst galaxies over a large range of s by adopting the formalism of those studies which probe deep enough to observe such a correlation.

In order to determine the fraction of starburst galaxies in each redshift bin we now turn to the s of the full SPIRE sample plotted in the right panel of Fig. 9. Again a Salpeter IMF is adopted in calculating stellar masses. Galaxies in our sample are observed over a large range of s, spanning nearly three orders of magnitude in all redshift bins. This effect is largely driven by the stellar mass of the galaxies, since the Herschel observations have a limited dynamic range for a given redshift range. We now adopt the definition of a starburst from a combination of those definitions in E11 and R11 for the remainder of this section. Applying these cuts to the full SPIRE sample at all redshift bins and imposing the stellar mass cuts from the previous section results in a fraction of starburst galaxies, so defined, of % over all redshift ranges. The fractional contribution of starburst galaxies to the global star formation rate density () in a given (stellar-mass and volume-limited) redshift bin varies immensely, ranging from a 28% contribution in the lowest redshift bin to 61% in the highest redshift bin. These numbers are, however, subject to the relatively shallow limits in each bin and should be viewed as strict upper limits in each redshift bin. This exercise shows that starburst galaxies make up a non-negligible fraction of the SPIRE population and contribute significantly to the overall number of stars being formed at all redshifts, at least to the depth of our data.

While the observed starburst fraction in the full SPIRE sample is considerably higher than that found in R11 for galaxies between (2-3%), except for the lowest redshifts probed in our sample, the data presented here does not go deep enough to probe the BzK selected samples of R11 that make up the bulk of the lower and lower stellar mass sample. Thus, due to the differing sample selections and various differential biases, these numbers can not be directly compared to those of R11 in anything other than a broad sense, and, given the limits of our data, are not in direct contradiction with their results. In an attempt to make a more rigorous comparison, a volume-limited sample was selected between from our data by imposing yr and a stellar mass cut of , which allows for a direct comparison with the sample of R11. In this volume-limited sample we find that 35.6% of the at for IR-bright galaxies with s in excess 200 yr is contributed by galaxies classified as starbursts. While this number is marginally inconsistent with the % contribution reported by R11, within the relative uncertainties of the two samples it is plausible that the two results are consistent.

However, a larger issue remains. Wide variations in the starburst fraction and fractional contributions to the by starbursts of up to nearly an order of magnitude in our own sample result if we impose slightly different stellar mass completeness criteria, redshift ranges, or vary the definition of a starburst by as little as 0.3 dex. The fraction of starburst galaxies and their contribution to the global appear highly subject to the stellar mass completeness limits and the definition of a starburst. The latter is intimately linked both with the precision at which one can constrain the MS at a given redshift, the method used and precision to which one can measure stellar masses and s (see Rodighiero et al. 2014 for a thorough disucssion of the latter), and the methodology used in determining galaxies are significantly offset from this sequence. The extreme level of variance observed in our data clearly shows the necessity of precisely and accurately determining the MS at a given redshift and to properly and precisely quantify the stellar mass and limits of the data in order to use such analysis to constrain galaxy models with any sort of reasonable accuracy. The problem of constraining the contribution of starburst galaxies to the overall star formation rate density of the universe using this framework appears, in the true sense of the term, a chaotic one.

While this discussion may appear to be a matter of semantics, the operational definition of the term starburst is an extremely important one in the study of galaxy evolution as such galaxies are connected to a physical phase of galaxy formation, which is then used in turn to contextualize the evolution of galaxies across a wide variety of environments (e.g., Goto 2006; Moran et al. 2007; Oemler et al. 2009; Dressler et al. 2013). Given the large variation in results seen both in our own sample and from contrasting our sample with those results from the literature, it is clear that significant work is necessary to understand the relationship between Herschel-bright sources and galaxies undergoing normal star formation, to better understand the nature of the star-forming MS, and to find proper methodologies to robustly define starburst populations. It is outside the scope of this paper to attempt a precise definition of the term starburst, and, indeed, our data generally preclude such possibility. Thus, we instead apply the term starburst throughout this paper as a proxy for “SPIRE-detected, strongly star-forming galaxy” in an agnostic view.

4.3 General properties of SPIRE galaxies at different TIR luminosities

In §4.1 the color, magnitude, and stellar mass properties of the full SPIRE sample were investigated without differentiating galaxies based on their TIR luminosities. While doing this allowed us to determine, broadly, the properties of our sample, marginalizing over TIR luminosity limits our ability to discern the more subtle properties of the full SPIRE sample, as it has been observed that dusty starburst galaxies of differing TIR luminosities exhibit different dust temperatures and SEDs (Symeonidis et al. 2013; Hwang et al. 2010b), rest-frame colors (Kartaltepe et al. 2010b), and are found to lie in different environments (Kocevski et al. 2011a). In this section we move from considering the full SPIRE sample as a whole to galaxies in individual TIR bins. We define here three samples: a “low-” sample (, low luminosity luminous infrared galaxy (LIRG) level emission), a “moderate-” sample (, high luminosity LIRG level emission), and a “high-” sample (, ultra-LIRG (ULIRG) level). We restrict the redshift range of each of these subsamples to redshifts , which is the redshift range containing the bulk of our sample. The TIR levels and redshift range were chosen in such a way to minimize redshift-dependent effects in each bin (i.e., sampling a representative population over the chosen redshift range) while still maximizing the number of SPIRE-detected galaxies in each bin. Furthermore, we attempted to construct the bins such that roughly an equal number of SPIRE-detected galaxies with spectral redshifts were placed in each bin. There is some concern that galaxies of the same have redshift-dependent properties, such that performing this analysis over such a large redshift bin will confuse galaxy populations which are fundamentally different. However, for galaxies which we could compare over the entire redshift range (i.e., the high- sample), the spectral, photometric, color, and stellar mass properties were all nearly identical from low () to high () redshift, which likely precludes the possibility that redshift-dependent properties will confuse the analysis which follows. The sample properties of the galaxies in each bin are summarized in Table 4.

Sample Range Equivalent Class 50
low- low-lum LIRG 175 42 11.35 -21.90 3.74
moderate- high-lum LIRG 266 52 11.70 -22.12 3.50
high- ULIRG 358 43 12.30 -22.67 3.14
Table 4: Properties of the luminosity subsamples of the full SPIRE sample

4.4 Spectral properties of SPIRE galaxies at different TIR luminosities

We begin by contrasting the VVDS rest-frame NUV/optical spectral properties of each of the three subclasses. The typical individual VVDS spectrum in our sample has an median signal-to-noise (S/N) ratio of order pix, which is not sufficient for much of the analysis presented in this section. Thus, rather than study the spectra of individual galaxies we instead chose to create an inverse-variance weighted composite spectrum (referred to hereafter as “coadded spectrum” or “coadd”) of the galaxies in each luminosity bin using the methodology described in Lemaux et al. (2009, 2012) adapted for the VVDS data. As was the case in Lemaux et al. (2012), the average flux density of the spectrum of each galaxy is normalized to unity prior to the stacking process. The normalization factor for each spectrum is determined over the range of wavelengths free of poorly subtracted airglow lines or other reduction artifacts. This stacking and weighting scheme has two main advantages: it allows to increase the S/N of our spectra significantly (the coadded spectrum in each of the three bins has S/N50), while allowing us to retain the ability to study the average galaxy in each luminosity bin (though subject to the caveats discussed in Lemaux et al. 2012). It is important to note that the subsample of galaxies in each bin with available high-quality VVDS spectra fully spans the stellar mass, color, and absolute magnitude phase space of their full photometric parent samples, and, further, are distributed through this phase space in an approximately representative manner relative to their parent samples. In addition, the distributions of the subsample of those galaxies with high-quality VVDS spectra have a (scaled) distribution and a median that are essentially identical to the parent photometric sample. The distributions of the two samples in an CMD also appear as nearly scaled versions of each other. The observed similarity in these two phase spaces precludes the possibility that a large population of TIR-bright, highly dust-extincted galaxies are lost in the spectroscopic sample. Thus, we feel confident in applying the results of the spectral analysis presented in this and the following two sections to the full galaxy populations of each bin.


[angle=0,width=1.0]Herschel.threelumbins.3sigmacut.zcut.overplot.eps \includegraphics[angle=0,width=1.0]Herschel.threelumbins.3sigmacut.zcut.differenced.eps

Figure 10: Left: Unit-weighted spectral coadditions of VVDS spectra of galaxies detected in at least one SPIRE band at for three different luminosity bins. Each coadded spectrum contains galaxies, with only VVDS galaxies being used as ORELSE galaxies comprise only % of any given sample. Important spectral lines are overplotted and spectra are offset vertically for clarity. Though all subsamples of the full SPIRE sample show strong [OII] emission and strong Balmer features, which indicate significant ongoing star formation, the subsamples also show marked differences. Right: Differenced coadded spectra of the three subsamples. The high- subsample appears to be comprised of galaxies which, on average, have stronger [OII] emission and much stronger Balmer absorption features than either of their lower- counterparts. Also clearly visible are weaker features associated with older stellar populations (CaII H&K and G-Band) in the average high- galaxy relative to the other two samples. The average galaxy in the moderate- subsample largely continues this trend relative to that of the low- subsample.

Plotted in Fig. 10 is the coadded spectra of those low-, moderate-, and high- with available VVDS spectra. Qualitatively, the average galaxy in each bin appears to exhibit considerably different spectral properties from the average galaxies in the other two bins. Considering the main features marked in Fig. 10, it is apparent that the average high- galaxy has stronger 3727Å [OII] emission, stronger Balmer features (i.e., 4340Å H, 4101Å H, and the higher order Balmer lines blueward of the CaII H&K features), and weaker spectral features typically associated with older stellar populations, i.e., CaII H&K, 4305Å -band, and the strength of the continuum break at 4000Å known as , than its lower TIR-luminosity counterparts. For the latter point, the one exception is the 3969Å CaII H line, which appears stronger in the average high- galaxy. However, the Balmer series absorption line 3970Å H is coincident with this feature at the resolution of the rest-frame VVDS spectrum. Given the strength of the Balmer features and the strong asymmetry in the strengths between the CaII H&K lines (see Rose et al. 1985 and Rumbaugh et al. 2012 for a detailed discussion of this asymmetry), it is almost certainly the case that absorption coming from younger stellar populations, i.e., B- and A-type stars, is the dominant contribution to the observed feature. The Balmer series lines, particularly that of H, also appear progressively wider at progressively higher due to Stark broadening, an effect attributable to a larger fractional contribution of earlier A-type stars in the higher coadds (Struve 1926; Mihalas 1964). These qualitative observations are reinforced by analyzing the spectra quantitatively. For each spectrum, bandpass equivalent width () measurements of the [OII] and H features were made using a methodology almost identical to Lemaux et al. (2012). A slight modification was made to this method to accommodate the decreased resolution of the VVDS data by widening the bandpass used to measure the H feature by 6Å on each side of the nominal Fisher et al. (1999) bandpasses. The strength of , typically used to measure luminosity-weighted stellar ages (e.g., Kauffmann et al. 2003), was also measured for each coadd using the ratio of continua blueward and redward of 4000Å, as defined by Balogh et al. (1999), using the methodology of Lemaux et al. (2012). These values are given in Table 5.

In Fig. 10 we also plot differenced spectra for each of the three combinations of TIR luminosity bins. The conclusions reached previously by a simple qualitatively visual inspection of the individual coadds are broadly confirmed here: the contribution from young stellar populations is dominant in the average high- galaxy relative to the average galaxies in the other two bins. Taking the extremes of the distribution, features associated with recently created stars in the high-/low- differenced spectrum all appear strongly in absorption, meaning the amplitude of the absorption is significantly higher in the average high- galaxy. Conversely, the 3933Å CaII K feature and the -band appear in emission in the differenced spectrum, implying a higher fractional contribution from late-type (i.e., older) stars in the average low- galaxy. Visual inspection of the high-/moderate- differenced spectrum implies that the luminosity-weighted stellar population of the average moderate- galaxy is intermediate to that of the other two samples, with significantly less contribution from early-type stars (proxied by the Balmer absorption) than that of the high- sample, but also significantly less contribution from late-type stars than that of the low- sample. The range plotted in Fig. 10 was determined by requiring that at least half of the spectra in each luminosity subsample went into the coadded spectrum over the entire wavelength range. Restricting the wavelength of interest in this manner is necessary so that differences between the various populations can be attributed to properties inherent to the galaxy population in each bin and not to artifacts from the stacking process. This point is particularly important here as we begin to to discuss the differences in average stellar continua and spectral fits to stellar population synthesis models, as the loss of spectra over a particular wavelength range can artificially induce a significant change in the color of the spectra.

Sample 52 53
low- -8.630.14Å 4.11 (4.30)0.11Å 1.2840.003 (1.2240.004) 0.73 0.6610.005
moderate- -14.470.13Å 3.80 (4.12)0.10Å 1.1370.002 (1.0810.003) 1.65 0.7030.004
high- -14.990.15Å 6.77 (7.09)0.14Å 1.1150.003 (1.0500.004) 9.00 0.8070.005
Table 5: Composite spectral properties of the luminosity subsamples of the full SPIRE sample

A/K Ratios and Continuum Measurements

Now that we are confident that the full spectral range can be used to probe properties inherent to each galaxy population, we move from measuring single indices to considering the coadded spectra as a whole. To further quantify the difference in the average galaxy of each bin, two component stellar population synthesis models from Bruzual (2007; hereafter CB07) were fit to each coadded spectrum to determine the average ratio between the contribution of late-type and early-type stars commonly referred to as an A/K ratio. The CB07 models are modified versions of the original Bruzual & Charlot (2003) models updated with a new prescription to account for the mass loss of thermally pulsing asymptotic giant branch stars (Charlot, priv. comm.). For the purposes of the analysis presented here, the results obtained using CB07 models are negligibly different from those obtained with those of Bruzual & Charlot (2003). The method used to calculate A/K ratios for each spectrum is similar to that used in Quintero et al. (2004) and Yan et al. (2006), though different models are used here that result in offsets in the absolute value of the A/K ratio. Since we are interested only in the relative difference between the average galaxy in each TIR luminosity bin, this absolute offset makes no difference to the results in presented in this study. The A component of the spectrum was represented by a single burst model with a Gyr, “observed” 300 Myr after the starburst, and smoothed to match the resolution of the VVDS coadded spectra. The time after the starburst was chosen to be consistent with those models used in the works mentioned above and is roughly the time at which the luminosity-weighted contribution of A-type stars to the synthetic spectrum is maximized (i.e., when H is maximized, see Fig. 12). The K component of the spectrum also used a smoothed single burst Gyr model but observed 5 Gyr after the initial starburst, a time when late-type stars dominate the synthetic spectrum. For each coadd, the [OII] emission feature was interpolated over and fit a linear combination of the two components to each coadded spectrum over the entire wavelength range presented in Fig. 10 was performed, adopting the best-fit as that linear combination which minimized the statistic. The A/K ratio was then defined as the quotient of the amplitudes of the A and K components of the best-fit composite model. The results of these fits, given in Table 5, add to the growing body of evidence that the luminosity-weighted stellar populations in the average galaxy of each of the three TIR luminosity bins is significantly different. The measured A/K ratio of the average high- galaxy shows that the continuum light of this galaxy is completely dominated by early-type stars (A/K=9.0). Conversely, the stellar continuum of the average low- galaxy has a K-star component which is more prominent than that of A-type stars (A/K=0.73), with the moderate- average galaxy again falling at an intermediate stage between the galaxies in the two other TIR luminosity bins (A/K=1.65).



Figure 11: Ratios of s determined from to those determined by [OII] emission for the average galaxies of the three subsamples plotted against . In the Herschel subsamples, [OII] emission underestimates the true by a factor of 4080. Increasing results in an increase in the ()/([OII]) ratio, presumably a consequence of the larger dust content in the average high- galaxy relative to its lower luminosity counterpart. Only random errors are used resulting in error bars the sizes of the symbols. Systematic errors related to conversions and extinction are quite large, but are unnecessary to include here as the comparison is relative and any systematic offset would affect each ratio in the same manner. The ()/([OII]) ratio is used to determine the nebular extinction, , for the average galaxy in each subsample. These values are given at the top of the plot along with the definitions of each subsample.

The color of the stellar continuum in the differenced coadded spectra also allows for insights on the average luminosity-weighted stellar population. Both the spectra which are differenced with respect to the low- sample appear quite blue, which, in the absence of dust effects, also implies a larger fractional contribution of late-type stars in the average low- galaxy relative to that of the high- and moderate- samples. Of course, given that our sample is selected, by definition, as galaxies containing a significant amount of obscuring dust, the interpretation of this color difference is not as straightforward. If the strength of the stellar extinction in the average low- galaxy (i.e., ) is higher than that of the high- or moderate- samples, the blueness observed in the differenced spectra could be attributable solely to dust. While it has been shown that the opposite trend is true, i.e., that galaxies with higher TIR luminosities have, on average, a higher value of , and by proxy, a higher value of , than those galaxies of lower TIR luminosities (Kocevski et al. 2011a), such statistical arguments may not hold when applied to the small number of galaxies in each spectral bin ().

Instead, we calculated for each coadded spectrum an average nebular extinction value, , which allows us to completely remove the effects of extinction when discussing average quantities of the various subsamples. In order to derive an for each sample it is necessary to compare star formation rates coming from two different indicators, one that is significantly affected by dust and one that is largely dust-independent. Because each of our coadded spectra has the [OII] emission feature within its usable range, the [OII]-derived can be used in concert with the derived from the average TIR luminosity to determine for each of our subsamples. A calculation involving this method has two inherent assumptions: that the derived from the TIR luminosity is, for our sample, a good representation of the true amount of star formation in a given galaxy, and the star-formation properties of the galaxy as a whole is well represented by the portion of the galaxy covered by the slit. As stated in §3.4, the TIR luminosity dominates the SED of these galaxies, meaning that we can safely ignore the first assumption. The second assumption is a bit more complex. In a recent study of galaxy clusters at , Oemler et al. (2013) found little difference between the strength of the EW(H) and EW([OII]) for galaxies in their sample when slit apertures of two different sizes were used. The covering fraction of these two slit sizes was disparate enough for them to conclude that a spectrum observed with a slit is sufficiently representative of the galaxy as a whole for a large fraction of their sample. The sample we consider here is bounded on the low-redshift end by the high-redshift galaxies in the sample of Oemler et al. (2013). As such, the covering fraction for the galaxies in our sample, given that the VVDS galaxies were observed with a , is equivalent or higher than that of Oemler et al. (2013) and approaches unity at . Thus, we ignore all considerations of slit effects for the remainder of this discussion and equate the properties observed in each coadded spectrum with the properties of the galaxies as a whole.

Because of the way the coaddition process was performed, absolute flux calibration is not preserved. Thus, rather than using an average [OII] line flux it is necessary to rely on the EW([OII]) value for each bin to derive the dust-dependent . For each bin, the average [OII]-derived was calculated by


where is the mean rest-frame absolute magnitude of a given bin, chosen because it provides a fair sampling of the rest-frame continuum emission surrounding the [OII] feature, is the effective wavelength of the filter curve cm, and c is the speed of light in cm s. The constant of proportionality is adopted from a combination of 4, where is the luminosity distance to each source, 10 pc in units of cm for the case of an absolute magnitude, conversion of the from Å to cm, and the [OII] formula of Kewley et al. (2004). The negative sign is a consequence of the convention adopted in this paper. The IMF is assumed to be Salpeter, the same that was used for the TIR derived in §3.4. The average nebular extinction, , for each coadd is then


where TIR is the mean TIR luminosity in a given bin and is the Calzetti et al. (2000) reddening curve evaluated at the wavelength of [OII]. The ratio of the indicators as a function of for the three bins are given in Fig. 11 and the values for each bin are given in Table 5 and Fig. 11. The values found for the first two subsamples are similar to those found in Choi et al. (2006), who applied almost identical methodologies to a near- (2.2m) and mid- (24m) infrared selected spectroscopic sample of ULIRGs spanning a slightly lower redshift range than the current sample. The mean -derived of the high- subsample is prohibitvely high to make a similar comparison. However, the value found for the high- subsample is similar to those values derived for a Herschel/PACS-selected ULIRG sample with roughly equivalent median and redshifts (Buat et al. 2011), albeit with an extinction law which is slightly modified from the one chosen here.



Figure 12: Measurements of EW(H) and of the coadded spectra of the high- (circles), moderate- (diamonds), and low- subsamples corrected (filled light gray symbols) and uncorrected (filled dark gray symbols) for the effects of extinction. These measurements are plotted against the backdrop of equivalently measured quantities from CB07 synthetic spectra for various time steps using various models. For models with secondary bursts, these bursts occur 2 Gyr after the initial burst. No extinction is applied to the CB07 models and the models employ a Salpeter (1955) IMF, identical to the IMF used for all calculations. The decreasing starburst strength relative to the already established stellar mass for galaxies of decreasing suggest that the LIRG population is experiencing rejuvenated star formation at moderate levels while the average ULIRG is undergoing either its initial starburst or an extremely vigorous (i.e., 100% by mass) rejuvenated star-formation event.

Mean Stellar and Starburst Ages as a Function of TIR Luminosity

As stated earlier, for the trend observed in the differenced continua observed in the previous section to be, at least in part, due to differential dust effects rather than to the difference in the mean stellar populations of each subsample, it would be necessary for the average low- galaxy to have a higher extinction than that of the high- subsample. What we found in the previous section is the opposite, namely that the nebular extinction of the average low- galaxy () is significantly lower than that of the average high- galaxy (). The average moderate- galaxy again falls in between that of the two other subsamples, with a . These results are consistent with those of Kocevski et al. (2011a), in which it was observed that, for Spitzer/MIPS selected galaxies in a large-scale structure at , the nebular extinction increased with increasing . These results are also broadly consistent with those found in the field by Roseboom et al. (2012) in a spectroscopic sample of MIPS- and SPIRE-selected galaxies at redshifts similar to the ones probed in our sample (). However, the difference in nebular extinction between the subsamples does not necessarily translate in a straightforward manner to a difference in the extinction of the stellar continuum. We assume here that the nebular extinction is related to the extinction of the stellar continuum, the quantity of interest in this discussion, via (Calzetti 1997; Calzetti et al. 2000). In principle it is possible or even likely that a variety of dust geometries for a given galaxy would change the constant of proportionality relating the two quantities, such that we could not draw a definite conclusion based on only nebular extinction values. However, in this analysis we are averaging over a large population of galaxies in each subsample and thus mitigating the effects of varying dust geometries when comparing the various subsamples. Given this and given the small variance in the constant of proportionality observed for individual galaxies in Calzetti et al. (2000), it is likely that the variation in the constant of proportionality is minimal for the average galaxy of the three subsamples. Therefore, the intrinsic properties of the stellar population of the mean galaxy of each subsample must be the cause of the observed color difference in the various subsample, as an extinction correction of the coadded spectra would only serve to increase the blue color of the differenced spectra. Thus, there remains two main effects that could cause a increasingly redder stellar continuum with decreasing TIR luminosity: metallicity and age of the mean stellar populations. While differing metallicities between the various subsamples could account for some of the observed color difference between the various subsamples, the effect of metallicity on the color of galaxies over the wavelength range considered in the coadded spectra is too weak to explain the observed difference (see, e.g., Cooper et al. 2008; Lemaux et al. 2012). Thus, these results, along with those obtained earlier on the analysis of the spectral lines and the A/K ratios, strongly implies that the mean luminosity-weighted age of the stellar populations of an average (U)LIRG decreases with increasing TIR luminosity.

This result has one of two implications for the galaxies studied here. One way to induce this correlation is if the average galaxy in the low- subsample is observed later in the development of the starburst than the average galaxy in either of the other subsamples, thereby decreasing the contribution of A-type stars to the observed spectral features and stellar continuum color. An alternative scenario is to appeal to a well-established older stellar population in the average low- galaxy, and to a lesser extent in the average moderate- galaxy, embedded into the galaxy from star-formation events that occurred prior to the onset of the currently observed burst. In the diagnostics used thus far, the two scenarios appear equally probable, as the two scenarios would manifest themselves in an identical manner. In order to disentangle these scenarios we revisit the EW(H) and measurements made of the coadded spectra of each subsample §4.4. To begin, in Fig. 12 we plot the measured EW(H) and values in dark filled symbols against the backdrop of measurements made of solar-metallicity synthetic spectra from CB07 using a variety of different models using the method described in Lemaux et al. (2012). Unlike the method of Lemaux et al. (2012) we do not, however, redden the models here as noted by the value in the legend of Fig. 12.

The values of EW(H) and plotted in dark symbols (and given outside the parentheses in Table 5) represent the measurement taken directly from each coadded spectrum. In order to make these measurements meaningful in the context of the CB07 models plotted in Fig. 12 corrections need to be made to both EW(H) and . The correction to EW(H) is straightforward. Because EW measurements are generally insensitive to the effects of dust (see discussion in Lemaux et al. 2010), it is only necessary to make a correction for the infill to the stellar H absorption from nebular emission. This is necessary because the CB07 models are only stellar and do not contain a component which models emission from HII regions. This correction was based on the method of Kocevski et al. (2011a) and uses the EW([OII]) of each coadded spectrum. The emission infill corrected EW(H) values are given in parentheses in Table 5. The correction to is also straightforward now that we have already derived the average extinction for each subsample. The average values for each coadd were used to apply a Calzetti reddening curve to each coadded spectrum and a de-reddened measurement was made of each coadded spectrum. These values are given inside the parentheses in Table 5. The corrected values of EW(H) and for each subsample should be completely independent of dust, nebular emission, and, given that metallicity effects are not strong in the age range considered here (Kauffmann et al. 2003), metallicity effects, and can thus be compared directly to the BC07 models to determine the age since the onset of the starburst for each subsample. The corrected values are plotted in Fig. 12 as lightly filled symbols.


[angle=90,width=1.0]Herschelvsnon.lumbinned.NUV-r.Mr.restframe.singlepanel.0.5z2.eps \includegraphics[angle=90,width=1.0]Herschelvsnon.lumbinned.NUV-r.stellarmass.restframe.singlepanel.0.5z2.eps

Figure 13: Similar to Fig. 7 except that only a single redshift range is plotted (), a redshift range which comprises the bulk of our sample. SPIRE-detected galaxies are broken into high- (red filled circles), moderate- (green filled squares), and low- (blue filled diamonds) subsamples. Again in both the CMD (left) and CSMD (right), galaxies which went undetected or were not detected at high enough significance in the SPIRE imaging are shown as small black points. Histograms are normalized in the same way as Fig. 7 and arrows again represent the movement in each phase space when an stellar continuum extinction of is applied. Only those galaxies with and stellar masses larger than the completeness limit () are plotted. While the LIRG subsamples show a high level of similarity in their colors and absolute magnitudes (left), the ULIRG population appears both bluer and brighter optically. The range of stellar masses spanned by the ULIRGs at any given color is generally larger than that of the LIRG samples (right) suggesting that the ULIRG population is inhomogenous, comprised of nascent galaxies undergoing their initial starbursts and older galaxies undergoing intense rejuvinated starburst events.

By comparing the corrected EW(H) and values of the three subsamples to the closest model track in Fig. 12 the first scenario is almost entirely ruled out, as the age of the starburst in the average galaxy of each of the three subsamples is nearly the same (between 50-140 Myr). Indeed, the populations that showed the most marked difference in the other diagnostics, the high- and low- subsamples, appear here to be nearly identical in terms of the starburst age of their average galaxies. The major difference between the three subsamples appears to be, rather, the strength of the burst. While this is not unexpected given that the populations were binned in terms of , which is directly proportional to the , the strength of the burst probed here is not the absolute strength (which is probed by ), but rather the strength of the burst relative to the already established stellar population in that galaxy. The results in Fig. 12 appear to directly confirm the second scenario proposed above. Both the average galaxy of the moderate- and low- are consistent with galaxies having already undergone significant star formation in the past, which are now undergoing a moderate (by mass) burst of star formation. Conversely, the average galaxy in the high- is inconsistent with either of these models. Instead, the measurements of the high- spectrum indicate a population of incipient galaxies undergoing their initial phase of star formation. However, caution must be used when interpreting this population, as a secondary-burst model that doubles the mass (a 100% burst in the language of Fig. 12) and occurs several Gyr after the initial star-formation event is essentially indistinguishable in this phase space. These two scenarios will be discussed further in the next section. Combining the results of all the spectral analysis here and in the previous sections clearly indicates that, amongst the full SPIRE sample, galaxies with progressively lower have a progressively higher abundance of established older stellar populations relative to the strength of the starburst. In other words, the LIRGs55 in the full SPIRE sample appear to be galaxies that have a high fractional abundance of older stellar populations born of previous star-formation events and are currently undergoing moderate rejuvenated star formation. The ULIRGs in the full SPIRE sample appear to be a decidedly different population, being largely comprised of galaxies that are either undergoing a primordial starburst or a massive secondary starburst event which is considerably increasing the amount of stellar mass in these galaxies. These results will have important consequences for the interpretation of SPIRE-detected galaxies hosting AGN discussed later in this paper.

Photometric and Stellar Mass Properties

Having established the nature of the LIRGs in the full SPIRE sample, we now turn from the spectral properties of the galaxies of different TIR luminosities to their color, magnitude, and stellar mass properties. This is done both to contextualize the results of the previous section and to investigate more fully the nature of the SPIRE-detected ULIRGs. As was mentioned in the beginning of §4.4, neither the color, absolute magnitude, or stellar mass distributions of the galaxies with spectra in each subsample used in the previous section are appreciably different than their parent population such that the conclusions drawn from the spectral analysis in the previous section can be directly combined with conclusions drawn from the full subsamples. With this in mind, plotted in the left panel of Fig. 13 is the CMD of the full galaxy populations of each of the subsamples along with SPIRE-undetected galaxies in the redshift range that satisfy the criteria given in §4.1. While the color histograms of all three subsamples peak at roughly the same color, the galaxies in the full LIRG subsample have significantly redder colors than the SPIRE-detected ULIRGs. This is confirmed by a KS test, which rejects the hypothesis that the galaxies of either of the LIRG subsamples are drawn from a parent galaxy population with the same color properties as the ULIRG subsample at % CL. As was established in the previous section, the redder colors of the LIRGs relative to that of the ULIRGs cannot be the result of differential extinction effects and is, in fact, only enhanced if the rest-frame colors are corrected by the of each subsample. A comparison of the rest-frame magnitudes of the three subsamples finds an even more marked difference between the populations. SPIRE-detected ULIRGs appear significantly brighter in the rest-frame band than SPIRE-detected LIRGs, with ULIRGs representing some of the brightest galaxies at any color in this redshift range. The magnitude distributions reflect this, as the magnitude distribution of the ULIRGs peaks 1.5 mags brighter than either of the LIRG populations, a difference that is also enhanced if the rest-frame magnitudes are corrected for the effects of extinction. A similar picture emerges from the CCD plotted in Fig. 14. Of the galaxies safely in the region of the diagram populated by galaxies with quiescent colors (i.e., the area where even large extinction effects could not move the galaxies out of the quiescent region), almost all are made up of low-luminosity LIRGs. The ULIRG galaxies that are in the quiescent region are easily moved out when the derived in the previous section is applied to the observed colors (i.e., one arrow length). These properties are broadly consistent with the results of the previous sections: the average stellar continuum of a LIRG in our sample has a larger fractional contribution from older stellar populations that that of a ULIRG. The blue colors and extremely bright rest-frame magnitudes of the ULIRGs show that this is a population largely undergoing a violent burst of star-formation in which its optical colors are completely dominated by the newly formed young stellar populations.

However, the ambiguity related to the ULIRG population revealed in the previous section still remains. In star-forming galaxies, the optical luminosity (i.e., absolute magnitude) is a confusing quantity to interpret, as a galaxy with a well-established stellar population undergoing a moderate strength starburst can appear equally luminous as a nascent galaxy undergoing a massive starburst event. While it is tempting to interpret those ULIRGs which are extremely blue, extremely (optically) luminous as being comprised of the latter population, it is necessary to examine the stellar masses of these galaxies in order to contextualize their luminosities. Plotted in the right panel of Fig. 13 is the CSMD of the three subsamples along with all SPIRE-undetected galaxies in the same redshift range. What is immediately noticeable is that the distributions of the three subsamples, which were extremely distinct in absolute magnitude space, have now become essentially indistinguishable in stellar mass space. A KS test of the various distributions confirms this, with only the low- and high- subsamples having even moderate levels of incompatibility (i.e., a rejection of the null hypothesis at the level).



Figure 14: Similar to Fig. 8 but plotting only those galaxies in the redshift range 0.52 broken into subsamples. The symbols are identical to those of Fig. 13. Again only those galaxies brighter than