A Spectroscopic Follow-up Program of Very Massive Galaxies at 3 < z < 4: Confirmation of Spectroscopic Redshifts, and a High Fraction of Powerful AGNs

X-shooter is a single-object, medium-resolution echelle spectrograph with simultaneous coverage of wavelength range 0.3–2.5 μm in three arms (UVB, VIS, NIR; D'Odorico et al. 2006). The observations of our targets were carried out in queue mode as part of the ESO program 087.A-0514 (PI: Brammer) in 2011 May, following an ABA'B' on-source dither pattern using the 11'' × 10 and 11'' × 09 slits for the UVB and VIS/NIR arms, respectively. This instrumental setup resulted in a spectral resolution of R = 4200, 8250, and 4000 for the NIR, VIS, and UVB arms, respectively. For calibration purposes, telluric and spectrophotometric standard stars were observed in the same setup as science observations. We refer the reader to Table 1 for the seeing conditions (FHWM, '') and exposure times for each object. The data reduction for the X-shooter observations used custom scripts based on the standard X-shooter reduction pipeline (Modigliani et al. 2010). The calibration steps (master darks, order prediction, flat fields, and the two-dimensional maps for later rectification of the spectra) were run with the default parameters in the pipeline (Goldoni 2011). After these five calibration steps, the echelle spectra were dark-subtracted, flat-fielded, and rectified, and the orders stitched (12, 15, and 16 for the UVB, VIS, and NIR arms, respectively). The sky was then subtracted using adjacent exposures. Standard star observations were reduced with the same calibration data as the science frames and used to correct for telluric absorption and detector response. A final spectrum was created by mean stacking all exposures. We refer to Geier et al. (2013) for a detailed description of reduction steps of X-shooter spectra.

Following Horne (1986), one-dimensional spectra were extracted by summing all adjacent lines (along the spatial direction) using weights corresponding to a Gaussian centered on the central row with a full width at half maximum equal to the slit width used in each observation. To correct for slit losses and obtain an absolute flux calibration, spectroscopic broad/medium-band fluxes were obtained by integrating over the corresponding filter curves, and a constant scaling was applied to each spectra individually. A binned, lower resolution spectrum with higher S/N was extracted for each spectra using optimal weighting, excluding parts of spectra contaminated by strong sky emission or strong atmospheric absorption. The resulting spectral resolutions of the binned X-shooter spectra were $R\approx 10\mbox{--}45,15\mbox{--}40$, and 25 − 50 for the UVB, VIS, and NIR arms, respectively, while the spectral resolutions of the binned NIRSPEC spectra were $R\approx 30\mbox{--}100$ and 30–300 for the H and K bands, respectively.

We found that five-sixths of our targets have significant (>3σ) MIPS 24 μm detections (four already studied in Marchesini et al. 2010), consistent with the high fraction of MIPS-detected sources in the sample of IRAC-selected massive galaxies at $z\gt 3.5$ over GOODS-North (Mancini et al. 2009). At the redshifts of our targets, the observed 24 μm probes the rest-frame ∼4.8–7.8 μm. Emission at these wavelengths can arise from warm/hot dust and polycyclic aromatic hydrocarbons (Draine & Li 2007), heated by dust-enshrouded star formation or AGNs.

Figure 6 shows the observed far-infrared (FIR) SEDs of our targets with significant MIPS 24 μm detections. We include the Herschel PACS 100 and 160 μm photometry from the PACS Evolutionary Probe (PEP) survey DR1 public data release⁹ (Lutz et al. 2011; Magnelli et al. 2013) and Herschel SPIRE 250, 350, and 500 μm photometry from the Herschel Multi-tiered Extragalactic Survey (HerMES) DR3/4 public data release¹⁰ (Roseboom et al. 2010, 2012; Oliver et al. 2012; Hurley et al. 2017). These fluxes are listed in Table 4. For C1-21316 (a.k.a. AzTEC 5), we also included the de-boosted photometry from SMA at 890 μm (${S}_{890\mu {\rm{m}}}=9.3\pm 1.3$ mJy; Younger et al. 2007), JCMT AzTEC at 1.1 mm (${S}_{1.1\mathrm{mm}}=6.5\pm 1.4$ mJy; Scott et al. 2008), ASTE AzTEC at 1.1 mm (${S}_{1.1\mathrm{mm}}=4.8\pm 1.1$ mJy; Aretxaga et al. 2011), and JCMT SCUBA2 at 450 and 850 μm (${S}_{450\mu {\rm{m}}}\,=25.35\pm 6.04$ mJy and ${S}_{850\mu {\rm{m}}}=11.42\pm 1.38$ mJy; Casey et al. 2013). Most of our sources are not detected in any of the Herschel bands. Only C1-21316, which is a submillimeter galaxy and a radio-loud AGN, is detected in SPIRE and at longer wavelengths; C1-2127 is detected in the blue and green SPIRE bands, but not at 500 μm. No source is detected in PACS despite all being robustly detected in the MIPS 24 μm band. Figure 6 shows the IR SEDs of the targeted sources, along with the 1, 2, and 3σ upper limits indicated for the bands without detections.

Zoom In Zoom Out Reset image size

Table 4.  FIR Flux Density Detections and Limits

ID	${S}_{\nu }$ (100 μm)	${S}_{\nu }$ (160 μm)	${S}_{\nu }$ (250 μm)	${S}_{\nu }$ (350 μm)	${S}_{\nu }$ (500 μm)
	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)
C1-2127	(1.8)	(3.4)	${13.4}_{-2.3}^{+2.2}$	${22.2}_{-4.7}^{+4.7}$	(10.9)
C1-15182	(1.7)	(3.4)	(3.6)	(2.1)	(2.9)
C1-19536	(1.7)	(3.4)	(2.7)	(3.6)	(5.1)
C1-19764	(1.7)	(3.4)	(2.7)	(3.6)	(5.1)
C1-21316	(2.0)	(3.4)	${29.1}_{-3.2}^{+2.5}$	${38.3}_{-3.4}^{+3.1}$	${31.6}_{-29.2}^{+6.6}$
A2-15753	(1.2)	(2.5)	(3.6)	(5.9)	(5.6)

Note. The Herschel PACS 100, 160 μm, SPIRE 250, 350, 500 μm detections and limits. Photometric points with detections are listed with the corresponding 1σ errors, whereas non-detections are represented with the 1σ error limits in parenthesis.

Download table as:  ASCIITypeset image

The 24 μm emission has been widely used in the literature to estimate total infrared luminosities (${L}_{\mathrm{IR}}={L}_{8\mbox{--}1000\mu {\rm{m}}};$ see, e.g., Papovich et al. 2007; Rigby et al. 2008), which can then be transformed into dust-enshrouded SFRs (Kennicutt 1998). In this section, we will show that this approach cannot be blindly adopted in ultra-massive galaxies at $z\gt 3$ given the almost ubiquitous presence of obscured AGNs in this population and the resulting AGN contamination of the 24 μm emission.

First, we derived dust-enshrouded SFRs from the MIPS 24 μm fluxes using the approach presented in Wuyts et al. (2008), which has become one of the most adopted approaches in the literature. Specifically, this method uses the mean log${L}_{\mathrm{IR},\alpha =1,\ldots ,2.5}$ of the infrared SEDs of star-forming galaxies provided by Dale & Helou (2002) to calculate ${\mathrm{SFR}}_{\mathrm{IR}}$. This method was validated out to $z\sim 3.5$ for lower mass (i.e., $\mathrm{log}({M}_{\star }/{M}_{\odot })\lt 11$) star-forming galaxies using Herschel data (Wuyts et al. 2011; Tomczak et al. 2016). The calculated total 8–1000 μm rest-frame luminosities and the corresponding ${\mathrm{SFR}}_{{IR}}$ adapted for a Kroupa (2001) IMF from Kennicutt (1998) are listed in Table 5, third and fourth columns. The corresponding errors are also listed, with and without the uncertainties due to random photometric redshift uncertainties when spectroscopic redshifts are not available. As shown in Table 5, the random uncertainties on ${L}_{\mathrm{IR}}$ and ${\mathrm{SFR}}_{\mathrm{IR}}$ as derived using the Wuyts et al. (2008) approach are dominated by the contribution from random photometric redshift errors. We note that a significant scatter of $\sim 0.2\mbox{--}0.3$ dex was found by Wuyts et al. (2011) and Tomczak et al. (2016) between the MIPS 24 μm derived SFRs and the SFRs derived including robust Herschel PACS detections. Therefore, we stress that the total error budget of ${L}_{\mathrm{IR}}$ and ${\mathrm{SFR}}_{\mathrm{IR}}$ listed in Table 5 as derived using the Wuyts et al. (2008) approach is dominated by this scatter and can be as large as a factor of a few. When the Wuyts et al. (2008) approach is adopted, the ${L}_{\mathrm{IR}}$ of our targets range from $\sim 5\times {10}^{12}{L}_{\odot }$ to $\sim 5\times {10}^{13}{L}_{\odot }$, typical of Ultra Luminous IR galaxies (ULIRGs) and Hyper Luminous IR galaxies (HLIRGs; Murphy et al. 2011), corresponding to SFR ∼ 500–5000 M_⊙ yr⁻¹ of obscured star formation. We used the rest-frame 2800 Å luminosity ($\nu {L}_{2800\mathring{\rm{A}} }$), determined via the best-fit templates in EAZY (using the same methodology for deriving rest-frame colors described in Brammer et al. 2011), as a proxy for ${L}_{\mathrm{UV}}$, reflecting the contribution of unobscured star formation. We used the approach detailed in Bell et al. (2005) to convert $\nu {L}_{2800\mathring{\rm{A}} }$ to ${\mathrm{SFR}}_{\mathrm{UV}}$. The calculated ${L}_{\mathrm{UV}}$ and ${\mathrm{SFR}}_{\mathrm{UV}}$ are listed in Table 5.

Table 5.  Spitzer-24 μm Fluxes and the Derived ${L}_{\mathrm{IR}}$ and SFR

ID	S24	${L}_{\mathrm{IR}}$ a	${\mathrm{SFR}}_{\mathrm{IR}}$ a	${L}_{\mathrm{IR}}$ b	${\mathrm{SFR}}_{\mathrm{IR}}$ b	$f{(\mathrm{AGN})}_{\mathrm{bol}}$ c	$f{(\mathrm{AGN})}_{\mathrm{MIR}}$ c	${L}_{\mathrm{UV}}$	${\mathrm{SFR}}_{\mathrm{UV}}$
	(μJy)	(${10}^{12}{L}_{\odot }$)	(${M}_{\odot }$ yr−1)	(1012 ${L}_{\odot }$)	(${M}_{\odot }$ yr−1)	%	%	(109 ${L}_{\odot }$)	(${M}_{\odot }$ yr−1)
C1-2127	110.6 ± 10.7	$9.6\pm 0.9{(}_{-1.5}^{+2.0}$)	$1049\pm 101{(}_{-168}^{+220}$)	$7.2\pm 0.7{(}_{-0.8}^{+0.8}$)	$761\pm 291{(}_{-294}^{+295}$)	2.7±1.0	∼15	5.3	1.7
C1-15182	78.6 ± 7.3	10.3 ± 1.0	1119 ± 104	<2.6 ± 0.24	<259 ± 96	>9.2 ± 3.2	$\gt 50$	5.6	1.8
C1-19536	61.6 ± 6.7	4.7 ± 0.5	512 ± 56	<3.6 ± 0.4	<384 ± 148	$\gt \,4\pm 1$	$\gt 15$	4.8	1.6
C1-21316	177.1 ± 7.8	46.0 ± 1.2	4993 ± 64	9.1 ± 0.4–14.9 ± 0.7	900 ± 325–1516 ± 540	4 − 12	∼35–50	5.4	1.8
A2-15753	165.7 ± 4.2	$13.3\pm 0.6{(}_{-2.1}^{+2.5}$)	$1357\pm 34{(}_{-227}^{+276}$)	$\lt 5.3\pm 0.1{(}_{-0.3}^{+0.3}$)	$\lt 533\pm 188{(}_{-190}^{+190}$)	$\gt \,10\pm 4$	$\gt 40$	8.5	2.7

Notes. The Spitzer 24 μm fluxes and derived ${L}_{\mathrm{IR}}$ and ${\mathrm{SFR}}_{\mathrm{IR}}$ assuming

^athe mean Dale & Helou (2002) IR template ^bthe best-fit mean Dale et al. (2014) IR template in agreement with Herschel 3σ detection limits; these values are strictly upper limits to the ${L}_{\mathrm{IR}}$ and SFRs allowed by the observed FIR fluxes and limits ^clists the fraction of AGN luminosity contribution to the total bolometric (5–1100 μm) and MIR (5–20 μm) luminosity (values in paranthesis) for the best-fit mean Dale et al. (2014) templates; these values are lower limits constrained by FIR fluxes and limits. For C1-21316, the two values provided in columns 5 to 8 correspond to values derived assuming the two Dale et al. (2014) templates constrained by either the SCUBA 450 μm or the $\lambda \gt 800\,\mu {\rm{m}}$ fluxes (solid blue curves in Figure 6). The errors listed for ${L}_{{IR}}$ and ${\mathrm{SFR}}_{\mathrm{IR}}$ were computed using just the 24 μm photometric errors (values not in parentheses) and the combination of the 24 μm photometric errors and the photometric redshift errors (values in parentheses).

Download table as:  ASCIITypeset image

The SFRs derived by assuming that all the flux at the observed 24 μm is associated with dust-enshrouded star formation are ∼2–3 orders of magnitude larger than the SFRs estimated from SED modeling (although values for C1-2127 and A2-15753 are consistent within 1σ uncertainties). We can see in Figure 6 that the median pure-SF IR SED template derived from the library of Dale & Helou (2002), and adopted in the Wuyts et al. (2008) approach, is not an adequate fit to the observed detections and limits. For each object, we overplot the median Dale & Helou (2002) template scaled to the observed MIPS 24 μm detections as dashed red curves. It can be seen that a model in which the IR emission is due only to obscured star formation cannot reproduce the observed IR SEDs. Specifically, the detections and upper limits in Herschel bands put strong constraints on the overall shape of the IR SEDs, rejecting the scenario in which the observed IR properties are due entirely to obscured star-formation activity, and pointing to the presence of a component with warmer dust temperature (arguably originating from the dusty torus of the AGN) not present in the average template from Dale & Helou (2002). This result was already presented in Marsan et al. (2015) for the Vega galaxy at z = 3.35, and we now show that it appears to be a common property of ultra-massive galaxies at $z\gt 3$.

To quantitatively assess the amount of AGN contribution to the IR SEDs, we used the IR template set of Dale et al. (2014), which includes fractional AGN contributions to the mid-infrared (MIR) radiation. We calculated the mean log(${L}_{\mathrm{IR}}$) for $\alpha =1$,..., 2.5 from this template set for varying AGN contributions. For visual ease in comparison, we only plot the templates with AGN contributions of 0%, 20%, 40%, 60%, 80%, and 100% scaled to the observed MIPS 24 μm photometry (dashed blue curves in Figure 6). The Dale et al. (2014) templates that minimize the ${\chi }^{2}$ to detections and are consistent with upper limits for each target are plotted as solid blue curves, with f(AGN)_MIR indicated in each panel. We note that the specified f(AGN)_MIR refers to the fraction of AGN light to the integrated mid-infrared (i.e., 5–20 μm) emission. A significant amount of AGN contribution to the MIR SEDs is required by all ultra-massive galaxies at $z\gt 3$ to explain the IR SEDs. If 3σ upper limits for the Herschel photometry are adopted, we find a lower limit on the fraction of AGNs to the MIR SEDs ranging from 15% to 50% depending on the galaxy. If the 1σ upper limits are adopted, we find a lower limit on the fraction of AGNs to the MIR SEDs ranging from 40% to 75%. Columns 5–8 of Table 5 list the values or upper limits of ${L}_{\mathrm{IR}}$ and ${\mathrm{SFR}}_{\mathrm{IR}}$, and the values or lower limits of ${f}_{\mathrm{Bol}}$(AGN) and ${f}_{\mathrm{MIR}}$(AGN) as obtained from adopting the best-fit templates from Dale et al. (2014). The upper limits were derived adopting the 3σ upper limits in the Herschel photometry. We stress that these newly derived ${L}_{\mathrm{IR}}$ and SFRs should be considered as upper limits, given that we only have lower limits on the AGN contribution to the IR. We note that the three galaxies with X-ray detections, namely A2-15753, C1-15182, and C1-19536 (see Section 4.4.3 below), are undetected in Herschel while being detected at several σ in the MIPS 24 μm band, f(AGN)_MIR > 40%(75%), 50%(80%), and 15%(60%), respectively, when adopting the 3σ (1σ) upper limits in Herschel. One target, C1-21316, which has f (AGN)_MIR ∼ 0.35–0.50 is detected in radio data, indicative of hosting a radio-loud AGN (Section 4.4.4).

We also used the templates from the observed FIR color-based library of Kirkpatrick et al. (2015) to investigate the AGN continuum contribution to the IR SEDs. In this study, libraries of empirical IR SED templates were created from a large sample of high-redshift ULIRGS with rest-frame mid-IR spectroscopy from the Spitzer Space Telescope Infrared Spectrograph (IRS). We display their color diagnostic in the bottom right panel of Figure 6, with the FIR colors of our targets calculated using photometry from Herschel SPIRE and Spitzer MIPS/IRAC plotted in solid colored circles (color scheme identical to Figures 2 and 5), with arrows indicating 2σ limits. Kirkpatrick et al. (2015) used mid-IR features to determine the strength of the AGN contribution to IR SEDs by jointly fitting the mid-IR spectrum of a prototypical starburst (M82) and a pure power law for the AGN component to the rest-frame mid-IR spectroscopy. FIR color-based templates were created by stacking the SEDs of galaxies that fall within each ${S}_{250}/{S}_{24}$ (observed) versus ${S}_{8}/{S}_{3.6}$ (observed) color box ($\langle {f}_{\mathrm{AGN}}\rangle$ of galaxies that fall into each color box is indicated in the lower right panel of Figure 6). The Kirkpatrick et al. (2015) color-based templates consistent with SPIRE, MIPS, and IRAC fluxes and detection limits of our targets are displayed in violet in Figure 6, scaled to the MIPS 24 μm flux of each galaxy. Using the appropriate color-based templates and the listed fraction of ${L}_{\mathrm{IR}}$ attributable to star formation for each template (Table 3 in Kirkpatrick et al. 2015), we re-calculated the dust-enshrouded SFRs. The SFRs calculated accounting for the AGN contribution to the FIR SEDs this way are broadly consistent with the observed UV-8 μm SED derived SFRs. Specifically, we calculated the dust-enshrouded SFRs for C1-15182, C1-19536, and A2-15753 to be $\lt 100\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, ∼100–200 M_⊙ yr⁻¹ for C1-21316, and ∼200 M_⊙ yr⁻¹ for C1-2127, in line with SED derived SFRs within uncertainties. It should be noted that the Kirkpatrick et al. (2015) library is more representative of our targets because the templates were created using a large sample of ULIRGs at higher redshifts ($z\sim 0.3\mbox{--}2.8$) than compared to Dale et al. (2014) templates, which used local galaxies.

To summarize, the investigation of the FIR SEDs of the massive galaxy sample reveals that none of them can be adequately described as purely star-forming systems. In fact, a MIR AGN contribution of at least $\sim 40 \%$, on average, is necessary to explain the marginal detections in the Herschel bands. It is promising that estimates of AGN contribution using the template sets of Dale et al. (2014) and Kirkpatrick et al. (2015) are consistent with each other. Although in both libraries the fractional AGN contribution is calculated over the rest-frame MIR spectrum, the quasar and starburst templates differ. Dale et al. (2014) uses a single quasar template (median PG quasar spectrum of Shi et al. 2013) to linearly mix with a suite of local normal star-forming galaxies spanning a range in α to characterize different heating levels ($\alpha =1$ and 2.5 for active and quiescent galaxies, respectively) over the 5–20 μm wavelength range. In Kirkpatrick et al. (2015), the mid-IR spectrum of the prototypical starburst M82 is used to represent the star-formation component, and the AGN component is characterized by a pure (free-slope) power law.

We use the luminosity-dependent O[iii] bolometric correction factor, ${C}_{{\rm{O}}{\rm{III}}}$ from Lamastra et al. (2009) for the two galaxies with [O iii] emission-line detections (namely C1-19536 and C1-15182), ${L}_{\mathrm{bol}}\approx 454\times {L}_{{\rm{O}}{\rm{III}}}$. Because both C1-19536 and C1-15182 have non-negligible best-fit ${A}_{{\rm{V}}}$ values from SED fitting, we assume that the observed [O iii] line luminosities are a lower limit to the intrinsic luminosity of [O iii]. Assuming all the observed [O iii] emission is due to AGN radiation, we find ${L}_{\mathrm{bol}}\gt 8\times {10}^{45}$ erg s⁻¹ and $\gt 3.5\times {10}^{45}$ erg s⁻¹ for C1-19536 and C1-15182, respectively, consistent with them hosting powerful hidden AGNs.

The stellar population parameters derived in Section 4.3 may be influenced by the presence of strong AGN continuum contamination to the observed photometry, most significantly biasing the stellar mass estimates. We investigated this by using the publicly available fully Bayesian Markov Chain Monte Carlo SED fitting algorithm of active galaxies, AGNfitter¹¹ (Calistro Rivera et al. 2016), to model the observed SEDs of the 24 μm detected targets. AGNfitter simultaneously fits the total active galaxy SED by decomposing it into four physically motivated components: the stellar population of the host galaxy, cold dust emission from star-forming regions, the accretion disk emission (Big Blue Bump, BBB) and hot dust emission surrounding the accretion disk (torus). We used the fiducial library of templates that are supplied with AGNfitter, which has been tested to perform well in classifying Type1 and Type 2 AGNs purely based on observed broadband photometry. We briefly describe the modeling assumptions (as adopted in Calistro Rivera et al. 2016) and the fiducial library of templates used in AGNfitter, and highlight important distinctions compared to the modeling assumptions employed in previous sections.

The BBB describing the thermal radiation emitted from the accretion disk surrounding the central supermassive black hole was modeled using a single modified template based on the composite spectrum of Sloan Digial Sky Survey Type 1 QSOs (Richards et al. 2006). Extinction to the emitted BBB spectrum was modeled assuming a Small Magellanic Cloud reddening law (Prevot et al. 1984). Emission from the dusty torus was modeled by using a library of SEDs with varying hydrogen column densities (${N}_{{\rm{H}}}=21\mbox{--}25$) created based on the set of empirical templates of Silva et al. (2004). The BC03 stellar population synthesis models are used to construct the library of templates to describe the host galaxy's stellar emission. Although the stellar population synthesis model assumptions are similar to those adopted in FAST as presented in Section 4.3, the explored parameter space is somewhat reduced. Specifically, a similar exponential declining SFH is adopted in AGNfitter, of the form $\mathrm{SFR}(t)\propto {e}^{-t/\tau }$, with a mildly more restricted range for the timescale τ from 0.1 to 10 Gyr. A model with constant SFR is also included in AGNfitter. The grid of stellar population ages is created from 0.2 Gyr to the age of the universe at the target's redshift in step of ${\rm{\Delta }}t\sim 0.1$ dex at fixed solar metallicity. AGNfitter adopts a Chabrier (2003) and the Calzetti et al. (2000) reddening law to model the extinction to the stellar light. The library of templates for modeling the emission from cold dust in star-forming regions includes the template set of Dale & Helou (2002) (also used in Section 4.4.1) and Chary & Elbaz (2001). For a more detailed description of the fitting algorithm, modeling assumptions, and template libraries, we refer the reader to Calistro Rivera et al. (2016).

Figure 7 illustrates the results when modeling the observed SEDs with the full parameter space range of the fiducial templates of AGNfitter and Table 6 lists the best-fit stellar masses converted to the Kroupa (2001). The fractions of AGN emission in the rest-frame UV–optical (0.1–1 μm) and MIR (1–30 μm) for each object, calculated by comparing the relative contribution of AGN luminosities in the considered wavelenghts are also listed in Table 6. The amount of AGN contribution to the rest-frame MIR estimated by AGNfitter is remarkably consistent with the values or lower limits derived in Section 4.4.1. The AGN contribution in the rest-frame UV–optical is found to range between a few percent (for A2-15753) to $\sim 40 \%$ (for C1-2127), whereas the AGN contribution in the rest-frame MIR is found to range between $\sim 10 \%$ (for C1-2127) to almost 100% (for C1-15182, C1-19536, and C1-21316).

Zoom In Zoom Out Reset image size

Table 6.  Best-fit Stellar Mass and AGN Contributions with AGNfitter

ID	$\mathrm{log}{M}_{* }$	${\rm{\Delta }}\mathrm{log}{M}_{* }$	${f}_{\mathrm{AGN}}$ (0.1–1 μm)	${f}_{\mathrm{AGN}}$ (1–30 μm)
C1-2127	${11.24}_{-0.24}^{+0.17}$	$-{0.11}_{-0.25}^{+0.17}$	${0.37}_{-0.07}^{+0.07}$	${0.09}_{-0.05}^{+0.18}$
C1-15182	${11.61}_{-0.09}^{+0.05}$	$+{0.05}_{-0.13}^{+0.17}$	${0.09}_{-0.07}^{+0.12}$	${0.96}_{-0.17}^{+0.04}$
C1-19536	${11.37}_{-0.04}^{+0.03}$	$-{0.08}_{-0.10}^{+0.10}$	${0.10}_{-0.07}^{+0.10}$	${0.98}_{-0.21}^{+0.02}$
C1-21316	${11.34}_{-0.68}^{+0.23}$	$-{0.24}_{-0.68}^{+0.57}$	${0.28}_{-0.23}^{+0.35}$	${0.39}_{-0.05}^{+0.04}$
A2-15753	${11.21}_{-0.03}^{+0.01}$	$+{0.04}_{-0.13}^{+0.15}$	${0.04}_{-0.02}^{+0.03}$	${0.99}_{-0.01}^{+0.01}$

Note. Results obtained from the modeling of the binned UV–NIR spectra, NMBS UV-8 μm photometry, and FIR detections and upper limits with AGNfitter. The best-fit stellar masses were adjusted to the Kroupa (2001) IMF in order to compare with results listed in Table 3. ${\rm{\Delta }}{M}_{* }$ lists the effect on the derived stellar mass when including AGN templates in the modeling of the panchromatic SEDs (${\rm{\Delta }}\mathrm{log}{M}_{* }=\mathrm{log}{({M}_{* })}_{(\mathrm{AGNfitter})}-\mathrm{log}{({M}_{* })}_{(\mathrm{FAST})}$) at fixed metallicity ($Z={Z}_{\odot }$). f_AGN is calculated as the relative contribution of the AGN templates in the considered wavelengths (i.e., the Big Blue Bump emission describing the accretion disk at 0.1–1 μm and the hot dusty torus emission at 1–30 μm).

Download table as:  ASCIITypeset image

The third column of Table 6 lists the difference between the stellar mass estimated by AGNfitter and the stellar mass derived using FAST (i.e., assuming no AGN contamination to the observed SEDs). It is interesting to note that, while the naive expectation would be that not accounting for AGN contamination to the observed photometry would lead to over-estimating the stellar masses, two galaxies (C1-15182 and A2-15753) have AGNfitter derived stellar masses marginally larger (by $\sim 0.04\mbox{--}0.05$ dex) than when no AGN contamination is assumed. This appears to be caused by slightly older best-fit ages of the stellar population (and hence slightly larger mass-to-light ratios) when the SEDs are modeled accounting for the AGN contribution. We note, however, that these differences are smaller than (or comparable to) the 1σ error on the stellar mass. For the other three galaxies, the stellar masses derived by accounting for the AGN contamination are found to be smaller (by $\sim 0.1\mbox{--}0.3$ dex on average) than the stellar masses derived when no AGN contamination is assumed. The best-fit stellar masses derived assuming AGN contamination ranges from $1.5\times {10}^{11}\,{M}_{\odot }$ to $4\times {10}^{11}\,{M}_{\odot }$. All galaxies except C1-21316 are found to be more massive than 10¹¹ M_⊙ within 1σ when modeled by AGNfitter. For all galaxies but C1-21316, the systematic effect to the derived stellar masses from AGN contamination is found to be smaller than a factor of ∼2 at 1σ. This is because the rest-frame UV–optical SEDs of these galaxies are dominated by the stellar light despite the rest-frame MIR SEDs being dominated by the obscured AGN. For C1-21316, the stellar mass is found to be overestimated by a factor of ∼2 on average, although the 1σ uncertainty allows for solutions that are smaller by as much as a factor of ∼8. Finally, we note that the stellar masses of all studied galaxies derived from AGNfitter are consistent with the stellar masses derived using FAST within the 1σ uncertainties.

As previously done in Marchesini et al. (2010), we used the publicly available Chandra X-ray catalogs (Laird et al. 2009 and Elvis et al. 2009 for the AEGIS and COSMOS fields, respectively) and the improved redshifts to revise the X-ray luminosities of the three detected sources (C1-15182, C1-19536, and A2-15753). Assuming a power-law photon index Γ = 1.9 (Nandra & Pounds 1994), the X-ray luminosities ${L}_{2\mbox{--}10\mathrm{keV}}$ are (6.4 ± 1.7) × 10⁴⁵ for C1-15182, (16.6 ± 2.6) × 10⁴⁵ for C1-19536 and (9.6 ± 1.5) × 10⁴⁵ erg s⁻¹ for A2-15753 (redshift uncertainties are included in quoted errors). We used the empirically derived observed ${L}_{[{\rm{O}}{\rm{III}}]}\mbox{--}{L}_{2\mbox{--}10\mathrm{keV}}$ relation of local AGNs from Ueda et al. (2015) to investigate the expected X-ray flux of the sources with detected [O iii] emission (C1-19536 and C1-15182). The estimated ${L}_{2\mbox{--}10\mathrm{keV}}$ for these galaxies using the observed ${L}_{[{\rm{O}}{\rm{III}}]}$ are an order of magnitude less than the luminosities calculated using X-ray detections (using ${\rm{\Gamma }}=1.4$ lowers this factor to $\sim 6\mbox{--}8$). We calculated upper limits to the rest-frame ${L}_{2\mbox{--}10\mathrm{keV}}$ for the targets not detected (C1-2127, C1-19764, and C1-21316) using the 3σ detection limits with Γ = 1.9. The 3σ upper limits are in the range of ${L}_{2\mbox{--}10\mathrm{keV}}\sim 5\mbox{--}20\ \times {10}^{44}$ erg s⁻¹, comparable to only the very powerful X-ray AGN, and therefore we cannot conclude or rule out whether these galaxies harbor AGNs.

From the sample of galaxies presented here, only C1-21316 (AzTEC-5) is detected at 1.4 GHz (Schinnerer et al. 2010), with a flux of 0.126 ± 0.015 mJy. We calculated the rest-frame flux densities assuming the canonical value of the radio spectral index of $\alpha =0.8$ (Condon 1992). The calculated rest-frame flux density is ${L}_{1.4\mathrm{GHz}}=10.6\pm 1.3\times {10}^{27}$ W Hz⁻¹, above the threshold to be classified as a radio-loud AGN (log(${L}_{1.4}\mathrm{GHz})\,\gt 25$; Schinnerer et al. 2007). The radio spectral index is not observed to evolve significantly over cosmic time (Magnelli et al. 2015), and is consistent with values derived for different samples of galaxies (Ibar et al. 2009, 2010; Ivison et al. 2010; Bourne et al. 2011).

The fraction of AGNs (${f}_{\mathrm{AGN}}$) within a galaxy population represents a probe of the level of AGN activity and of the growth of supermassive black holes at the center of galaxies. The energetic output from AGNs is considered to be one of the dominant feedback processes that can regulate, and even halt altogether, the infalling of gas and star formation in massive galaxies. The measurement of the level of AGN activity in massive galaxies as a function of redshift is therefore an indirect probe of the impact of AGN feedback at a given cosmic time.

Figure 8 shows the evolution of ${f}_{\mathrm{AGN}}$ as a function of redshift for galaxies with stellar mass ${M}_{\star }\gt {10}^{11}$ M_⊙. In our sample of very massive galaxies at $3\lt z\lt 4$, AGNs appear almost ubiquitous, with a very large AGN fraction of ${f}_{\mathrm{AGN}}\approx 0.8\pm 0.2$ (in which we assumed that C1-2127 does not host an AGN). Figure 8 also shows the fraction of AGN in massive galaxies at $0.6\lt z\lt 3.0$ from Cowley et al. (2016), at $1.7\lt z\lt 2.7$ from Kriek et al. (2007), and in the local universe from Kauffmann et al. (2003). Our measurement is consistent with the very high AGN fraction of ${f}_{\mathrm{AGN}}={0.86}_{-0.22}^{+0.14}$ at $2.6\lt z\lt 3.0$ as measured by Cowley et al. (2016) using the zFOURGE data set. On the contrary, the AGN fraction in massive galaxies at $z\sim 0$ is much smaller, ${f}_{\mathrm{AGN}}\approx 0.3$, and it appears to remain relatively constant all the way out to $z\sim 2.5$. Figure 8 clearly shows a dramatic transition in the fraction of AGNs in galaxies with ${M}_{\star }\gt {10}^{11}$ M_⊙ happening at $z\sim 2.5$, i.e., when the universe was ∼3 Gyr old. We quantified the significance of the transition in the AGN fraction at $z\sim 2.5$ by assuming that the fraction of AGN is constant with redshift, and equal to the bi-weight mean of the observed fraction at $z\lt 2.5$ (i.e., ${f}_{\mathrm{AGN}}\approx 0.29$) and testing this hypothesis with the chi-squared statistics of the two points at $z\gt 2.5$. If each point at $z\gt 2.5$ is considered by itself, we find a probability $p\lt 0.0006$ of obtaining each $z\gt 2.5$ measurement at least this discrepant from the no evolution assumption. If the two $z\gt 2.5$ measurements of the AGN fraction are considered together, we find a probability $p\lt {10}^{-7}$ of obtaining both $z\gt 2.5$ points at least this discrepant from the no evolution assumption. At earlier times, most (all) supermassive black holes at the center of massive galaxies appear to be actively accreting, whereas at later times most (∼2/3) of them are found dormant. Although the small sample does not allow us to robustly quantify a potential enhancement of the AGN fraction in our spectroscopically confirmed sample of very massive galaxies at $z\gt 3$ due to [O iii] emission-line contamination to the stellar continuum. However, we note that the [O iii] line emissions detected in the spectra are found to contribute at most 15%, and more generally $\leqslant 5 \%$, to the observed K-band fluxes. The limited observed [O iii] emission-line contamination suggests that our sample is not significantly biased toward strong [O iii] emitters, i.e., toward massive galaxies hosting AGNs, and that the derived AGN fraction in our sample is not considerably enhanced. Future wide area surveys will be able to constrain more robustly the fraction of AGN in massive galaxies at $z\gt 2.5$, whereas deeper surveys will be able to investigate the evolution of the AGN fraction in lower-mass galaxies.

Zoom In Zoom Out Reset image size