The MOSDEF Survey: Metallicity Dependence of PAH Emission at High Redshift and Implications for 24 μm Inferred IR Luminosities and Star Formation Rates at z ∼ 2

During the MOSDEF survey, we obtained rest-frame optical spectra of ∼1500 galaxies with the MOSFIRE spectrograph on the Keck I telescope (McLean et al. 2012). The parent sample was selected based on H-band magnitude in three redshift ranges: z = 1.37–1.70, 2.09–2.61, and 2.95–3.80, down to H = 24.0, 24.5, and 25.0 mag, respectively. These redshift ranges were adopted to ensure coverage of the strong optical emission lines ([O ii], [O iii], ${\rm{H}}\beta$, ${\rm{H}}\alpha$) in the JHK bands. The MOSDEF survey was conducted in the five CANDELS fields: AEGIS, COSMOS, GOODS-N, GOODS-S, and UDS (Grogin et al. 2011; Koekemoer et al. 2011). These fields are also covered by the Hubble Space Telescope (HST)/WFC3 grism observations of the 3D-HST survey (Skelton et al. 2014; Momcheva et al. 2016). The details of the MOSDEF survey, observing strategy, and data reduction are reported in Kriek et al. (2015).

We limit our study to the first two redshift bins (1.37 ≤ z ≤ 1.70 and 2.09 ≤ z ≤ 2.61), where the MIPS 24 μm filter covers the rest-frame 7.7 μm PAH feature (Figure 1). Objects with active galactic nucleus (AGN) contamination are removed from our study, based on their X-ray emission, IRAC colors (using the Donley et al. 2012 criteria), and/or [N ii]/${\rm{H}}\alpha$ line ratios ([N ii]/${\rm{H}}\alpha$ > 0.5; Coil et al. 2015; Azadi et al. 2017).

Zoom In Zoom Out Reset image size

Emission line fluxes were estimated by fitting Gaussian functions to the line profiles. Uncertainties in the line fluxes were derived by perturbing the spectrum of each object according to its error spectrum and measuring the line fluxes in these perturbed spectra (see Kriek et al. 2015; Reddy et al. 2015). The ${\rm{H}}\alpha$ line and [N ii] doublet were fit with three Gaussian functions, and the [O ii] doublet was fit with two Gaussians. To account for the loss of flux outside of the spectroscopic slits, corrections were applied by normalizing the spectrum of a slit star in each observing mask to match the 3D-HST total photometric flux (Skelton et al. 2014). Additionally, we used the HST images of our resolved galaxy targets to estimate and correct for differential flux loss relative to that of the slit star (Kriek et al. 2015). ${\rm{H}}\alpha$ and ${\rm{H}}\beta$ fluxes were further corrected for underlying Balmer absorption as determined from the best-fit stellar population models to the broadband photometry (Section 2.3).

We use the [N ii]λ6585/${\rm{H}}\alpha$ (N2) and ([O iii]λ5008/${\rm{H}}\beta$)/([N ii]λ6585/${\rm{H}}\alpha$) (O3N2) indicators to derive gas-phase oxygen abundances (metallicities) based on the empirical calibrations of Pettini & Pagel (2004) (see Sanders et al. 2015; Shapley et al. 2015). As both N2 and O3N2 indicators include ratios of emission lines that are close in wavelength space, no dust correction is applied. Although there are concerns regarding biases in the metallicity indicators that use nitrogen (Shapley et al. 2015; Sanders et al. 2016), the N2 and O3N2 indicators still distinguish the lower- and higher-metallicity galaxies, which is sufficient for the purpose of this study. Furthermore, we use the [O iii]λλ4960, 5008/[O ii]λλ3727, 3730 ratio (${{\rm{O}}}_{32}$) as a proxy for ionization parameter. As [O iii]λ5008 has a higher signal-to-noise ratio (S/N) than [O iii]λ4960, we assume a fixed [O iii]λ5008/[O iii]λ4960 line ratio of 2.98 (Storey & Zeippen 2000) to calculate the sum of [O iii]λ5008 and [O iii]λ4960. We correct [O iii] and [O ii] lines for dust extinction by using the Balmer decrement (${\rm{H}}\alpha$/${\rm{H}}\beta$) and assuming the Cardelli et al. (1989) extinction curve (Reddy et al. 2015; Shivaei et al. 2015a). Where necessary, we calculate ${{\rm{O}}}_{32}$ gas-phase metallicity based on the calibrations of Jones et al. (2015).

Stellar masses and ages are derived by fitting rest-UV to near-IR photometry from 3D-HST (Skelton et al. 2014; Momcheva et al. 2016) with Bruzual & Charlot (2003) models through a minimum χ² method. The photometry is corrected for emission line contamination according to the MOSDEF spectra. For the SED fitting we assume a solar metallicity, a Chabrier (2003) IMF, and an exponentially rising star formation history (Reddy et al. 2015). The latter is assumed because it has been shown that rising star formation histories best reproduce the observed SFRs at z ∼ 2 (e.g., Wuyts et al. 2011a; Reddy et al. 2012b). As described in Reddy et al. (2012b), ages that are inferred from exponentially rising star formation histories are ambiguous and often larger than those derived from constant star formation histories. However, a majority (86%) of our galaxies have very large characteristic timescales (τ = 5000 Myr), which makes their star formation histories very similar to a constant one, along with well-defined ages.

We convert ${\rm{H}}\alpha$ luminosities to SFRs by adopting the Kennicutt (1998) relation, modified for the Chabrier (2003) IMF. The line luminosities are corrected for dust attenuation using the Balmer decrement and assuming the Cardelli et al. (1989) curve (Shivaei et al. 2015a). Hereafter, we refer to the dust-corrected ${\rm{H}}\alpha$ SFR as ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$.

We perform scaled point-spread function (PSF) photometry on the Spitzer/MIPS and Herschel/PACS images in COSMOS, GOODS-N, GOODS-S, and AEGIS fields from multiple observing programs (PI: M. Dickinson; Dickinson & FIDEL Team 2007; Elbaz et al. 2011; Magnelli et al. 2013). The accuracy of the measured fluxes is determined by adding simulated sources to the science images and recovering their fluxes in the same way as for real sources. Details of the photometry and simulations are described in Shivaei et al. (2016). In brief, we use the IRAC and MIPS priors down to a 3σ limit for MIPS and PACS images, respectively. A 40 × 40 pixel subimage is made around each target. In each subimage, the PSFs are scaled to match the prior objects and the target simultaneously. In this process, a covariance matrix is defined to determine the robustness of the fit, which is affected by confusion with nearby sources as described in Reddy et al. (2010). Based on the total covariance value of the MIPS 24 μm fits, we remove objects that have nearby companions from our analysis. To determine the noise and background level in each subimage, we fit PSFs to 20 random positions, ensuring that these positions are more than 1 FWHM away from any real sources. The average and standard deviation of the background fluxes are adopted as the background level and noise in the image, respectively. The typical measurement uncertainty for 24 μm detected objects is σ₂₄ ∼ 15%. The majority (∼80%) of these 24 μm detected objects are not detected in 100 and 160 μm images. Out of 406 non-AGN objects without any nearby sources (see above), 128 of them (32%) are detected at 24 μm with S/N > 3, and only 13 (3%) are detected with S/N > 3 at both 24 and 100 μm.

Due to the relatively shallow depth of the available far-IR data, we stack the PACS images to obtain a sufficiently high S/N, so that we can measure average IR luminosities. We also stack the Spitzer/MIPS 24 μm images because of the high fraction of 24 μm undetected objects (S/N < 3). In this analysis, all the stacks include both the detected and undetected objects. The stacking technique is described below.

We extract 40 × 40 pixel subimages centered on each target and subtract all prior sources brighter than the detection limit of S/N = 3, using the scaled PSF method described in Section 2.4. We then combine these residual images by making an inverse-SFR-weighted average stack as follows. We normalize each 24 μm image by its ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ and divide the sum of the normalized images by the sum of the weights (i.e., 1/SFR). We adopt inverse-SFR-weighted average stacks because ultimately we are interested in investigating the average of f₂₄/SFR and ${L}_{7.7}$/SFR ratios (Section 3.1). As we do not have direct far-IR detections for a majority of our galaxies, we cannot construct inverse-${L}_{\mathrm{IR}}$-weighted average stacks. As a result, we use a 3σ-clipped mean or a median stack of 24, 100, and 160 μm images for the f₂₄/${L}_{\mathrm{IR}}$ and ${L}_{7.7}$/${L}_{\mathrm{IR}}$ analyses (Section 3.2). A comparison of these different stacking methods is presented in the Appendix.

We calculate the stacked flux by performing aperture photometry on the stacked image. Based on the measured S/N of various aperture sizes, we adopt a 4 pixel radius aperture. The amount of light falling outside of the 4 pixel aperture is calculated using the PSF, and the fluxes are corrected accordingly. Similar to the PSF-fitting technique, we measure background fluxes in 20 random positions to determine the noise (i.e., the uncertainty in the stacked 24 μm flux). There is little variation in the aperture corrections from field to field (≲5%), which is negligible compared to the measurement errors.

The k-correction required for converting 24 μm flux densities to rest-frame 7.7 μm luminosities is computed using the Chary & Elbaz (2001, hereafter CE01) templates. First, we shift the templates to the redshift of each object (or median/weighted-average redshift of the stack), and then we convolve the models with the 24 μm filter transmission curve. The best-fit model is determined through a least-${\chi }^{2}$ method, and $\nu {L}_{\nu }(7.7\,\mu {\rm{m}})$ (${L}_{7.7}$, in units of ${L}_{\odot }$) is extracted from the best-fit model.

To quantify the systematic biases from using different IR templates, we measure ${L}_{7.7}$ from the best-fit Dale & Helou (2002) and Rieke et al. (2009) models and compare them with those of CE01 (Figure 2). Adopting the Dale & Helou (2002) and Rieke et al. (2009) models would systematically increase the calculated ${L}_{7.7}$ by ∼0.06 and ∼0.10 dex, respectively. For a small fraction (4%) of galaxies, the Rieke et al. (2009) ${L}_{7.7}$ is larger than the CE01 ${L}_{7.7}$ by >0.3 dex. However, these galaxies do not have systematically different masses or metallicities compared to that of the rest of the sample. Therefore, their ${L}_{7.7}$ offsets do not affect the mass and metallicity trends in this study.

Zoom In Zoom Out Reset image size

Total IR luminosities (${L}_{\mathrm{IR}}$, integrated from 8 to 1000 μm) are measured from the best-fit CE01 models to PACS 100 and 160 μm stacks. We measure ${L}_{\mathrm{IR}}$ errors by perturbing the 100 and 160 μm flux densities of each stack by their measurement uncertainties 10,000 times and calculating the 68% confidence intervals from these realizations. Furthermore, we verify that the best-fit models of Dale & Helou (2002) alter the inferred IR luminosities by only ∼5%, which is negligible compared to the typical ${L}_{\mathrm{IR}}$ measurement uncertainties of ∼10%–20%.

Due to the relatively shallow depth and poor spatial resolution of the far-IR data, obtaining robust individual far-IR measurements for a representative number of galaxies in our sample is not possible. Instead, we rely on our individual measurements of ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ (Section 2.3). It was shown in Shivaei et al. (2016) that ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ accurately traces SFRs up to $\sim 300\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, when compared with those from the IR.

We calculate inverse-${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$-weighted average stacks of 24 μm images (see Section 2.5) in four bins of metallicity with a roughly equal number of galaxies in each bin. The stacked 24 μm flux density is converted to ${L}_{7.7}$ using the CE01 IR templates (Section 2.6) and then divided by the weighted-average ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$. Dividing the inverse-${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$-weighted average ${L}_{7.7}$ to weighted average ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ is mathematically equivalent to calculating an average of the L_7.7/${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ ratios in each bin (Appendix). Properties of the stacks are listed in Table 1. Figures 3(a) and (b) show the L_7.7/${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ ratio as a function of N2 and O3N2 metallicities⁷ (Section 2.2). The uncertainties in L_7.7/SFR ratios are calculated by adding the L_7.7 and SFR fractional uncertainties in quadrature. There is a clear increase of L_7.7/SFR with increasing metallicity. Considering the small fraction of 24 μm detected objects in the two lowest-metallicity bins, it is expected that the average L_7.7/${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ ratios, which include all the detected and undetected 24 μm objects in the corresponding bins, will lie below the few individually 24 μm detected sources. We consider separately the trends for the full redshift range, and for z > 2 alone. The L_7.7/SFR ratio increases by a factor of ∼10 (∼15) from the lowest- to the highest-metallicity bin in stacks of galaxies at 1.37 ≤ z ≤ 2.61 (2.0 ≤ z ≤ 2.61).

Zoom In Zoom Out Reset image size

Table 1.  Properties of L_7.7/${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ Stacks

Parameter	Parameter Range	N	$\tilde{z}$	$\langle {f}_{24}\rangle \ (\mu \mathrm{Jy})$	$\langle {L}_{7.7}\rangle$ $({10}^{9}\,{L}_{\odot }$)	$\langle {\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }\rangle \ ({M}_{\odot }\,{\mathrm{yr}}^{-1})$	$\langle {L}_{7.7}/{\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }\rangle ({10}^{8}\,{L}_{\odot }/{M}_{\odot }\,{\mathrm{yr}}^{-1})$
$12+\mathrm{log}{({\rm{O}}/{\rm{H}})}_{{\rm{N}}2}$	8.04–8.37	45	2.22	16.4 ± 0.8	31.9 ± 1.5	37 ± 3	6 ± 1
(${N}_{\mathrm{tot}}=185$)	8.38–8.48	47	2.23	13.8 ± 0.8	27.2 ± 1.5	30 ± 4	5 ± 2
	8.48–8.56	47	2.21	39.9 ± 0.9	74.6 ± 1.6	37 ± 4	29 ± 11
	8.56–8.73	46	2.17	52.0 ± 0.9	90.2 ± 1.5	29 ± 4	41 ± 15
$12+\mathrm{log}{({\rm{O}}/{\rm{H}})}_{{\rm{O}}3{\rm{N}}2}$	8.03–8.26	42	2.22	10.8 ± 0.7	21.2 ± 1.3	35 ± 2	2.9 ± 0.7
(${N}_{\mathrm{tot}}=171$)	8.27–8.35	43	2.21	18.8 ± 0.8	35.9 ± 1.5	30 ± 5	7 ± 3
	8.35–8.44	43	2.26	28.4 ± 1.0	57.5 ± 2.1	31 ± 4	30 ± 9
	8.44–8.74	43	2.13	52.0 ± 1.0	86.4 ± 1.6	39 ± 5	35 ± 12
${{\rm{O}}}_{32}$	0.23–1.46	42	2.29	34.4 ± 1.1	72.7 ± 2.3	45 ± 5	26 ± 8
(${N}_{\mathrm{tot}}=170$)	1.46–2.21	43	2.29	17.3 ± 0.7	36.6 ± 1.5	26 ± 4	11 ± 3
	2.22–3.28	43	2.27	4.7 ± 0.7	9.9 ± 1.5	26 ± 3	4.8a
	3.28–13.52	42	2.29	4.1 ± 0.6	8.9 ± 1.2	20 ± 2	4.4a
$\mathrm{log}({M}_{* }/{M}_{\odot }$)	9.00–9.84	98	2.22	0.9 ± 0.4	2.5a	9 ± 1	7a
(${N}_{\mathrm{tot}}=296$)	9.84–10.04	50	2.29	4.9 ± 0.6	10.8 ± 1	19 ± 2	11 ± 2
	10.04–10.28	50	2.23	16.0 ± 0.7	10.8 ± 1.3	31.3 ± 1	25 ± 4
	10.29–10.52	50	2.27	30.3 ± 0.8	31.3 ± 1.3	62.2 ± 2	32 ± 5
	10.53–11.37	48	2.24	69.5 ± 1.0	62.2 ± 1.8	130.9 ± 2	22 ± 3
Age (Myr)	71–453	49	2.27	6.5 ± 0.5	13.8 ± 1.1	30 ± 3	6 ± 1
(${N}_{\mathrm{tot}}=296$)	453–570	50	2.23	9.9 ± 0.6	19.4 ± 1.3	24 ± 2	10 ± 2
	570–904	50	2.29	12.6 ± 0.5	26.9 ± 1.1	17 ± 2	9 ± 3
	904–1278	50	2.16	18.1 ± 0.8	32.4 ± 1.4	24 ± 3	27 ± 11
	1278–2600	49	2.18	22.3 ± 1.0	41.0 ± 1.8	20 ± 3	22 ± 9
	2600–4250	48	2.24	34.9 ± 1.3	67.9 ± 2.5	29 ± 5	29 ± 12

Notes. Entries show properties of the stacks shown with orange symbols (1.37 ≤ z ≤ 2.60) in Figures 3 and 6. Each entry includes the parameter range, number of objects in each stack, median redshift, 24 μm stacked flux density and its measurement uncertainty, rest-frame 7.7 μm luminosity and its measurement uncertainty, 3σ-clipped mean ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ and its error of the mean, and the average of L_7.7/${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ ratios and its error. The latter should not be confused with the ratio of averages (refer to the text and the Appendix).

^aThese objects are undetected in the 24 μm image stacks (before the aperture correction), and the values are 3σ upper limits.

Download table as:  ASCIITypeset image

We investigate the potential incompleteness of our samples due to the requirement of detecting the [N ii] and [O iii] lines for N2 and O3N2 metallicities. The fractions of galaxies detected in [N ii] ([O iii]) for the four bins of stellar mass indicated in Table 2 are 0.28, 0.56, 0.72, and 0.93 (0.28, 0.55, 0.68, and 0.83) proceeding from low to high stellar masses. As expected, more galaxies with undetected emission lines are missed from the N2 and O3N2 samples as we go to lower masses. However, in each mass bin, the distributions of the 24 μm S/N and f₂₄ of the detected-line subsamples are very similar to that of the full sample (including detected and undetected [N ii] and [O iii] lines). Additionally, in Section 5.1 we show the L_7.7 and ${L}_{\mathrm{IR}}$ stacks in bins of stellar mass regardless of the emission-line detection. Both the L_7.7/SFR and ${L}_{7.7}/{L}_{\mathrm{IR}}$ ratios are suppressed at low masses, which correspond to low metallicities according to the mass–metallicity relation (Tremonti et al. 2004). We conclude that the sample incompleteness at low N2 and O3N2 metallicities is unlikely to be the dominant factor in driving the L_7.7/SFR and ${L}_{7.7}/{L}_{\mathrm{IR}}$ trends with metallicity.

Table 2.  Properties of ${L}_{7.7}/{L}_{\mathrm{IR}}$ Stacks

Parameter	Parameter Range	N	$\tilde{z}$	$\langle {f}_{24}\rangle \ (\mu \mathrm{Jy})$	$\langle {f}_{100}\rangle (\mu \mathrm{Jy})$	$\langle {f}_{160}\rangle (\mu \mathrm{Jy})$	$\langle {L}_{7.7}\rangle \ ({10}^{9}\,{L}_{\odot })$	$\langle {L}_{\mathrm{IR}}\rangle \ ({10}^{10}\,{L}_{\odot })$ a	$\langle {L}_{7.7}\rangle /\langle {L}_{\mathrm{IR}}\rangle$
$12+\mathrm{log}{({\rm{O}}/{\rm{H}})}_{{\rm{N}}2}$ b	8.02–8.47	80	2.27	12.4±0.5	59 ± 30	472 ± 88	26 ± 1	12 ± 3	0.22 ± 0.05
(${N}_{\mathrm{tot}}=160$)	8.47–8.55	40	2.27	27.0 ± 0.9	347 ± 41	413 ± 126	56 ± 2	25 ± 4	0.22 ± 0.03
	8.55–8.72	40	2.30	34.0 ± 0.9	411 ± 51	486 ± 108	72 ± 2	26 ± 4	0.28 ± 0.04
$12+\mathrm{log}{({\rm{O}}/{\rm{H}})}_{{\rm{O}}3{\rm{N}}2}$	7.98–8.35	86	2.21	11.5 ± 0.4	199 ± 30	374 ± 98	22.3 ± 0.9	16 ± 3	0.14 ± 0.03
(${N}_{\mathrm{tot}}=172$)	8.35–8.44	43	2.26	26.6 ± 0.8	550 ± 50	454 ± 109	54 ± 2	30 ± 3	0.18 ± 0.02
	8.44–8.74	43	2.13	46.3 ± 0.9	616 ± 54	875 ± 110	77 ± 2	33 ± 3	0.24 ± 0.02
${{\rm{O}}}_{32}$	0.25–1.46	42	2.29	37.6 ± 0.8	329 ± 42	607 ± 105	79 ± 2	27 ± 3	0.29 ± 0.04
(${N}_{\mathrm{tot}}=171$)	1.46–2.71	65	2.27	10.8 ± 0.5	52 ± 36	441 ± 101	23 ± 1	11 ± 3	0.21 ± 0.06
	2.71–18.71	64	2.29	3.7 ± 0.4	286 ± 40	209 ± 81	8 ± 1	15 ± 2	0.05 ± 0.01
$\mathrm{log}({M}_{* }/{M}_{\odot })$	8.26–9.6	92	2.15	−3.2 ± 0.4	−125 ± 27	−281 ± 101	2.2c	8.9c	⋯
(${N}_{\mathrm{tot}}=476$)	9.6–10	137	2.20	4.7 ± 0.4	127 ± 27	247 ± 61	9.0 ± 0.7	10 ± 1	0.09 ± 0.02
	10–10.6	173	2.19	21.3 ± 0.5	261 ± 22	306 ± 63	39.5 ± 0.9	16 ± 2	0.24 ± 0.02
	10.6–11.6	74	2.19	69.4 ± 0.9	854 ± 38	1625 ± 102	121.7 ± 1.6	62 ± 3	0.20 ± 0.01
Age (Myr)	71–570	100	2.26	7.6 ± 0.4	172 ± 31	515 ± 84	15.8 ± 0.9	17 ± 3	0.09 ± 0.02
(${N}_{\mathrm{tot}}=299$)	570–1278	100	2.22	15.8 ± 0.5	291 ± 36	−166 ± 97	30.8 ± 1.0	10 ± 3	0.31 ± 0.11
	1278–2600	50	2.18	22.3 ± 1.0	318 ± 44	−88 ± 148	41.0 ± 1.8	15 ± 4	0.27 ± 0.07
	2600–4250	49	2.23	34.3 ± 1.3	263 ± 55	793 ± 131	65.8 ± 2.5	25 ± 5	0.26 ± 0.06

Notes. Entries show properties of stacks in Figures 4 and 6. Each entry includes the parameter range, number of objects in each bin, median redshift, 24, 100, and 160 μm stacked flux densities and their measurement uncertainties, rest-frame 7.7 μm luminosity and its measurement uncertainty, ${L}_{\mathrm{IR}}$ and its uncertainty, and ratio of stacked L_7.7 to stacked ${L}_{\mathrm{IR}}$ and its error.

^a ${L}_{\mathrm{IR}}$ is derived from f₁₀₀ and f₁₆₀ fit to CE01 templates. Its error is the standard deviation of 10,000 ${L}_{\mathrm{IR}}$ realizations that are calculated by fitting IR templates to the perturbed f₁₀₀ and f₁₆₀. If both f₁₀₀ and f₁₆₀ are undetected, we used upper limits, but if only one is undetected, we used the flux and its error to find the best-fit model and the associated ${L}_{\mathrm{IR}}$. ^bN2 stacks are limited to objects at 2.0 ≤ z ≤ 2.6, because otherwise the highest-metallicity bin is dominated by low-z galaxies with $\hat{z}\sim 1.6$. The rest of the stacks (O3N2, ${{\rm{O}}}_{32}$, M_*) include all objects at 1.37 ≤ z ≤ 2.60. The median redshifts in these bins are consistent and are above 2. We should note that the ${{\rm{O}}}_{32}$ sample has only eight galaxies with z < 2.0. ^cStacks are undetected in 24, 100, and 160 μm; hence, L_7.7 and ${L}_{\mathrm{IR}}$ are 3σ upper limits. ${L}_{7.7}/{L}_{\mathrm{IR}}$ is meaningless in this case.

Download table as:  ASCIITypeset image

Another important trend is the anticorrelation between the PAH intensity and ${{\rm{O}}}_{32}$ ratio. ${{\rm{O}}}_{32}$ is used as an empirical proxy for ionization parameter, which is in turn anticorrelated with metallicity (e.g., Pérez-Montero 2014; see also Figure 12 in Sanders et al. 2016 for the relation of ${{\rm{O}}}_{32}$ with O3N2 and N2 in MOSDEF). The weighted-average stacks of 24 μm images in the two highest ${{\rm{O}}}_{32}$ bins yield nondetections. In addition, objects with detected and undetected 24 μm emission are distinctly separated toward lower and higher ${{\rm{O}}}_{32}$, respectively. The fraction of 24 μm undetected objects increases significantly from 38% in the lowest ${{\rm{O}}}_{32}$ bin (median ${{\rm{O}}}_{32}$ = 0.9) to 91% in the highest bin (median ${{\rm{O}}}_{32}$ = 4.5), while the 24 μm nondetection fraction changes from ∼30% at 12 +log(O/H) ∼ 8.6 (8.5) to ∼60%–70% at 12 + log(O/H) ∼ 8.3 (8.2) for N2 (O3N2) metallicities (refer to the bottom panels of Figure 3). These trends suggest that the PAH intensity is strongly affected by the ionization parameter. We will return to this point in Section 6.

Excluding the 24 μm undetected objects from the analysis leads to a significantly biased sample. However, to explore the variations of 24 μm bright sources with ${{\rm{O}}}_{32}$, we stack only the 24 μm detected objects in the same four bins of ${{\rm{O}}}_{32}$ as before. The results indicate that the L_7.7/SFR ratios of 24 μm bright sources are almost independent of ${{\rm{O}}}_{32}$ within the uncertainties (L_7.7/SFR $\sim 2.5\times {10}^{9}\,{L}_{\odot }/{M}_{\odot }\,{\mathrm{yr}}^{-1}$ in all ${{\rm{O}}}_{32}$ bins.)

The PAH–metallicity correlation found at z ∼ 2 is consistent with the trend observed at z ∼ 0 (Engelbracht et al. 2005; Draine et al. 2007; Galliano et al. 2008; Hunt et al. 2010, among others). In local galaxies, there is a paucity of PAH emission at 12+log(O/H) ≲ 8.1–8.2 (Engelbracht et al. 2005; Wu et al. 2006; Draine et al. 2007). In Figure 3(c) we see the same threshold in the ${{\rm{O}}}_{32}$ plot at z ∼ 2: there is a sharp difference in the L_7.7/${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ ratio below and above ${{\rm{O}}}_{32}\sim 2$, which corresponds to $12+\mathrm{log}{({\rm{O}}/{\rm{H}})}_{{{\rm{O}}}_{32}}\sim 8.2$ based on the metallicity calibrations of Jones et al. (2015). The threshold is at $12+\mathrm{log}{({\rm{O}}/{\rm{H}})}_{{\rm{N}}2}\sim 8.5$ for the N2 metallicity calibration. We will discuss this observed threshold in more detail in Section 6.

Additionally, we examine possible variations of the L_7.7/${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ (L_7.7/${L}_{\mathrm{IR}}$) ratio in the BPT diagram ([O iii]/${\rm{H}}\beta$ vs. [N ii]/${\rm{H}}\alpha$; Baldwin et al. 1981), as the position of galaxies below and above the BPT locus is known to be sensitive to the hardness of ionizing radiation (Kewley et al. 2013). We use the locus of z ∼ 2 star-forming galaxies in the BPT diagram from Shapley et al. (2015) and stack 24 μm images of galaxies below and above the locus. Surprisingly, we do not find any evidence for variation of L_7.7 intensity in this parameter space, as galaxies at z ≥ 2 below and above the BPT locus have the same ${L}_{7.7}/{L}_{\mathrm{IR}}$ ratio of 0.20 ± 0.03. It is not clear why we do not see a significant change of the PAH intensity across the BPT diagram; it may be due to the small sample size.

Our results directly show the correlation between the PAH emission and intensity of the radiation field as was speculated by Elbaz et al. (2011, hereafter E11). In Elbaz et al. (2011), the reduced L₈/${L}_{\mathrm{IR}}$ ratio was attributed to star formation "compactness" and "starburstiness" characterized by the SFR surface density and specific SFR, respectively. These authors concluded that the weak PAH emission in compact starbursts compared to their typical star-forming counterparts is a consequence of an increased radiation field intensity in these galaxies, which is directly shown by our results of decreasing PAH intensity with ${{\rm{O}}}_{32}$. In other words, based on the assumption of a constant electron density (as electron density appears to be almost independent from other galaxy properties in MOSDEF; Sanders et al. 2016) and the same ionizing spectrum, there is a tight correspondence between ${{\rm{O}}}_{32}$ and ionization parameter, such that increasing ${{\rm{O}}}_{32}$ reflects a more intense radiation field. This conclusion still holds even if we assume a different (harder) ionizing spectrum for the higher ${{\rm{O}}}_{32}$ bins, which would lower the inferred ionization parameter for a given ${{\rm{O}}}_{32}$ value (see Figure 10 in Sanders et al. 2016).

Estimating the total IR luminosity independent of f₂₄ requires access to longer-wavelength IR data. As mentioned before, in our sample, only very dusty star-forming galaxies with bright dust continua are detected in Herschel/PACS at 100 and 160 μm bands. Therefore, we rely on stacks of images to obtain a sufficiently high S/N to calculate robust IR luminosities.

Figure 4 shows ${L}_{7.7}/{L}_{\mathrm{IR}}$ stacks in bins of metallicity and ${{\rm{O}}}_{32}$. We combine the two lowest-metallicity bins and the two highest-${{\rm{O}}}_{32}$ bins in Figure 3 and Table 1 to gain higher S/N in the PACS stacks. The N2 stacks are adopted only for galaxies at z > 2, because otherwise the highest N2 metallicity bin is dominated by lower-redshift (z ∼ 1.5) galaxies. In O3N2 and ${{\rm{O}}}_{32}$ plots, the median redshifts of all bins are similar and all above z = 2. We note that there are only eight galaxies at z < 2 in the ${{\rm{O}}}_{32}$ sample, as [O ii] is not covered by the MOSFIRE Y filter at z < 1.6 (see Figure 3(c)). Properties of ${L}_{7.7}/{L}_{\mathrm{IR}}$ stacks are listed in Table 2.

Zoom In Zoom Out Reset image size

To place our analysis in the context of other studies, we compare our results with E11, Reddy et al. (2012a, R12), and Wuyts et al. (2008, W08). E11 found that typical main-sequence galaxies between z = 0 and 2.5 follow a Gaussian distribution of ${L}_{\mathrm{IR}}/{L}_{8}$ ($\nu {L}_{\nu }(8\,\mu {\rm{m}})\equiv {L}_{8}$) centered at ${L}_{\mathrm{IR}}/{L}_{8}=4$ ($\sigma =1.6$) with a tail of starburst galaxies with larger ${L}_{\mathrm{IR}}/{L}_{8}$ ratios. The median ${L}_{\mathrm{IR}}/{L}_{8}$ ratio of all galaxies in their sample was ${4.9}_{-2.2}^{+2.9}$. Furthermore, the R12 study found ${L}_{\mathrm{IR}}/{L}_{8}=7.7\pm 1.6$ for a sample of UV-selected galaxies at 1.5 ≤ z < 2.6 (with mean redshift of 2.08). R12 attributed their higher ${L}_{\mathrm{IR}}/{L}_{8}$ to redshift evolution and larger IR luminosity surface densities compared to those of the local galaxies. For a consistent comparison of these ratios with our results, we converted ${L}_{8}/{L}_{\mathrm{IR}}$ ratios of R12 and E11 to ${L}_{7.7}/{L}_{\mathrm{IR}}$ ratios. The E11 and R12 samples have ${L}_{\mathrm{IR}}=5\times {10}^{9}\mbox{--}3\times {10}^{12}\,{L}_{\odot }$ and ${10}^{10}\mbox{--}5\times {10}^{12}\,{L}_{\odot }$, respectively. According to the CE01 templates, $\nu {L}_{\nu }$ at 7.7 μm is higher than $\nu {L}_{\nu }$ at 8 μm by 29% for a template with ${L}_{\mathrm{IR}}=5\times {10}^{9}\,{L}_{\odot }$ and 12% for ${L}_{\mathrm{IR}}=5\times {10}^{12}\,{L}_{\odot }$. The average ${L}_{\mathrm{IR}}$ in the R12 and E11 high-z samples is $\sim 2\times {10}^{11}\,{L}_{\odot }$, which corresponds to an ${L}_{7.7}/{L}_{8}$ ratio of 1.25. Therefore, we adopt a correction factor of 1.25 to convert the E11 and R12 ${L}_{8}/{L}_{\mathrm{IR}}$ ratios to ${L}_{7.7}/{L}_{\mathrm{IR}}$ ratios. The variations of the ${L}_{7.7}/{L}_{\mathrm{IR}}$ ratios that resulted from adopting other conversion factors (from 1.12 to 1.29) are well within the reported 1σ measurement uncertainty of ${L}_{7.7}/{L}_{\mathrm{IR}}$. After the conversion, we obtain ${L}_{7.7}/{L}_{\mathrm{IR}}={0.25}_{-0.11}^{+0.15}$ and 0.16 ± 0.03 for the E11 and R12 studies, respectively. These ratios and their confidence intervals are shown with lines and shaded regions in Figure 4. The ${L}_{7.7}/{L}_{\mathrm{IR}}$ of E11 is consistent with our highest-metallicity and lowest-${{\rm{O}}}_{32}$ stacks. This is expected as at the same stellar mass, low-redshift galaxies have typically higher metallicities and less intense ionizing radiation fields. The ${L}_{7.7}/{L}_{\mathrm{IR}}$ of R12 is in agreement with our lower-metallicity stacks and the middle-${{\rm{O}}}_{32}$ stack, as well as with the stacks of the full sample (yellow stars).

We also compare our results with those of Wuyts et al. (2008, 2011a). In these studies, the authors derived an empirical IR template and luminosity-independent monochromatic conversions from 24 μm flux density to ${L}_{\mathrm{IR}}$ as a function of redshift. We used the median redshifts of our three ${{\rm{O}}}_{32}$ stacks (from the lowest to the highest bin: z = 2.29, 2.27, and 2.29; Table 2) to adopt appropriate f₂₄-to-${L}_{\mathrm{IR}}$ conversions and compare them with our results in Figure 4(d). The width of the purple line indicates the range of the W08 ${f}_{24}/{L}_{\mathrm{IR}}$ ratios for different redshifts of our three stacks. As expected, f₂₄/${L}_{\mathrm{IR}}$ of W08 is consistent with our middle-${{\rm{O}}}_{32}$ bin, as well as the stack of all objects, but there is discrepancy with the low- and high-${{\rm{O}}}_{32}$ stacks. The same conclusion holds for the metallicity stacks. In Section 5, we will discuss the implications of the inconsistency between the L_7.7-to-${L}_{\mathrm{IR}}$ conversions of W08, E11, and R12 and our observed values at low stellar masses.

In Shivaei et al. (2015a) we discussed how SFR–M_* studies at z ∼ 2 that rely on ${\rm{H}}\alpha$ SFRs predict shallower log(SFR)–log(M_*) slopes (∼0.5–0.7; Erb et al. 2006; Zahid et al. 2012; Koyama et al. 2013; Atek et al. 2014; Darvish et al. 2014; Shivaei et al. 2015a), compared to those that utilize IR SFRs (slope of ∼1; Daddi et al. 2007; Santini et al. 2009; Wuyts et al. 2011b; Reddy et al. 2012b; Whitaker et al. 2014). In order to understand the cause of this discrepancy, in Shivaei et al. (2016) we compared ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ with SFRs inferred from panchromatic SED models that were fit to rest-frame UV to far-IR photometry. The comparison indicated that ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ does not underpredict total SFR even for the most star-forming and dusty galaxies in our sample (up to $\sim 300\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$). In this section, we examine whether the discrepancy in the slope of the SFR–M_* relation arises from the 24 μm inferred IR SFRs that are not corrected for the PAH–metallicity effect.

Figure 7 shows SFR in bins of stellar mass for galaxies with ${M}_{* }\geqslant {10}^{9.6}\,{M}_{\odot }$. In Figure 7, orange circles denote mean ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$, and other symbols are mean ${\mathrm{SFR}}_{\mathrm{UV}}$ added to mean ${\mathrm{SFR}}_{\mathrm{IR}}$, where the latter is derived from various conversions of f₂₄ to ${L}_{\mathrm{IR}}$ adopted from the literature (diamonds) and from this analysis (blue stars), which is described as follows. We derive a mass-dependent conversion of L_7.7 to ${L}_{\mathrm{IR}}$ based on our stacks presented in Figures 6(b) and (c) and Table 2. The ${L}_{7.7}/{L}_{\mathrm{IR}}$ ratio is 0.09, 0.24, and 0.20 for $\mathrm{log}({M}_{* }/{{\rm{M}}}_{\odot })=9.6\mbox{--}10.0$, 10.0–10.6, and 10.6–11.6, respectively. We adopt ${L}_{7.7}/{L}_{\mathrm{IR}}$ = 0.09 for ${M}_{* }\lt {10}^{10}\,{M}_{\odot }$ and an average of the two highest-mass bins, 0.22, for ${M}_{* }\geqslant {10}^{10}\,{M}_{\odot }$.⁸ The result is shown with blue stars in Figure 7. The diamonds in Figure 7 reflect results of the mass-independent 24 μm-to-IR conversions of R12, W08, and E11. The 24 μm stack below ${10}^{9.6}\,{M}_{\odot }$ (shaded gray region in Figure 7) is not robustly detected, owing to the reduced dust emission at low masses and low SFRs.

Zoom In Zoom Out Reset image size

The ${\mathrm{SFR}}_{\mathrm{UV}}+{\mathrm{SFR}}_{\mathrm{IR}}$ derived from the mass-dependent 24 μm-to-IR conversion presented in this study is in a very good agreement with ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$, particularly considering that these two SFRs are derived completely independently from each other. As expected, the other IR SFRs that are inferred from a single 24 μm-to-IR conversion are underestimated for low-mass (${M}_{* }\lt {10}^{10}\,{M}_{\odot }$) galaxies. Figure 7 clearly demonstrates how using a single 24 μm-to-IR conversion can underestimate SFRs at lower masses and lead to a steeper log(SFR)–log(M_*) slope compared to that derived in this study. A simple linear least-squares regression calculation gives a slope of 1.0 ± 0.1 for the stacks of ${\mathrm{SFR}}_{\mathrm{UV}}+{\mathrm{SFR}}_{\mathrm{IR}}$ that assume the W08, E11, and R12 conversions and a slope of 0.7 ± 0.1 for the ${\mathrm{SFR}}_{{\rm{H}}\alpha ,{\rm{H}}\beta }$ and ${\mathrm{SFR}}_{\mathrm{UV}}+{\mathrm{SFR}}_{\mathrm{IR}}$ stacks that adopt our mass-dependent conversion.

The slope of 0.7 found in this study is consistent with observations at lower redshifts (slope ∼0.7 at z ∼ 0.1 and 1; Salim et al. 2007; Noeske et al. 2007, respectively), suggesting no evolution in the slope of the SFR–M_* relation from z ∼ 0 to 2. Moreover, the increased ${\mathrm{SFR}}_{\mathrm{UV}}+{\mathrm{SFR}}_{\mathrm{IR}}$ of the low-mass galaxies estimated from our mass-dependent 24 μm-to-IR conversion removes evidence for a "flattening" at the high-mass end of the SFR–M_* relation, as suggested in some previous studies (e.g., Whitaker et al. 2014; Schreiber et al. 2015). However, it is plausible that the slope of the relation steepens at masses lower than those accessible in this study—a mass-independent shallow slope of 0.7 would lead to a too rapid evolution in the stellar mass function at low masses (Weinmann et al. 2012; Leja et al. 2015).

Our results indicate that specific SFR (sSFR ≡ SFR/M_*) is a decreasing function of M_*. If we assume that the gas and stars occupy the same regions in the galaxies, then the Schmidt relation (Kennicutt 1998) implies that $\mathrm{SFR}\propto {M}_{\mathrm{gas}}^{1.4}$. Consequently, sSFR is related to the cold gas fraction defined as ${M}_{\mathrm{gas}}/({M}_{* }+{M}_{\mathrm{gas}})$ (Reddy et al. 2006a). The negative slope of the log(sSFR)–log(M_*) relation indicates that, as expected, there is a correspondence between the cold gas fraction and stellar mass.

There is tension between our results and those from simulations. According to simulations, the slope of the SFR–M_* relation is close to unity and the sSFR is almost a constant function of mass (Dutton et al. 2010a; Davé et al. 2011; Kannan et al. 2014; Torrey et al. 2014; Sparre et al. 2015, among many others). Furthermore, our analysis shows an increased sSFR at ${M}_{* }=5\times {10}^{9}$ at z ∼ 2 by a factor of 1.8 ($\mathrm{log}(\mathrm{sSFR}/\mathrm{Gyr})=0.4$), compared to the previous sSFR measurements that used mass-independent 24 μm-to-IR conversions. The increased sSFR found here also aggravates the inconsistency between observations and models, where the models predict lower sSFRs at z ∼ 2 than those reported in previous studies (Furlong et al. 2015; Sparre et al. 2015; Davé et al. 2016). From the observational viewpoint, we showed in Shivaei et al. (2015b) that our MOSDEF ${\rm{H}}\alpha$-detected sample is complete down to ${M}_{* }={10}^{9.5}\,{M}_{\odot }$. Therefore, the shallow slope of 0.7 in this study is immune from the common sample selection biases at low masses. On the other hand, Davé (2008) and Torrey et al. (2014) demonstrated that changing feedback models in simulations does not affect the slope of the relation, as it moves galaxies along the SFR–M_* relation. More investigations are needed to resolve these issues and improve the agreement between simulations and observations.

Luminous infrared galaxies with ${L}_{\mathrm{IR}}={10}^{11}\mbox{--}{10}^{12}\,{L}_{\odot }$ are known to contribute significantly to the total IR luminosity density at z ∼ 2 (Reddy et al. 2008; Rodighiero et al. 2010; Magnelli et al. 2011). Current measurements of the IR luminosity density (and the obscured SFR density) at z ∼ 2 either are based solely on ultraluminous infrared galaxies (ULIRGs) with ${L}_{\mathrm{IR}}\gt {10}^{12}\,{L}_{\odot }$ (e.g., Gruppioni et al. 2013) or have included lower IR luminosities by relying on different methods of converting 24 μm flux density to ${L}_{\mathrm{IR}}$ (e.g., Pérez-González et al. 2005; Caputi et al. 2007; Reddy et al. 2008, 2010; Rodighiero et al. 2010; Riguccini et al. 2011). The most common method uses empirical IR templates that are calibrated with local galaxies (e.g., CE01; Dale & Helou 2002; Rieke et al. 2009). These IR templates are shown to overestimate ${L}_{\mathrm{IR}}$ for ULIRGs at z ∼ 2 (Nordon et al. 2010; Elbaz et al. 2011; Reddy et al. 2012a). To overcome this so-called "mid-IR excess," other 24 μm-to-IR conversions that were calibrated with samples of high-redshift galaxies were introduced (Bavouzet et al. 2008; Elbaz et al. 2011, among many others). However, the dependence of ${L}_{7.7}/{L}_{\mathrm{IR}}$ on metallicity and stellar mass was neglected in these calibrations. As we showed in previous sections, these conversions underpredict ${L}_{\mathrm{IR}}$ at low and moderate masses (${M}_{* }\lesssim {10}^{10}\,{M}_{\odot }$) and IR luminosities (${L}_{\mathrm{IR}}\lesssim {10}^{11}\,{L}_{\odot }$).

The increased ${L}_{\mathrm{IR}}$ of the low and moderately IR-luminous galaxies leads to a higher IR luminosity density at z ∼ 2. To quantify this effect, we simulate 1000 galaxies drawn from the double-power-law luminosity function of Magnelli et al. (2011), increase the IR luminosities of galaxies below ${L}_{\mathrm{IR}}={10}^{11}\,{L}_{\odot }$ in accordance with the relations found in this study (as discussed below), bin the galaxies (we used bins with 0.3 dex width, similar to Magnelli et al. 2011), and refit the bins with the double-power-law function. We use two alternative methods to include the increased ${L}_{\mathrm{IR}}$ of low-luminosity galaxies. In the first method, we increase ${L}_{\mathrm{IR}}$ by a factor of 2 for all simulated galaxies with ${L}_{\mathrm{IR}}\leqslant {10}^{11}\,{L}_{\odot }$. This assumption is motivated by our finding for the low-mass galaxies in the previous section (Section 5.1). However, from our data it is not clear whether ${L}_{\mathrm{IR}}$ of faint IR galaxies should be increased by a greater factor or an equal factor compared to that of the moderately IR-luminous galaxies. Therefore, as a sanity check we adopt a second approach, where ${L}_{\mathrm{IR}}$ is increased by a factor of 2 for galaxies with ${L}_{\mathrm{IR}}\leqslant {10}^{10}\,{L}_{\odot }$ and by a factor of 1.5 for those with ${10}^{10}\lt {L}_{\mathrm{IR}}\leqslant {10}^{11}\,{L}_{\odot }$. In both cases, the number density of IR bins below ${10}^{11}\,{L}_{\odot }$ increases by ≳0.3 dex, and the faint-end slope steepens from −0.6 to ∼−0.7. The increase in the slope is within its 1σ error of 0.1 (Sanders et al. 2003; Magnelli et al. 2011). A stronger variation in ${L}_{\mathrm{IR}}$ would cause a steeper faint-end slope, resulting in a higher fractional contribution of galaxies with ${L}_{\mathrm{IR}}\lesssim {10}^{11}\,{L}_{\odot }$ to the total IR budget and obscured SFR density at z ∼ 2 than predicted before (Caputi et al. 2007; Reddy et al. 2010).

If we assume the same fractional contribution of the faint and moderately luminous IR galaxies to the bolometric (i.e., IR+UV) luminosity density at z ∼ 2, as derived in previous studies (∼50%–80%; Caputi et al. 2007; Reddy et al. 2008; Rodighiero et al. 2010), our estimate of a factor of ∼2 increase in the IR luminosity of galaxies with ${L}_{\mathrm{IR}}\lesssim {10}^{11}\,{L}_{\odot }$ would increase the IR luminosity density by 30%. Applying the 30% increase to the IR luminosity density measurement of Reddy et al. (2008) at z ∼ 2 and converting it to IR SFR density using the Kennicutt (1998) relation, modified for a Chabrier (2003) IMF,⁹ yields an IR SFR density of 0.15 ${M}_{\odot }$ yr⁻¹ Mpc⁻³. Adding this value to the UV SFR density measured in the same study (Reddy et al. 2008) results in a total SFR density of 0.18 ${M}_{\odot }$ yr⁻¹ Mpc⁻³ (assuming a Chabrier 2003 IMF), which is ∼30% higher than the initial value calculated from the unchanged IR luminosity density. This finding exacerbates the discrepancy between the stellar mass density at z ∼ 2 measured from observations and that obtained by integrating the SFR density over time (Hopkins & Beacom 2006; Madau & Dickinson 2014).